text:
struct Regularize{N,F}

  "x values of points, normalized to grid index space"
  x :: Vector{Float64}

  "y values of points, normalized to grid index space"
  y :: Vector{Float64}

  "1/dV factor"
  overdv :: Float64

  "weights for each point (e.g. arclengths), divided by dV"
  wgt :: Vector{Float64}

  "Discrete Delta function"
  ddf :: AbstractDDF

  "Radius of the DDF"
  ddf_radius :: Float64

  "Symmetry flag"
  _issymmetric :: Bool

end

"""
    Regularize(x,y,dx,[ddftype=Yang3],[graddir=0],[I0=(1,1)], [weights=1.0], [filter=false],
                       [issymmetric=false])

Constructor to set up an operator for regularizing and interpolating data from/to
points immersed in the grid to/from fields on the grid itself. The supplied
`x` and `y` represent physical coordinates of the immersed points, and `dx`
denotes a uniform physical cell size of the grid. The separate arguments `x` and
`y` can be replaced by a single argument `X` of type `VectorData` holding the
coordinates.

The operations of regularization and interpolation are carried out with a discrete
delta function (ddf), which defaults to the type `Yang3`. Others are also possible,
such as `Roma`, `Goza` or `M3`. The optional argument `graddir`, if set to 1 or 2, will
generate an interpolation operator that evaluates the negative of the
respective component of the gradient of a grid field at the immersed points. The
default value of this argument is 0, which simply interpolates. Note that the
regularization form of this gradient type is also possible.

The optional tuple
`I0` represents the indices of the primary node that coincides with `(x,y) = (0,0)`.
This defaults to `(1,1)`, which leaves one layer of ghost (dual) cells and sets
the physical origin in the lower left corner of the grid of interior dual cells.

Another optional parameter, `weights`, sets the weight of each point in the
regularization. This would generally be set with, say, the differential arc
length for regularization of data on a curve. It can be a vector (of the same length
as x and y) or a scalar if uniform. It defaults to 1.0.

The optional Boolean parameter `filter` can be set to `true` if it is desired to
apply filtering (see Goza et al, J Comput Phys 2016) to the grid data before
interpolating. This is generally only used in the context of preconditioning
the solution for forces on the immersed points.

If the optional Boolean parameter `issymmetric` is set to `true`, then the
regularization and interpolation are constructed to be transposes of each other.
Note that this option overrides any supplied weights. The default of this
parameter is `false`.

The resulting operator can be used in either direction, regularization and
interpolation, with the first argument representing the *target* (the entity
to regularize/interpolate to), and the second argument
the *source* (the entity to regularize/interpolate from). The regularization
does not use the filtering option.

# Example

In the example below, we set up a 12 x 12 grid. Using the default value for `I0`
and setting `dx = 0.1`, the physical dimensions of the non-ghost part of the grid
are 1.0 x 1.0. Three points are set up in the interior, and a vector field is assigned
to them, with the x component of each of them set to 1.0. These data are regularized
to a field of primal edges on the grid, using the Roma DDF kernel.

```jldoctest
julia> x = [0.25,0.75,0.25]; y = [0.75,0.25,0.25];

julia> X = VectorData(x,y);

julia> q = Edges(Primal,(12,12));

julia> dx = 0.1;

julia> H = Regularize(x,y,dx;ddftype=Roma)
Regularization/interpolation operator with non-filtered interpolation
  DDF type CartesianGrids.Roma
  3 points in grid with cell area 0.01

julia> f = VectorData(X);

julia> fill!(f.u,1.0);

julia> H(q,f)
Edges{Primal,12,12,Float64} data
u (in grid orientation)
11×12 Array{Float64,2}:
 0.0  0.0  0.0       0.0     0.0      …  0.0       0.0     0.0      0.0  0.0
 0.0  0.0  0.0       0.0     0.0         0.0       0.0     0.0      0.0  0.0
 0.0  0.0  8.33333  33.3333  8.33333     0.0       0.0     0.0      0.0  0.0
 0.0  0.0  8.33333  33.3333  8.33333     0.0       0.0     0.0      0.0  0.0
 0.0  0.0  0.0       0.0     0.0         0.0       0.0     0.0      0.0  0.0
 0.0  0.0  0.0       0.0     0.0      …  0.0       0.0     0.0      0.0  0.0
 0.0  0.0  0.0       0.0     0.0         0.0       0.0     0.0      0.0  0.0
 0.0  0.0  8.33333  33.3333  8.33333     8.33333  33.3333  8.33333  0.0  0.0
 0.0  0.0  8.33333  33.3333  8.33333     8.33333  33.3333  8.33333  0.0  0.0
 0.0  0.0  0.0       0.0     0.0         0.0       0.0     0.0      0.0  0.0
 0.0  0.0  0.0       0.0     0.0      …  0.0       0.0     0.0      0.0  0.0
v (in grid orientation)
12×11 Array{Float64,2}:
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
```
"""
function Regularize(x::AbstractVector{T},y::AbstractVector{T},dx::T;
                    ddftype=Yang3,graddir::Int=0,
                    I0::Tuple{Int,Int}=(1,1),
                    weights::Union{T,Vector{T}}=1.0,
                    filter::Bool = false,
                    issymmetric::Bool = false) where {T<:Real}

  _issymmetric = (filter ? false : issymmetric)

  n = length(x)
  @assert length(y)==n
  if !_issymmetric
    if typeof(weights) == T
      wtvec = similar(x)
      fill!(wtvec,weights/(dx*dx))
    else
      @assert length(weights)==n
      wtvec = deepcopy(weights)./(dx*dx)
    end
  else
    # if the regularization and interpolation are symmetric, then the
    # weights are automatically set to be the cell area in order to cancel it
    # in the denominator of the regularization operator.
    wtvec = similar(x)
    fill!(wtvec,1.0)
  end

  baseddf = DDF(ddftype=ddftype,dx=1.0)
  if graddir == 0
    ddf = baseddf
  else
    ddf = GradDDF(graddir,ddftype=ddftype,dx=1.0)
  end

  Regularize{length(x),filter}(x/dx.+I0[1],y/dx.+I0[2],1.0/(dx*dx),
                      wtvec,ddf,_get_regularization_radius(baseddf),_issymmetric)
end

Regularize(x::T,y::T,a...;b...) where {T<:Real} = Regularize([x],[y],a...;b...)

Regularize(x::VectorData,a...;b...) = Regularize(x.u,x.v,a...;b...)

function Base.show(io::IO, H::Regularize{N,F}) where {N,F}
    filter = F ? "filtered" : "non-filtered"
    ddftype,_ = typeof(H.ddf).parameters
    op = H._issymmetric ? "Symmetric regularization/interpolation" : "Regularization/interpolation"
    println(io, "$op operator with $filter interpolation")
    println(io, "  DDF type $ddftype")
    println(io, "  $N points in grid with cell area $(sprint(show,1.0/H.overdv;context=:compact => true))")
end

function _get_regularization_radius(ddf::DDF)
  v = 1.0
  r = 0.0
  dr = 0.01
  while (v = abs(ddf(r))) > eps()
    r += dr
  end
  return r
end

"""
    RegularizationMatrix(H::Regularize,f::PointData,u::CellData) -> Hmat

Construct and store a matrix representation of regularization associated with `H`
for data of type `f` to data of type `u`. The resulting matrix `Hmat` can then be
used to apply on point data of type `f` to regularize it to grid data of type `u`,
using `mul!(u,Hmat,f)`. It can also be used as just `Hmat*f`.

If `H` is a symmetric regularization and interpolation operator, then this
actually returns a tuple `Hmat, Emat`, where `Emat` is the interpolation matrix.
"""
struct RegularizationMatrix{TU,TF} <: AbstractMatrix{Float64}
  M :: SparseMatrixCSC{Float64,Int64}
end

"""
    InterpolationMatrix(H::Regularize,u::CellData,f::PointData) -> Emat

Construct and store a matrix representation of interpolation associated with `H`
for data of type `u` to data of type `f`. The resulting matrix `Emat` can then be
used to apply on grid data of type `u` to interpolate it to point data of type `f`,
using `mul!(f,Emat,u)`. It can also be used as just `Emat*u`.
"""
struct InterpolationMatrix{TU,TF} <: AbstractMatrix{Float64}
  M :: SparseMatrixCSC{Float64,Int64}
end

@wraparray RegularizationMatrix M 2
@wraparray InterpolationMatrix M 2

# ======  Regularization and interpolation operators of scalar types ======== #

pointtype = :ScalarData
for (gridtype,ctype,dnx,dny,shiftx,shifty) in @generate_scalarlist(SCALARLIST)

# Regularization
  @eval function (H::Regularize{N,F})(target::$gridtype{$ctype,NX,NY,T,DDT},source::$pointtype{N,S,DT}) where {N,F,NX,NY,S,T,DT,DDT}
        radius = H.ddf_radius
        fill!(target,0.0)
        xmin, xmax = -radius+$shiftx, radius+$shiftx
        ymin, ymax = -radius+$shifty, radius+$shifty
        @inbounds for pt in 1:N
            prangex = max(1,ceil(Int,H.x[pt]+xmin)):min(NX-$dnx,floor(Int,H.x[pt]+xmax))
            prangey = max(1,ceil(Int,H.y[pt]+ymin)):min(NY-$dny,floor(Int,H.y[pt]+ymax))
            for y in prangey, x in prangex
                target[x,y] += source[pt]*H.wgt[pt]*H.ddf(x-$shiftx-H.x[pt],y-$shifty-H.y[pt])
            end
        end
        target
  end

# Interpolation
  @eval function (H::Regularize{N,false})(target::$pointtype{N,S,DT},source::$gridtype{$ctype,NX,NY,T,DDT}) where {N,NX,NY,S,T,DT,DDT}
        radius = H.ddf_radius
        fill!(target,0.0)
        xmin, xmax = -radius+$shiftx, radius+$shiftx
        ymin, ymax = -radius+$shifty, radius+$shifty
        @inbounds for pt in 1:N
            prangex = max(1,ceil(Int,H.x[pt]+xmin)):min(NX-$dnx,floor(Int,H.x[pt]+xmax))
            prangey = max(1,ceil(Int,H.y[pt]+ymin)):min(NY-$dny,floor(Int,H.y[pt]+ymax))
            for y in prangey, x in prangex
                target[pt] += source[x,y]*H.ddf(x-$shiftx-H.x[pt],y-$shifty-H.y[pt])
            end
        end
        target
  end

# Interpolation with filtering
  @eval function (H::Regularize{N,true})(target::$pointtype{N,S,DT},source::$gridtype{$ctype,NX,NY,T,DDT}) where {N,NX,NY,S,T,DT,DDT}
        tmp = typeof(source)()
        radius = H.ddf_radius
        fill!(target,0.0)
        xmin, xmax = -radius+$shiftx, radius+$shiftx
        ymin, ymax = -radius+$shifty, radius+$shifty
        @inbounds for pt in 1:N
            prangex = max(1,ceil(Int,H.x[pt]+xmin)):min(NX-$dnx,floor(Int,H.x[pt]+xmax))
            prangey = max(1,ceil(Int,H.y[pt]+ymin)):min(NY-$dny,floor(Int,H.y[pt]+ymax))
            for y in prangey, x in prangex
                tmp[x,y] += H.wgt[pt]*H.ddf(x-$shiftx-H.x[pt],y-$shifty-H.y[pt])
            end
        end
        @inbounds for pt in 1:N
            prangex = max(1,ceil(Int,H.x[pt]+xmin)):min(NX-$dnx,floor(Int,H.x[pt]+xmax))
            prangey = max(1,ceil(Int,H.y[pt]+ymin)):min(NY-$dny,floor(Int,H.y[pt]+ymax))
            @inbounds for y in prangey, x in prangex
                w = tmp[x,y]
                w = w > eps() ? source[x,y]/w : 0.0
                target[pt]  += w*H.ddf(x-$shiftx-H.x[pt],y-$shifty-H.y[pt])
            end
        end
        target
  end

  # Construct regularization matrix
  @eval function RegularizationMatrix(H::Regularize{N,F},
    f::$pointtype{N,S,DT},
    u::$gridtype{$ctype,NX,NY,T,DDT}) where {N,F,NX,NY,S,T,DT,DDT}

    Hmat = spzeros(length(u),length(f))
    g = deepcopy(f)
    v = deepcopy(u)
    fill!(g,0.0)
    for i = 1:N
      g[i] = 1.0
      Hmat[:,i] = sparsevec(H(v,g))
      g[i] = 0.0
    end
    if H._issymmetric
      # In symmetric case, these matrices are identical. (Interpolation is stored
      # as its transpose.)
      return RegularizationMatrix{typeof(u),typeof(f)}(Hmat),InterpolationMatrix{typeof(u),typeof(f)}(Hmat)
    else
      return RegularizationMatrix{typeof(u),typeof(f)}(Hmat)
    end
  end

  # Construct interpolation matrix
  @eval function InterpolationMatrix(H::Regularize{N,false},
    u::$gridtype{$ctype,NX,NY,T,DDT},
    f::$pointtype{N,S,DT}) where {N,NX,NY,S,T,DT,DDT}

    Emat = spzeros(length(u),length(f))
    g = deepcopy(f)
    v = deepcopy(u)
    fill!(g,0.0)
    for i = 1:N
      g[i] = 1.0/H.wgt[i]  # unscale for interpolation
      Emat[:,i] = sparsevec(H(v,g))
      g[i] = 0.0
    end
    InterpolationMatrix{typeof(u),typeof(f)}(Emat)
  end

  # Construct interpolation matrix with filtering. Because we construct
  # the matrix by supplying unit vectors from the points, it is easier
  # to carry out the filtering here rather than call the filtered operator
  # we defined above
  @eval function InterpolationMatrix(H::Regularize{N,true},
    u::$gridtype{$ctype,NX,NY,T,DDT},
    f::$pointtype{N,S,DT}) where {N,NX,NY,S,T,DT,DDT}

    Emat = spzeros(length(u),length(f))
    g = deepcopy(f)
    v = deepcopy(u)
    fill!(g,1.0)
    wt = sparsevec(H(v,g))
    wt.nzval .= 1 ./ wt.nzval
    fill!(g,0.0)
    for i = 1:N
      g[i] = 1.0/H.wgt[i]  # unscale for interpolation
      Emat[:,i] = wt.*sparsevec(H(v,g))
      g[i] = 0.0
    end
    InterpolationMatrix{typeof(u),typeof(f)}(Emat)
  end

