text:
function Projection(args::Vector{<: Union{Pair, Tuple{String, Any}}})

    str = ""

    for (arg, val) in args
        if val === nothing
            str *= "+" * string(arg)
        else
            str *= "+$arg=$val"
        end
        str *= " "
    end

    return Projection(str)
end

Projection(args::Pair...) = Projection(args)

function Proj4.transform(
        src::Projection, dest::Projection,
        x::AbstractVecOrMat{<: AbstractFloat},
        y::AbstractVecOrMat{<: AbstractFloat},
        z::AbstractVecOrMat{<: AbstractFloat} = zeros(size(y))
    )

    @assert size(x) == size(y) == size(z)

    xv, yv, zv = vec(Float64.(x)), vec(Float64.(y)), vec(Float64.(z))

    Proj4.transform!(src, dest, xv, yv, zv)

    return xv, yv, zv
end

function Proj4.transform!(
    src::Projection, dest::Projection,
    x::VecOrMat{Float64},
    y::VecOrMat{Float64},
    z::VecOrMat{Float64} = zeros(size(y))
    )

    return Proj4.transform!(src, dest, vec(x), vec(y), vec(z))
end

function Proj4.transform!(src::Projection, dest::Projection, x::Vector{Float64}, y::Vector{Float64}, z::Vector{Float64} = zeros(length(y)); radians = false)
    @assert length(x) == length(y) == length(z) "Input vectors must have the same length!"
    if !radians && is_latlong(src)
        x .= deg2rad.(x)
        y .= deg2rad.(y)
    end
    Proj4._transform!(src.rep, dest.rep, length(x), 1, x, y, z)
    if !radians && is_latlong(dest)
        x .= deg2rad.(x)
        y .= deg2rad.(y)
    end
    position
end

Proj4.transform(src::Projection, dest::Projection, x::Point2f0) = Point2f0(Proj4.transform(src, dest, [x[1], x[2]]))

# Proj4.transform!(src::Projection, dest::Projection, points::T) where T <: StructVector{Point2f0} = Proj4.transform!(src, dest, points.data)

const PtrOrVecCdouble = Union{Ptr{Cdouble}, Vector{Cdouble}}
"Low level interface to libproj transform, C_NULL can be passed in for z, if it's 2-dimensional"
function Proj4._transform!(src_ptr::Ptr{Cvoid}, dest_ptr::Ptr{Cvoid}, point_count::Integer, point_stride::Integer,
                     x::PtrOrVecCdouble, y::PtrOrVecCdouble, z::PtrOrVecCdouble)
    @assert src_ptr != C_NULL && dest_ptr != C_NULL
    err = ccall((:pj_transform, Proj4.libproj), Cint, (Ptr{Cvoid}, Ptr{Cvoid}, Clong, Cint, Ptr{Cdouble}, Ptr{Cdouble},
                Ptr{Cdouble}), src_ptr, dest_ptr, point_count, point_stride, x, y, z)
    err != 0 && error("transform error: $(_strerrno(err))")
end

"""
    LonLat()

Standard latitude-longitude grid.
"""
LonLat() = Projection("+proj=longlat +datum=WGS84")
const PlateCarree = LonLat

"""
    TransverseMercator(; central_longitude=0.0, central_latitude=0.0, false_easting=0.0, false_northing=0.0, scale_factor=1.0, units = "m")

## Parameters

* `central_longitude`:
    The true longitude of the central meridian in degrees.
    Defaults to 0.
* `central_latitude`:
    The true latitude of the planar origin in degrees. Defaults to 0.
* `false_easting`:
    X offset from the planar origin in metres. Defaults to 0.
* `false_northing`:
    Y offset from the planar origin in metres. Defaults to 0.
* `scale_factor`:
    Scale factor at the central meridian. Defaults to 1.
"""
function TransverseMercator(; central_longitude=0.0, central_latitude=0.0, false_easting=0.0, false_northing=0.0, scale_factor=1.0, units = "m")
    proj4_params = [
        ("proj", "tmerc"),
        ("lon_0", central_longitude),
        ("lat_0", central_latitude),
        ("k", scale_factor),
        ("x_0", false_easting),
        ("y_0", false_northing),
        ("units", units)
    ]

    return Projection(proj4_params)
end

"""
UTM(zone, southern_hemisphere=false; units = "m")

## Parameters
zone:
    The numeric zone of the UTM required.
`southern_hemisphere`:
    Set to True if the zone is in the southern hemisphere. Defaults to
    False.
"""
function UTM(zone, southern_hemisphere=false; units = "m")
    proj4_params = [
        ("proj", "utm"),
        ("units", units),
        ("zone", zone)
    ]

    if southern_hemisphere
        push!(proj4_params, ("south", nothing))
    end

    return Projection(proj4_params)
end

"""
    Mercator(; central_longitude=0.0,
            latitude_true_scale=nothing, false_easting=0.0, false_northing=0.0,
            scale_factor=nothing, units = "m")

## Parameters
*  `central_longitude`:
   The central longitude. Defaults to 0.
*  `latitude_true_scale`:
   The latitude where the scale is 1. Defaults to 0 degrees.
*  `false_easting`:
   X offset from the planar origin in metres. Defaults to 0.
*  `false_northing`:
   Y offset from the planar origin in metres. Defaults to 0.
*  `scale_factor`:
   Scale factor at natural origin. Defaults to unused.

!!! note
    Only one of ``latitude_true_scale`` and ``scale_factor`` should
    be included.
"""
function Mercator(; central_longitude=0.0,
                    latitude_true_scale=nothing, false_easting=0.0, false_northing=0.0,
                    scale_factor=nothing, units = "m"
        )

    @assert isnothing(latitude_true_scale) && isnothing(scale_factor) "Only one of latitude_true_scale and scale_factor can be included!"

    proj4_params = [
        ("proj", "merc"),
        ("lon_0", central_longitude),
        ("x_0", false_easting),
        ("y_0", false_northing),
        ("units", units)
    ]

    if latitude_true_scale !== nothing
        push!(proj4_params, ("lat_true_scale", latitude_true_scale))
    elseif scale_factor !== nothing
        push!(proj4_params, ("k_0", scale_factor))
    end

    return Projection(proj4_params)

end

"""
AlbersEqualArea(; central_longitude=0.0, central_latitude=0.0,
                 false_easting=0.0, false_northing=0.0,
                 standard_parallels=(20.0, 50.0)
        )

This projection is conic and equal-area, and is commonly used for maps of
the conterminous United States.

## Parameters
* `central_longitude`:
    The central longitude. Defaults to 0.
* `central_latitude`:
    The central latitude. Defaults to 0.
* `false_easting`:
    X offset from planar origin in metres. Defaults to 0.
* `false_northing`:
    Y offset from planar origin in metres. Defaults to 0.
* `standard_parallels`:
    The one or two latitudes of correct scale. Defaults to (20, 50).
"""
function AlbersEqualArea(; central_longitude=0.0, central_latitude=0.0,
                 false_easting=0.0, false_northing=0.0,
                 standard_parallels=(20.0, 50.0), units = "m")

    proj4_params = [
        ("proj", "aea"),
        ("lon_0", central_longitude),
        ("lat_0", central_latitude),
        ("x_0", false_easting),
        ("y_0", false_northing),
        ("lat_1", standard_parallels[1]),
        ("lat_2", standard_parallels[2]),
        ("units", units)
    ]

    return Projection(proj4_params)
end

"""
    Robinson(; central_longitude=0.0, false_easting=0.0, false_northing=0.0, units = "m")

The Robinson projection is the successor of the van der Grinten series of projections,
and was created specifically for the purpose of world mapping.

## Parameters
* `central_longitude`:
    The central longitude. Defaults to 0.
* `false_easting`:
    X offset from planar origin in metres. Defaults to 0.
* `false_northing`:
    Y offset from planar origin in metres. Defaults to 0.

"""
function Robinson(;
    central_longitude=0.0,
    radius = nothing,
    false_easting=0.0,
    false_northing=0.0,
    units = "m"
    )

    proj4_params = [
        ("proj", "robin"),
        ("lon_0", central_longitude),
        ("x_0", false_easting),
        ("y_0", false_northing),
    ]

    if radius !== nothing
        push!(proj4_params, ("R", radius))
    end

    return Projection(proj4_params)

end

"""
WinkelTripel(; central_longitude = 0, first_parallel = 0, false_easting=0,
            false_northing=0, units = "m")

The Winkel tripel projection was the successor of the Robinson projection, and
minimizes the triplet (_tripel_) of distortions: area, direction, and distance.
"""
function WinkelTripel(;
            central_longitude = 0,
            first_parallel = 0,
            false_easting=0,
            false_northing=0,
            units = "m"
        )

    proj4_params = [
        ("proj", "wintri"),
        ("lon_0", central_longitude),
        ("lat_1", first_parallel),
        ("x_0", false_easting),
        ("y_0", false_northing),
        ("units", units)
    ]

    return Projection(proj4_params)
end