end

# ===== Regularization and interpolation operators of vector data to edges ===== #
pointtype = :VectorData
for (gridtype,ctype,dunx,duny,dvnx,dvny,shiftux,shiftuy,shiftvx,shiftvy) in @generate_collectionlist(VECTORLIST)

# Regularization
  @eval function (H::Regularize{N,F})(target::$gridtype{$(ctype...),NX,NY,T,DDT},source::$pointtype{N,S,DT}) where {N,F,NX,NY,S,T,DT,DDT}
        H(target.u,source.u)
        H(target.v,source.v)
        target
  end

# Interpolation
  @eval function (H::Regularize{N,F})(target::$pointtype{N,S,DT},source::$gridtype{$(ctype...),NX,NY,T,DDT}) where {N,F,NX,NY,S,T,DT,DDT}
        H(target.u,source.u)
        H(target.v,source.v)
        target
  end

  # Construct regularization matrix
  @eval function RegularizationMatrix(H::Regularize{N,F},src::$pointtype{N,S,DT},target::$gridtype{$(ctype...),NX,NY,T,DDT}) where {N,F,NX,NY,S,T,DT,DDT}

    lenu = length(target.u)
    lenv = length(target.v)
    Hmat = spzeros(lenu+lenv,2N)

    if H._issymmetric
      # In symmetric case, these matrices are identical. (Interpolation is stored
      # as its transpose.)
      Hmat[1:lenu,          1:N]    = RegularizationMatrix(H,src.u,target.u)[1].M
      Hmat[lenu+1:lenu+lenv,N+1:2N] = RegularizationMatrix(H,src.v,target.v)[1].M
      return RegularizationMatrix{typeof(target),typeof(src)}(Hmat),InterpolationMatrix{typeof(target),typeof(src)}(Hmat)
    else
      Hmat[1:lenu,          1:N]    = RegularizationMatrix(H,src.u,target.u).M
      Hmat[lenu+1:lenu+lenv,N+1:2N] = RegularizationMatrix(H,src.v,target.v).M
      return RegularizationMatrix{typeof(target),typeof(src)}(Hmat)
    end
  end

  # Construct interpolation matrix
  @eval function InterpolationMatrix(H::Regularize{N,F},src::$gridtype{$(ctype...),NX,NY,T,DDT},target::$pointtype{N,S,DT}) where {N,F,NX,NY,S,T,DT,DDT}

    # note that we store interpolation matrices in the same shape as regularization matrices
    lenu = length(src.u)
    lenv = length(src.v)
    Emat = spzeros(lenu+lenv,2N)
    Emat[1:lenu,          1:N]    = InterpolationMatrix(H,src.u,target.u).M
    Emat[lenu+1:lenu+lenv,N+1:2N] = InterpolationMatrix(H,src.v,target.v).M
    InterpolationMatrix{typeof(src),typeof(target)}(Emat)
  end

end

# ======  Regularization and interpolation operators of tensor data to edge gradients ==== #
# Here, u describes both diagonal components and v the off diagonal
# We do not use the scalar-wise operations here because there is double use of
# the same ddf evaluation, so this saves us quite a bit of time.

pointtype = :TensorData
for (gridtype,ctype,dunx,duny,dvnx,dvny,shiftux,shiftuy,shiftvx,shiftvy) in @generate_collectionlist(TENSORLIST)

# Regularization
  @eval function (H::Regularize{N,F})(target::$gridtype{$(ctype...),NX,NY,T,DDT},source::$pointtype{N,S,DT}) where {N,F,NX,NY,S,T,DT,DDT}
        radius = H.ddf_radius
        fill!(target.dudx,0.0)
        fill!(target.dvdy,0.0)
        xmin, xmax = -radius+$shiftux, radius+$shiftux
        ymin, ymax = -radius+$shiftuy, radius+$shiftuy
        @inbounds for pt in 1:N
            prangex = max(1,ceil(Int,H.x[pt]+xmin)):min(NX-$dunx,floor(Int,H.x[pt]+xmax))
            prangey = max(1,ceil(Int,H.y[pt]+ymin)):min(NY-$duny,floor(Int,H.y[pt]+ymax))
            for y in prangey, x in prangex
                d = H.wgt[pt]*H.ddf(x-$shiftux-H.x[pt],y-$shiftuy-H.y[pt])
                target.dudx[x,y] += source.dudx[pt]*d
                target.dvdy[x,y] += source.dvdy[pt]*d
            end
        end
        fill!(target.dudy,0.0)
        fill!(target.dvdx,0.0)
        xmin, xmax = -radius+$shiftvx, radius+$shiftvx
        ymin, ymax = -radius+$shiftvy, radius+$shiftvy
        @inbounds for pt in 1:N
            prangex = max(1,ceil(Int,H.x[pt]+xmin)):min(NX-$dvnx,floor(Int,H.x[pt]+xmax))
            prangey = max(1,ceil(Int,H.y[pt]+ymin)):min(NY-$dvny,floor(Int,H.y[pt]+ymax))
            for y in prangey, x in prangex
                d = H.wgt[pt]*H.ddf(x-$shiftvx-H.x[pt],y-$shiftvy-H.y[pt])
                target.dudy[x,y] += source.dudy[pt]*d
                target.dvdx[x,y] += source.dvdx[pt]*d
            end
        end
        target
  end

  # @eval function (H::Regularize{N,F})(target::$ctype,source::$ftype) where {N,F,NX,NY,S,T}
  #   H(target.dudx,source.dudx)
  #   H(target.dudy,source.dudy)
  #   H(target.dvdx,source.dvdx)
  #   H(target.dvdy,source.dvdy)
  #   target
  # end

# Interpolation
  @eval function (H::Regularize{N,false})(target::$pointtype{N,S,DT},source::$gridtype{$(ctype...),NX,NY,T,DDT}) where {N,NX,NY,S,T,DT,DDT}
        radius = H.ddf_radius
        fill!(target.dudx,0.0)
        fill!(target.dvdy,0.0)
        xmin, xmax = -radius+$shiftux, radius+$shiftux
        ymin, ymax = -radius+$shiftuy, radius+$shiftuy
        @inbounds for pt in 1:N
            prangex = max(1,ceil(Int,H.x[pt]+xmin)):min(NX-$dunx,floor(Int,H.x[pt]+xmax))
            prangey = max(1,ceil(Int,H.y[pt]+ymin)):min(NY-$duny,floor(Int,H.y[pt]+ymax))
            for y in prangey, x in prangex
                d = H.ddf(x-$shiftux-H.x[pt],y-$shiftuy-H.y[pt])
                target.dudx[pt] += source.dudx[x,y]*d
                target.dvdy[pt] += source.dvdy[x,y]*d
            end
        end
        fill!(target.dudy,0.0)
        fill!(target.dvdx,0.0)
        xmin, xmax = -radius+$shiftvx, radius+$shiftvx
        ymin, ymax = -radius+$shiftvy, radius+$shiftvy
        @inbounds for pt in 1:N
            prangex = max(1,ceil(Int,H.x[pt]+xmin)):min(NX-$dvnx,floor(Int,H.x[pt]+xmax))
            prangey = max(1,ceil(Int,H.y[pt]+ymin)):min(NY-$dvny,floor(Int,H.y[pt]+ymax))
            for y in prangey, x in prangex
                d = H.ddf(x-$shiftvx-H.x[pt],y-$shiftvy-H.y[pt])
                target.dudy[pt] += source.dudy[x,y]*d
                target.dvdx[pt] += source.dvdx[x,y]*d
            end
        end
        target
  end

  # @eval function (H::Regularize{N,F})(target::$ftype,source::$ctype) where {N,F,NX,NY,S,T}
  #   H(target.dudx,source.dudx)
  #   H(target.dudy,source.dudy)
  #   H(target.dvdx,source.dvdx)
  #   H(target.dvdy,source.dvdy)
  #   target
  # end

# Interpolation with filtering -- need to speed up
  @eval function (H::Regularize{N,true})(target::$pointtype{N,S,DT},source::$gridtype{$(ctype...),NX,NY,T,DDT}) where {N,NX,NY,S,T,DT,DDT}
        tmp = typeof(source)()
        radius = H.ddf_radius
        fill!(target.dudx,0.0)
        fill!(target.dvdy,0.0)
        xmin, xmax = -radius+$shiftux, radius+$shiftux
        ymin, ymax = -radius+$shiftuy, radius+$shiftuy
        @inbounds for pt in 1:N
            prangex = max(1,ceil(Int,H.x[pt]+xmin)):min(NX-$dunx,floor(Int,H.x[pt]+xmax))
            prangey = max(1,ceil(Int,H.y[pt]+ymin)):min(NY-$duny,floor(Int,H.y[pt]+ymax))
            for y in prangey, x in prangex
                tmp.dudx[x,y] += H.wgt[pt]*H.ddf(x-$shiftux-H.x[pt],y-$shiftuy-H.y[pt])
            end
        end
        tmp.dvdy .= tmp.dudx
        @inbounds for pt in 1:N
            prangex = max(1,ceil(Int,H.x[pt]+xmin)):min(NX-$dunx,floor(Int,H.x[pt]+xmax))
            prangey = max(1,ceil(Int,H.y[pt]+ymin)):min(NY-$duny,floor(Int,H.y[pt]+ymax))
            @inbounds for y in prangey, x in prangex
                w = tmp.dudx[x,y]
                w = w > eps() ? source.dudx[x,y]/w : 0.0
                target.dudx[pt]  += w*H.ddf(x-$shiftux-H.x[pt],y-$shiftuy-H.y[pt])
                w = tmp.dvdy[x,y]
                w = w > eps() ? source.dvdy[x,y]/w : 0.0
                target.dvdy[pt]  += w*H.ddf(x-$shiftux-H.x[pt],y-$shiftuy-H.y[pt])
            end
        end
        fill!(target.dudy,0.0)
        fill!(target.dvdx,0.0)
        xmin, xmax = -radius+$shiftvx, radius+$shiftvx
        ymin, ymax = -radius+$shiftvy, radius+$shiftvy
        @inbounds for pt in 1:N
            prangex = max(1,ceil(Int,H.x[pt]+xmin)):min(NX-$dvnx,floor(Int,H.x[pt]+xmax))
            prangey = max(1,ceil(Int,H.y[pt]+ymin)):min(NY-$dvny,floor(Int,H.y[pt]+ymax))
            for y in prangey, x in prangex
                tmp.dudy[x,y] += H.wgt[pt]*H.ddf(x-$shiftvx-H.x[pt],y-$shiftvy-H.y[pt])
            end
        end
        tmp.dvdx .= tmp.dudy
        @inbounds for pt in 1:N
            prangex = max(1,ceil(Int,H.x[pt]+xmin)):min(NX-$dvnx,floor(Int,H.x[pt]+xmax))
            prangey = max(1,ceil(Int,H.y[pt]+ymin)):min(NY-$dvny,floor(Int,H.y[pt]+ymax))
            @inbounds for y in prangey, x in prangex
                w = tmp.dudy[x,y]
                w = w > eps() ? source.dudy[x,y]/w : 0.0
                target.dudy[pt] += w*H.ddf(x-$shiftvx-H.x[pt],y-$shiftvy-H.y[pt])
                w = tmp.dvdx[x,y]
                w = w > eps() ? source.dvdx[x,y]/w : 0.0
                target.dvdx[pt] += w*H.ddf(x-$shiftvx-H.x[pt],y-$shiftvy-H.y[pt])
            end
        end
        target
  end

  # Construct regularization matrix
  @eval function RegularizationMatrix(H::Regularize{N,F},src::$pointtype{N,S,DT},target::$gridtype{$(ctype...),NX,NY,T,DDT}) where {N,F,NX,NY,S,T,DT,DDT}

    # note that we only need to compute two distinct matrices, since there are
    # only two types of cell data in this tensor

    #Hmat = (spzeros(length(target.u),length(src.u)),spzeros(length(target.v),length(src.v)))
    lenu = length(target.dudx)
    lenv = length(target.dudy)
    Hmat = spzeros(2lenu+2lenv,4N)
    g = deepcopy(src)
    v = deepcopy(target)
    g.dudx .= g.dudy .= g.dvdx .= g.dvdy .= zeros(Float64,N)
    for i = 1:N
      g.dudx[i] = 1.0  # these two are sufficient
      g.dudy[i] = 1.0
      H(v,g)
      Hmat[1:lenu,i]           = sparsevec(v.dudx)
      Hmat[lenu+1:lenu+lenv,i+N] = sparsevec(v.dudy)
      Hmat[lenu+lenv+1:lenu+2lenv,i+2N] = sparsevec(v.dudy)
      Hmat[lenu+2lenv+1:2lenu+2lenv,i+3N] = sparsevec(v.dudx)
      g.dudx[i] = 0.0
      g.dudy[i] = 0.0
    end
    if H._issymmetric
      # In symmetric case, these matrices are identical. (Interpolation is stored
      # as its transpose.)
      return RegularizationMatrix{typeof(target),typeof(src)}(Hmat),InterpolationMatrix{typeof(target),typeof(src)}(Hmat)
    else
      return RegularizationMatrix{typeof(target),typeof(src)}(Hmat)
    end
  end

  # Construct interpolation matrix
  @eval function InterpolationMatrix(H::Regularize{N,false},src::$gridtype{$(ctype...),NX,NY,T,DDT},target::$pointtype{N,S,DT}) where {N,NX,NY,S,T,DT,DDT}

    # note that we store interpolation matrices in the same shape as regularization matrices
    #Emat = (spzeros(length(src.u),length(target.u)),spzeros(length(src.v),length(target.v)))
    lenu = length(src.dudx)
    lenv = length(src.dudy)
    Emat = spzeros(2lenu+2lenv,4N)
    g = deepcopy(target)
    v = deepcopy(src)
    g.dudx .= g.dudy .= g.dvdx .= g.dvdy .= zeros(Float64,N)
    for i = 1:N
      g.dudx[i] = 1.0/H.wgt[i]  # unscale for interpolation
      g.dudy[i] = 1.0/H.wgt[i]  # unscale for interpolation
      H(v,g)
      Emat[1:lenu,i]           = sparsevec(v.dudx)
      Emat[lenu+1:lenu+lenv,i+N] = sparsevec(v.dudy)
      Emat[lenu+lenv+1:lenu+2lenv,i+2N] = sparsevec(v.dudy)
      Emat[lenu+2lenv+1:2lenu+2lenv,i+3N] = sparsevec(v.dudx)
      g.dudx[i] = 0.0
      g.dudy[i] = 0.0
    end
    InterpolationMatrix{typeof(src),typeof(target)}(Emat)
  end

  # Construct interpolation matrix with filtering
  @eval function InterpolationMatrix(H::Regularize{N,true},src::$gridtype{$(ctype...),NX,NY,T,DDT},target::$pointtype{N,S,DT}) where {N,NX,NY,S,T,DT,DDT}

    # note that we store interpolation matrices in the same shape as regularization matrices
    #Emat = (spzeros(length(src.u),length(target.u)),spzeros(length(src.v),length(target.v)))
    lenu = length(src.dudx)
    lenv = length(src.dudy)
    Emat = spzeros(2lenu+2lenv,4N)
    g = deepcopy(target)
    v = deepcopy(src)
    fill!(g,1.0)
    H(v,g)
    wtu = sparsevec(v.dudx)
    wtu.nzval .= 1 ./ wtu.nzval
    wtv = sparsevec(v.dudy)
    wtv.nzval .= 1 ./ wtv.nzval
    fill!(g,0.0)
    for i = 1:N
      g.dudx[i] = 1.0/H.wgt[i]  # unscale for interpolation
      g.dudy[i] = 1.0/H.wgt[i]  # unscale for interpolation
      H(v,g)
      Emat[1:lenu,i]             = wtu.*sparsevec(v.dudx)
      Emat[lenu+1:lenu+lenv,i+N] = wtv.*sparsevec(v.dudy)
      Emat[lenu+lenv+1:lenu+2lenv,i+2N] = wtv.*sparsevec(v.dudy)
      Emat[lenu+2lenv+1:2lenu+2lenv,i+3N] = wtu.*sparsevec(v.dudx)
      g.dudx[i] = 0.0
      g.dudy[i] = 0.0
    end
    InterpolationMatrix{typeof(src),typeof(target)}(Emat)
  end

end

###### Matrix multiplication extended to these Regularization/Interpolation matrices ########

# Regularize
function mul!(u::C,Hmat::RegularizationMatrix{C,F},f::F) where {C<:GridData,F<:PointData}
  fill!(u,0.0)
  nzv = Hmat.M.nzval
  rv = Hmat.M.rowval
  @inbounds for col = 1:Hmat.M.n
    fj = f[col]
    for j = Hmat.M.colptr[col]:(Hmat.M.colptr[col + 1] - 1)
        u[rv[j]] += nzv[j]*fj
    end
  end
  u
end

# Interpolate
function mul!(f::F,Emat::InterpolationMatrix{C,F},u::C) where {C<:GridData,F<:PointData}
  fill!(f,0.0)
  nzv = Emat.M.nzval
  rv = Emat.M.rowval
  @inbounds for col = 1:Emat.M.n
      tmp = zero(eltype(f))
      for j = Emat.M.colptr[col]:(Emat.M.colptr[col + 1] - 1)
          tmp += transpose(nzv[j])*u[rv[j]]
      end
      f[col] += tmp
  end
  f
end

# Interpolation of regularization, used for developing the filtering matrix
function mul!(C::Array{Float64},Emat::InterpolationMatrix{G,F},
                              Hmat::RegularizationMatrix{G,F}) where {G<:GridData,F<:PointData}
  fill!(C,0.0)
  Enzv = Emat.M.nzval
  Erv = Emat.M.rowval
  @inbounds for row = 1:Emat.M.n, col = 1:Hmat.M.n
      tmp = zero(eltype(C))
      for j = Emat.M.colptr[row]:(Emat.M.colptr[row + 1] - 1)
          tmp += transpose(Enzv[j])*Hmat[Erv[j],col]
      end
      C[row,col] += tmp
  end
  return C
end

(*)(Hmat::RegularizationMatrix{TU,TF},src::TF) where {TU,TF<:PointData} =
        mul!(TU(),Hmat,src)

(*)(Emat::InterpolationMatrix{TU,TF},src::TU) where {TU<:Union{Nodes,Edges,EdgeGradient,NodePair},TF<:PointData} =
                mul!(TF(),Emat,src)

(*)(Emat::InterpolationMatrix,Hmat::RegularizationMatrix) =
        mul!(Array{eltype(Emat),2}(undef,Emat.M.n,Hmat.M.n),Emat,Hmat)

function Base.summary(io::IO, H::RegularizationMatrix{TU,TF}) where {TU,TF}
    print(io, "Regularization matrix acting on type $TF and returning type $TU")
end

function Base.summary(io::IO, H::InterpolationMatrix{TU,TF}) where {TU,TF}
    print(io, "Interpolation matrix acting on type $TU and returning type $TF")
end


text:
using NLIDatasets, Test

@testset "NLIDatasets.jl" begin

    @testset "SNLI" begin
        using NLIDatasets: SNLI
        @test isfile(SNLI.train_tsv())
        @test isfile(SNLI.train_jsonl())
        @test isfile(SNLI.dev_tsv())
        @test isfile(SNLI.dev_jsonl())
        @test isfile(SNLI.test_tsv())
        @test isfile(SNLI.test_jsonl())
    end

    @testset "MultiNLI" begin
        using NLIDatasets: MultiNLI
        @test isfile(MultiNLI.train_tsv())
        @test isfile(MultiNLI.train_jsonl())
        @test isfile(MultiNLI.dev_matched_tsv())
        @test isfile(MultiNLI.dev_matched_jsonl())
        @test isfile(MultiNLI.dev_mismatched_tsv())
        @test isfile(MultiNLI.dev_mismatched_jsonl())
    end

    @testset "XNLI" begin
        using NLIDatasets: XNLI
        @test isfile(XNLI.dev_tsv())
        @test isfile(XNLI.dev_jsonl())
        @test isfile(XNLI.test_tsv())
        @test isfile(XNLI.test_jsonl())
    end

    @testset "SciTail" begin
        using NLIDatasets: SciTail
        @test isfile(SciTail.train_tsv())
        @test isfile(SciTail.train_jsonl())
        @test isfile(SciTail.dev_tsv())
        @test isfile(SciTail.dev_jsonl())
        @test isfile(SciTail.test_tsv())
        @test isfile(SciTail.test_jsonl())
    end

    @testset "HANS" begin
        using NLIDatasets: HANS
        @test isfile(HANS.test_tsv())
    end

    @testset "Breaking_NLI" begin
        using NLIDatasets: BreakingNLI
        @test isfile(BreakingNLI.test_jsonl())
    end
end

ext: jl


text:
# dtft.jl

import MIRT: dtft_init, dtft, dtft_adj
import MIRT: dtft_loop_n
import MIRT: dtft_loop_m
import MIRT: dtft_dist_m
import MIRT: dtft_pmap_m
import MIRT: dtft_matvec

using FFTW: fft
#using LinearAlgebra: norm
using Test: @test, @testset
using LinearMapsAA: LinearMapAA, LinearMapAO

# 1D small test data to verify correctness
function dtft_test1( ; N::Int=2^10)
	x = randn(ComplexF64, N)
	w = (2*pi) * (0:N-1) / N

	o0 = fft(x)

	o1 = dtft_loop_n(w, x)
	o2 = dtft_loop_m(w, x)
	o3 = dtft_dist_m(w, x)
	o4 = dtft_pmap_m(w, x)
	o5 = dtft_matvec(w, x)
#	@show norm(o1 - o0) / norm(o0)

	@test isapprox(o1, o0)
	@test isapprox(o2, o0)
	@test isapprox(o3, o0)
	@test isapprox(o4, o0)
	@test isapprox(o5, o0)

	true
end

# 2D small test
function dtft_test2( ; N::Dims = (2^3,2^2))
	x = randn(ComplexF64, N)
	w1 = (2*pi) * (0:N[1]-1) / N[1]
	w2 = (2*pi) * (0:N[2]-1) / N[2]
	w1 = repeat(w1, 1, N[2])
	w2 = repeat(w2', N[1], 1)
	w = [vec(w1) vec(w2)]
	d = dtft_init(w, N)

	o0 = fft(x)
#	o1 = dtft_loop_n(w, x)
	o1 = d.dtft(x)
	o1 = reshape(o1, N)
	@test isapprox(o1, o0)
	true
end

# 1D verify consistency
function dtft_test1c( ; N::Int=2^10, M::Int=2^11, n_shift::Real=7)
	x = randn(ComplexF64, N)
	w = randn(M)

	o1 = dtft_loop_n(w, x ; n_shift=n_shift)
	o2 = dtft_loop_m(w, x ; n_shift=n_shift)
	o3 = dtft_dist_m(w, x ; n_shift=n_shift)
	o4 = dtft_pmap_m(w, x ; n_shift=n_shift)
	o5 = dtft_matvec(w, x ; n_shift=n_shift)

	@test isapprox(o2, o1)
	@test isapprox(o3, o1)
	@test isapprox(o4, o1)
	@test isapprox(o5, o1)

	d = dtft_init(w, N ; n_shift=n_shift)
	o6 = d.dtft(x)
	@test isapprox(o6, o1)

	o7 = d.A * x
	@test isapprox(o7, o1)

	b1 = dtft_adj(w, o1, N ; n_shift=n_shift)
	b2 = d.adjoint(o1)
	@test isapprox(b2, b1)
	b3 = d.A' * o1
	@test isapprox(b3, b1)

	A = d.A
	@test A.name == "dtft1"
	@test A.N == (N,)

	true
end

# 2D verify consistency
function dtft_test2c( ;
	M::Int=2^9,
	N::Dims=(2^6,2^4),
	n_shift::AbstractVector{<:Real} = zeros(Int, length(N)),
)

	x = randn(ComplexF64, N)
	w = (rand(M,2) .- 0.5) * 2 * pi

	o1 = dtft_loop_n(w, x ; n_shift=n_shift)
	sd = dtft_init(w, N ; n_shift=n_shift)
	o2 = sd.dtft(x)
	@test isequal(o2, o1)
	o3 = sd.A * x
	@test isequal(o3, o1)

	b1 = dtft_adj(w, o1, N ; n_shift=n_shift)
	b2 = sd.adjoint(o1)
	@test isequal(b2, b1)

	b3 = sd.A' * o1
	@test isequal(b3, b1)

	A = sd.A
	@test A.name == "dtft2"
	@test A.N == N
	@test A isa LinearMapAO
	true
end

"""
dtft_test1_adj()
test adjoint
"""
function dtft_test1_adj( ; N::Int=20, M::Int=30, n_shift::Int=5)
	w = randn(M)
	_, _, A = dtft_init(w, N ; n_shift = n_shift)
	@test Matrix(A)' ≈ Matrix(A')
	true
end

"""
dtft_test2_adj()
test adjoint
"""
function dtft_test2_adj( ;
	M::Int=2^4,
	N::Dims=(2^3,2^2),
	n_shift::AbstractVector{<:Real} = [6,1.7],
)
	w = (rand(M,2) .- 0.5) * 2 * pi
	sd = dtft_init(w, N ; n_shift=n_shift)
	A = sd.A
	@test A isa LinearMapAO
	@test Matrix(A)' ≈ Matrix(A')
	true
end

@testset "1d" begin
	@test dtft_test1()
	@test dtft_test1c()
	@test dtft_test1_adj()
end

@testset "2d" begin
	@test dtft_test2()
	@test dtft_test2c()
	@test dtft_test2_adj()
end


text:
using Pkg
using SafeTestsets
const LONGER_TESTS = false

const GROUP = get(ENV, "GROUP", "All")
const is_APPVEYOR = Sys.iswindows() && haskey(ENV,"APPVEYOR")

function activate_downstream_env()
    Pkg.activate("downstream")
    Pkg.develop(PackageSpec(path=dirname(@__DIR__)))
    Pkg.instantiate()
end

if GROUP == "All" || GROUP == "Core"
  @time @safetestset "Basic Tests" begin include("basictests.jl") end
end

if !is_APPVEYOR && GROUP == "GPU"
  activate_downstream_env()
  @time @safetestset "CUDA" begin include("cuda.jl") end
end


text:
#Ost - "E6" - восток, Nord - "N8" - север, West - "W4" - запад, Sud - "S2" - юг
function main_function11!(::Robot)
    A=[]
    B=[]
    help_main_2(r,Ost)
    help_main_1(r, Nord)
    help_main_1(r, West)
    help_main_2(r, Sud)
end

function obhod(::Robot,side1::HorizonSide,side2::HorizonSide,A,B)
    while isborder(r,side1)
        move!(r,side2),push!(A,side2)
    end
    while length(B)>0
        side=pop!(B)
        move!(r,inverse(side)),push!(A,inverse(side))
    end
        putmarker!(r)