text:
export iauPmp
"""
P-vector subtraction.

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  vector/matrix support function.

Given:
   a        double[3]      first p-vector
   b        double[3]      second p-vector

Returned:
   amb      double[3]      a - b

Note:
   It is permissible to re-use the same array for any of the
   arguments.

This revision:  2013 June 18

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# void iauPmp(double a[3], double b[3], double amb[3])

function iauPmp(a::AbstractVector{<:Real}, b::AbstractVector{<:Real})
   amb = zeros(Float64, 3)

   ccall((:iauPmp, libsofa_c), Cdouble,
       (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
       convert(Array{Float64, 1}, a),
       convert(Array{Float64, 1}, b),
       amb)

   return SVector{3}(amb)
end

ext: jl


text:
include("libavcodec_h.jl")

#include("ac3_parser.jl")
#include("adts_parser.jl")
include("avcodec.jl")
#include("avdct.jl")
#include("d3d11va.jl")
#include("dirac.jl")
#include("dv_profile.jl")
#include("jni.jl")
#include("mediacodec.jl")
#include("qsv.jl")
#include("vaapi.jl")
#include("vdpau.jl")
#include("version.jl")
#include("videotoolbox.jl")
#include("vorbis_parser.jl")
#include("xvmc.jl")


text:
mutable struct SparseIterate{T, N} <: AbstractArray{T, N}
    nzval::Vector{T}            # stored nonzero values (1:nnz)
    nzval2ind::Vector{Int}      # Mapping from nzval to indices
    ind2nzval::Array{Int, N}    # Mapping from indices to nzval
    nnz::Int                    # number of stored zeros
end

SparseIterate(n::Int) = SparseIterate(Float64, n)
SparseIterate(::Type{T}, n::Int) where {T} = SparseIterate{T, 1}(zeros(T, n), zeros(Int, n), zeros(Int, n), 0)
SparseIterate(n::Int, m::Int) = SparseIterate(Float64, n, m)
SparseIterate(::Type{T}, n::Int, m::Int) where {T} = SparseIterate{T, 2}(zeros(T, n*m), zeros(Int, n*m), zeros(Int, n, m), 0)

SparseArrays.issparse(x::SparseIterate) = true

Base.length(x::SparseIterate) = length(x.ind2nzval)
Base.size(x::SparseIterate) = size(x.ind2nzval)
SparseArrays.nnz(x::SparseIterate) = x.nnz
numCoordinates(x::SparseIterate) = length(x.ind2nzval)

## similar
#
Base.similar(A::SparseIterate{T}) where {T} = SparseIterate(T, size(A)...)
Base.similar(A::SparseIterate, ::Type{S}) where {S} = SparseIterate(S, size(A)...)

### Conversion
function SparseIterate(x::SparseMatrixCSC{T}) where {T}
  n, m = size(x)
  out = SparseIterate(T, n, m)

  I, J, V = findnz(x)
  for i=1:length(I)
    @inbounds out[I[i], J[i]] = V[i]
  end

  out
end

function SparseIterate(x::SparseVector{T}) where {T}
  p = length(x)
  out = SparseIterate(T, p)
  nzval = SparseArrays.nonzeros(x)
  rowval = SparseArrays.nonzeroinds(x)
  for i=1:length(nzval)
    out[rowval[i]] = nzval[i]
  end
  out
end

function SparseIterate(x::Array{T}) where {T}
  out = SparseIterate(T, size(x)...)

  @inbounds for i in eachindex(x)
      out[i] = x[i]
  end

  out
end

function Base.Vector(x::SparseIterate{Tv, 1}) where Tv
    n = length(x)
    n == 0 && return Vector{Tv}(0)
    r = zeros(Tv, n)
    for k in 1:nnz(x)
        i = x.nzval2ind[k]
        v = x.nzval[k]
        r[i] = v
    end
    return r
end

function Base.Matrix(x::SparseIterate{Tv, 2}) where Tv
  n, m = size(x)
  A = zeros(Tv, n, m)
  for k=1:nnz(x)
    i = x.nzval2ind[k]
    v = x.nzval[k]
    A[i] = v
  end
  return A
end

# Base.convert(::Type{Array}, x::SparseIterate{T, 1}) where {T} = convert(Vector, x)
# Base.convert(::Type{Array}, x::SparseIterate{T, 2}) where {T} = convert(Matrix, x)
# Base.full(x::SparseIterate) = convert(Array, x)

Base.IndexStyle(::Type{<:SparseIterate}) = IndexLinear()

Base.getindex(x::SparseIterate{T}, i::Int) where {T} = (checkbounds(x, i); x.ind2nzval[i] == 0 ? zero(T) : x.nzval[x.ind2nzval[i]])
function Base.setindex!(x::SparseIterate{T}, v::T, i::Int) where {T}
  checkbounds(x, i)
  if x.ind2nzval[i] == 0
    if v != zero(T)
      x.nnz += 1
      x.nzval[x.nnz] = v
      x.nzval2ind[x.nnz] = i
      x.ind2nzval[i] = x.nnz
    end
  else
    x.nzval[x.ind2nzval[i]] = v
  end
  x
end

function Base.copyto!(x::SparseIterate, y::SparseIterate)
    size(x) == size(y) || throw(DimensionMismatch())
    copyto!(x.nzval, y.nzval)
    copyto!(x.nzval2ind, y.nzval2ind)
    copyto!(x.ind2nzval, y.ind2nzval)
    x.nnz = y.nnz
    x
end

### show and friends

Base.show(io::IO, x::SparseIterate) = show(convert(IOContext, io), x)

function Base.show(io::IO, ::MIME"text/plain", x::SparseIterate{T, 1}) where T
    xnnz = nnz(x)
    print(io, length(x), "-element ", typeof(x), " with ", xnnz,
           " stored ", xnnz == 1 ? "entry" : "entries")
    if xnnz != 0
        println(io, ":")
        show(io, x)
    end
end

function Base.show(io::IOContext, x::SparseIterate{T, 1}) where T
    n = length(x)
    xnnz = nnz(x)
    if xnnz == 0
        return show(io, MIME("text/plain"), x)
    end
    limit::Bool = get(io, :limit, false)
    half_screen_rows = limit ? div(displaysize(io)[1] - 8, 2) : typemax(Int)
    pad = ndigits(n)
    if !haskey(io, :compact)
        io = IOContext(io, :compact => true)
    end
    for k = 1:n
        if k < half_screen_rows || k > xnnz - half_screen_rows
          if x.ind2nzval[k] != 0
            print(io, "  ", '[', rpad(k, pad), "]  =  ")
            show(io, x[k])
            k != n && println(io)
          end
        elseif k == half_screen_rows
            println(io, "   ", " "^pad, "   \u22ee")
        end
    end
end

function Base.show(io::IO, ::MIME"text/plain", S::SparseIterate{T,2}) where T
    xnnz = nnz(S)
    print(io, size(S, 1), "×", size(S, 2), " ", typeof(S), " with ", xnnz, " stored ",
              xnnz == 1 ? "entry" : "entries")
    if xnnz != 0
        print(io, ":")
        show(io, S)
    end
end

function Base.show(io::IOContext, S::SparseIterate{T, 2}) where T
  if nnz(S) == 0
    return show(io, MIME("text/plain"), S)
  end

  limit::Bool = get(io, :limit, false)
  if limit
    rows = displaysize(io)[1]
    half_screen_rows = div(rows - 8, 2)
  else
    half_screen_rows = typemax(Int)
  end
  pad = ndigits(maximum( size(S) ))
  sep = "\n  "
  if !haskey(io, :compact)
    io = IOContext(io, :compact => true)
  end
  ind = 0
  _ci = CartesianIndices(size(S))
  for ifull in eachindex(S)
    if S.ind2nzval[ifull] != 0
      ind += 1
      if ind < half_screen_rows || ind > nnz(S)-half_screen_rows
        row, col = Tuple(_ci[ifull])
        print(io, sep, '[', rpad(row, pad), ", ", lpad(col, pad), "]  =  ")
        if isassigned(S.nzval, S.ind2nzval[ifull])
          show(io, S.nzval[S.ind2nzval[ifull]])
        else
          print(io, Base.undef_ref_str)
        end
      elseif ind == half_screen_rows
        print(io, sep, '\u22ee')
      end
    end
  end
end

##

function SparseArrays.dropzeros!(x::SparseIterate{T}) where T
  i = 1
  while i <= x.nnz
    if x.nzval[i] == zero(T)
      x.ind2nzval[x.nzval2ind[i]] = 0
      if i != x.nnz
        x.nzval[i] = x.nzval[x.nnz]
        x.ind2nzval[x.nzval2ind[x.nnz]] = i
        x.nzval2ind[i] = x.nzval2ind[x.nnz]
      end
      x.nnz -= 1
      i -= 1
    end
    i += 1
  end
  x
end

## multiplication

function *(A::StridedMatrix{Ta}, x::SparseIterate{Tx, 1}) where {Ta,Tx}
    m, n = size(A)
    length(x) == n || throw(DimensionMismatch())
    Ty = promote_type(Ta, Tx)
    y = Vector{Ty}(undef, m)
    mul!(y, A, x)
end

function LinearAlgebra.mul!(out::Vector, A::StridedMatrix, coef::SparseIterate{T, 1}) where {T}
  length(coef) == size(A, 2) || throw(DimensionMismatch("A has second dimension $(size(A,2)), length of coef is $(length(coef))"))
  length(out) == size(A, 1) || throw(DimensionMismatch("A has first dimension $(size(A,1)), length of out is $(length(out))"))

  fill!(out, zero(eltype(out)))
  @inbounds for inz = 1:nnz(coef)
    indColumn = coef.nzval2ind[inz]
    c = coef.nzval[inz]
    @simd for i = 1:size(A, 1)
      out[i] += c*A[i, indColumn]
    end
  end
  out
end

function *(transA::Transpose{<:Any,<:StridedMatrix{Ta}}, x::SparseIterate{Tx, 1}) where {Ta,Tx}
    A = transA.parent
    m, n = size(A)
    length(x) == m || throw(DimensionMismatch())
    Ty = promote_type(Ta, Tx)
    y = Vector{Ty}(undef, n)
    mul!(y, transpose(A), x)
end

function LinearAlgebra.mul!(out::AbstractVector, transA::Transpose{<:Any,<:StridedMatrix}, coef::SparseIterate{T, 1}) where {T}
    A = transA.parent
    length(coef) == size(A, 1) || throw(DimensionMismatch("A has first dimension $(size(A,1)), length of coef is $(length(coef))"))
    length(out) == size(A, 2) || throw(DimensionMismatch("A has second dimension $(size(A,2)), length of out is $(length(out))"))

    fill!(out, zero(eltype(out)))
    @inbounds for j = 1:size(A, 2)
        s = zero(eltype(out))
        for inz = 1:nnz(coef)
            s += A[coef.nzval2ind[inz], j] * coef.nzval[inz]
        end
        out[j] = s
    end
    return out
end

LinearAlgebra.dot(coef::SparseIterate{T, 1}, x::Vector{T}) where {T} = dot(x, coef)

function LinearAlgebra.dot(x::Vector{T}, coef::SparseIterate{T, 1}) where {T}
    v = zero(T)
    @inbounds @simd for icoef = 1:nnz(coef)
        v += x[coef.nzval2ind[icoef]]*coef.nzval[icoef]
    end
    v
end

####################################
#
#  Symmetric Sparse Iterate
#
####################################

# data are stored in the lower triangle
mutable struct SymmetricSparseIterate{T} <: AbstractMatrix{T}
    nzval::Vector{T}            # stored nonzero values (1:nnz)
    nzval2ind::Vector{Int}      # Mapping from nzval to indices
    ind2nzval::Vector{Int}      # Mapping from indices to nzval
    nnz::Int                    # number of stored zeros
    p::Int                      # dimension of matrix
end

SymmetricSparseIterate(n::Int) = SymmetricSparseIterate(Float64, n)
SymmetricSparseIterate(::Type{T}, n::Int) where {T} =
  (_t = div(n*(n+1), 2); SymmetricSparseIterate{T}(zeros(T, _t), zeros(Int, _t), zeros(Int, _t), 0, n))

Base.length(x::SymmetricSparseIterate) = div(x.p*(x.p+1), 2)
Base.size(x::SymmetricSparseIterate) = (x.p, x.p)
SparseArrays.nnz(x::SymmetricSparseIterate) = x.nnz
function numCoordinates(x::SymmetricSparseIterate)
  p = x.p
  div(p*(p+1), 2)
end

Base.IndexStyle(::Type{<:SymmetricSparseIterate}) = IndexCartesian()
function Base.getindex(x::SymmetricSparseIterate, r::Int, c::Int)
  linInd = r >= c ? sub2indLowerTriangular(x.p, r, c) : sub2indLowerTriangular(x.p, c, r)
  getindex(x, linInd)
end
Base.getindex(x::SymmetricSparseIterate{T}, i::Int) where {T} = x.ind2nzval[i] == 0 ? zero(T) : x.nzval[x.ind2nzval[i]]
function Base.setindex!(x::SymmetricSparseIterate{T}, v::T, r::Int, c::Int) where {T}
  linInd = r >= c ? sub2indLowerTriangular(x.p, r, c) : sub2indLowerTriangular(x.p, c, r)
  setindex!(x, v, linInd)
end
function Base.setindex!(x::SymmetricSparseIterate{T}, v::T, i::Int) where {T}
  if x.ind2nzval[i] == 0
    if v != zero(T)
      x.nnz += 1
      x.nzval[x.nnz] = v
      x.nzval2ind[x.nnz] = i
      x.ind2nzval[i] = x.nnz
    end
  else
    x.nzval[x.ind2nzval[i]] = v
  end
  x
end

function SymmetricSparseIterate(X::SparseMatrixCSC{T}) where {T}
  n, m = size(X)
  n == m || throw(ArgumentError("X needs to be square matrix"))
  out = SymmetricSparseIterate(T, n)

  I, J, V = findnz(X)
  @inbounds for i=1:length(I)
    if I[i] >= J[i]
      out[I[i], J[i]] = V[i]
    end
  end

  out
end

function Base.Matrix(S::SymmetricSparseIterate{Tv}) where Tv
  p = S.p
  A = zeros(Tv, p, p)
  for k=1:nnz(S)
    i = S.nzval2ind[k]
    v = S.nzval[k]
    r, c = ind2subLowerTriangular(p, i)
    A[r, c] = v
    if r != c
      A[c, r] = v
    end
  end
  return A
end

Base.similar(A::SymmetricSparseIterate{T}) where {T} = SymmetricSparseIterate(T, size(A, 1))
Base.similar(A::SymmetricSparseIterate, ::Type{S}) where {S} = SymmetricSparseIterate(S, size(A, 1))

function Base.copyto!(x::SymmetricSparseIterate, y::SymmetricSparseIterate)
    size(x) == size(y) || throw(DimensionMismatch())
    copyto!(x.nzval, y.nzval)
    copyto!(x.nzval2ind, y.nzval2ind)
    copyto!(x.ind2nzval, y.ind2nzval)
    x.nnz = y.nnz
    x.p = y.p
    x
end

function SparseArrays.dropzeros!(x::SymmetricSparseIterate{T}) where T
    i = 1
    while i <= x.nnz
        if x.nzval[i] == zero(T)
            x.ind2nzval[x.nzval2ind[i]] = 0
            if i != x.nnz
                x.nzval[i] = x.nzval[x.nnz]
                x.ind2nzval[x.nzval2ind[x.nnz]] = i
                x.nzval2ind[i] = x.nzval2ind[x.nnz]
            end
            x.nnz -= 1
            i -= 1
        end
        i += 1
    end
    x
end

function LinearAlgebra.mul!(C::AbstractMatrix{T}, A::SymmetricSparseIterate{T}, B::AbstractMatrix{T}) where T
  size(A, 2) == size(B, 1) || throw(DimensionMismatch())
  size(A, 1) == size(C, 1) || throw(DimensionMismatch())
  size(B, 2) == size(C, 2) || throw(DimensionMismatch())

  fill!(C, zero(T))
  i = 1
  @inbounds while i <= A.nnz
    v = A.nzval[i]
    ind = A.nzval2ind[i]
    r, c = ind2subLowerTriangular(size(A, 1), ind)
    if r == c
      for k=1:size(C, 2)
        @inbounds C[r, k] += B[c, k] * v
      end
    else
      for k=1:size(C, 2)
        @inbounds C[r, k] += B[c, k] * v
        @inbounds C[c, k] += B[r, k] * v
      end
    end
    i += 1
  end
  C
end

function LinearAlgebra.mul!(C::AbstractVector{T}, A::SymmetricSparseIterate{T}, B::AbstractVector{T}) where T
  size(A, 2) == size(B, 1) || throw(DimensionMismatch())
  size(A, 1) == size(C, 1) || throw(DimensionMismatch())

  fill!(C, zero(T))
  i = 1
  @inbounds while i <= A.nnz
    v = A.nzval[i]
    ind = A.nzval2ind[i]
    r, c = ind2subLowerTriangular(size(A, 1), ind)
    if r == c
        @inbounds C[r] += B[c] * v
    else
        @inbounds C[r] += B[c] * v
        @inbounds C[c] += B[r] * v
    end
    i += 1
  end
  C
end

####################################
#
#  Group Sparse Iterate
#
####################################

struct AtomIterate{T, N, V} <: AbstractArray{T, N}
  storage::Array{T, N}
  atoms::Vector{V}                                 # each atom is a view of the storage
end

function AtomIterate(n::Int, numAtom::Int)
  mod(n, numAtom) == 0 || throw(ArgumentError("numAtom must divide n"))
  storage = zeros(n)
  lenAtom = div(n, numAtom)
  b, e = 1, lenAtom
  a1 = view(storage, b:e)
  V = typeof(a1)
  atoms = Vector{V}(undef, numAtom)
  atoms[1] = a1
  for i=2:numAtom
    b += lenAtom
    e += lenAtom
    atoms[i] = view(storage, b:e)
  end
  AtomIterate{Float64, 1, V}(storage, atoms)
end

# if byRow == true then each atom is a group of the form [b:e, :]
# otherwise each atom is a group of the form [:, b:e]
function AtomIterate(dims::Dims{2}, numAtom::Int, byRow::Bool)
  nr = dims[1]
  nc = dims[2]
  if byRow
    mod(nr, numAtom) == 0 || throw(ArgumentError("numAtom must divide dims[1]"))
    lenAtom = div(nr, numAtom)
  else
    mod(nc, numAtom) == 0 || throw(ArgumentError("numAtom must divide dims[1]"))
    lenAtom = div(nc, numAtom)
  end
  storage = zeros(nr, nc)
  b, e = 1, lenAtom
  if byRow
    a1 = view(storage, b:e, :)
  else
    a1 = view(storage, :, b:e)
  end
  V = typeof(a1)
  atoms = Vector{V}(undef, numAtom)
  atoms[1] = a1
  for i=2:numAtom
    b += lenAtom
    e += lenAtom
    if byRow
      atoms[i] = view(storage, b:e, :)
    else
      atoms[i] = view(storage, :, b:e)
    end
  end
  AtomIterate{Float64, 2, V}(storage, atoms)
end

#

function Base.similar(a::AtomIterate{T, N, V}) where {T,N,V}
  storage = zero(a.storage)
  numAtoms = length(a.atoms)
  atoms = Vector{V}(undef, numAtoms)
  for i=1:numAtoms
    atoms[i] = SubArray(storage, a.atoms[i].indices)
  end
  AtomIterate{T,N,V}(storage, atoms)
end

Base.copy(a::AtomIterate) = copyto!(similar(a), a)
function Base.copyto!(dest::AtomIterate, src::AtomIterate)
  copyto!(dest.storage, src.storage)
  dest
end

#

Base.length(x::AtomIterate) = length(x.storage)
Base.size(x::AtomIterate) = size(x.storage)
Base.IndexStyle(::Type{<:AtomIterate}) = IndexLinear()

Base.getindex(x::AtomIterate, i::Int) = getindex(x.storage, i)
Base.getindex(x::AtomIterate, I...) = getindex(x.storage, I...)
Base.setindex!(x::AtomIterate, v, i::Int) = setindex!(x.storage, v, i)
Base.setindex!(x::AtomIterate, v, I...) = setindex!(x.storage, v, I...)


text:
using Rainsimulator
using Test

@testset "Rainsimulator.jl" begin
    # Write your tests here.
end


text:
"""
Module with graph theoretical methods for maximum matching, Pantelides index reduction and Tarjan's strongly connected components.