end
function go_back_position(::Robot,A)
    while length(A)>0
        side=pop!(A)
        move!(r,inverse(side))
    end
end

function help_main_1(r, side)
    A = []
    B = []
    if !isborder(r,side)
        move!(r,side)
        putmarker!(r)
        move!(r,inverse(side))
    else
        if isborder(r,hourly(side))
            side3=inverse(hourly(side))
        else
            side3=hourly(side)
        end
        while isborder(r,side) && !isborder(r,side3)
            move!(r,side3),push!(A,side3),push!(B,side3)
        end
        if !isborder(r,side3) && !isborder(r,side)
            move!(r,side),push!(A,side)
            obhod(r,inverse(side3),side,A,B)
            putmarker!(r)
            go_back_position(r,A)
        else
            go_back_position(r,A)
            B=[]
        end
    end
end

function help_main_2(r, side)
    A = []
    B = []
    if !isborder(r,side)
        move!(r,side)
        putmarker!(r)
        move!(r,inverse(side))
    else
        if isborder(r,counterclockwise(side))
            side4=inverse(counterclockwise(side))
        else
            side4=counterclockwise(side)
        end
        while isborder(r,side) && !isborder(r,side4)
            move!(r,side4),push!(A,side4),push!(B,side4)
        end
        if !isborder(r,side4) && !isborder(r,side)
            move!(r,side),push!(A,side)
            obhod(r,inverse(side4),side,A,B)
            putmarker!(r)
            go_back_position(r,A)
        else
            go_back_position(r,A)
            B=[]
        end
    end
end

counterclockwise(side) = HorizonSide(mod(Int(side)+1, 4))
hourly(side) = HorizonSide(mod(Int(side)- 1 , 4))
inverse(side::HorizonSide) = HorizonSide(mod(Int(side)+2,4))


text:
@testset "1486.xor-operation-in-an-array.jl" begin
    @test xor_operation(5, 0) == 8
    @test xor_operation(4, 3) == 8
    @test xor_operation(1, 7) == 7
    @test xor_operation(10, 5) == 2
end


text:
module Data

using Random
using Photon: makeArray

include("dataset.jl")
include("transformer.jl")

end


text:
##########
# TArray #
##########

"""
    TArray{T}(dims, ...)

Implementation of data structures that automatically perform deepcopy during task copying.

More specifically, we store the underlying arrays in the field
`TArray.data`, and this filed will be deepcopied while we are copying
a task who manipulates it.

Usage:

```julia
TArray(dim)
```

Example:

```julia
ta = TArray(4)              # init
for i in 1:4 ta[i] = i end  # assign
Array(ta)                   # convert to 4-element Array{Int64,1}: [1, 2, 3, 4]
```
"""
mutable struct TArray{T, N, A <: AbstractArray{T, N}} <: AbstractArray{T, N}
    data::A
end

TArray{T}(d::Integer...) where T = TArray(T, d)
TArray{T}(::UndefInitializer, d::Integer...) where T = TArray(T, d)
TArray{T}(::UndefInitializer, dim::NTuple{N,Int}) where {T,N} = TArray(T, dim)
TArray{T,N}(d::Vararg{<:Integer,N}) where {T,N} = TArray(T, d)
TArray{T,N}(::UndefInitializer, d::Vararg{<:Integer,N}) where {T,N} = TArray{T,N}(d)
TArray{T,N}(dim::NTuple{N,Int}) where {T,N} = TArray(T, dim)
TArray(T::Type, dim) = TArray(Array{T}(undef, dim))

localize(x) = x
localize(x::AbstractArray) = TArray(x)
getdata(x) = x
getdata(x::TArray) = x.data
tape_copy(x::TArray) = TArray(deepcopy(x.data))

# Constructors
"""
     tzeros(dims, ...)

Construct a TArray of zeros.
Trailing arguments are the same as those accepted by `TArray`.

```julia
tzeros(dim)
```

Example:

```julia
tz = tzeros(4)              # construct
Array(tz)                   # convert to 4-element Array{Int64,1}: [0, 0, 0, 0]
```
"""
function tzeros(T::Type, dim)
    d = zeros(T, dim)
    TArray(d)
end

tzeros(::Type{T}, d1::Integer, drest::Integer...) where T = tzeros(T, convert(Dims, tuple(d1, drest...)))
tzeros(d1::Integer, drest::Integer...) = tzeros(Float64, convert(Dims, tuple(d1, drest...)))
tzeros(d::Dims) = tzeros(Float64, d)

"""
     tfill(val, dim, ...)

Construct a TArray of a specified value.

```julia
tfill(val, dim)
```

Example:

```julia
tz = tfill(9.0, 4)            # construct
Array(tz)                     # convert to 4-element Array{Float64,1}:  [9.0  9.0  9.0  9.0]
```
"""
function tfill(val::Real, dim)
    d = fill(val, dim)
    TArray(d)
end

#
# Conversion between TArray and Array
#

function Base.convert(::Type{Array}, x::TArray)
    convert(Array{eltype(x), ndims(x)}, x)
end
function Base.convert(::Type{Array{T, N}}, x::TArray{T, N}) where {T, N}
    convert(Array{T, N}, getdata(x))
end

function Base.convert(::Type{TArray}, x::AbstractArray)
    convert(TArray{eltype(x), ndims(x)}, x)
end
function Base.convert(::Type{TArray{T, N}}, x::AbstractArray{T, N}) where {T, N}
    TArray(x)
end

#
# Representation
#
function Base.show(io::IO, ::MIME"text/plain", x::TArray)
    show(io,  MIME("text/plain"), getdata(x))
end

Base.show(io::IO, x::TArray) = show(io, getdata(x))

function Base.summary(io::IO, x::TArray)
    print(io, "TArray: ")
    summary(io, getdata(x))
end

#
# Forward many methods to the underlying array
#
for F in (:size,
          :iterate,
          :firstindex, :lastindex, :axes)
    @eval Base.$F(a::TArray, args...) = $F(getdata(a), args...)
end

#
# Similarity implementation
#

Base.similar(x::TArray, ::Type{T}, dims::Dims) where T = TArray(similar(getdata(x), T, dims))

for op in [:(==), :≈]
    @eval Base.$op(x::TArray, y::AbstractArray) = Base.$op(getdata(x), y)
    @eval Base.$op(x::AbstractArray, y::TArray) = Base.$op(x, getdata(y))
    @eval Base.$op(x::TArray, y::TArray) = Base.$op(getdata(x), getdata(y))
end

#
# Array Stdlib
#

# Indexing Interface
Base.@propagate_inbounds function Base.getindex(x::TArray{T, N}, I::Vararg{Int,N}) where {T, N}
    return getdata(x)[I...]
end

Base.@propagate_inbounds function Base.setindex!(x::TArray{T, N}, e, I::Vararg{Int,N}) where {T, N}
    getdata(x)[I...] = e
end

function Base.push!(x::TArray, e)
    push!(getdata(x), e)
end

function Base.pop!(x::TArray)
    pop!(getdata(x))
end

# Other methods from stdlib

Base.view(x::TArray, inds...; kwargs...) =
    Base.view(getdata(x), inds...; kwargs...) |> localize
Base.:-(x::TArray) = (-getdata(x)) |> localize
Base.transpose(x::TArray) = transpose(getdata(x)) |> localize
Base.adjoint(x::TArray) = adjoint(getdata(x)) |> localize
Base.repeat(x::TArray; kw...) = repeat(getdata(x); kw...) |> localize

Base.hcat(xs::Union{TArray{T,1}, TArray{T,2}}...) where T =
    hcat(getdata.(xs)...) |> localize
Base.vcat(xs::Union{TArray{T,1}, TArray{T,2}}...) where T =
    vcat(getdata.(xs)...) |> localize
Base.cat(xs::Union{TArray{T,1}, TArray{T,2}}...; dims) where T =
    cat(getdata.(xs)...; dims = dims) |> localize

Base.reshape(x::TArray, dims::Union{Colon,Int}...) = reshape(getdata(x), dims) |> localize
Base.reshape(x::TArray, dims::Tuple{Vararg{Union{Int,Colon}}}) =
    reshape(getdata(x), Base._reshape_uncolon(getdata(x), dims)) |> localize
Base.reshape(x::TArray, dims::Tuple{Vararg{Int}}) = reshape(getdata(x), dims) |> localize

Base.permutedims(x::TArray, perm) = permutedims(getdata(x), perm) |> localize
Base.PermutedDimsArray(x::TArray, perm) = PermutedDimsArray(getdata(x), perm) |> localize
Base.reverse(x::TArray; dims) = reverse(getdata(x), dims = dims) |> localize

Base.sum(x::TArray; dims = :) = sum(getdata(x), dims = dims) |> localize
Base.sum(f::Union{Function,Type},x::TArray) = sum(f.(getdata(x))) |> localize
Base.prod(x::TArray; dims=:) = prod(getdata(x); dims=dims) |> localize
Base.prod(f::Union{Function, Type}, x::TArray) = prod(f.(getdata(x))) |> localize

Base.findfirst(x::TArray, args...) = findfirst(getdata(x), args...) |> localize
Base.maximum(x::TArray; dims = :) = maximum(getdata(x), dims = dims) |> localize
Base.minimum(x::TArray; dims = :) = minimum(getdata(x), dims = dims) |> localize

Base.:/(x::TArray, y::TArray) = getdata(x) / getdata(y) |> localize
Base.:/(x::AbstractArray, y::TArray) = x / getdata(y) |> localize
Base.:/(x::TArray, y::AbstractArray) = getdata(x) / y |> localize
Base.:\(x::TArray, y::TArray) = getdata(x) \ getdata(y) |> localize
Base.:\(x::AbstractArray, y::TArray) = x \ getdata(y) |> localize
Base.:\(x::TArray, y::AbstractArray) = getdata(x) \ y |> localize
Base.:*(x::TArray, y::TArray) = getdata(x) * getdata(y) |> localize
Base.:*(x::AbstractArray, y::TArray) = x * getdata(y) |> localize
Base.:*(x::TArray, y::AbstractArray) = getdata(x) * y |> localize

# broadcast
Base.BroadcastStyle(::Type{<:TArray}) = Broadcast.ArrayStyle{TArray}()
Broadcast.broadcasted(::Broadcast.ArrayStyle{TArray}, f, args...) = f.(getdata.(args)...) |> localize

import LinearAlgebra
import LinearAlgebra:  \, /, inv, det, logdet, logabsdet, norm

LinearAlgebra.inv(x::TArray) = inv(getdata(x)) |> localize
LinearAlgebra.det(x::TArray) = det(getdata(x)) |> localize
LinearAlgebra.logdet(x::TArray) = logdet(getdata(x)) |> localize
LinearAlgebra.logabsdet(x::TArray) = logabsdet(getdata(x)) |> localize
LinearAlgebra.norm(x::TArray, p::Real = 2) =
    LinearAlgebra.norm(getdata(x), p) |> localize

import LinearAlgebra: dot
dot(x::TArray, ys::TArray) = dot(getdata(x), getdata(ys)) |> localize
dot(x::AbstractArray, ys::TArray) = dot(x, getdata(ys)) |> localize
dot(x::TArray, ys::AbstractArray) = dot(getdata(x), ys) |> localize

using Statistics
Statistics.mean(x::TArray; dims = :) = mean(getdata(x), dims = dims) |> localize
Statistics.std(x::TArray; kw...) = std(getdata(x), kw...) |> localize

# TODO
# * NNlib

max_stars_repo_name: TuringLang/Libtask

max_issues_repo_name: TuringLang/Libtask

max_forks_repo_name: TuringLang/Libtask


text:
### A Pluto.jl notebook ###
# v0.12.16

using Markdown
using InteractiveUtils

# ╔═╡ e077ab16-22e5-11eb-15ed-5f0654bf7328
using Pkg, DrWatson

# ╔═╡ e75a9a06-22e5-11eb-2e36-8d4f62b830ed
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using Statistics, Distributions, DataFrames, CSV
	using HypothesisTests
end;

# ╔═╡ ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
md"## Listing6.10"

# ╔═╡ fa68607e-22e5-11eb-0558-c9a4d9f77426
begin
	data = CSV.read(datadir("machine1.csv"), DataFrame; delim=",")[:, 1]
end;

# ╔═╡ 1aa41898-3afa-11eb-2f49-570853ce6435
begin
	n, s, alpha = length(data), std(data), 0.1
	ci = (  (n-1)*s^2/quantile(Chisq(n-1),1-alpha/2),
			(n-1)*s^2/quantile(Chisq(n-1),alpha/2)  )

	Text("Point estimate for the variance: $(s^2)")
end

# ╔═╡ d0a56572-3b2a-11eb-383a-f94fe2f45a5b
Text("Confidence interval for the variance: $ci")

# ╔═╡ 475ff888-22e6-11eb-2354-09f8f40a8e12
md"## End of listing6.10"

# ╔═╡ Cell order:
# ╟─ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
# ╠═e077ab16-22e5-11eb-15ed-5f0654bf7328
# ╠═e75a9a06-22e5-11eb-2e36-8d4f62b830ed
# ╠═fa68607e-22e5-11eb-0558-c9a4d9f77426
# ╠═1aa41898-3afa-11eb-2f49-570853ce6435
# ╠═d0a56572-3b2a-11eb-383a-f94fe2f45a5b
# ╟─475ff888-22e6-11eb-2354-09f8f40a8e12


text:
function cpairsum(a::TwoFloat{T}, b::TwoFloat{T}) where {T} =
    TwoFloat{T}(cpairsum(a[1],a[2],b[1],b[2]))

function cpairdiff(a::TwoFloat{T}, b::TwoFloat{T}) where {T} =
    TwoFloat{T}(cpairdiff(a[1],a[2],b[1],b[2]))

function cpairprod(a::TwoFloat{T}, b::TwoFloat{T}) where {T} =
    TwoFloat{T}(cpairprod(a[1],a[2],b[1],b[2]))

function cpairdiv(a::TwoFloat{T}, b::TwoFloat{T}) where {T} =
    TwoFloat{T}(cpairdiv(a[1],a[2],b[1],b[2]))

function cpairinv(a::TwoFloat{T}) where {T} =
    TwoFloat{T}(cpairinv(a[1],a[2]))

function cpairinv(a::TwoFloat{T}) where {T} =
    TwoFloat{T}(cpairinv(a[1],a[2]))

@inline function cpairsum(ahi::T, alo::T, bhi::T, blo::T) where {T}
    hihi, hilo = two_sum_ordered(ahi, bhi)
    t  = ahi + bhi - hihi
    lo = (t + hilo) + (alo + blo)
    return hi, lo
end

@inline function cpairsum(ahi::T, alo::T, bhi::T, blo::T) where {T}
    hi = ahi + bhi
    t  = ahi + bhi - hi
    lo = t + (alo + blo)
    return hi, lo
end

@inline function cpairdiff(ahi::T, alo::T, bhi::T, blo::T) where {T}
    hi = ahi - bhi
    t  = ahi - bhi - hi
    lo = t + (alo - blo)
    return hi, lo
end

function cpairprod(ahi::T, alo::T, bhi::T, blo::T) where {T}
    hi = a * b
    t  = fma(ahi, bhi, -hi)
    lo = t + (ahi*blo + bhi*alo)
    return hi, lo
end

function cpairinv(bhi::T, blo::T) where {T}
    hi = inv(bhi)
    t  = fma(-hi, bhi, one(T))
    lo = fma(-hi, blo, t) / (bhi + blo)
    return hi, lo
end

function cpairdiv(ahi::T, alo::T, bhi::T, blo::T) where {T}
    hi = ahi / bhi
    t  = fma(-hi, bhi, ahi) + alo
    lo = fma(-hi, blo, t) / (bhi + blo)
    return hi, lo
end

function cpairsqrt(ahi::T, alo::T) where {T}
    hi = sqrt(ahi)
    t  = -fma(hi, hi, -ahi)
    lo = (t + alo) / (2hi)
    return hi, lo
end

Base.:(-)(a::TwoFloat{T}) where {T} = (-a[1], -a[2])
Base.signbit(a::TwoFloat{T}) where {T} = signbit(a[1])
Base.abs(a::TwoFloat{T}) where {T} = signbit(a[1]) ? (abs(a[1]), -a[2]) : a
Base.copysign(a::TwoFloat, b) = signbit(b) ? (signbit(a[1]) ? a : -a) : a
Base.flipsign(a::TwoFloat, b) = signbit(b) ? -a : a

function Base.:(+)(a::TwoFloat{T}, b::TwoFloat{T}) where {T}
    if abs(b[1]) <= abs(a[1])
         ahi, alo = a
         bhi, blo = b
    else
         ahi, alo = b
         bhi, blo = a
    end
    return Base.:(+)(ahi, alo, bhi, blo)
end

function Base.:(+)(ahi::T, alo::T, bhi::T, blo::T)  where {T<:FastFloat}
    t0, t1 = two_hilo_sum(ahi, alo)
    t2, t3 = two_hilo_sum(bhi, blo)
    hi, t4 = two_hilo_sum(t0, t2)
    lo = t4 + (t1 + t3)
    hi, lo = two_hilo_sum(hi, lo)

    return TwoFloat{T}((hi, lo))
end

function Base.:(+)(a::TwoFloat{T}, b::TwoFloat{T}) where {T}
    if abs(b[1]) <= abs(a[1])
         ahi, alo = a
         bhi, blo = b
    else
         ahi, alo = b
         bhi, blo = a
    end
    t0, t1 = two_hilo_sum(ahi, alo)
    t2, t3 = two_hilo_sum(bhi, blo)
    hi, t4 = two_hilo_sum(t0, t2)
    lo = t4 + (t1+t3)
    hi, lo = two_hilo_sum(hi, lo)
    return TwoFloat{T}((hi, lo))
end

function Base.:(-)(a::TwoFloat{T}, b::TwoFloat{T}) where {T}
    ahi, alo = a
    bhi, blo = b
    ahi, alo, bhi, blo = mxxmnn_abs(ahi, alo, -bhi, -blo)
    t0, t1 = two_hilo_sum(ahi, alo)
    t2, t3 = two_hilo_sum(bhi, blo)
    hi, t4 = two_hilo_sum(t0, t2)
    lo = t4 + (t1+t3)
    hi, lo = two_hilo_sum(hi, lo)
    return TwoFloat{T}((hi, lo))
end

function Base.:(*)(a::TwoFloat{T}, b::TwoFloat{T}) where {T}
    ahi, alo = a
    bhi, blo = b
    if abs(alo) > abs(bhi)
        alo, bhi = bhi, alo
        alo, blo = maxmin_abs(alo,blo)
    end
    hi = ahi * bhi
    t = fma(ahi, bhi, -hi)
    lo = t + (ahi*blo + bhi*alo)
    return TwoFloat{T}((hi, lo))
end

function Base.:(/)(a::TwoFloat{T}, b::TwoFloat{T}) where {T}
    ahi, alo = a
    bhi, blo = b
    hi = ahi / bhi
    t = fma(-bhi, hi, ahi)
    lo = ((t + alo) - hi*blo) / (bhi + blo)
    return TwoFloat{T}((hi, lo))
end

function Base.inv(b::TwoFloat{T}) where {T}
    bhi, blo = b
    hi = inv(bhi)
    t = fma(-bhi, hi, one(T))
    lo = -(hi * blo) / (bhi + blo)
    return TwoFloat{T}((hi, lo))
end

function Base.sqrt(a::TwoFloat{T}) where {T}
    ahi, alo = a
    hi = sqrt(ahi)
    t = fma(hi, -hi, ahi)
    lo = (t + alo) / (2 * hi)
    return TwoFloat{T}((hi, lo))
end


text:
module DarkTracks

export main

using FFmpegPipe, OnlineStats, Statistics, IdentityRanges, Dates, Images, Unitful, JSON, DelimitedFiles#, Printf, , ImageDraw

function lastframe(videofile)
    str = read(`ffprobe -show_streams -of json -v quiet -i $videofile`, String)
    jsn = JSON.parse(str)
    stream = findfirst(x -> x["codec_type"] == "video", jsn["streams"])
    txt = jsn["streams"][stream]["nb_frames"]
    parse(Int, txt)
end

function getdimensions(videofile, start_frame, stop_frame)
    str = read(`ffprobe -show_streams -of json -v quiet -i $videofile`, String)
    jsn = JSON.parse(str)
    stream = findfirst(x -> x["codec_type"] == "video", jsn["streams"])
    txt = jsn["streams"][stream]["r_frame_rate"]
    m = match(r"(\d*)/(\d*)", txt)
    f, s = parse.(Int, m.captures)
    fps₀ = f/(s*Unitful.s)
    w₀ = jsn["streams"][stream]["width"]
    h₀ = jsn["streams"][stream]["height"]
    start_second = ustrip(start_frame/fps₀)
    duration_second = ustrip((stop_frame - start_frame)/fps₀)
    w₀, h₀, fps₀, start_second, duration_second
end

function warpdimensions(w₀, h₀, fps₀, spatialscale, temporalscale)
    w, h = round.(Int, [w₀, h₀]./spatialscale)
    fps = round(Int, ustrip(fps₀)/temporalscale)
    w, h, fps
end

function detectnext(r, img, noise, ind)
    rows = IdentityRange((-r:r) .+ ind[1])
    cols = IdentityRange((-r:r) .+ ind[2])
    pimg = padarray(img, Fill(zero(img[1]), (r, r)))
    roi = @view pimg[rows, cols]
    M, ind = findmax(roi)
    M > noise[ind] ? ind : missing
end

function detectfirst(vi, noise, sz, n)
    ind₀ = round.(Int, sz./2)
    r₀ = min(sz...) ÷ 11
    for i in 1:n
        img = readframe(vi)
        ind = detectnext(r₀, img, noise, ind₀)
        if !ismissing(ind)
            return i, ind
        end
    end
    error("failed to detect the beetle's initial starting point!")
end

function track(r₀, vi, noise, n, sz)
    inds = Vector{Union{Missing, CartesianIndex{2}}}(missing, n)
    i1, ind = detectfirst(vi, noise, sz, n)
    inds[i1] = ind
    r = copy(r₀)
    for i in i1+1:n
        img = readframe(vi)
        inds[i] = detectnext(r, img, noise, inds[findlast(!ismissing, inds)])
        if ismissing(inds[i])
            r += r₀
        else
            r = copy(r₀)
        end
    end
    inds
end

function savedata(videofile, ts, inds)
    path, nameext = splitdir(videofile)
    name, _ = splitext(nameext)
    csvname = joinpath(path, "$name.csv")
    xs = ["x"; getindex.(inds, 1)]
    ys = ["y"; getindex.(inds, 2)]
    writedlm(csvname, zip(["seconds"; ts], xs, ys))
end

#=function drawit(videofile, n, inds, imgs)
    path, nameext = splitdir(videofile)
    name, _ = splitext(nameext)
    imgs2 = joinpath(path, name)
    mkpath(imgs2)
    rm.(joinpath.(imgs2, readdir(imgs2)))
    for (i, ind, img) in zip(1:n, inds, imgs)
        img2 = RGB.(img)
        if !ismissing(ind)
            draw!(img2, CirclePointRadius(ind, 2), RGB(1,0,0))
        end
        save(joinpath(imgs2, @sprintf "%06i.png" i), img2)
    end
    run(`ffmpeg -y -framerate 25 -i $(joinpath(imgs2, "%06d.png")) $(joinpath(imgs2, "0.mp4"))`)
end=#

function getnoise(sz, videofile, start_second, duration_second)
    h, w = sz
    o = [Mean() for i in 1:h, j in 1:w]
    nframes = 9
    vi = openvideo(videofile, "r", ss_in = start_second, t = duration_second, r = nframes/duration_second, s = string(w, "x", h), pix_fmt = "gray")
    for i in 1:nframes
        img = readframe(vi)
        for j in eachindex(img)
            fit!(o[j], Float64(img[j]))
        end
    end
    close(vi)
    2value.(o)
end

getradius(sscale, tscale) = 10exp(-1.020819402836599log(sscale) + 0.7668247162044156log(tscale) + 0.8641290371601322)

"""
    seconds, coordinates = main(videofile; start_frame = 1, stop_frame = lastframe(videofile),  spatialscale = 10, temporalscale = 10)

Track a bright spot in a dim video file, `videofile`. Optionally give a frame number to start at, `start_frame`, and stop, `stop_frame`. These default to the beginning and ending of the video file. `spatialscale` and `temporalscale` help subsample the video for faster processing. When these equal 1 then there is no sub-sampling and the video is processed in original spatial and temporal resolutions. For example, for scaling the frame by half its original size and sampling every second frame use `spatialscale = 2, temporalscale = 2`. A default value of 10 for both works well.

# Return variables
`seconds` is a `Vector{Float64}` of times in second for each of the `coordinates` (itself a `Vector{Tuple{Int, Int}}`).

# Examples

```jldoctest; setup = :(using DarkTracks; videofile = abspath(joinpath(dirname(pathof(DarkTracks)), "..", "test", "testvideo.mp4")))
julia> main(videofile, stop_frame = 100, spatialscale = 10)
([1.32, 1.54, 1.76, 1.98], Tuple{Int64,Int64}[(510, 730), (520, 730), (520, 730), (520, 730)])

```

"""
function main(videofile; start_frame = 1, stop_frame = lastframe(videofile),  spatialscale = 10, temporalscale = 10)
    w₀, h₀, fps₀, start_second, duration_second = getdimensions(videofile, start_frame, stop_frame)
    w, h, fps = warpdimensions(w₀, h₀, fps₀, spatialscale, temporalscale)
    sz = (h, w)
    n = round(Int, duration_second*fps)
    ts = range(0, stop = duration_second, length = n)
    r₀ = round(Int, getradius(spatialscale, temporalscale))
    noise = getnoise(sz, videofile, start_second, duration_second)
    vi = openvideo(videofile, "r", ss_in = start_second, t = duration_second, r = fps, s = string(w, "x", h), pix_fmt = "gray")
    inds = track(r₀, vi, noise, n, sz)
    close(vi)

    tss = [t for (t, ind) in zip(ts, inds) if !ismissing(ind)]
    indss = [spatialscale.*Tuple(ind) for ind in inds if !ismissing(ind)]
    savedata(videofile, tss, indss)
    tss, indss
    # drawit(videofile, n, inds, imgs)
end

end # module


text:
"""
    momentum(mechanism, storage)

    mechanism's linear and angular momentum

    mechanism: Mechanism
    storage: Storage
"""
function momentum(mechanism::Mechanism, storage::Storage{T,N}) where {T,N}
    m = [szeros(T,6) for i = 1:N]
    for i = 1:N
        m[i] = momentum(mechanism, storage, i)
    end
    return m # in world frame
end

function momentum(mechanism::Mechanism{T}, body::Body{T}) where T
    timestep= mechanism.timestep
    state = body.state
    mass = body.mass
    inertia = body.inertia
    # x1, q1 = previous_configuration(state)
    x2, q2 = current_configuration(state)
    x3, q3 = next_configuration(state, timestep)

    v15 = body.state.vsol[2] # v1.5
    ω15 = body.state.ϕsol[2] # ω1.5

    D2x = 1 / timestep * mass * (x3 - x2) - 0.5 * timestep * mass * mechanism.gravity
    D2q = -2.0 / timestep * LVᵀmat(q2)' * Tmat() * Rmat(q3)' * Vᵀmat() * inertia * Vmat() * Lmat(q2)' * vector(q3)
    p_linear_body = D2x - 0.5 * state.F2
    p_angular_body = D2q - 0.5 * state.τ2

    for joint in mechanism.joints
        if body.id ∈ [joint.parent_id, joint.child_id]
            f_joint = impulse_map(mechanism, joint, body) * joint.impulses[2] # computed at t = 1.5

            joint.spring && (f_joint += spring_impulses(mechanism, joint, body)) # computed at t = 1.5
            joint.damper && (f_joint += damper_impulses(mechanism, joint, body)) # computed at t = 1.5

            p_linear_body -= 0.5 * f_joint[1:3]
            p_angular_body -= 0.5 * f_joint[4:6]
        end
    end

    return [p_linear_body; rotation_matrix(q2) * p_angular_body]