* Author: Hilding Elmqvist, Mogram AB
* Date: July-August 2016 (rewritten).
* License: MIT

A bipartite graph (E, V) is defined by an array of arrays. Each E-entry has an integer array of indices to V-entries.

Example bipartite graph:

    G = [
      [3, 5],
      [4, 6],
      [1, 7, 9],
      [2, 8, 9],
      [1, 2]
    ]
"""
module BLTandPantelides

export matching, pantelides!, BLT, checkAssign

using ..BLTandPantelidesUtilities

"Controls logging"
const log = false

"""
    function augmentPath!(G, i, assign, vColour, eColour, vPassive)
Construction of augmenting path

Reference:
Pantelides, C.: The consistent initialization of differential-algebraic systems. SIAM Journal
of Scientific and Statistical Computing, 9(2), pp. 213–231 (1988).
"""
function augmentPath!(G, i, assign, vColour, eColour, vPassive)
    # returns pathFound
    # assign: assign[j] contains the E-node to which V-node j is assigned or 0 if V-node j not assigned
    # i: E-node
    # vPassive: set to != 0 has the same effect as deleting V-node and corresponding edges
    # j: V-node

    if log
        println("augmentPath: equation $i")
    end

    pathFound = false
    eColour[i] = true

    # If a V-node j exists such that edge (i-j) exists and assign[j] == 0
    for j in G[i]
        if vPassive[j] == 0 && assign[j] == 0
            pathFound = true
            assign[j] = i
            return pathFound
        end
    end

    # For every j such that edge (i-j) exists and j is uncoloured
    for j in G[i]
        if vPassive[j] == 0 && !vColour[j]
            vColour[j] = true
            k = assign[j]
            pathFound = augmentPath!(G, k, assign, vColour, eColour, vPassive)

            if pathFound
                assign[j] = i
                return pathFound
            end
        end
    end
    return pathFound
end

function checkAssign(assign, VSizes, VTypes, ESizes, ETypes, equationsInfix, variableNames, A, vPassive=A)
    println("Checking assignment")
    assignmentOK = true
    for j in 1:length(assign)
        if vPassive[j] == 0
            i = assign[j]
            if i > 0 && VSizes[j] != ESizes[i]
                assignmentOK = false
                print("Error: Variable ")
                printList(variableNames, [j], A, newLine=false)
                println(" (($j)) with size=$(VSizes[j]) is assigned in equation (($i)) with size $(ESizes[i])")
            end
        end
    end

    if assignmentOK
        println("Assignment is OK")
    else
        # error("Assignment not OK")
  end
end

"""
    function matching(G, M, vActive=fill(true, M))
Find maximum matching in bipartite graph

* `G`: bipartite graph
* `M`: number of V-nodes
* `vActive`: set to false has the same effect as deleting V-node and corresponding edges
* `return assign`: assign[j] contains the E-node to which V-node j is assigned or 0 if V-node j not assigned

Reference:
Pantelides, C.: The consistent initialization of differential-algebraic systems. SIAM Journal
of Scientific and Statistical Computing, 9(2), pp. 213–231 (1988).
"""
function matching(G, M, vActive=fill(true, M))
    assign::Array{Int,1} = fill(0, M)
    eColour::Array{Bool,1} = fill(false, length(G))
    vColour::Array{Bool,1} = fill(false, M)
    vPassive::Array{Int,1} = [if va; 0 else 1 end for va in vActive]
    for i in 1:length(G)
        fill!(eColour, false)
        fill!(vColour, false)
        pathFound = augmentPath!(G, i, assign, vColour, eColour, vPassive)
    end
    return assign
end

# -------------------------------------------------------

"""
    function pantelides!(G, M, A)
Perform index reduction with Pantelides algorithm.

* `G`: bipartite graph (updated)
* `M`: number of V-nodes
* `A`: A[j] = if V[k] = der(V[j]) then k else 0
* `return assign`: assign[j] contains the E-node to which V-node j is assigned or 0 if V-node j not assigned
* `return A`: A[j] = if V[k] = der(V[j]) then k else 0
* `return B`: B[i] = if E[l] = der(E[i]) then l else 0

Reference:
Pantelides, C.: The consistent initialization of differential-algebraic systems. SIAM Journal
of Scientific and Statistical Computing, 9(2), pp. 213–231 (1988).
"""
function pantelides!(G, M, A)
    assign::Array{Int,1} = fill(0, M)
    B::Array{Int,1} = fill(0, length(G))
    eColour::Array{Bool,1} = fill(false, length(G))
    vColour::Array{Bool,1} = fill(false, M)
    N = length(G)
    N2 = N
    for k in 1:N2
        pathFound = false
        i = k
        while !pathFound
            # Delete all V-nodes with A[.] != 0 and all their incidence edges from the graph
            # Designate all nodes as "uncoloured"
            if length(eColour) == length(G)
                fill!(eColour, false)
            else
                eColour = fill(false, length(G))
            end
            if length(vColour) == length(M)
                fill!(vColour, false)
            else
                vColour = fill(false, M)
            end

            pathFound = augmentPath!(G, i, assign, vColour, eColour, A)
            if !pathFound
                if log
                    println("\nDifferentiate:")
                end

                # For every coloured V-node j do
                for j in 1:length(vColour)
                    if vColour[j]
                        M += 1
                        if log
                            println("New variable derivative: var($M) = der($j)")
                        end
                        push!(A, 0)
                        A[j] = M
                        push!(assign, 0)
                    end
                end

                # For every coloured E-node l do
                for l in 1:N
                    if eColour[l]
                        N += 1
                        if log
                            println("New equation derivative:  equ($N) = DER($l)")
                        end

                        # Create new E-node N
                        push!(G, copy(G[l]))

                        # Create edges from E-node N to all V-nodes j and A[j] such that edge (l-j) exists
                        for m in 1:length(G[l])
                            j = G[l][m]
                            if !(A[j] in G[N])
                                push!(G[N], A[j])
                            end
                        end
                        push!(B, 0)

                        # Set B[l] = N
                        B[l] = N
                    end
                end

                # For every coloured V-node j
                for j in 1:length(vColour)
                    if vColour[j]
                        if log
                            println("Assigning derivative of variable $(A[j]) to derivative of equation: $(B[assign[j]])")
                        end
                        assign[A[j]] = B[assign[j]]
                    end
                end
                i = B[i]
            end
        end
    end

    return assign, A, B
end

const notOnStack = 1000000000

"""
Find minimal systems of equations that have to be solved simultaneously.
Reference:
Tarjan, R. E. (1972), "Depth-first search and linear graph algorithms", SIAM Journal on Computing 1 (2): 146–160, doi:10.1137/0201010
"""
function strongConnect!(G, assign, v, nextnode, stack, components, lowlink, number)
    # println("strongConnect: ", v)

    if v == 0
        return nextnode
    end

    nextnode += 1
    lowlink[v] = number[v] = nextnode
    push!(stack, v)

    for w in [assign[j] for j in G[v]] # for w in the adjacency list of v
        if w > 0   # Is assigned
            if number[w] == 0 # if not yet numbered
                nextnode = strongConnect!(G, assign, w, nextnode, stack, components, lowlink, number)
                lowlink[v] = min(lowlink[v], lowlink[w])
            else
                if number[w] < number[v]
                    # (v, w) is a frond or cross-link
                    # if w is on the stack of points. Always valid since otherwise number[w]=notOnStack (a big number)
                    lowlink[v] = min(lowlink[v], number[w])
                end
            end
        end
    end

    if lowlink[v] == number[v]
        # v is the root of a component
        # start a new strongly connected component
        comp = []
        repeat = true
        while repeat
            # delete w from point stack and put w in the current component
            # println("delete w from point stack and put w in the current component")
            w = pop!(stack)
            number[w] = notOnStack
            push!(comp, w)
            repeat = w != v
        end
        push!(components, comp)
    end
    return nextnode
end

"""
    function BLT(G, assign)

Find Block Lower Triangular structure for a bipartite graph `G` with assignment `assign`

* `G`: bipartite graph
* `assign`: assign[j] contains the E-node to which V-node j is assigned or 0 if V-node j not assigned
* `return components`: cell array of components. Each component is a list of indices to E-nodes
"""
function BLT(G, assign)
    nextnode::Int = 0
    stack = []
    components = []
    lowlink = fill(0, length(G))
    number = fill(0, length(G))

    for v in 1:length(G)
        if number[v] == 0
            nextnode = strongConnect!(G, assign, v, nextnode, stack, components, lowlink, number)
        end
    end
    return components
end

end


text:
abstract type AbstractMachine end

function filtrations end

function (m::AbstractMachine)(y₀)
    forward, backward = filtrations(m, y₀)
    return solve(y₀, m.W, m.σ, (forward, backward))
end

sum_dims(dims::Tuple) = prod(dims[1:end-2]) * sum(dims[end-1:end])

# glorot initialization, from https://github.com/FluxML/Flux.jl
glorot_normal(dims...) = randn(Float32, dims...) .* sqrt(2.0f0 / sum_dims(dims))
glorot_uniform(dims...) = (rand(Float32, dims...) .- 0.5f0) .* sqrt(24.0f0 / sum_dims(dims))

function consecutive_ranges(dims::Vector{Int})
    return map(cumsum(dims), dims) do cumdim, dim
        return (cumdim - dim + 1):cumdim
    end
end

function dense_filtrations(y, dims::Vector{Int})
    channels, minibatch = sum(dims), size(y, 2)
    shapes = map(rg -> (rg, 1:minibatch), consecutive_ranges(dims))
    input(_) = (1:channels, 1:minibatch)
    function project(input_shape, output_shape)
        return (first(input_shape), first(output_shape))
    end
    return (
        forward=Filtration(shapes, input, project),
        backward=Filtration(reverse(shapes), input, project)
    )
end

"""
    DenseMachine(W, σ, dims::Vector{Int})

Create a `DenseMachine` object from a square weight matrix `W`, pointwise nonlinearity `σ`
and filtration sequence `dims`.
The values of `dims` specify how to split the input space into a sequence of subspaces.
In particular, it is required that `size(W) == (sum(dims), sum(dims))`.
"""
struct DenseMachine{T, F} <: AbstractMachine
    W::T
    σ::F
    dims::Vector{Int}
end

@functor DenseMachine (W, σ)

DenseMachine(σ, dims::Vector{Int}; init=glorot_uniform) = DenseMachine(dims, σ; init)

"""
    DenseMachine(dims::Vector{Int}, σ; init=glorot_uniform)

Return a `DenseMachine(W, σ, dims)` object, where `W = init(sum(dims), sum(dims))`.
Default to Glorot uniform initialization.
"""
function DenseMachine(dims::Vector{Int}, σ; init=glorot_uniform)
    W = init(sum(dims), sum(dims))
    return DenseMachine(W, σ, dims)
end

filtrations(m::DenseMachine, y₀) = dense_filtrations(y₀, m.dims)

nest_pads(::Tuple{}) = ()

function nest_pads(pad::Tuple)
    a, b, tail... = pad
    return ((a, b), nest_pads(tail)...)
end

function conv_range(output_range, (pad0, pad1)::Tuple{Int,Int})
    i0, i1 = first(output_range), last(output_range)
    return (i0 - pad0):(i1 + pad1)
end

function conv_filtrations(y::AbstractArray{T, N}, dims::Vector{Int}; pad::NTuple{M, Int}) where {T, N, M}
    M == 2 * (N - 2) || throw(ArgumentError("padding must have twice the length of image dimensions"))
    pads = nest_pads(pad)

    ranges = consecutive_ranges(dims)
    channels, minibatch = sum(dims), size(y, N)
    shapes = [ntuple(n -> ifelse(n == N - 1, rg, UnitRange(axes(y, n))), N) for rg in ranges]

    function input(shape)
        data = ntuple(n -> shape[n], N - 2)
        zdata = map(conv_range, data, pads)
        return (zdata..., 1:channels, 1:minibatch)
    end

    function input′(shape)
        data = ntuple(n -> shape[n], N - 2)
        zdata = map(conv_range, data, map(reverse, pads))
        return (zdata..., 1:channels, 1:minibatch)
    end

    function project(input_shape, output_shape)
        return ntuple(N) do n
            n == N && return output_shape[end - 1]
            n == N - 1 && return input_shape[end - 1]
            p0, p1 = pads[n]
            return 1:(p0 + p1 + 1)
        end
    end

    return (
        forward=Filtration(shapes, input, project),
        backward=Filtration(reverse(shapes), input′, project)
    )
end

"""
    ConvMachine(W, σ, dims::Vector{Int}, pad::Dims)

Create a `ConvMachine` object from a weight array `W`, pointwise nonlinearity `σ`,
filtration sequence `dims` and padding `pad`.
The values of `dims` specify how to split the input space into a sequence of subspaces.
In particular, it is required that `size(W) == (kernelsize..., sum(dims), sum(dims))`,
where `kernelsize` is such that convolution by `W` with padding `pad` preserves
input dimension.
Padding `pad` has length twice the number of kernel dimensions (for example, it takes
four values for image convolutions).
"""
struct ConvMachine{T, F, M} <: AbstractMachine
    W::T
    σ::F
    dims::Vector{Int}
    pad::NTuple{M, Int}
end

@functor ConvMachine (W, σ)

ConvMachine(W, σ, dims; pad) = ConvMachine(W, σ, dims, pad)

ConvMachine(σ, dims::Vector{Int}; pad, init=glorot_uniform) = ConvMachine(dims, σ; pad, init)

"""
    ConvMachine(dims::Vector{Int}, σ; pad, init=glorot_uniform)

Return a `ConvMachine(W, σ, dims, pad)` object, where `W = init(kernelsize..., sum(dims), sum(dims))`.
Default to Glorot uniform initialization.
Here `kernelsize` is such that convolution with a kernel of size `kernelsize`
and padding `pad` preserves input dimension.
Padding `pad` has length twice the number of kernel dimensions (for example, it takes
four values for image convolutions).
"""
function ConvMachine(dims::Vector{Int}, σ; pad, init=glorot_uniform)
    kernelsize = map(((p0, p1),) -> p0 + p1 + 1, nest_pads(pad))
    W = init(kernelsize..., sum(dims), sum(dims))
    return ConvMachine(W, σ, dims; pad)
end

filtrations(m::ConvMachine, y₀) = conv_filtrations(y₀, m.dims; m.pad)

function recur_filtrations(y, dims::Vector{Int}; pad, timeblock)
    ranges = consecutive_ranges(dims)
    datalength, channels, minibatch = size(y, 1), sum(dims), size(y, 3)
    N = ceil(Int, datalength / timeblock)
    shapes = [(timeblock  * (n  - 1) + 1 : timeblock * n, range, 1:minibatch) for n in 1:N for range in ranges]

    input((timerange, _, mb)) =  (conv_range(timerange, (pad, 0)), 1:channels, mb)
    input′((timerange, _, mb)) =  (conv_range(timerange, (0, pad)), 1:channels, mb)
    function project(input_shape, output_shape)
        input_range, output_range = first(input_shape), first(output_shape)
        w0 = 1 + clamp(first(output_range) - last(input_range), 0:pad)
        w1 = 1 + clamp(last(output_range) - first(input_range), 0:pad)
        return (w0:w1, input_shape[end-1], output_shape[end-1])
    end
    flip = (pad+1):-1:1
    function project′(input_shape, output_shape)
        input_range, output_range = first(input_shape), first(output_shape)
        w0 = 1 + clamp(first(input_range) - last(output_range), 0:pad)
        w1 = 1 + clamp(last(input_range) - first(output_range), 0:pad)
        return (flip[w1]:flip[w0], input_shape[end-1], output_shape[end-1])
    end
    return (
        forward=Filtration(shapes, input, project),
        backward=Filtration(reverse(shapes), input′, project′)
    )
end

"""
    RecurMachine(W, σ, dims::Vector{Int}, pad::Int, timeblock::Int)

Create a `RecurMachine` object from a weight array `W`, pointwise nonlinearity `σ`,
filtration sequence `dims`, padding `pad` and time block `timeblock`.
The values of `dims` specify how to split the input space into a sequence of subspaces.
In particular, it is required that `size(W) == (pad + 1, sum(dims), sum(dims))`.
"""
struct RecurMachine{T, F} <: AbstractMachine
    W::T
    σ::F
    dims::Vector{Int}
    pad::Int
    timeblock::Int
end

@functor RecurMachine (W, σ)

RecurMachine(W, σ, dims; pad, timeblock) = RecurMachine(W, σ, dims, pad, timeblock)

RecurMachine(σ, dims::Vector{Int}; pad, timeblock, init=glorot_uniform) = RecurMachine(dims, σ; pad, timeblock, init)

"""
    RecurMachine(dims::Vector{Int}, σ; pad, timeblock, init=glorot_uniform)

Return a `RecurMachine(W, σ, dims, pad, timeblock)` object, where `W = init(pad + 1, sum(dims), sum(dims))`.
Default to Glorot uniform initialization.
"""
function RecurMachine(dims::Vector{Int}, σ; pad, timeblock, init=glorot_uniform)
    W = init(pad + 1, sum(dims), sum(dims))
    return RecurMachine(W, σ, dims; pad, timeblock)
end

filtrations(m::RecurMachine, y₀) = recur_filtrations(y₀, m.dims; m.pad, m.timeblock)

function ChainRulesCore.rrule(::typeof(filtrations), m::AbstractMachine, y)
    res = filtrations(m, y)
    function filtrations_pullback(_, _)
        return NoTangent(), NoTangent(), NoTangent()
    end
    return res, filtrations_pullback
end


text:
using JET, InteractiveUtils

const CC = JET.CC

import .CC:
    widenconst, ⊑

import JET:
    JETInterpreter,
    VirtualGlobalVariable,
    profile_call,
    get_result,
    virtual_process!,
    gen_virtual_module,
    ToplevelErrorReport,
    InferenceErrorReport,
    ExceptionReport,
    print_reports

for sym in Symbol.(last.(Base.Fix2(split, '.').(string.(vcat(subtypes(JET, ToplevelErrorReport),
                                                             subtypes(JET, InferenceErrorReport),
                                                             subtypes(JET, ExceptionReport),
                                                             )))))
    Core.eval(@__MODULE__, :(import JET: $(sym)))
end

import Base:
    Fix1,
    Fix2

import Base.Meta:
    isexpr

# define virtual module and setup fixtures
macro def(ex)
    @assert isexpr(ex, :block)
    return quote let
        vmod = $(gen_virtual_module)()
        Core.eval(vmod, $(QuoteNode(ex)))
        vmod # return virtual module
    end end
end

# profile from file name
profile_file′(filename, args...; kwargs...) =
    return profile_text′(read(filename, String), args...; filename, kwargs...)

# profile from string
function profile_text′(s,
                       virtualmod = gen_virtual_module(@__MODULE__);
                       filename = "top-level",
                       actualmodsym = Symbol(parentmodule(virtualmod)),
                       interp = JETInterpreter(),
                       )
    return virtual_process!(s,
                            filename,
                            virtualmod,
                            actualmodsym,
                            interp)
end

# profile from toplevel expression
macro profile_toplevel(virtualmod, ex)
    return profile_toplevel(virtualmod, ex, __source__)
end

macro profile_toplevel(ex)
    virtualmod = gen_virtual_module(__module__)
    return profile_toplevel(virtualmod, ex, __source__)
end

function profile_toplevel(virtualmod, ex, lnn)
    toplevelex = (isexpr(ex, :block) ?
                  Expr(:toplevel, lnn, ex.args...) : # flatten here
                  Expr(:toplevel, lnn, ex)
                  ) |> QuoteNode
    return quote let
        interp = $(JETInterpreter)()
        ret = $(JET.gen_virtual_process_result)()
        virtualmod = $(esc(virtualmod))
        actualmodsym = Symbol(parentmodule(virtualmod))
        $(virtual_process!)($(toplevelex), $(string(lnn.file)), virtualmod, actualmodsym, interp, ret)
    end end
end

is_concrete(mod, sym) = isdefined(mod, sym) && !isa(getfield(mod, sym), VirtualGlobalVariable)
is_abstract(mod, sym) = isdefined(mod, sym) && isa(getfield(mod, sym), VirtualGlobalVariable)
isa_abstract(vgv, @nospecialize(typ)) = isa(vgv, VirtualGlobalVariable) && vgv.t ⊑ typ


text:
"""
# References
- Controller (it may be modified in this implementation)
[1] G. P. Falconi and F. Holzapfel,
“Adaptive Fault Tolerant Control Allocation for a Hexacopter System,”
Proc. Am. Control Conf., vol. 2016-July, pp. 6760–6766, 2016.
- Reference model (e.g., xd, vd, ad, ad_dot, ad_ddot)
[2] S. J. Su, Y. Y. Zhu, H. R. Wang, and C. Yun, “A Method to Construct a Reference Model for Model Reference Adaptive Control,” Adv. Mech. Eng., vol. 11, no. 11, pp. 1–9, 2019.

# NOTICE
- (Tracking control) Add the argument `pos_cmd_func` (a function of time such that (t) -> position(t) ∈ ℝ^3) of constructor,
please do not add the argument `pos_cmd` in the inner function of `Dynamics!(env::BacksteppingPositionController)`.
- (Set-point control) Do not add the argument `pos_cmd_func` (a function of time such that (t) -> position(t) ∈ ℝ^3) of constructor,
please add the argument `pos_cmd` in the inner function of `Dynamics!(env::BacksteppingPositionController)`.
"""
struct BacksteppingPositionController <: AbstractEnv
    Ref_model
    Ap
    Bp
    P
    Kp
    Kt
    Kω
    function BacksteppingPositionController(m::Real; pos_cmd_func=nothing)
        @assert m > 0
        # from (s+5)^5 = s^5 + 25s^4 + ... (for **faster response**)
        Ks = []
        push!(Ks, Diagonal(3125*ones(3)))
        push!(Ks, Diagonal(3125*ones(3)))
        push!(Ks, Diagonal(1250*ones(3)))
        push!(Ks, Diagonal(250*ones(3)))
        push!(Ks, Diagonal(25*ones(3)))
        auto_diff = pos_cmd_func == nothing ? false : true
        Ref_model = ReferenceModel(length(Ks)-1, Ks, auto_diff, pos_cmd_func)
        # position
        Kx = m*1*Matrix(I, 3, 3)
        Kv = m*1*1.82*Matrix(I, 3, 3)
        Kp = hcat(Kx, Kv)
        # Q = Diagonal(1*ones(6))
        Q = diagm(1*ones(6))
        # thrust
        Kt = Diagonal(4*ones(3))
        # angular
        Kω = Diagonal(20*ones(3))
        # reference model
        Ap = [zeros(3, 3) Matrix(I, 3, 3);
              -(1/m)*Kx -(1/m)*Kv]
        Bp = [zeros(3, 3);
              (1/m)*Matrix(I, 3, 3)]
        P = MatrixEquations.lyapc(Ap, Q)
        new(Ref_model, Ap, Bp, P, Kp, Kt, Kω)
    end
end

"""
# Notes
ref_model.x_0 = xd
ref_model.x_1 = vd
ref_model.x_2 = ad
ref_model.x_3 = ȧd
ref_model.x_4 = äd
"""
function State(controller::BacksteppingPositionController)
    @unpack Ref_model = controller
    return function (pos0, m, g)
        @assert m > 0
        ref_model = State(Ref_model)(pos0)
        Td = m*g
        ComponentArray(ref_model=ref_model, Td=Td)
    end
end

function Dynamics!(controller::BacksteppingPositionController)
    @unpack Ref_model = controller
    @Loggable function dynamics!(dX, X, p, t; pos_cmd=nothing, Ṫd)
        @log state = X
        Dynamics!(Ref_model)(dX.ref_model, X.ref_model, (), t; x_cmd=pos_cmd)  # be careful; parameter = ()
        dX.Td = Ṫd
        nothing
    end
end

"""
Several functions are exported from `utils.jl`, e.g., T_u_inv(T).
"""
function Command(controller::BacksteppingPositionController)
    @unpack Ap, Bp, P, Kp, Kt, Kω = controller
    T_u_inv(T) = [   0 1/T  0;
                  -1/T   0  0;
                     0   0 -1]
    T_u_inv_dot(T, Ṫ) = [    0 -Ṫ/T^2 0;
                         Ṫ/T^2      0 0;
                             0      0 0]
    T_ω(T) = [0 -T  0;
              T  0  0;
              0  0  0]
    skew(x) = [    0  -x[3]   x[2];
                x[3]      0  -x[1];
               -x[2]   x[1]     0]
    return function(p, v, R, ω,
                    xd, vd, ad, ȧd, äd, Td,
                    m::Real, J, g::Real,)
        ex = xd - p
        ev = vd - v
        ep = vcat(ex, ev)
        # u1
        g_vec = [0, 0, g]
        u1 = m*(ad - g_vec) + Kp*ep
        zB = R' * [0, 0, 1]
        td = -Td * zB
        et = u1 - td
        ėp = Ap*ep + Bp*et
        u̇1 = m*ȧd + Kp*ėp
        T = Td  # TODO: no lag
        # u2
        u2 = T_u_inv(T) * R * (2*Bp' * P * ep + u̇1 + Kt*et)
        Ṫd = u2[end]  # third element
        Ṫ = Ṫd  # TODO: no lag
        żB = -R' * T_ω(1.0) * ω
        ėt = u̇1 + u2[end]*zB + Td*żB
        ëp = Ap*ėp + Bp*ėt
        ü1 = m*äd + Kp*ëp
        Ṙ = -skew(ω) * R
        u̇2 = (
              (T_u_inv_dot(T, Ṫ)*R + T_u_inv(T)*Ṙ) * (2*Bp'*P*ep + u̇1 + Kt*et)
              + T_u_inv(T) * R * (2*Bp'*P*ėp + ü1 + Kt*ėt)
             )
        ω̇d = [1 0 0;
              0 1 0;
              0 0 0] * u̇2
        ωd = [u2[1:2]..., 0]
        eω = ωd - ω
        Md = cross(ω, J*ω) + J*(T_ω(T)'*R*et + ω̇d + Kω*eω)
        νd = vcat(Td, Md)
        e = [ep;
             et;
             eω]
        νd, Ṫd, e, zB, T
    end
end


text:
module EndowmentScripts

const Nothing = Void #NOTE: 0.7 compatabiltiy hack, delete on ugrade

if  pwd() ∉ LOAD_PATH
  push!(LOAD_PATH,pwd())
end

using Endowment, CSV, CodecZlib, DataFrames

const WORKERS=2
const CLEAN_WORKERS=false

#set up for parallel operations
initializePar(WORKERS, reset=CLEAN_WORKERS)

#a list of high-level functions to run at the script level
function endowmentTasks()::Nothing
  preProcessData(update=true, validateDownloads=false, refreshMetadata = false, refreshData = false)

  return nothing
end

@time endowmentTasks()

#cleanup the extra processes if desired
if CLEAN_WORKERS
  cleanup()
end
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
module ModelKit

export @var,
    @unique_var,
    Variable,
    Expression,
    variables,
    nvariables,
    subs,
    evaluate,
    differentiate,
    monomials,
    expand,
    System,
    Homotopy,
    evaluate!,
    jacobian!,
    evaluate_and_jacobian!,
    diff_t!,
    evaluate,
    jacobian,
    evaluate_and_jacobian,
    diff_t

import LinearAlgebra

include("symengine.jl")
include("symbolic.jl")
include("instructions.jl")
include("codegen.jl")

end # module


text:
""" IN&OUT INTERFACE FOR SOLVER DEVELOPING \n
    {AnalysisDataStruct} = UI_dev.parsing_xml_file(xmlFileDir::string)
          nothing        = UI_dev.Writing_ascii_file(inSimDataStruct::SimulationDataStruct,  fileName::string) \n
    __________________________________________________________________________________________________________ \n