end

function momentum(mechanism::Mechanism{T,Nn,Ne,Nb,Ni}, storage::Storage{T,Ns}, t::Int) where {T,Nn,Ne,Nb,Ni,Ns}
    p = zeros(T, 6)
    com = center_of_mass(mechanism, storage, t)
    mass = total_mass(mechanism)
    p_linear_body = [storage.px[i][t] for i = 1:Nb] # in world frame
    p_angular_body = [storage.pq[i][t] for i = 1:Nb] # in world frame

    p_linear = sum(p_linear_body)
    p_angular = zeros(T, 3)
    v_com = p_linear ./ mass
    for (i, body) in enumerate(mechanism.bodies)
        r = storage.x[i][t] - com
        v_body = p_linear_body[i] ./ body.mass
        p_angular += p_angular_body[i]
        p_angular += cross(r, body.mass * (v_body - v_com)) #TODO maybe there is cleaner way to handle the factor 2
    end

    return [p_linear; p_angular] # in world frame
end

function center_of_mass(mechanism::Mechanism{T}, storage::Storage{T,N}, t::Int) where {T,N}
    r = zeros(T, 3)
    for (i,body) in enumerate(mechanism.bodies)
        r += body.mass * storage.x[i][t]
    end
    return r ./ total_mass(mechanism)
end

function total_mass(mechanism::Mechanism{T}) where T
    w = 0.0
    for body in mechanism.bodies
        w += body.mass
    end
    return w
end


text:
export has_dice_binary

function has_dice_binary()
    !isnothing(Sys.which("dice"))
end

function rundice(code::String;
            showinternal=false, skiptable=false,
            determinism=true, showsize=false,
            printstatebdd=false)
    dice = Sys.which("dice")
    isnothing(dice) && error("Ocaml dice binary not found in path")
    mktemp() do path, io
        write(io, code)
        close(io)
        flags = String[]
        if showinternal
            push!(flags, "-show-internal")
        end
        if skiptable
            push!(flags, "-skip-table")
        end
        if determinism
            push!(flags, "-determinism")
        end
        if showsize
            push!(flags, "-show-size")
        end
        if printstatebdd
            push!(flags, "-print-state-bdd")
        end
        cmd = `$dice $path $flags`
        Base.read(cmd, String)
    end
end

function num_nodes_ocml(code)
    out = rundice(code; skiptable=true, showsize=true)
    regex = r"================.*================\n(.+)\n"
    size_str = match(regex, out)
    @assert size_str !== nothing "$out did not match expected pattern"
    Base.parse(Int, size_str[1])
end

function infer_ocml(code)
    out = rundice(code)
    regex = r"\ntrue\t(.+)\t\n"
    size_str = match(regex, out)
    @assert size_str !== nothing "$out did not match expected pattern"
    Base.parse(Float64, size_str[1])
end


text:
using Test
using ParserCombinator
using GraphIO.DOT

@testset "DOT" begin
    g = CompleteGraph(6)
    dg = DiGraph(4)
    for e in [Edge(1,2), Edge(1,3), Edge(2,2), Edge(2,3), Edge(4,1), Edge(4,3)]
        add_edge!(dg, e)
    end
    fname = joinpath(testdir, "testdata", "twographs.dot")
    read_test(DOTFormat(), g, "g1", fname, testfail=true)
    read_test(DOTFormat(), dg, "g2", fname)
    read_test_mult(DOTFormat(), Dict{String,AbstractGraph}("g1"=>g, "g2"=>dg), fname)
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiExtensionsV1beta1NetworkPolicyPeer <: SwaggerModel
    ipBlock::Any # spec type: Union{ Nothing, IoK8sApiExtensionsV1beta1IPBlock } # spec name: ipBlock
    namespaceSelector::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1LabelSelector } # spec name: namespaceSelector
    podSelector::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1LabelSelector } # spec name: podSelector

    function IoK8sApiExtensionsV1beta1NetworkPolicyPeer(;ipBlock=nothing, namespaceSelector=nothing, podSelector=nothing)
        o = new()
        validate_property(IoK8sApiExtensionsV1beta1NetworkPolicyPeer, Symbol("ipBlock"), ipBlock)
        setfield!(o, Symbol("ipBlock"), ipBlock)
        validate_property(IoK8sApiExtensionsV1beta1NetworkPolicyPeer, Symbol("namespaceSelector"), namespaceSelector)
        setfield!(o, Symbol("namespaceSelector"), namespaceSelector)
        validate_property(IoK8sApiExtensionsV1beta1NetworkPolicyPeer, Symbol("podSelector"), podSelector)
        setfield!(o, Symbol("podSelector"), podSelector)
        o
    end
end # type IoK8sApiExtensionsV1beta1NetworkPolicyPeer

const _property_map_IoK8sApiExtensionsV1beta1NetworkPolicyPeer = Dict{Symbol,Symbol}(Symbol("ipBlock")=>Symbol("ipBlock"), Symbol("namespaceSelector")=>Symbol("namespaceSelector"), Symbol("podSelector")=>Symbol("podSelector"))
const _property_types_IoK8sApiExtensionsV1beta1NetworkPolicyPeer = Dict{Symbol,String}(Symbol("ipBlock")=>"IoK8sApiExtensionsV1beta1IPBlock", Symbol("namespaceSelector")=>"IoK8sApimachineryPkgApisMetaV1LabelSelector", Symbol("podSelector")=>"IoK8sApimachineryPkgApisMetaV1LabelSelector")
Base.propertynames(::Type{ IoK8sApiExtensionsV1beta1NetworkPolicyPeer }) = collect(keys(_property_map_IoK8sApiExtensionsV1beta1NetworkPolicyPeer))
Swagger.property_type(::Type{ IoK8sApiExtensionsV1beta1NetworkPolicyPeer }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiExtensionsV1beta1NetworkPolicyPeer[name]))}
Swagger.field_name(::Type{ IoK8sApiExtensionsV1beta1NetworkPolicyPeer }, property_name::Symbol) =  _property_map_IoK8sApiExtensionsV1beta1NetworkPolicyPeer[property_name]

function check_required(o::IoK8sApiExtensionsV1beta1NetworkPolicyPeer)
    true
end

function validate_property(::Type{ IoK8sApiExtensionsV1beta1NetworkPolicyPeer }, name::Symbol, val)
end


text:
using Distributions

function evsi_an(μ = [0, 0], σ = [1, 1], nsamples = [10000, 10000])
    m   = μ[1] - μ[2]
    s   = σ .* σ
    ps  = nsamples .* s
    θ   = sqrt(s[1] * (ps[1] / (ps[1] + 1)) + s[2] * (ps[2] / (ps[2] + 1)))
    (θ * sqrt(2 / π) * exp((-m^2) / (2 * θ^2)) + m * erf(m / (θ * sqrt(2))) - abs(m)) / 2
end

function evsi_sim(μ = [0, 0], σ = [1, 1], nsamples = [10000, 10000]; nsims = 1000000)

    ndists        = length(μ)
    true_values   = Array(Float64, nsims, ndists)
    prior_weight  = Array(Float64, ndists)
    sample_weight = Array(Float64, ndists)
    sample_sd     = Array(Float64, ndists)

    for dist in 1:ndists

        true_values[:, dist] = rand(Normal(μ[dist], σ[dist]), nsims)
        prior_weight[dist]   = (1 / σ[dist]) / (nsamples[dist] + (1 / σ[dist]))
        sample_weight[dist]  = nsamples[dist] / (nsamples[dist] + (1 / σ[dist]))
        sample_sd[dist]      = sqrt(1 / nsamples[dist])

    end

    benefit = Array(Float64, nsims)

    for sim in 1:nsims

        posterior_mean = Array(Float64, ndists)

        for dist in 1:ndists
            posterior_mean[dist] = μ[dist] * prior_weight[dist]
            if nsamples[dist] > 0
                posterior_mean[dist] +=
                    rand(Normal(true_values[sim, dist], sample_sd[dist])) *
                        sample_weight[dist]
            end
        end

        benefit[sim] = true_values[sim, indmax(posterior_mean)]

    end

    mean(benefit) - maximum(μ)

end


text:
function moi_mimo(optimizer, seed, n; verbose = false, test = false, scalar = false)

    MOI.empty!(optimizer)
    if test
        @test MOI.is_empty(optimizer)
    end

    m = 10n
    s, H, y, L = mimo_data(seed, m, n)

    nvars = ProxSDP.sympackedlen(n + 1)

    X = MOI.add_variables(optimizer, nvars)
    if scalar
        for i in 1:nvars
            MOI.add_constraint(optimizer, MOI.SingleVariable(X[i]), MOI.LessThan(1.0))
            MOI.add_constraint(optimizer, MOI.SingleVariable(X[i]), MOI.GreaterThan(-1.0))
        end
    else
        MOI.add_constraint(optimizer,
                            MOI.VectorAffineFunction(
                                MOI.VectorAffineTerm.(
                                    collect(1:nvars), MOI.ScalarAffineTerm.(1.0, X)),
                                -ones(nvars)),
                            MOI.Nonpositives(nvars))
        MOI.add_constraint(optimizer,
                            MOI.VectorAffineFunction(
                                MOI.VectorAffineTerm.(
                                    collect(1:nvars), MOI.ScalarAffineTerm.(-1.0, X)),
                                -ones(nvars)),
                            MOI.Nonpositives(nvars))
    end

    Xsq = Matrix{MOI.VariableIndex}(undef, n+1,n+1)
    ivech!(Xsq, X)
    Xsq = Matrix(Symmetric(Xsq,:U))

    vov = MOI.VectorOfVariables(X)
    cX = MOI.add_constraint(optimizer, vov, MOI.PositiveSemidefiniteConeTriangle(n+1))
    if scalar
        for i in 1:n+1
            MOI.add_constraint(optimizer, MOI.SingleVariable(Xsq[i,i]), MOI.EqualTo(1.0))
        end
    else
        MOI.add_constraint(optimizer,
                            MOI.VectorAffineFunction(
                                MOI.VectorAffineTerm.(
                                    collect(1:n+1), MOI.ScalarAffineTerm.(1.0, [Xsq[i,i] for i in 1:n+1])),
                                -ones(n+1)),
                            MOI.Zeros(n+1))
    end

    objf_t = vec([MOI.ScalarAffineTerm(L[i,j], Xsq[i,j]) for i in 1:n+1, j in 1:n+1])
    MOI.set(optimizer, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), MOI.ScalarAffineFunction(objf_t, 0.0))

    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    MOI.optimize!(optimizer)

    objval = MOI.get(optimizer, MOI.ObjectiveValue())

    stime = -1.0
    try
        stime = MOI.get(optimizer, MOI.SolveTime())
    catch
        println("could not query time")
    end

    Xsq_s = MOI.get.(optimizer, MOI.VariablePrimal(), Xsq)

    if test
        for i in 1:n+1, j in 1:n+1
            @test 1.01 > abs(Xsq_s[i,j]) > 0.99
        end
    end

    verbose && mimo_eval(s, H, y, L, Xsq_s)

    rank = -1
    status = 0
    if MOI.get(optimizer, MOI.TerminationStatus()) == MOI.OPTIMAL
        status = 1
    end
    return (objval, stime, rank, status)
end


text:
module BenchDataType

using BenchmarkTools
using Statistics
using Gadfly
using MLStyle
using DataFrames
import Match
import Rematch
using ..ArbitrarySampler
using ..Utils

@data NormalData begin
    Normal1(a::Int, b::Int, c::Any)
    Normal2(a::Int)
end

@data GenericData{A} begin
    Generic1(a::A, b::Int)
    Generic2(A)
end

mod′(n) = x -> mod(x, n)

specs = [
    :s1 => @spec(Normal1(::Int, ::Int, ::Real || ::Complex || _)),
    :s2 => @spec(Normal1(::Int, ::Int, Normal2(1 || ::Int))),
    :s3 => @spec(Generic1(3 || _, 4 || ::Int)),
    :s4 => @spec(Generic2((1 || _, 2 || _))),
    :s5 => @spec(Generic1(::String || _, 0) || Generic2(::String || _)),
    :s6 => @spec(
        Generic1(_, ::Int) || Generic2(_) || Normal2(::Int) || Normal1(::Int, ::Int, _)
    ),
    :_ => @spec(_),
]

implementations = [
    :MLStyle => (@λ begin
        Normal1(c = ::Number) -> 1
        Normal1(c = Normal2(1)) -> 2
        Generic1(a = 3, b = 3) -> 3
        Generic2((1, 2)) -> 4
        ::GenericData{String} -> 5
        _ -> 0
    end),
    :Rematch => function (x)
        Rematch.@match x begin
            Normal1(_, _, _::Number) => 1
            Normal1(_, _, Normal2(1)) => 2
            Generic1(3, 3) => 3
            Generic2((1, 2)) => 4
            _::GenericData{String} => 5
            _ => 0
        end
    end,
    Symbol("Match.jl") => function (x)
        Match.@match x begin
            Normal1(_, _, _::Number) => 1
            Normal1(_, _, Normal2(1)) => 2
            Generic1(3, 3) => 3
            Generic2((1, 2)) => 4
            _::GenericData{String} => 5
            _ => 0
        end
    end,
]