"""
module IO_dev
    export
        # PreProcessing
        Parsing_xml_file,
        # PostProcessing
        Writing_ascii_file

    include("PreProcessing.jl")
    include("PostProcessing.jl")
end

max_stars_repo_name: Hipgineer/JuliaPBF

max_issues_repo_name: Hipgineer/JuliaPBF

max_forks_repo_name: Hipgineer/JuliaPBF


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: forecastquery
using AWS.Compat
using AWS.UUIDs

"""
    query_forecast(filters, forecast_arn)
    query_forecast(filters, forecast_arn, params::Dict{String,<:Any})

Retrieves a forecast for a single item, filtered by the supplied criteria. The criteria is
a key-value pair. The key is either item_id (or the equivalent non-timestamp, non-target
field) from the TARGET_TIME_SERIES dataset, or one of the forecast dimensions specified as
part of the FeaturizationConfig object. By default, QueryForecast returns the complete date
range for the filtered forecast. You can request a specific date range. To get the full
forecast, use the CreateForecastExportJob operation.  The forecasts generated by Amazon
Forecast are in the same timezone as the dataset that was used to create the predictor.

# Arguments
- `filters`: The filtering criteria to apply when retrieving the forecast. For example, to
  get the forecast for client_21 in the electricity usage dataset, specify the following:
  {\"item_id\" : \"client_21\"}  To get the full forecast, use the CreateForecastExportJob
  operation.
- `forecast_arn`: The Amazon Resource Name (ARN) of the forecast to query.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"EndDate"`: The end date for the forecast. Specify the date using this format:
  yyyy-MM-dd'T'HH:mm:ss (ISO 8601 format). For example, 2015-01-01T20:00:00.
- `"NextToken"`: If the result of the previous request was truncated, the response includes
  a NextToken. To retrieve the next set of results, use the token in the next request. Tokens
  expire after 24 hours.
- `"StartDate"`: The start date for the forecast. Specify the date using this format:
  yyyy-MM-dd'T'HH:mm:ss (ISO 8601 format). For example, 2015-01-01T08:00:00.
"""
function query_forecast(
    Filters, ForecastArn; aws_config::AbstractAWSConfig=global_aws_config()
)
    return forecastquery(
        "QueryForecast",
        Dict{String,Any}("Filters" => Filters, "ForecastArn" => ForecastArn);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function query_forecast(
    Filters,
    ForecastArn,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return forecastquery(
        "QueryForecast",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}("Filters" => Filters, "ForecastArn" => ForecastArn),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end


text:
##########################################################################
#
#    Misc expression manipulation function
#
##########################################################################

#### translates ~ into regular syntax
function translate(ex::Expr)
	if ex.head == :block
		return Expr(:block, translate(ex.args)...)

	#  handles ~ for julia 0.2
	elseif :(a ~ b).head == :call &&  # test if < julia 0.3 (cf. #4882)
		ex.head == :call && ex.args[1] == :~

		length(ex.args) == 3 || error("Syntax error in ($ex)")

		ex2 = ex.args[3]
		if isa(ex2, Expr) && length(ex2.args)==2 && ex2.args[1] == :+   #  ~+  (right censoring) statement
			return :( $ACC_SYM += logccdf( $(ex2.args[2]), $(ex.args[2]) ) )

		elseif isa(ex2, Expr) && length(ex2.args)==2 && ex2.args[1] == :-  #  ~-  (left censoring) statement
			return :( $ACC_SYM += logcdf( $(ex2.args[2]), $(ex.args[2]) ) )

		elseif isa(ex2, Expr) || isa(ex2, Symbol)   # ~ statement
			return :( $ACC_SYM += logpdf( $(ex2), $(ex.args[2]) ) )

		else
			error("Syntax error in ($ex)")

		end

	#  handles ~ for julia 0.3+
	elseif :(a ~ b).head == :macrocall &&  # test if >= julia 0.3 (cf. #4882)
		ex.head == :macrocall && ex.args[1] == symbol("@~")

		length(ex.args) == 3 || error("Syntax error in ($ex)")

		ex2 = ex.args[3]
		if isa(ex2, Expr) && length(ex2.args)==2 && ex2.args[1] == :+   #  ~+  (right censoring) statement
			return :( $ACC_SYM += logccdf( $(ex2.args[2]), $(ex.args[2]) ) )

		elseif isa(ex2, Expr) && length(ex2.args)==2 && ex2.args[1] == :-  #  ~-  (left censoring) statement
			return :( $ACC_SYM += logcdf( $(ex2.args[2]), $(ex.args[2]) ) )

		elseif isa(ex2, Expr) || isa(ex2, Symbol)   # ~ statement
			return :( $ACC_SYM += logpdf( $(ex2), $(ex.args[2]) ) )

		else
			error("Syntax error in ($ex)")

		end

	else
		return ex
	end
end
translate(ex::Vector) = map(translate, ex)
translate(ex::Any) = ex

### creates mapping statements from Vector{Float64} to model parameter variables
function vec2var(;init...)
	ex = Expr[]
	pos = 1
	for (v,i) in init
		sz = size(i)
		if length(sz) == 0  # scalar
			push!(ex, :($v = $PARAM_SYM[ $pos ]) )
			pos += 1
		elseif length(sz) == 1  # vector
			r = pos:(pos+sz[1]-1)
			push!(ex, :($v = $PARAM_SYM[ $(Expr(:(:), pos, pos+sz[1]-1)) ]) )
			pos += sz[1]
		else # matrix case  (needs a reshape)
			r = pos:(pos+prod(sz)-1)
			push!(ex, :($v = reshape($PARAM_SYM[ $(Expr(:(:), pos, pos+prod(sz)-1)) ], $(sz[1]), $(sz[2]))) )
			pos += prod(sz)
		end
	end
	ex
end

### creates mapping statements from model parameter variables to Vector{Float64}
# FIXME : using undocumented dprefix function of ReverseDiffSource (should be replaced)
function var2vec(;init...)
	ex = {}
	for (v,i) in init
		sz = size(i)
		if in(length(sz), [0,1]) # scalar or vector
			push!(ex, ReverseDiffSource.dprefix(v))
		else # matrix case  (needs a reshape)
			push!(ex, :( vec($(ReverseDiffSource.dprefix(v))) ) )
		end
	end
	Expr(:vcat, ex...)
end

### returns parameter info : total size, vector <-> model parameter map, inital values vector
function modelVars(;init...)
	# init = [(:x, 3.)]
    pars = Dict{Symbol, NTuple{2}}()
    pos = 1
    vi = Float64[]

    for (par, def) in init  # par, def = init[1]
    	isa(def, Real) || isa(def, AbstractVector) || isa(def, AbstractMatrix) ||
    		error("unsupported type for parameter $(par)")

        pars[par] = (pos, size(def))
        pos += length(def)
        vi = [vi, float64([def...])]
    end
    (pos-1, pars, vi)
end


text:
module TestMeta
    using Test
    using CBSOData4

    include("mock_getodata.jl")

    meta_880072ned = CBSOData4.get_meta("80072ned")
    @test "Properties" in keys(meta_880072ned)
    @test Set(["PeriodenGroups", "BedrijfskenmerkenSBI2008Codes", "BedrijfskenmerkenSBI2008Groups", "Dimensions", "MeasureCodes", "PeriodenCodes", "Properties"]) == Set(keys(meta_880072ned))

end


text:
module SimpleElastixRegistration

using PyCall, DelimitedFiles
export prepare_elastix_registration, prepare_elastix_transformix, load_points

# Write your package code here.
""" Prepare for elastix registration
`Params` is a tuple of parameters or of files that includes parameters.
"""
function prepare_elastix_registration(sitk, outputdir, fixed, moving, params::Tuple)
  elastixImageFilter = sitk.ElastixImageFilter()
  elastixImageFilter.SetFixedImage(fixed)
  elastixImageFilter.SetMovingImage(moving)
  parameterMapVector = sitk.VectorOfParameterMap()
  for param in params
    if isfile(param)
      parameterMapVector.append(sitk.ReadParameterFile(param))
    else
      parameterMapVector.append(sitk.GetDefaultParameterMap(param))
    end
  end
  elastixImageFilter.SetParameterMap(parameterMapVector)
  isdir(outputdir) ? nothing : mkdir(outputdir); elastixImageFilter.SetOutputDirectory(outputdir)
  return(elastixImageFilter)
end
prepare_elastix_registration(sitk, outputdir, fixed, moving, param::String) = prepare_elastix_registration(sitk, outputdir, fixed, moving, (param,))

""" Add manual points to help registration
'fixedpts_fn' and 'mvpts_fn' follow the point format of Simpleelastix.
"""
function prepare_elastix_registration(sitk, outputdir, fixed, moving, params::Tuple, metric0, metric1, fixedpts_fn::String, mvpts_fn::String)
  elastixImageFilter = sitk.ElastixImageFilter()
  elastixImageFilter.SetFixedImage(fixed)
  elastixImageFilter.SetMovingImage(moving)
  parameterMapVector = sitk.VectorOfParameterMap()
  for param in params
    if isfile(param)
      parameterMapVector.append(sitk.ReadParameterFile(param))
    else
      parameterMapVector.append(sitk.GetDefaultParameterMap(param))
    end
  end
  elastixImageFilter.SetParameterMap(parameterMapVector)
  elastixImageFilter.AddParameter("Metric", ("NormalizedMutualInformation", "CorrespondingPointsEuclideanDistanceMetric" ))
  elastixImageFilter.SetParameter("Metric0Weight", metric0)
  elastixImageFilter.SetParameter("Metric1Weight", metric1)
  elastixImageFilter.SetFixedPointSetFileName(fixedpts_fn)
  elastixImageFilter.SetMovingPointSetFileName(mvpts_fn)
  isdir(outputdir) ? nothing : mkdir(outputdir); elastixImageFilter.SetOutputDirectory(outputdir)
  return(elastixImageFilter)
end

""" Prepare for elastix transformix: This is useful when transformation has been acquired by elastix registration.
`tformfns` is a tuple of files that includes parameters.
"""
function prepare_elastix_transformix(sitk, outputdir, moving, tformfns::Tuple)
  transformixImageFilter = sitk.TransformixImageFilter()
  for (i, tformfn) in enumerate(tformfns)
    parameterMap = sitk.ReadParameterFile(tformfn)
    if i == 1
      transformixImageFilter.SetTransformParameterMap(parameterMap)
    else
      transformixImageFilter.AddTransformParameterMap(parameterMap)
    end
  end
  transformixImageFilter.SetMovingImage(moving)
  isdir(outputdir) ? nothing : mkdir(outputdir); transformixImageFilter.SetOutputDirectory(outputdir)
  return(transformixImageFilter)
end
prepare_elastix_transformix(sitk, outputdir, moving, tformfns::String) = prepare_elastix_transformix(sitk, outputdir, moving, (tformfn,))

""" load points from elastix outputpoint file"""
function load_points(fn)
  outputpts = readdlm(fn, '\t')
  strs = split.(outputpts[:,5])
  pts  = [(parse(Int, pt[5]), parse(Int, pt[6]), parse(Int, pt[7])) for pt in strs]
  return(pts)
end

end


text:
@doc """
    powertostring(base, number)
Converts a number b^p into a string "bᵖ".
#Example
```julia-repl
julia> powertostring(2, 2^BigInt(1345))
"2¹³⁴⁵"
```
""" powertostring!
function powertostring(base, number)
    power = Int(round(log(base,number))) #scale = base^power
    s = "$base"*superscript(power)
    return s
end

@doc """
    photon_dens_integrated(sim; species = "photon", direction = direction, plot_title = title)

Plot of photon number density along a given direction (direction::Symbol argument; default is :x).
The distribution is normalized to the number density of the electrons at t=0.
""" photon_dens_integrated!
function photon_dens_integrated(sim; species = "photon", direction = :x, plot_title = "Photon number density integrated along x-axis")
    file = sim[1]
    n₀ = file["number_density/electron"]

    file = sim[end]
    nᵧ = file["number_density/$species"]

    kwargs = (color=:rainbow,
              colorbar_title="nᵧ/n₀",
              tick_direction = :out,
              minorticks = true,
              title = plot_title
    )
    dim = dir_to_idx(direction)
    nᵧ_int = dropdims(sum(nᵧ*cell_volume(file), dims=dim), dims=dim)
    n₀_int = dropdims(sum(n₀*cell_volume(file), dims=dim), dims=dim)

    if direction == :x
        a1 = cell_length(file, :y)
        a2 = cell_length(file, :z)
        plt = Plots.heatmap([1:1:size(nᵧ_int,1)]*a1, [1:1:size(nᵧ_int,2)]*a2, transpose(nᵧ_int ./ n₀_int); kwargs...) #damn heatmap is dumbly implemented
        Plots.plot!(plt,
              xlabel = "y (μm)",
              ylabel = "z (μm)",
              size = (600,500)
        )
        return plt
    elseif direction == :y
        a1 = cell_length(file, :x)
        a2 = cell_length(file, :z)
        plt = Plots.heatmap([1:1:size(nᵧ_int,1)]*a1, [1:1:size(nᵧ_int,2)]*a2, transpose(nᵧ_int ./ n₀_int); kwargs...)
        Plots.plot!(plt,
              xlabel = "x (μm)",
              ylabel = "z (μm)",
              size = (665,500)
        )
        return plt
    elseif direction == :z
        a1 = cell_length(file, :x)
        a2 = cell_length(file, :y)
        plt = Plots.heatmap([1:1:size(nᵧ_int,1)]*a1, [1:1:size(nᵧ_int,2)]*a2, transpose(nᵧ_int ./ n₀_int); kwargs...)
        Plots.plot!(plt,
              xlabel = "x (μm)",
              ylabel = "y (μm)",
              size = (665,500)
        )
        return plt
    end
end

@doc """
    photon_en_dens_integrated(sim; species = "photon", direction = direction, endens_scale = scale , plot_title = title)

Plot of photon energy density along a given direction (direction::Symbol argument; default is :x) in MeV/μm².
One can choose to add a 10ˣ scale to the units with the endens_scale argument (pro tip: use BigInt for large powers).
""" photon_en_dens_integrated!
function photon_en_dens_integrated(sim; species = "photon", direction = :x, endens_scale = 10^12 , plot_title = "Photon energy density integrated along x-axis")
    file = sim[end]
    nᵧ, enᵧ= file["number_density/$species", "ekbar/$species"]
    enᵧ = ustrip.(u"MeV",enᵧ)

    scale = powertostring(10, endens_scale)
    kwargs = (color=:rainbow,
             colorbar_title="ϵᵧ ($scale MeV/μm²)",
             size = (665,600),
             tick_direction = :out,
             minorticks = true,
             title = plot_title
    )
    dim = dir_to_idx(direction)
    n = uconvert(NoUnits, nᵧ*cell_volume(file).*enᵧ)

    if direction == :x
        a1 = ustrip(u"μm", cell_length(file, :y))
        a2 = ustrip(u"μm", cell_length(file, :z))
        nᵧ_int = dropdims(sum(n, dims=dim), dims=dim) / (a1*a2) / endens_scale
        plt = Plots.heatmap([1:1:size(nᵧ_int,1)]*a1, [1:1:size(nᵧ_int,2)]*a2, transpose(nᵧ_int); kwargs...)
        Plots.plot!(plt,
              xlabel = "y (μm)",
              ylabel = "z (μm)",
        )
        return plt
    elseif direction == :y
        a1 = ustrip(u"μm", cell_length(file, :x))
        a2 = ustrip(u"μm", cell_length(file, :z))
        nᵧ_int = dropdims(sum(n, dims=dim), dims=dim) / (a1*a2) / endens_scale
        plt = Plots.heatmap([1:1:size(nᵧ_int,1)]*a1, [1:1:size(nᵧ_int,2)]*a2, transpose(nᵧ_int); kwargs...)
        Plots.plot!(plt,
              xlabel = "x (μm)",
              ylabel = "z (μm)",
        )
        return plt
    elseif direction == :z
        a1 = ustrip(u"μm",cell_length(file, :x))
        a2 = ustrip(u"μm",cell_length(file, :y))
        nᵧ_int = dropdims(sum(n, dims=dim), dims=dim) / (a1*a2) / endens_scale
        plt = Plots.heatmap([1:1:size(nᵧ_int,1)]*a1, [1:1:size(nᵧ_int,2)]*a2, transpose(nᵧ_int); kwargs...)
        Plots.plot!(plt,
              xlabel = "x (μm)",
              ylabel = "y (μm)",
        )
        return plt
    end
end

@doc """
    photon_energy_spectrum(sim; species = photon_species, plot_title = title, angle = pi)

Photon energy spectra sᵧ² dN/dsᵧ = f(sᵧ = log₁₀(E/keV)). The photon species can be selected.
One can filter photons in cone. default angle is π (180⁰), i.e. all photons going to the right
""" photon_energy_spectrum!
function photon_energy_spectrum(sim; label = "label", mark_max = false, angle = pi, species = "photon", plot_title = "Compton scattering spectrum")
    file = sim[end]
    px, py, pz = file["px/$species", "py/$species", "pz/$species"]
    # filtering data: keep only photons moving in the direction of the laser
    idx = findall(p -> p > zero(eltype(px)), px)
    px = px[idx]
    py = py[idx]
    pz = pz[idx]

    # θ is measured from the x axis
    # filtering data: keep photons that propagate within 2θ
    pₚₑᵣₚ = sqrt.(py.^2 + pz.^2)   #transversal momentum of each photon
    θ = atan.(pₚₑᵣₚ./px)           #emission angles for each photon
    idx = findall(θ₀ -> θ₀ < angle/2, θ)
    px = px[idx]
    py = py[idx]
    pz = pz[idx]

    # make plot
    Eᵧ = sqrt.(px.^2 + py.^2 + pz.^2) * c_0
    sᵧ = log.(ustrip(u"keV", Eᵧ))

    h = fit(Histogram, sᵧ; nbins = 200)
    s = midpoints(collect(h.edges[1]))
    N = h.weights
    dNds = N./abs(s[1]-s[2])

    idx = findall(x -> x > 0, s)
    s = s[idx]
    dNds = dNds[idx]

    kwargs = (xlabel = "sᵧ = log₁₀(Eᵧ/keV)",
              ylabel = "sᵧ² dN/dsᵧ",
              label = label,
              tickfontsize = 10,
              guidefontsize = 15,
              title = plot_title,
              legendfontsize = 15,
              legendtitlefont = 15,
              titlefont = 16
    )

    plt = Plots.plot(s, dNds.*s.*s, linewidth = 4; kwargs...)

    if mark_max
        dNds = dNds.*s.*s
        p0 = findfirst(smax -> smax == maximum(dNds), dNds)
        Plots.plot!(plt, [s[p0]], [dNds[p0]], seriestype = :scatter,
                                    markersize = 7,
                                    markercolor = :black,
                                    label = false
        )
    end

    return plt
end

# an alternative version of photon_energy_spectrum() that plots multiple spectra on the same plot
function photon_energy_spectrum(sims, labels; legend_title = "labels", mark_max = false, angle = pi, species = "photon", plot_title = "Compton scattering spectrum")
    kwargs = (xlabel = "sᵧ = log₁₀(Eᵧ/keV)",
            ylabel = "sᵧ² dN/dsᵧ",
            tickfontsize = 10,
            guidefontsize = 15,
            title = plot_title,
            legendfontsize = 15,
            legendtitlefont = 15,
            titlefont = 16,
            legendtitle = legend_title
    )
    plt = Plots.plot(;kwargs...)
    for (i, sim) in enumerate(sims)
        file = sim[end]
        px, py, pz = file["px/$species", "py/$species", "pz/$species"]
        # filtering data: keep only photons moving in the direction of the laser
        idx = findall(p -> p > zero(eltype(px)), px)
        px = px[idx]
        py = py[idx]
        pz = pz[idx]

        # θ is measured from the x axis
        # filtering data: keep photons that propagate within 2θ
        pₚₑᵣₚ = sqrt.(py.^2 + pz.^2)   #transversal momentum of each photon
        θ = atan.(pₚₑᵣₚ./px)           #emission angles for each photon
        idx = findall(θ₀ -> θ₀ < angle/2, θ)
        px = px[idx]
        py = py[idx]
        pz = pz[idx]

        # make plot
        Eᵧ = sqrt.(px.^2 + py.^2 + pz.^2) * c_0
        sᵧ = log.(ustrip(u"keV", Eᵧ))

        h = fit(Histogram, sᵧ; nbins = 200)
        s = midpoints(collect(h.edges[1]))
        N = h.weights
        dNds = N./abs(s[1]-s[2])

        idx = findall(x -> x > 0, s)
        s = s[idx]
        dNds = dNds[idx]
        plt = Plots.plot!(plt, s, dNds.*s.*s, linewidth = 4; label = labels[i])

        if mark_max
            dNds = dNds.*s.*s
            p0 = findfirst(smax -> smax == maximum(dNds), dNds)
            Plots.plot!(plt, [s[p0]], [dNds[p0]], seriestype = :scatter,
                                    markersize = 7,
                                    markercolor = :black,
                                    label = false
            )
        end
    end
    return plt
end

@doc """
    energy_spectrum(sim; mass = m_e, species = "electron", plot_title = title)

Energy spectra s² dN/ds = f(s = log₁₀(E/MeV)) for particles with mass.
The mass of the particle has to be added by hand with the mass argument (default is the electron mass).
""" energy_spectrum!
function energy_spectrum(sim; mass = m_e, species = "electron", plot_title = "Electron energy spectrum")
    file = sim[end]
    px, py, pz = file["px/$species", "py/$species", "pz/$species"]

    #should be implemented: read species mass from input.deck
    E₀ = sqrt.(mass^2*c_0^4 .+ (px.^2 + py.^2 + pz.^2) * c_0^2)
    s₀ = log.(ustrip(u"MeV", E₀))

    h = fit(Histogram, s₀; nbins = 200)
    s = midpoints(collect(h.edges[1]))
    N = h.weights
    dNds = N./abs(s[1]-s[2])

    idx = findall(x -> x > 0, s)
    s = s[idx]
    dNds = dNds[idx]

    kwargs = (xlabel = "s = log₁₀(E/MeV)",
              ylabel = "s² dN/ds",
              tickfontsize = 10,
              guidefontsize = 15,
              title = plot_title,
              legendfontsize = 15,
              legendtitlefont = 15,
              titlefont = 16
    )

    plt = Plots.plot(s, dNds.*s.*s, linewidth = 4; kwargs...)

    return plt
end

# an alternative version of energy_spectrum() that plots multiple spectra on the same plot
function energy_spectrum(sims, labels; legend_title = "labels", mass = m_e, species = "electron", plot_title = "Electron energy spectrum")
    kwargs = (xlabel = "s = log₁₀(E/MeV)",
              ylabel = "s² dN/ds",
              tickfontsize = 10,
              guidefontsize = 15,
              title = plot_title,
              legendfontsize = 15,
              legendtitlefont = 15,
              titlefont = 16,
              legendtitle = legend_title
    )
    plt = Plots.plot(;kwargs...)
    for (i, sim) in enumerate(sims)
        file = sim[end]
        px, py, pz = file["px/$species", "py/$species", "pz/$species"]

        #should be implemented: read species mass from input.deck
        E₀ = sqrt.(mass^2*c_0^4 .+ (px.^2 + py.^2 + pz.^2) * c_0^2)
        s₀ = log.(ustrip(u"MeV", E₀))

        h = fit(Histogram, s₀; nbins = 200)
        s = midpoints(collect(h.edges[1]))
        N = h.weights
        dNds = N./abs(s[1]-s[2])

        idx = findall(x -> x > 0, s)
        s = s[idx]
        dNds = dNds[idx]

        plt = Plots.plot!(plt, s, dNds.*s.*s, linewidth = 4; label = labels[i])
    end
    return plt
end

@doc """
    photon_solid_angle_emission(sim;  angle = pi, species = photon_species, plot_title = title)

Angular distribution of gamma energy emission (MeV/srad). The angles are given by the direction of the momentum
(the reference axis is Ox). One can filter photons in cone. default angle is π (180⁰), i.e. all photons going
to the right.
""" photon_solid_angle_emission!
function photon_solid_angle_emission(sim;  angle = pi, species = "photon", plot_title = "Solid angle emission")
    file = sim[end]

    px, py, pz= file["px/$species",
                     "py/$species",
                     "pz/$species"
    ]

    idx = findall(p -> p > zero(eltype(px)), px)
    px = px[idx]
    py = py[idx]
    pz = pz[idx]

    pₚₑᵣₚ = sqrt.(py.^2 + pz.^2)   #transversal momentum of each photon
    θ = atan.(pₚₑᵣₚ./px)           #emission angles for each photon
    ϕ = atan.(py, pz)

    # filtering data: keep photons that propagate within 2θ
    idx = findall(θ₀ -> θ₀ < angle/2, θ)
    px = px[idx]
    py = py[idx]
    pz = pz[idx]
    θ = θ[idx]
    ϕ = ϕ[idx]

    Eᵧ = sqrt.(px.^2 + py.^2 + pz.^2) * c_0

    data = (ϕ .+ π, θ)
    w = weights(ustrip.(u"MeV", Eᵧ))

    h = fit(Histogram, data, w, nbins = (360,90))
    xedges = collect(h.edges[1])
    yedges = rad2deg.(h.edges[2])
    H = h.weights

    ax = PyPlot.Axes3D(PyPlot.figure())
    PyPlot.subplot(projection="polar")
    PyPlot.pcolormesh(xedges[2:end], yedges[2:end], transpose(H))
    PyPlot.colorbar(label = "dEᵧ/dΩ (MeV/srad)")
    PyPlot.title(plot_title)
    ax[:grid](true)
    # PyPlot.display(PyPlot.gcf())
    return PyPlot.gcf()
end


text:
using OrdinaryDiffEq
prob = ODEProblem((u,p,t) -> -u,BigFloat(1.0),(0.0,1.0))
solve(prob,Tsit5())
solve(prob,KenCarp4())

# Function initial condition
f(u,p,t) = u
u0(p,t0) = ones(2)
prob = ODEProblem(f,u0,(0.0,1.0))
sol = solve(prob,Tsit5())


text:
@testset "Computations: find_root!" begin
  iterator = PolynomialIterator([1], 1)
  # => 1, 1 + x, only one root -1

  image = RootsImage(-1.5-1.0im, 1.5+1.5im, 3, 1)

  find_roots!(iterator, image)

  @test image.data[1,1] == 1.0
  @test image.data[2,1] == 0.0
  @test image.data[3,1] == 0.0
end


text:
function SimpleVectorValueChecker(arg0::jdouble, arg1::jdouble)
    return SimpleVectorValueChecker((jdouble, jdouble), arg0, arg1)
end

function SimpleVectorValueChecker(arg0::jdouble, arg1::jdouble, arg2::jint)
    return SimpleVectorValueChecker((jdouble, jdouble, jint), arg0, arg1, arg2)
end

function converged(obj::SimpleVectorValueChecker, arg0::jint, arg1::PointVectorValuePair, arg2::PointVectorValuePair)
    return jcall(obj, "converged", jboolean, (jint, PointVectorValuePair, PointVectorValuePair), arg0, arg1, arg2)
end


text:
# module Sorting
#
# using LightGraphs
# using Random

export
    topological_sort,
    find_index_in_sorted_array,
    insert_to_sorted_array!,
    linear_interp,
    get_interp_interval,
    resample_array

"""
    `topological_sort(G)`

Returns a topological sort of the vertices of a graph, with the property that
v1 < v2 iff there is not forward path through the graph from v2 to v1.
G must be non-cyclic. Can handle disjoint graphs.
"""
function topological_sort(G)
    @assert(!is_cyclic(G))
    frontier = Set{Int}(rand(1:nv(G)))
    # explored = Set{Int}()
    explored = fill!(Vector{Bool}(undef,nv(G)), false)
    ordering = []
    while any(explored .!= true)
        while length(frontier) > 0
            v = pop!(frontier)
            tip = true # check whether all downstream nodes are explored
            for v2 in inneighbors(G,v)
                if !explored[v2]
                    tip = false
                    push!(frontier,v2)
                    break
                end
            end
            if tip == true
                # push!(explored,v)
                explored[v] = true
                push!(ordering, v)
                for v2 in outneighbors(G,v)
                    push!(frontier,v2)
                end
            end
        end
        for v in 1:nv(G)
            if !explored[v]
                push!(frontier, v)
            end
        end
    end
    ordering
end

"""
    `find_index_in_sorted_array(array, x)`

    Assumes that array is already sorted. Returns index at which x would need to
    be inserted in order to maintain ordering of array. Chooses the smallest
    index in the case of a tie.

    Looks like  "Base.searchsorted" does the same thing as this.
"""
function find_index_in_sorted_array(array, x)
    A = 0
    C = length(array)+1
    B = max(1,Int(round((A+C) / 2)))
    while C-A > 1
        if x < array[B] || ( !(array[B] < x) && !(x < array[B]))
            A = A
            C = B
            B = Int(ceil((A+C) / 2))
        else
            A = B
            C = C
            B = Int(ceil((A+C) / 2))
        end
    end
    return B
end

"""
    `insert_to_sorted_array!(array, x)`

    Assumes that array is already sorted. Inserts new element x so that
    array remains sorted. Requires that Base.isless(a::C,b::C) where
    C==typeof(x) be implemented.
"""
function insert_to_sorted_array!(array, x)
    B = find_index_in_sorted_array(array, x)
    insert!(array, B, x)
    array
end

"""
    linear_interpolation
"""
function linear_interp(a,b,t)
    return a*(1-t) + b*t
end
function linear_interp(a,b,t1,t2,t)
    lin_interp(a,b,(t-t1)/(t-t2))
    return a*(1-t) + b*t
end
function get_interp_interval(array,x)
    idx = find_index_in_sorted_array(array, x)
    if idx == 1
        return idx,idx+1
    elseif idx > length(array)
        return max(1,idx-2),idx-1
    else
        return idx-1,idx
    end
end

"""
    resample array
"""
function resample_array(array::Matrix{Float64},old_t_vec,new_t_vec)
    new_array = Vector{Vector{Float64}}()
    for t in new_t_vec
        idx1,idx2 = get_interp_interval(old_t_vec,t)
        push!(new_array, linear_interp(array[idx1,:],array[idx2,:],(t-old_t_vec[idx1])/(old_t_vec[idx2]-old_t_vec[idx1])))
    end
    return collect(hcat(new_array...)')
end
function resample_array(array::A,old_t_vec,new_t_vec) where {A}
    new_array = Vector{Vector{Float64}}()
    for t in new_t_vec
        idx1,idx2 = get_interp_interval(old_t_vec,t)
        push!(new_array, linear_interp(array[idx1],array[idx2],(t-old_t_vec[idx1])/(old_t_vec[idx2]-old_t_vec[idx1])))
    end
    return new_array
end

export
    clip
"""
    `clip(a,b,c)`

    Returns the closest element to `a` on the interval `[b,c]`
"""
clip(a,b,c) = max.(min.(a,c),b)

export
    draw_random_uniform

"""
    draw_random_elements(vec,n)

Draw `n` elements uniformly at random (without replacement) from vec.
"""
function draw_random_uniform(vec,n)
    N = length(vec)
    @assert N >= n "Cannot draw $n elements from vector of length $N"
    idxs = collect(1:N)
    selected_elements = Vector{eltype(vec)}()
    for i in 1:n
        k = rand(1:length(idxs))
        push!(selected_elements, vec[idxs[k]])
        deleteat!(idxs,k)
    end
    return selected_elements
end

# end # End of module GraphSorting


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 Random, Distributions, Plots, LaTeXStrings
	Random.seed!(0)
end

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

# ╔═╡ fa68607e-22e5-11eb-0558-c9a4d9f77426
begin
	prior(lam) = pdf(TriangularDist(0, 10, 3), lam)
	data = [2,1,0,0,1,0,2,2,5,2,4,0,3,2,5,0]

	like(lam) = *([pdf(Poisson(lam),x) for x in data]...)
	posteriorUpToK(lam) = like(lam)*prior(lam)
end

# ╔═╡ 2bd08a88-3afb-11eb-03e1-c3af32635947
begin
	delta = 10^-4.
	lamRange = 0:delta:10
	K = sum([posteriorUpToK(lam)*delta for lam in lamRange])
	posterior(lam) = posteriorUpToK(lam)/K
end

# ╔═╡ 2bd0c412-3afb-11eb-1e28-27323df7bbe6
begin
	bayesEstimate = sum([lam*posterior(lam)*delta for lam in lamRange])
	Text("Bayes estimate: $(bayesEstimate)")
end

# ╔═╡ 2bd13636-3afb-11eb-3c12-2d5a0c0e4b13
begin
	plot(lamRange, prior.(lamRange),
		c=:blue, label="Prior distribution")
	plot!(lamRange, posterior.(lamRange),
		c=:red, label="Posterior distribution",
		xlims=(0, 10), ylims=(0, 1.2),
		xlabel=L"\lambda",ylabel="Density")
end

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

# ╔═╡ Cell order:
# ╟─ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
# ╠═e077ab16-22e5-11eb-15ed-5f0654bf7328
# ╠═e75a9a06-22e5-11eb-2e36-8d4f62b830ed
# ╠═fa68607e-22e5-11eb-0558-c9a4d9f77426
# ╠═2bd08a88-3afb-11eb-03e1-c3af32635947
# ╠═2bd0c412-3afb-11eb-1e28-27323df7bbe6
# ╠═2bd13636-3afb-11eb-3c12-2d5a0c0e4b13
# ╟─475ff888-22e6-11eb-2354-09f8f40a8e12


text:
using FFTW, Statistics

function KuramotoSivashinsky(Q, N, tmax, h; M=16, tr=1000)
    #@assert N/λ % 1 == 0 "N/λ % 1 != 0"
    nmax = round(Int,tmax/h)
    x = N*(1:Q) / Q
    #x = N*(-(Q-1)/2:(Q-1)/2) / Q
    u = cos.(π*x/N).*(1 .+sin.(π*x/N))
    v = fft(u)
    #cc= μ* fft(cos.(2π*x/λ))
    #Precompute various ETDRK4 scalar quantities
    #Why all the mean operations?
    k = 2π/N*(0.0:Q-1)
    #k = 2π/N*[(0:(Q/2-1))..., 0, (-Q/2+1:-1)...]
    L = (k.^2 - k.^4) # Fourier Multiples
    E =  exp.(h*L)
    E2 = exp.(h*L/2)
    M = 16 # No. of points for complex means
    r = exp.(im*π*((1:M) .-0.5)/M)
    LR = h*L*ones(M)' + ones(Q)*r'
    QQ = h*real.(mean((exp.(LR/2).-1)./LR, dims=2))[:]
    f1 = h*real(mean((-4 .-LR+exp.(LR).*(4 .-3*LR+ LR.^2))./LR.^3,dims=2))[:]
    f2 = h*real(mean((2 .+LR+exp.(LR).*(-2 .+LR))./LR.^3,dims=2))[:]
    f3 = h*real(mean((-4 .-3*LR-LR.^2+exp.(LR).*(4 .-LR))./LR.^3,dims=2))[:]
    g = -0.5im*k
    #g = 0.5im*k*Q
    uu = [u]
    tt = [0.]

    T = plan_fft(v)
    Ti = plan_ifft(v)
    T! = plan_fft!(v)
    Ti! = plan_ifft!(v)
    a = Complex.(zeros(Q))
    b = Complex.(zeros(Q))
    c = Complex.(zeros(Q))
    Nv = Complex.(zeros(Q))
    Na = Complex.(zeros(Q))
    Nb = Complex.(zeros(Q))
    Nc = Complex.(zeros(Q))
    tmp = Complex.(zeros(Q))

    for n = 1:nmax+tr
                            Nv .= g .* (T*real(Ti*v).^2) #.+ cc
        @.  a  =  E2*v + QQ*Nv
                           Na .= g .* (T!*real(Ti!*a).^2) #.+ cc
        @. b  =  E2*v + QQ*Na
                            Nb .= g.* (T!*real(Ti!*b).^2) #.+ cc
        @. c  =  E2*a + QQ*(2Nb-Nv)
        Nc .= g.* (T!*real(Ti!*c).^2) #.+ cc
        @. v =  E*v + Nv*f1 + 2*(Na+Nb)*f2 + Nc*f3
        if n > tr
            t = n*h
            u = real(Ti*v)
            push!(uu,u)
            push!(tt,t)
        end
    end
    return uu,tt
end


text:
function myequal(a,b)
    println(typeof(a)==typeof(b) ? (:type,true) : (:type,false))
    for n in union(fieldnames(a), fieldnames(b))
        if isdefined(a,n) && isdefined(b,n)
            println(isequal(a.(n), b.(n)) ? (n,true) : (n,false))
        elseif isdefined(a,n) || isdefined(b,n)
            println((n,false))
        else
            println((n,:undef))
        end
    end
end


text:
@testset "sequences.jl" begin

    @testset "PGSE" begin
        δ, Δ = 500.0, 2000.0
        f = PGSE(δ, Δ)
        @test f.δ ≈ δ
        @test f.Δ ≈ Δ
        @test f(-100.0) ≈ 0.0
        @test f(prevfloat(0.0)) ≈ 0.0
        @test f(0.0) ≈ 1.0
        @test f(prevfloat(δ)) ≈ 1.0
        @test f(δ) ≈ 0.0
        @test f(nextfloat(δ)) ≈ 0.0
        @test f(prevfloat(Δ)) ≈ 0.0
        @test f(Δ) ≈ -1.0
        @test f(nextfloat(Δ)) ≈ -1.0
        @test f(Δ + δ + 100.0) ≈ 0.0

        # Integral
        T = LinRange(-δ, Δ + 2δ, 10)
        for t ∈ T
            i = intervals(f)
            i = i[0.0 .< i .< t]
            @test integral(f, t) ≈ quadgk(f, 0.0, i..., t)[1] rtol = 1e-8
        end

        # B-value
        b = quadgk(t -> integral(f, t)^2, 0.0, Δ + δ)[1]
        @test int_F²(f) ≈ b rtol = 1e-8
    end

    @testset "DoublePGSE" begin
        δ, Δ, p = 500.0, 2000.0, 1000.0
        tmid = p + δ + Δ
        f = DoublePGSE(δ, Δ, p)

        # Properties
        @test f.δ ≈ δ
        @test f.Δ ≈ Δ
        @test f.tpause ≈ p

        # Function call
        @test f(-100.0) ≈ 0.0
        @test f(prevfloat(0.0)) ≈ 0.0
        @test f(0.0) ≈ 1.0
        @test f(nextfloat(0.0)) ≈ 1.0
        @test f(prevfloat(δ)) ≈ 1.0
        @test f(δ) ≈ 0.0
        @test f(nextfloat(δ)) ≈ 0.0
        @test f(prevfloat(Δ)) ≈ 0.0
        @test f(Δ) ≈ -1.0
        @test f(nextfloat(Δ)) ≈ -1.0
        @test f(prevfloat(Δ + δ)) ≈ -1.0
        @test f(nextfloat(Δ + δ)) ≈ 0.0
        @test f(-prevfloat(tmid)) ≈ 0.0
        @test f(tmid) ≈ 1.0
        @test f(nextfloat(tmid)) ≈ 1.0
        @test f(prevfloat(tmid + δ)) ≈ 1.0
        @test f(nextfloat(tmid + δ)) ≈ 0.0
        @test f(tmid + δ) ≈ 0.0
        @test f(prevfloat(tmid + Δ)) ≈ 0.0
        @test f(tmid + Δ) ≈ -1.0
        @test f(nextfloat(tmid + Δ)) ≈ -1.0
        @test f(prevfloat(tmid + Δ + δ)) ≈ -1.0
        @test f(tmid + Δ + δ) ≈ 0.0
        @test f(nextfloat(tmid + Δ + δ)) ≈ 0.0
        @test f(tmid + Δ + δ + 100.0) ≈ 0.0

        # Integral
        T = LinRange(-δ, p + 2Δ + 3δ, 10)
        for t ∈ T
            i = intervals(f)
            i = i[0.0 .< i .< t]
            @test integral(f, t) ≈ quadgk(f, 0.0, i..., t)[1] rtol = 1e-8
        end

        # B-value
        b = quadgk(t -> integral(f, t)^2, 0.0, tmid + Δ + δ)[1]
        @test int_F²(f) ≈ b rtol = 1e-8
    end

    @testset "CosOGSE" begin

    end

    @testset "SinOGSE" begin

    end

    @testset "UnknownSequence" begin
        struct UnknownSequence{T} <: TimeProfile{T} end

        f = UnknownSequence{Float64}()

        @test_throws Exception f(0.0)
        @test_throws Exception integral(f, 0.0)
        @test_throws Exception int_F²(f)
        @test_throws Exception intervals(f)
        @test_throws Exception echotime(f)
        @test_throws Exception get_interval(f, 0.0)
        @test_throws Exception constant_intervals(f)
        @test_throws Exception isconstant(f, 0.0)
    end
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 629d929c-6a3d-11eb-0413-5185929928cb
begin
	using DataFrames, Random
	using PlutoUI
end

# ╔═╡ 0b08257e-6a3d-11eb-3da0-ed75be40d90e
md"""
## Introduction to DataFrames

ref. [Introduction to DataFrames, Bogumił Kamiński](https://github.com/bkamins/Julia-DataFrames-Tutorial/blob/master/01_constructors.ipynb)

$(html"
© Pascal, April 2021
")
"""

# ╔═╡ 2c93b1b1-39f7-49c3-8d19-060fac1dba52
PlutoUI.TableOfContents(indent=true, depth=4, aside=true)

# ╔═╡ 81665132-6a3d-11eb-37b7-991da3953909
md"""
### Constructors and Conversions

#### Constructors
"""

# ╔═╡ 82e2153c-6a3d-11eb-0397-e78654843e7e
DataFrame();  ## create an empty DataFrame

# ╔═╡ b037eb42-6a3d-11eb-0443-c3fef979ac19
begin
	n = 5
	DataFrame(A=1:n, B=rand(n), C=randstring.([n for _ ∈ 1:n]), fixed=42)
	## a 5 × 4 DataFrame -
	## Remark on fixed column, the scalar gets automatically broadcasted
end

# ╔═╡ 18ff9300-6a3e-11eb-158b-5bec35b58c84
begin
	## Create a DataFrame from a dictionary
	hsh = Dict(:A => collect(1:n), :B => [x % 2 == 0 for x ∈ 1:n],
		:C => [x % 2 == 0 ? 'a' : 'b' for x ∈ 1:n],
		:fixed => Ref([1 for _ ∈ 1:n]))

	DataFrame(hsh)
end

# ╔═╡ 6c87db40-6a3e-11eb-3a1c-7d3c95ffc793
md"""
**Remarks**

This time we used `Ref` to protect a vector from being treated as a column and forcing broadcasting it into every row of :fixed column (note that the [1, 1, 1, 1, 1] vector is aliased in each row).

Rather than explicitly creating a dictionary first, as above, we could pass `DataFrame` arguments with the syntax of dictionary key-value pairs.
"""

# ╔═╡ 9c069eba-6a3e-11eb-0b99-69c73c8db3f8
DataFrame(:A => collect(1:n), :B => [x % 2 == 0 for x ∈ 1:n],
	:C => [x % 2 == 0 ? 'a' : 'b' for x ∈ 1:n],
	:fixed => Ref([1 for _ ∈ 1:n]))

# ╔═╡ c7eb6274-6a3e-11eb-3846-e337d215f871
## or pass a Vector of pairs
DataFrame([:A => collect(1:n),
		:B => [x % 2 == 0 for x ∈ 1:n],
		:C => [x % 2 == 0 ? 'a' : 'b' for x ∈ 1:n],
    	:fixed => "const"])

# ╔═╡ eafc3138-6a3e-11eb-11e8-c333d5bf13c1
## Create a DataFrame from a vector of vectors, and each vector becomes a column.
DataFrame([rand(Float32, n) for _ ∈ 1:n])

# ╔═╡ d53123d6-6a3e-11eb-1304-c938d3ffc423
## As above with explicit column naming
DataFrame([rand(n) for _ ∈ 1:n], [Symbol(String("y$(ix)")) for ix ∈ 1:n])

# ╔═╡ 5f83c5b6-6a3f-11eb-09ad-2d01fa5e3a6e
## It is not allowed to pass a vector of scalars to the DataFrame constructor
## instead use a transposed vector
DataFrame(permutedims(collect(1:n)))

# ╔═╡ a3dd623a-6a3f-11eb-392d-5d188b60f5b3
## We can also pass a vector of NamedTuples to construct a DataFrame
DataFrame([(a=1, b=2, c=10), (a=3, b=4, c=6), (a=3, b=7, c=8)])

# ╔═╡ c7a51708-6a3f-11eb-0a86-ad52d4366887
## DataFrame from a matrix
DataFrame(rand(3, n), Symbol.('a':'e'))

# ╔═╡ dc5bc4ee-6a3f-11eb-07fa-21d227880d8e
## We can also create a DataFrame with no rows, but with predefined columns and their types:

DataFrame(A=Int[], B=Float64[], C=String[], D=Bool[], E=Symbol[])

# ╔═╡ fc951166-6a3f-11eb-0cce-71a959aed87a
## Finally, we can create a DataFrame by copying an existing DataFrame.
## Note that copy also copies the vectors.
begin
	xdf = DataFrame(a=1:n, b='a':'e')
	ydf = copy(xdf)

	## ===: same reference, isequal 2 different references, same content
	(xdf === ydf), isequal(xdf, ydf), (xdf.a == ydf.a), (xdf.a === ydf.a)
end

# ╔═╡ 2dc885e0-6a40-11eb-2bae-395fd78e1e4f
## Calling DataFrame on a DataFrame object works like copy.

begin
	df1 = DataFrame(a=1:n, b='a':'e')
	dfc1 = DataFrame(df1)
	(df1 === dfc1), isequal(df1, dfc1), (df1.a == dfc1.a), (df1.a === dfc1.a)
end

# ╔═╡ 65e148ba-6a40-11eb-3119-5faa428a40e9
## WE can avoid copying of columns of a data frame (if it is possible) by passing copycols=false keyword argument:

begin
	df2 = DataFrame(a=1:n, b='a':'e')
	dfc2 = DataFrame(df2, copycols=false)
	(df2 === dfc2), isequal(df2, dfc2), (df2.a == dfc2.a), (df2.a === dfc2.a)
end

# ╔═╡ aa41798c-6a40-11eb-32af-57cabd65363a
## Create a similar uninitialized DataFrame based on an original one:
ndf = DataFrame(a=1, b=1.0)

# ╔═╡ c8182d3e-6a40-11eb-014e-a344fba2bd93
similar(ndf)

# ╔═╡ d17d06d6-6a40-11eb-297d-7b94161b5187
## With number of rows explicitly stated
similar(ndf, 3)

# ╔═╡ f01bbc40-6a40-11eb-3c31-0fc5b5ea28da
## We can also create a new DataFrame from SubDataFrame or DataFrameRow
sdf = view(ndf, [1, 1], :)

# ╔═╡ 06efd550-6a41-11eb-0a51-6b7d64d82815
typeof(sdf)

# ╔═╡ 126e028a-6a41-11eb-1466-8f3c8c02143e
dfr = ndf[1, :]

# ╔═╡ 1b067fd0-6a41-11eb-3c17-777d1204480b
DataFrame(dfr)

# ╔═╡ 82c5cb20-6a3d-11eb-34bb-8fa036f13b8b
md"""
#### Conversion to matrix
"""

# ╔═╡ 82aca774-6a3d-11eb-3477-f91467210f50
x_df = DataFrame(x=1:2, y=["A", "B"])

# ╔═╡ 82954234-6a3d-11eb-138a-f3e2372d13ed
Matrix(x_df)

# ╔═╡ 827557d0-6a3d-11eb-1d97-857d6de9f26a
Array(x_df)

# ╔═╡ 82178042-6a3d-11eb-346e-0f1886c3856f
## This would work even if the DataFrame had some missings:

y_df = DataFrame(x=1:2, y=[missing,"B"])

# ╔═╡ 6161d65a-6a41-11eb-09fc-0b70e80f7ba1
Matrix(y_df)

# ╔═╡ 77dc3074-6a41-11eb-2a05-474caf947205
md"""
#### Conversion to `NamedTuple` related tabular structures
"""

# ╔═╡ 817775f8-6a41-11eb-0e93-1d1cfc09a00e
nt_df = DataFrame(x=1:n, y=["A", "B", "C", "D", "E"])

# ╔═╡ 8fd83650-6a41-11eb-2a0a-8727e31109a8
## Conversion into a NamedTuple of vectors
ct = Tables.columntable(nt_df)

# ╔═╡ a943ea9e-6a41-11eb-022e-8dc8a75562b3
## Conversion into a vector of NamedTuples
rt = Tables.rowtable(nt_df)

# ╔═╡ c3e077a0-6a41-11eb-3f3a-cd35f67ac77b
## Conversion back to DataFrame

DataFrame(ct)

# ╔═╡ d62e7e8c-6a41-11eb-226c-75489e440d63
DataFrame(rt)

# ╔═╡ e8c4ebdc-6a41-11eb-26c9-cf0060b276da
md"""
### Iterating data frame by rows or columns
"""

# ╔═╡ f53608f6-6a41-11eb-2bcf-fb9dbc15c8a0
ec = eachcol(xdf)

# ╔═╡ 00596ffa-6a42-11eb-36a0-f7ca0c1616ce
ec isa AbstractVector

# ╔═╡ 1ed33062-6a42-11eb-279a-a7de8bfd1436
ec[1]

# ╔═╡ 24d7a5ec-6a42-11eb-0743-997a55fa7e31
## Can also use index
ec[:a]

# ╔═╡ 4ea7ed00-6a42-11eb-3039-c904e26fdfd5
er = eachrow(xdf)

# ╔═╡ 5654c37a-6a42-11eb-1bd4-db7aacba9907
er isa AbstractVector

# ╔═╡ 58be76ba-6a42-11eb-040c-4d424ef55092
er[end]

# ╔═╡ 187c1633-6aa4-4e27-8707-cb8c24500273
html"""

"""

# ╔═╡ Cell order:
# ╟─0b08257e-6a3d-11eb-3da0-ed75be40d90e
# ╠═629d929c-6a3d-11eb-0413-5185929928cb
# ╟─2c93b1b1-39f7-49c3-8d19-060fac1dba52
# ╟─81665132-6a3d-11eb-37b7-991da3953909
# ╠═82e2153c-6a3d-11eb-0397-e78654843e7e
# ╠═b037eb42-6a3d-11eb-0443-c3fef979ac19
# ╠═18ff9300-6a3e-11eb-158b-5bec35b58c84
# ╟─6c87db40-6a3e-11eb-3a1c-7d3c95ffc793
# ╠═9c069eba-6a3e-11eb-0b99-69c73c8db3f8
# ╠═c7eb6274-6a3e-11eb-3846-e337d215f871
# ╠═eafc3138-6a3e-11eb-11e8-c333d5bf13c1
# ╠═d53123d6-6a3e-11eb-1304-c938d3ffc423
# ╠═5f83c5b6-6a3f-11eb-09ad-2d01fa5e3a6e
# ╠═a3dd623a-6a3f-11eb-392d-5d188b60f5b3
# ╠═c7a51708-6a3f-11eb-0a86-ad52d4366887
# ╠═dc5bc4ee-6a3f-11eb-07fa-21d227880d8e
# ╠═fc951166-6a3f-11eb-0cce-71a959aed87a
# ╠═2dc885e0-6a40-11eb-2bae-395fd78e1e4f
# ╠═65e148ba-6a40-11eb-3119-5faa428a40e9
# ╠═aa41798c-6a40-11eb-32af-57cabd65363a
# ╠═c8182d3e-6a40-11eb-014e-a344fba2bd93
# ╠═d17d06d6-6a40-11eb-297d-7b94161b5187
# ╠═f01bbc40-6a40-11eb-3c31-0fc5b5ea28da
# ╠═06efd550-6a41-11eb-0a51-6b7d64d82815
# ╠═126e028a-6a41-11eb-1466-8f3c8c02143e
# ╠═1b067fd0-6a41-11eb-3c17-777d1204480b
# ╟─82c5cb20-6a3d-11eb-34bb-8fa036f13b8b
# ╠═82aca774-6a3d-11eb-3477-f91467210f50
# ╠═82954234-6a3d-11eb-138a-f3e2372d13ed
# ╠═827557d0-6a3d-11eb-1d97-857d6de9f26a
# ╠═82178042-6a3d-11eb-346e-0f1886c3856f
# ╠═6161d65a-6a41-11eb-09fc-0b70e80f7ba1
# ╟─77dc3074-6a41-11eb-2a05-474caf947205
# ╠═817775f8-6a41-11eb-0e93-1d1cfc09a00e
# ╠═8fd83650-6a41-11eb-2a0a-8727e31109a8
# ╠═a943ea9e-6a41-11eb-022e-8dc8a75562b3
# ╠═c3e077a0-6a41-11eb-3f3a-cd35f67ac77b
# ╠═d62e7e8c-6a41-11eb-226c-75489e440d63
# ╟─e8c4ebdc-6a41-11eb-26c9-cf0060b276da
# ╠═f53608f6-6a41-11eb-2bcf-fb9dbc15c8a0
# ╠═00596ffa-6a42-11eb-36a0-f7ca0c1616ce
# ╠═1ed33062-6a42-11eb-279a-a7de8bfd1436
# ╠═24d7a5ec-6a42-11eb-0743-997a55fa7e31
# ╠═4ea7ed00-6a42-11eb-3039-c904e26fdfd5
# ╠═5654c37a-6a42-11eb-1bd4-db7aacba9907
# ╠═58be76ba-6a42-11eb-040c-4d424ef55092
# ╟─187c1633-6aa4-4e27-8707-cb8c24500273

max_stars_repo_name: pascal-p/julia-notebooks

max_issues_repo_name: pascal-p/julia-notebooks

max_forks_repo_name: pascal-p/julia-notebooks