records = NamedTuple{(:time_mean, :implementation, :case)}[]
for (spec_id, spec) in specs
    # group_key = string(spec_id)
    # suite[group_key] = BenchmarkGroup()
    for (impl_id, impl_fn) in implementations
        bench′ =
            @benchmark $impl_fn(sample) setup = (sample = $generate($spec)) samples = 2000
        time′ = mean(bench′.times)
        @info :bench (spec_id, impl_id, time′)
        push!(records, (time_mean = time′, implementation = impl_id, case = spec_id))
    end
end

df = DataFrame(records)

@info df

report_meantime, df_time = report(df, Guide.title("Data Types"); benchfield = :time_mean)

open("stats/bench-datatype.txt", "w") do f
    write(f, string(df))
end

draw(SVG("stats/bench-datatype.svg", 14inch, 6inch), report_meantime)
end


text:
#=

Author: Shunsuke Hori

=#

using Polynomials

struct LQFilter{TR<:Real, TI<:Integer, TF<:AbstractFloat}
    d::Vector{TF}
    h::TR
    y_m::Vector{TF}
    m::TI
    ϕ::Vector{TF}
    β::TR
    ϕ_r::Union{Vector{TF},Void}
    k::Union{TI,Void}
end

"""
Parameters
----------
d : Vector
        The order of the coefficients: [d_0, d_1, ..., d_m]
h : Real
        Parameter of the objective function (corresponding to the
        quadratic term)
y_m : Vector
        Initial conditions for y
r : Vector
        The order of the coefficients: [r_0, r_1, ..., r_k]
        (optional, if not defined -> deterministic problem)
β : Real or nothing
        Discount factor (optional, default value is one)
h_eps :
"""
function LQFilter{TR<:Real}(d::Vector{TR},
                   h::TR,
                   y_m::Vector{TR};
                   r::Union{Vector{TR},Void}=nothing,
                   β::Union{TR,Void}=nothing,
                   h_eps::Union{TR,Void}=nothing,
                   )

    m = length(d) - 1

    m == length(y_m) ||
        throw(ArgumentError("y_m and d must be of same length = $m"))

    #---------------------------------------------
    # Define the coefficients of ϕ up front
    #---------------------------------------------

    ϕ = Vector{TR}(2m + 1)
    for i in -m:m
        ϕ[m-i+1] = sum(diag(d*d', -i))
    end
    ϕ[m+1] = ϕ[m+1] + h

    #-----------------------------------------------------
    # If r is given calculate the vector ϕ_r
    #-----------------------------------------------------

    if r == nothing
        k=nothing
        ϕ_r = nothing
    else
        k = size(r, 1) - 1
        ϕ_r = Vector{TR}(2k + 1)

        for i = -k:k
            ϕ_r[k-i+1] = sum(diag(r*r', -i))
        end

        if h_eps != nothing
            ϕ_r[k+1] = ϕ_r[k+1] + h_eps
        end
    end

    #-----------------------------------------------------
    # If β is given, define the transformed variables
    #-----------------------------------------------------
    if β == nothing
        β = 1.0
    else
        d = β.^(collect(0:m)/2) * d
        y_m = y_m * β.^(- collect(1:m)/2)
    end

    return LQFilter(d, h, y_m, m, ϕ, β, ϕ_r, k)
end

"""
This constructs the matrices W and W_m for a given number of periods N
"""
function construct_W_and_Wm(lqf::LQFilter, N::Integer)

    d, m = lqf.d, lqf.m

    W = zeros(N + 1, N + 1)
    W_m = zeros(N + 1, m)

    #---------------------------------------
    # Terminal conditions
    #---------------------------------------

    D_m1 = zeros(m + 1, m + 1)
    M = zeros(m + 1, m)

    # (1) Constuct the D_{m+1} matrix using the formula

    for j in 1:(m+1)
        for k in j:(m+1)
            D_m1[j, k] = dot(d[1:j, 1], d[k-j+1:k, 1])
        end
    end

    # Make the matrix symmetric
    D_m1 = D_m1 + D_m1' - diagm(diag(D_m1))

    # (2) Construct the M matrix using the entries of D_m1

    for j in 1:m
        for i in (j + 1):(m + 1)
            M[i, j] = D_m1[i-j, m+1]
        end
    end
    M

    #----------------------------------------------
    # Euler equations for t = 0, 1, ..., N-(m+1)
    #----------------------------------------------
    ϕ, h = lqf.ϕ, lqf.h

    W[1:(m + 1), 1:(m + 1)] = D_m1 + h * eye(m + 1)
    W[1:(m + 1), (m + 2):(2m + 1)] = M

    for (i, row) in enumerate((m + 2):(N + 1 - m))
        W[row, (i + 1):(2m + 1 + i)] = ϕ'
    end

    for i in 1:m
        W[N - m + i + 1 , end-(2m + 1 - i)+1:end] = ϕ[1:end-i]
    end

    for i in 1:m
        W_m[N - i + 2, 1:(m - i)+1] = ϕ[(m + 1 + i):end]
    end

    return W, W_m
end

"""
This function calculates z_0 and the 2m roots of the characteristic equation
associated with the Euler equation (1.7)
Note:
------
`poly(roots)` from Polynomial.jll defines a polynomial using its roots that can be
evaluated at any point by `polyval(Poly,x)`. If x_1, x_2, ... , x_m are the roots then
    polyval(poly(roots),x) = (x - x_1)(x - x_2)...(x - x_m)
"""
function roots_of_characteristic(lqf::LQFilter)
    m, ϕ = lqf.m, lqf.ϕ

    # Calculate the roots of the 2m-polynomial
    ϕ_poly=Poly(ϕ[end:-1:1])
    proots = roots(ϕ_poly)
    # sort the roots according to their length (in descending order)
    roots_sorted = sort(proots, by=abs)[end:-1:1]
    z_0 = sum(ϕ) / polyval(poly(proots), 1.0)
    z_1_to_m = roots_sorted[1:m]     # we need only those outside the unit circle
    λ = 1 ./ z_1_to_m
    return z_1_to_m, z_0, λ
end

"""
This function computes the coefficients {c_j, j = 0, 1, ..., m} for
    c(z) = sum_{j = 0}^{m} c_j z^j
Based on the expression (1.9). The order is
    c_coeffs = [c_0, c_1, ..., c_{m-1}, c_m]
"""
function coeffs_of_c(lqf::LQFilter)
    m = lqf.m
    z_1_to_m, z_0, λ = roots_of_characteristic(lqf)
    c_0 = (z_0 * prod(z_1_to_m) * (-1.0)^m)^(0.5)
    c_coeffs = coeffs(poly(z_1_to_m)) * z_0 / c_0
    return c_coeffs
end

"""
This function calculates {λ_j, j=1,...,m} and {A_j, j=1,...,m}
of the expression (1.15)
"""
function solution(lqf::LQFilter)
    z_1_to_m, z_0, λ = roots_of_characteristic(lqf)
    c_0 = coeffs_of_c(lqf)[end]
    A = zeros(lqf.m)
    for j in 1:m
        denom = 1 - λ/λ[j]
        A[j] = c_0^(-2) / prod(denom[1:m .!= j])
    end
    return λ, A
end

"""
This function constructs the covariance matrix for x^N (see section 6)
for a given period N
"""
function construct_V(lqf::LQFilter; N::Integer=nothing)
    if N == nothing
        error("N must be provided!!")
    end
    if !(typeof(N) <: Integer)
        throw(ArgumentError("N must be Integer!"))
    end

    ϕ_r, k = lqf.ϕ_r, lqf.k
    V = zeros(N, N)
    for i in 1:N
        for j in 1:N
            if abs(i-j) <= k
                V[i, j] = ϕ_r[k + abs(i-j)+1]
            end
        end
    end
    return V
end

"""
Assuming that the u's are normal, this method draws a random path
for x^N
"""
function simulate_a(lqf::LQFilter, N::Integer)
    V = construct_V(N + 1)
    d = MVNSampler(zeros(N + 1), V)
    return rand(d)
end

"""
This function implements the prediction formula discussed is section 6 (1.59)
It takes a realization for a^N, and the period in which the prediciton is formed

Output:  E[abar | a_t, a_{t-1}, ..., a_1, a_0]
"""
function predict(lqf::LQFilter, a_hist::Vector, t::Integer)
    N = length(a_hist) - 1
    V = construct_V(N + 1)

    aux_matrix = zeros(N + 1, N + 1)
    aux_matrix[1:t+1 , 1:t+1 ] = eye(t + 1)
    L = chol(V)'
    Ea_hist = inv(L) * aux_matrix * L * a_hist

    return Ea_hist
end

"""
- if t is NOT given it takes a_hist (Vector or Array) as a deterministic a_t
- if t is given, it solves the combined control prediction problem (section 7)
　(by default, t == nothing -> deterministic)

for a given sequence of a_t (either determinstic or a particular realization),
it calculates the optimal y_t sequence using the method of the lecture

Note:
------
lufact normalizes L, U so that L has unit diagonal elements
To make things cosistent with the lecture, we need an auxiliary diagonal
matrix D which renormalizes L and U
"""

function optimal_y(lqf::LQFilter, a_hist::Vector, t = nothing)
    β, y_m, m = lqf.β, lqf.y_m, lqf.m

    N = length(a_hist) - 1
    W, W_m = construct_W_and_Wm(lqf, N)

    F = lufact(W, Val{true})

    L, U = F[:L], F[:U]
    D = diagm(1.0./diag(U))
    U = D * U
    L = L * diagm(1.0./diag(D))

    J = flipdim(eye(N + 1), 2)

    if t == nothing                      # if the problem is deterministic
        a_hist = J * a_hist

        #--------------------------------------------
        # Transform the a sequence if β is given
        #--------------------------------------------
        if β != 1
            a_hist =  reshape(a_hist * (β^(collect(N:0)/ 2)), N + 1, 1)
        end

        a_bar = a_hist - W_m * y_m        # a_bar from the lecutre
        Uy = \(L, a_bar)                  # U @ y_bar = L^{-1}a_bar from the lecture
        y_bar = \(U, Uy)                  # y_bar = U^{-1}L^{-1}a_bar
        # Reverse the order of y_bar with the matrix J
        J = flipdim(eye(N + m + 1), 2)
        y_hist = J * vcat(y_bar, y_m)     # y_hist : concatenated y_m and y_bar
        #--------------------------------------------
        # Transform the optimal sequence back if β is given
        #--------------------------------------------
        if β != 1
            y_hist = y_hist .* β.^(- collect(-m:N)/2)
        end

    else                                  # if the problem is stochastic and we look at it
        Ea_hist = reshape(predict(a_hist, t), N + 1, 1)
        Ea_hist = J * Ea_hist

        a_bar = Ea_hist - W_m * y_m       # a_bar from the lecutre
        Uy = \(L, a_bar)                  # U @ y_bar = L^{-1}a_bar from the lecture
        y_bar = \(U, Uy)                  # y_bar = U^{-1}L^{-1}a_bar

        # Reverse the order of y_bar with the matrix J
        J = flipdim(eye(N + m + 1), 2)
        y_hist = J * vcat(y_bar, y_m)     # y_hist : concatenated y_m and y_bar
    end
    return y_hist, L, U, y_bar
end


text:
# # Chapter 2 - Python..ではなくJulia入門

## #src
# Juliaでは、改行あり出力は`println`
println("hello")

## #src
# ## 2.2 変数と変数出力
# 変数や計算の埋め込みは`$()`で行う。
# 型を調べるのは`type`ではなく`typeof`を用いる。
num = 10
string = "Hello World"
println("num=$(num), typeof(num)=$(typeof(num))")
println("string=$(string), typeof(string)=$(typeof(string))")

## #src
# フォーマット書式ありの出力は標準ライブラリに含まれている`Printf.@sprintf`を使う。
# Pythonの`.format`に近い表記になる。
using Printf

num = 10.5678
println(@sprintf("小数2位:%.2f, 小数5位:%.5f, 小数2位の指数表記:%.2e", num, num, num))

## #src
# ## 2.3 データ構造
# ### 2.3.1 リスト
# Juliaではベクトル型が組み込まれているので、Pythonのようにリストとnumpy.ndarrayの相互変換
# を意識するする必要はない
price1 = [120, 180]
price2 = [300, 150]
all = [price1, price2]
all

## #src
# ### 2.3.2 リストの大きさ
length(price1)

## #src
# ### 2.3.3 リストの各要素へのアクセス
# インデックスが0から始まるPythonと違い、Juliaのインデックスは1から始まる
## 最初(1番目)の要素: 120
price1[1]

## #src
# 最初の要素はfirstでもアクセスできる
first(price1)

## #src
#-
## 2番目の要素
all[2]

## #src
#-
## 2行 2番目の要素: 150
all[2][2]

## #src
# ### 2.3.4 リストの拡張
# リストの拡張は`push!`を使う
## 空のリスト
price = []

## #src
#-
## 要素の追加: [120, 180]
push!(price, [120, 180])

## #src
#-
## 要素の追加: [[120, 180], [300, 150]]
push!(price, [300, 150])

## #src
# ### 2.3.5 Vector, Array
# Juliaではリストと行列、ベクトルの区別を意識しなくて良い
price1 = [120, 180]

## #src
#-
## Juliaはcolumn oriented
price2 = [200 300; 50 150.0]

## #src
# 要素の型は`eltype`を使って確認できる
eltype(price1)

## #src
#-
eltype(price2)

## #src
# ### 2.3.6 Arrayの大きさ
# `Array`のサイズは`size`で調べる
size(price1)

## #src
#-
size(price2)

## #src
# ### 2.3.7 Arrayのスライシング
price = [120, 180, 300, 150]

## 最初(1番目)の要素: 120
price[1]

## #src
#-
## 2番目の要素: 180
price[2]

## #src
#-
## 最後の要素: 150
price[end]

## #src
#-
## 最後から2番目の要素: 300
price[end-1]

## #src
#-
## 2番目から3番目の要素: [180, 300]
price[2:3]

## #src
#-
## 最初から2番目までの要素: [120, 180]
price[1:2]

## #src
#-
## 最初から3番目までの要素: [120. 180, 300]
price[1:end-1]

## #src
#-
## 最初から最後までの要素
price[:]

## #src
# ### 2.3.8 Arrayの簡単な計算
price = [120, 180, 300, 150]

## 和: 750
sum(price)

## #src
#-
## 積: 972000000
prod(price)

## #src
#-
## 二乗: [14400, 32400, 90000, 22500]
## `.`でbroadcastする
price .^ 2

## #src
#-
## 平方根: [10.9544515, 13.41640786, 17.32050808, 12.247448711]
sqrt.(price)

## #src
# 統計関連の計算には標準ライブラリであるStatisticsを使う
#-
using Statistics

## 平均値: 187.5
mean(price)

## #src
#-
## 分散: 4668.75
## Statistics.varはデフォルトで不偏分散を計算するので、標本分散は`corrected=false`とする
var(price, corrected=false)

## #src
#-
var(price)

## #src
#-
## 標準偏差
std(price, corrected=false)

## #src
#-
## 最小値: 120
minimum(price)

## #src
#-
## 最大値: 300
maximum(price)

## #src
#-
## 最小のインデックス: 1
argmin(price)

## #src
#-
## 最大のインデックス: 3
argmax(price)

## #src
# ### 2.3.9 Vectorのソート
price = [120, 180, 300, 150]

## 昇順に並び替えた値: [120, 150, 180, 300]
sort(price)

## #src
#-
## 降順に並び替えた値: [300, 180, 150, 120]
sort(price, rev=true)

## #src
#-
## 昇順に並び替えた値のインデックス: [1, 4, 2, 3]
sortperm(price)

## #src
#-
## 降順に並び替えた値のインデックス: [3, 2, 4, 1]
sortperm(price, rev=true)

## #src
# ### 2.3.10 Vectorのユニーク
price = [120, 180, 300, 150, 150]

## unique後: [120, 180, 300, 150]
## 順番は保存される
unique(price)

## #src
# ### 2.3.11 Vectorの条件を満たすインデックス
price = [120, 180, 300, 150]

## 条件を満たすインデックスを取得するにはfindall
findall(x -> x > 150, price)

## #src
#-
## マスクを作成するのはbroadcastを使えば簡単
indices = price .> 150

## #src
# ### 2.3.12 Arrayの結合
price1 = [120, 180]
price2 = [300, 150]

## horizontal方向への結合
hcat(price1, price2)

## #src
#-
## catを使っても書ける
cat(price1, price2; dims=2)

## #src
#-
## vertical方向への結合
vcat(price1, price2)

## #src
#-
## catを使っても書ける
cat(price1, price2; dims=1)

## #src
#-
price1 = reshape([120, 180], 1, 1, :)
price2 = reshape([300, 150], 1, 1, :)

## #src
#-
hcat(price1, price2)

## #src
#-
vcat(price1, price2)

## #src
#-
cat(price1, price2; dims=3)

## #src
# ### 2.3.13 Arrayの変形
price = [120, 180, 300, 150]

## #src
#-
## 2行 2列に変形
reshape(price, 2, 2)

## #src
#-
## 2行 2列に変形(自動)
reshape(price, 2, :)

## #src
#-
## 4行 1列に変形
reshape(price, 4, 1)

## #src
# ### 2.3.14 Arrayの繰り返し
price = [120, 180, 300, 150]

## 行方向に2回、列方向に1回繰り返す
repeat(price, 2, 1)

## #src
#-
## 行方向に1回、列方向に2回繰り返す
repeat(price, 1, 2)

## #src
#-
## 奥方向に2回、行・列方向に1回繰り返す
repeat(price, 1, 1, 2)

## #src
# ### 2.3.15 Arrayの次元拡張
price = [120, 180, 300, 150]

size(price)

## #src
#-
## 行方向に増やした後のsize
size(reshape(price, 1, :))

## #src
#-
## 列方向に増やした後のsize
size(reshape(price, :, 1))

## #src
# ### 2.3.16 Arrayの(軸方向の)回転
price = [120 180; 300 150]

## #src
#-
## 行方向に1つ回転する
circshift(price, (1, 0))

## #src
#-
## 列方向に1つ回転する
circshift(price, (0, 1))

## #src
#-
## 回転
rotr90(price, 1)

## #src
#-
rotl90(price, 2)

## #src
# ### 2.3.17 Arrayのインデックスの生成
price = [120, 180, 300, 150]

## 1:4
1:length(price)

## #src
#-
## [1, 2, 3, 4]
collect(1:length(price))

## #src
# ### 2.3.18 Arrayのインデックスのシャッフル
using Random

## シードを固定しておく
Random.seed!(1234)

price = [120, 180, 300, 150]

randperm(length(price))

## #src
# ### 2.3.19 単位行列の生成
using LinearAlgebra

Matrix(I, 5, 5)

## #src
# ### 2.3.20 CSV.jl & DataFrames.jl
# https://www.kaggle.com/c/house-prices-advanced-regression-techniques/data?select=train.csv
# からファイルをダウンロードしてdataフォルダに保存しているとする

## CSVファイルの読み込み
# pandas.DataFrameのようなDataFrameの構造を扱うにはDataFrames.jlを使う。
using CSV
using DataFrames

#-
data = CSV.read("../data/house-prices-advanced-regression-techniques/train.csv", DataFrame)

## #src
#-

## 5行目まで
first(data, 5)

## #src
#-
## カラムの一覧
names(data)

## #src
#-
## 建物の特徴(MSSubClass)
unique(data[!, :MSSubClass])

## #src
# ### 2.3.21 DataFrameからArrayへの変換
Array(data[!, [:GrLivArea, :SalePrice]])

## #src
# ### 2.3.22 DafaFrameのカラムの削除
length(names(data))

## #src
#-
## PoolAreaカラムを削除
select!(data, Not(:PoolArea))

## #src
#-
length(names(data))

## #src
# ## 2.4 グラフのプロット
# 別ファイル

## #src
# ### 2.5 if文とfor文
#-
values = [10, 3, 1, 5, 8, 6]

## 通常のfor文
## 空のリスト
passed_values = Int64[]
for v in values
    if v > 5
        push!(passed_values, v)
    end
end

## 結果を標準出力
println("5以上の値 $(passed_values)")

## #src
#-
## リスト内表記のforとif
passed_values = [v for v in values if v > 5]

## 結果を標準出力
println("5以上の値 $(passed_values)")

## #src
#-
## filterを使った表記
passed_values = filter(v -> v > 5, values)

## 結果を標準出力
println("5以上の値 $(passed_values)")

## #src
# ### 2.5.3 例外処理
#-
## Juliaでは1/0はInfになる
for i in [-1, 0, 1]
    println("$(1/i)")
end

## #src
#-
## 例外を発生させてみる
## for i in [-1, 0, 1]
##     println("$(i//i)")
##∂ end

## #src
#-
for i in [-1, 0, 1]
    try
        println("$(i//i)")
    catch
        continue
    end
end


text:
"Parameters for right-hand side function of Variational Partitioned Runge-Kutta methods."
const ParametersVPRKdegenerate = AbstractParametersVPRK{:vprk_degenerate}

"Variational partitioned Runge-Kutta integrator."
mutable struct IntegratorVPRKdegenerate{DT, TT, D, S,
                PT  <: ParametersVPRK{DT,TT},
                PPT <: ParametersVPRKdegenerate{DT,TT},
                ST  <: NonlinearSolver{DT},
                PST <: NonlinearSolver{DT},
                IT  <: InitialGuessIODE{TT}} <: AbstractIntegratorVPRK{DT,TT,D,S}

    params::PT
    pparams::PPT
    solver::ST
    projector::PST
    iguess::IT
    caches::CacheDict{PPT}

    function IntegratorVPRKdegenerate(params::ParametersVPRK{DT,TT,D,S},
                    pparams::ParametersVPRKdegenerate{DT,TT,D,S},
                    solver::ST, projector::PST, iguess::IT, caches) where {DT,TT,D,S,ST,PST,IT}
        new{DT, TT, D, S, typeof(params), typeof(pparams), ST, PST, IT}(params, pparams, solver, projector, iguess, caches)
    end

    function IntegratorVPRKdegenerate{DT,D}(equations::NamedTuple, tableau::TableauVPRK{TT}, Δt::TT; R=[1,1]) where {DT, TT, D}
        # get number of stages
        S = tableau.s

        # create solver params
        sparams = ParametersVPRK{DT,D}(equations, tableau, Δt)

        # create projector params
        pparams = ParametersVPRKdegenerate{DT,D}(equations, tableau, Δt, NamedTuple{(:V,:F)}((create_internal_stage_vector(DT,D,S), create_internal_stage_vector(DT,D,S))))

        # create cache dict
        caches = CacheDict(pparams)

        # create nonlinear solver
        solver = create_nonlinear_solver(DT, D*S, sparams, caches)

        # create projector
        projector = create_nonlinear_solver(DT, D, pparams, caches)

        # create initial guess
        iguess = InitialGuessIODE(get_config(:ig_extrapolation), equations[:v̄], equations[:f̄], Δt)

        # create integrator
        IntegratorVPRKdegenerate(sparams, pparams, solver, projector, iguess, caches)
    end

    function IntegratorVPRKdegenerate(equation::Union{IODE{DT},LODE{DT}}, tableau, Δt; kwargs...) where {DT}
        IntegratorVPRKdegenerate{DT, ndims(equation)}(get_functions(equation), tableau, Δt; kwargs...)
    end
end

function initial_guess!(int::IntegratorVPRKdegenerate{DT}, sol::AtomicSolutionPODE{DT},
                        cache::IntegratorCacheVPRK{DT}=int.caches[DT]) where {DT}
    for i in eachstage(int)
        evaluate!(int.iguess, sol.q̄, sol.p̄, sol.v̄, sol.f̄,
                              sol.q, sol.p, sol.v, sol.f,
                              cache.q̃, cache.ṽ,
                              tableau(int).q.c[i])
        for k in eachdim(int)
            int.solver.x[ndims(int)*(i-1)+k] = cache.ṽ[k]
        end
    end
end

function initial_guess_projection!(int::IntegratorVPRKdegenerate{DT}, sol::AtomicSolutionPODE{DT},
                                   cache::IntegratorCacheVPRK{DT}=int.caches[DT]) where {DT}
    for k in eachdim(int)
        int.projector.x[k] = sol.q[k]
    end
end

function compute_solution!(x::Vector{ST}, q̄::Vector{ST}, v̄::Vector{ST}, p̄::Vector{ST},
                           params::ParametersVPRKdegenerate{DT,TT,D,S}) where {ST,DT,TT,D,S}

    @assert length(x) == length(q̄) == length(p̄)

    # copy x to q
    q̄ .= x

    # compute p̄ = ϑ(q)
    params.equ[:ϑ](params.t̄ + params.Δt, q̄, v̄, p̄)
end

"Compute solution of degenerate symplectic partitioned Runge-Kutta methods."
function Integrators.function_stages!(x::Vector{ST}, b::Vector{ST},
                params::ParametersVPRKdegenerate{DT,TT,D,S},
                caches::CacheDict) where {ST,DT,TT,D,S}

    @assert length(x) == length(b)

    # get cache for internal stages
    cache = caches[ST]

    compute_solution!(x, cache.q̃, cache.ṽ, cache.p̃, params)

    # compute b = - [q̄-q-BV]
    for k in 1:div(D,2)
        b[k] = params.q̄[k] - cache.q̃[k]
        for i in 1:S
            b[k] += params.Δt * params.tab.q.b[i] * params.pparams[:V][i][k]
        end
    end

    # compute b = - [p̄-p-BF]
    for k in 1:div(D,2)
        b[div(D,2)+k] = params.p̄[k] - cache.p̃[k]
        for i in 1:S
            b[div(D,2)+k] += params.Δt * params.tab.p.b[i] * params.pparams[:F][i][k]
        end
    end
end

function Integrators.integrate_step!(int::IntegratorVPRKdegenerate{DT,TT}, sol::AtomicSolutionPODE{DT,TT},
                                     cache::IntegratorCacheVPRK{DT}=int.caches[DT]) where {DT,TT}
    # update nonlinear solver parameters from cache
    update_params!(int.params, sol)
    update_params!(int.pparams, sol)

    # compute initial guess
    initial_guess!(int, sol, cache)

    # reset solution
    reset!(sol, timestep(int))

    # call nonlinear solver
    solve!(int.solver)

    # print solver status
    print_solver_status(int.solver.status, int.solver.params)

    # check if solution contains NaNs or error bounds are violated
    check_solver_status(int.solver.status, int.solver.params)

    # compute vector fields at internal stages
    compute_stages!(int.solver.x, cache.Q, cache.V, cache.P, cache.F, int.params)

    # compute unprojected solution
    update_solution!(int, sol, cache)

    # copy vector fields
    for i in eachstage(int)
        int.pparams.pparams[:V][i] .= cache.V[i]
        int.pparams.pparams[:F][i] .= cache.F[i]
    end

    # set initial guess for projection
    initial_guess_projection!(int, sol)

    # call solution solver
    solve!(int.projector)

    # print solver status
    print_solver_status(int.projector.status, int.projector.params)

    # check if solution contains NaNs or error bounds are violated
    check_solver_status(int.projector.status, int.projector.params)

    # compute projection vector fields
    compute_solution!(int.projector.x, sol.q, cache.ṽ, sol.p, int.pparams)

    # copy solution to initial guess
    update_vector_fields!(int.iguess, sol.t, sol.q, sol.p, sol.v, sol.f)
end