text:
"""
Read/write a region of image data from multiple bands.

This method allows reading a region of one or more `GDALRasterBands` from this
dataset into a buffer, or writing data from a buffer into a region of the
`GDALRasterBands`. It automatically takes care of data type translation if the
data type (`eBufType`) of the buffer is different than that of the
`GDALRasterBand`. The method also takes care of image decimation / replication
if the buffer size (`nBufXSize x nBufYSize`) is different than the size of the
region being accessed (`nXSize x nYSize`).

The `pxspace`, `linespace` and `bandspace` parameters allow reading into
or writing from various organization of buffers.

For highest performance full resolution data access, read and write on \"block
boundaries\" as returned by `GetBlockSize()`, or use the `ReadBlock()` and
`WriteBlock()` methods.

### Parameters
* `access`      Either `GF_Read` to read a region of data, or `GF_Write` to
                write a region of data.
* `xoffset`     The pixel offset to the top left corner of the band to be
                accessed. This would be 0 to start from the left side.
* `yoffset`     The line offset to the top left corner of the region of the
                band to be accessed. This would be zero to start from the top.
* `xsize`       The width of the region of the band to be accessed in pixels.
* `ysize`       The height of the region of the band to be accessed in lines.
* `buffer`      The buffer into which the data should be read, or from which it
                should be written. It must contain
                    ≥`nBufXSize * nBufYSize * nBandCount`
                words of type `eBufType`. It is organized in left to right,
                top to bottom pixel order. Spacing is controlled by the
                `nPixelSpace`, and `nLineSpace` parameters
* `xsz`         The width of the buffer image into which the desired region is
                to be read, or from which it is to be written.
* `ysz`         The height of the buffer image into which the desired region is
                to be read, or from which it is to be written.
* `bands`       The list of bands (1 based) to be read/written.
* `pxspace`     The byte offset from the start of one pixel value in `pBuffer`
                to the start of the next pixel value within a scanline.
                If defaulted (0) the size of the datatype `eBufType` is used.
* `linespace`   The byte offset from the start of one scanline in pBuffer to
                the start of the next. If defaulted (0) the size of the datatype
                `eBufType * nBufXSize` is used.
* `bandspace`   The byte offset from the start of one bands data to the start
                of the next. If defaulted (0) the value will be
                    `nlinespace * nBufYSize`
                implying band sequential organization of the data buffer.
* `psExtraArg`  (new in GDAL 2.0) pointer to a GDALRasterIOExtraArg structure
with additional arguments to specify resampling and progress callback, or
`NULL` for default behaviour. The `GDAL_RASTERIO_RESAMPLING` configuration
option can also be defined to override the default resampling to one of
`BILINEAR`, `CUBIC`, `CUBICSPLINE`, `LANCZOS`, `AVERAGE` or `MODE`.
"""
function rasterio!(
        dataset::Dataset,
        buffer::Array{T, 3},
        bands::Vector{Cint},
        access::GDALRWFlag  = GDAL.GF_Read,
        pxspace::Integer    = 0,
        linespace::Integer  = 0,
        bandspace::Integer  = 0
    ) where T <: Real
    rasterio!(dataset, buffer, bands, 0, 0, width(dataset), height(dataset),
        access, pxspace, linespace, bandspace
    )
end

function rasterio!(
        dataset::Dataset,
        buffer::Array{T, 3},
        bands::Vector{Cint},
        rows::UnitRange{U},
        cols::UnitRange{U},
        access::GDALRWFlag  = GDAL.GF_Read,
        pxspace::Integer    = 0,
        linespace::Integer  = 0,
        bandspace::Integer  = 0
    ) where {T <: Real, U <: Integer}
    xsize = cols[end] - cols[1] + 1; xsize < 0 && error("invalid window width")
    ysize = rows[end] - rows[1] + 1; ysize < 0 && error("invalid window height")
    rasterio!(dataset, buffer, bands, cols[1], rows[1], xsize, ysize, access,
        pxspace, linespace, bandspace
    )
end

"""
Read/write a region of image data for this band.

This method allows reading a region of a `GDALRasterBand` into a buffer, or
writing data from a buffer into a region of a `GDALRasterBand`. It
automatically takes care of data type translation if the data type (`eBufType`)
of the buffer is different than that of the `GDALRasterBand`. The method also
takes care of image decimation / replication if the buffer size
`(nBufXSize x nBufYSize)` is different than the size of the region being
accessed `(nXSize x nYSize)`.

The `nPixelSpace` and `nLineSpace` parameters allow reading into or writing
from unusually organized buffers. This is primarily used for buffers containing
more than one bands raster data in interleaved format.

Some formats may efficiently implement decimation into a buffer by reading from
lower resolution overview images.

For highest performance full resolution data access, read and write on "block
boundaries" returned by `GetBlockSize()`, or use the `ReadBlock()` and
`WriteBlock()` methods.

### Parameters
* `eRWFlag`     Either GF_Read to read a region of data, or GF_Write to write a
region of data.
* `nXOff`       The pixel offset to the top left corner of the region of the
band to be accessed. This would be zero to start from the left side.
* `nYOff`       The line offset to the top left corner of the region of the
band to be accessed. This would be zero to start from the top.
* `nXSize`      The width of the region of the band to be accessed in pixels.
* `nYSize`      The height of the region of the band to be accessed in lines.
* `pData`       The buffer into which the data should be read, or from which it
should be written. This buffer must contain at least `(nBufXSize * nBufYSize)`
words of type `eBufType`. It is organized in left to right, top to bottom pixel
order. Spacing is controlled by the `nPixelSpace`, and `nLineSpace` parameters.
* `nBXSize`     The width of the buffer image into which the desired region is
to be read, or from which it is to be written.
* `nBYSize`     The height of the buffer image into which the desired region is
to be read, or from which it is to be written.
* `eBufType`    The type of the pixel values in the `buffer`. The pixel values
will be auto-translated to/from the `GDALRasterBand` data type as needed.
* `nPixelSpace` The byte offset from the start of one pixel value in `buffer`
to the start of the next pixel value within a scanline. If defaulted (0) the
size of the datatype `eBufType` is used.
* `nLineSpace`  The byte offset from the start of one scanline in `buffer` to
the start of the next. If defaulted (0) the size of the datatype
`(eBufType * nBufXSize)` is used.
* `psExtraArg`  (new in GDAL 2.0) pointer to a GDALRasterIOExtraArg structure
with additional arguments to specify resampling and progress callback, or
`NULL` for default behaviour. The `GDAL_RASTERIO_RESAMPLING` configuration
option can also be defined to override the default resampling to one of
`BILINEAR`, `CUBIC`, `CUBICSPLINE`, `LANCZOS`, `AVERAGE` or `MODE`.

### Returns
`CE_Failure` if the access fails, otherwise `CE_None`.
"""
function rasterio!(
        rasterband::RasterBand,
        buffer::Array{T,2},
        access::GDALRWFlag  = GDAL.GF_Read,
        pxspace::Integer    = 0,
        linespace::Integer  = 0
    ) where T <: Real
    rasterio!(rasterband, buffer, 0, 0, width(rasterband), height(rasterband),
        access, pxspace, linespace
    )
end

function rasterio!(
        rasterband::RasterBand,
        buffer::Array{T,2},
        rows::UnitRange{U},
        cols::UnitRange{U},
        access::GDALRWFlag  = GDAL.GF_Read,
        pxspace::Integer    = 0,
        linespace::Integer  = 0
    ) where {T <: Real, U <: Integer}
    xsize = length(cols); xsize < 1 && error("invalid window width")
    ysize = length(rows); ysize < 1 && error("invalid window height")
    rasterio!(rasterband, buffer, cols[1]-1, rows[1]-1, xsize, ysize,
        access, pxspace, linespace
    )
end

read!(rb::RasterBand, buffer::Array{T,2}) where {T <: Real} =
    rasterio!(rb, buffer, GDAL.GF_Read)

function read!(
        rb::RasterBand,
        buffer::Array{T,2},
        xoffset::Integer,
        yoffset::Integer,
        xsize::Integer,
        ysize::Integer
    ) where T <: Real
    rasterio!(rb, buffer, xoffset, yoffset, xsize, ysize)
end

function read!(
        rb::RasterBand,
        buffer::Array{T,2},
        rows::UnitRange{U},
        cols::UnitRange{U}
    ) where {T <: Real, U <: Integer}
    rasterio!(rb, buffer, rows, cols)
end

read(rb::RasterBand) =
    rasterio!(rb, Array{getdatatype(rb)}(undef, width(rb), height(rb)))

function read(
        rb::RasterBand,
        xoffset::Integer,
        yoffset::Integer,
        xsize::Integer,
        ysize::Integer
    )
    buffer = Array{getdatatype(rb)}(undef, width(rb), height(rb))
    rasterio!(rb, buffer, xoffset, yoffset, xsize, ysize)
end

function read(
        rb::RasterBand,
        rows::UnitRange{U},
        cols::UnitRange{U}
    ) where U <: Integer
    rasterio!(rb,
        Array{getdatatype(rb)}(undef, length(cols), length(rows)),
        rows, cols
    )
end

write!(rb::RasterBand, buffer::Array{T,2}) where {T <: Real} =
    rasterio!(rb, buffer, GDAL.GF_Write)

function write!(
        rb::RasterBand,
        buffer::Array{T, 2},
        xoffset::Integer,
        yoffset::Integer,
        xsize::Integer,
        ysize::Integer
    ) where T <: Real
    rasterio!(rb, buffer, xoffset, yoffset, xsize, ysize, GDAL.GF_Write)
end

function write!(
        rb::RasterBand,
        buffer::Array{T, 2},
        rows::UnitRange{U},
        cols::UnitRange{U}
    ) where {T <: Real, U <: Integer}
    rasterio!(rb, buffer, rows, cols, GDAL.GF_Write)
end

read!(dataset::Dataset, buffer::Array{T,2}, i::Integer) where {T <: Real} =
    read!(getband(dataset, i), buffer)

read!(dataset::Dataset, buffer::Array{T,3}, indices::Vector{Cint}) where {T <: Real} =
    rasterio!(dataset, buffer, indices, GDAL.GF_Read)

function read!(dataset::Dataset, buffer::Array{T,3}) where T <: Real
    nband = nraster(dataset); @assert size(buffer, 3) == nband
    rasterio!(dataset, buffer, collect(Cint, 1:nband), GDAL.GF_Read)
end

function read!(
        dataset::Dataset,
        buffer::Array{T, 2},
        i::Integer,
        xoffset::Integer,
        yoffset::Integer,
        xsize::Integer,
        ysize::Integer
    ) where T <: Real
    read!(getband(dataset, i), buffer, xoffset, yoffset, xsize, ysize)
end

function read!(
        dataset::Dataset,
        buffer::Array{T, 3},
        indices::Vector{Cint},
        xoffset::Integer,
        yoffset::Integer,
        xsize::Integer,
        ysize::Integer
    ) where T <: Real
    rasterio!(dataset, buffer, indices, xoffset, yoffset, xsize, ysize)
end

function read!(
        dataset::Dataset,
        buffer::Array{T, 2},
        i::Integer,
        rows::UnitRange{U},
        cols::UnitRange{U}
    ) where {T <: Real, U <: Integer}
    read!(getband(dataset, i), buffer, rows, cols)
end

function read!(
        dataset::Dataset,
        buffer::Array{T, 3},
        indices::Vector{Cint},
        rows::UnitRange{U},
        cols::UnitRange{U}
    ) where {T <: Real, U <: Integer}
    rasterio!(dataset, buffer, indices, rows, cols)
end

read(dataset::Dataset, i::Integer) = read(getband(dataset, i))

function read(dataset::Dataset, indices::Vector{Cint})
    buffer = Array{getdatatype(getband(dataset, indices[1]))}(
        undef, width(dataset), height(dataset), length(indices)
    )
    rasterio!(dataset, buffer, indices)
end

function read(dataset::Dataset)
    read!(dataset, Array{getdatatype(getband(dataset, 1))}(
        undef, width(dataset), height(dataset), nraster(dataset)
    ))
end

function read(
        dataset::Dataset,
        i::Integer,
        xoffset::Integer,
        yoffset::Integer,
        xsize::Integer,
        ysize::Integer
    )
    read(getband(dataset, i), xoffset, yoffset, xsize, ysize)
end

function read(
        dataset::Dataset,
        indices::Vector{T},
        xoffset::Integer,
        yoffset::Integer,
        xsize::Integer,
        ysize::Integer
    ) where T <: Integer
    buffer = Array{getdatatype(getband(dataset, indices[1]))}(
        undef, width(dataset), height(dataset), length(indices)
    )
    rasterio!(dataset, buffer, indices, xsize, ysize, xoffset, yoffset)
end

function read(
        dataset::Dataset,
        i::Integer,
        rows::UnitRange{U},
        cols::UnitRange{U}
    ) where U <: Integer
    read(getband(dataset, i), rows, cols)
end

function read(
        dataset::Dataset,
        indices::Vector{Cint},
        rows::UnitRange{U},
        cols::UnitRange{U}
    ) where U <: Integer
    buffer = Array{getdatatype(getband(dataset, indices[1]))}(
        undef, width(dataset), height(dataset), length(indices)
    )
    rasterio!(dataset, buffer, indices, rows, cols)
end

write!(dataset::Dataset, buffer::Array{T,2}, i::Integer) where {T <: Real} =
    write!(getband(dataset, i), buffer)

write!(dataset::Dataset, buffer::Array{T,3}, indices::Vector{Cint}) where {T <: Real} =
    rasterio!(dataset, buffer, indices, GDAL.GF_Write)

function write!(
        dataset::Dataset,
        buffer::Array{T, 2},
        i::Integer,
        xoffset::Integer,
        yoffset::Integer,
        xsize::Integer,
        ysize::Integer
    ) where T <: Real
    write!(getband(dataset, i), buffer, xoffset, yoffset, xsize, ysize)
end

function write!(
        dataset::Dataset,
        buffer::Array{T, 3},
        indices::Vector{Cint},
        xoffset::Integer,
        yoffset::Integer,
        xsize::Integer,
        ysize::Integer
    ) where T <: Real
    rasterio!(dataset, buffer, indices, xoffset, yoffset,
        xsize, ysize, GDAL.GF_Write
    )
end

function write!(
        dataset::Dataset,
        buffer::Array{T, 2},
        i::Integer,
        rows::UnitRange{U},
        cols::UnitRange{U}
    ) where {T <: Real, U <: Integer}
    write!(getband(dataset, i), buffer, rows, cols)
end

function write!(
        dataset::Dataset,
        buffer::Array{T, 3},
        indices::Vector{Cint},
        rows::UnitRange{U},
        cols::UnitRange{U}
    ) where {T <: Real, U <: Integer}
    rasterio!(dataset, buffer, indices, rows, cols, GDAL.GF_Write)
end

for (T,GT) in _GDALTYPE
    eval(quote
        function rasterio!(dataset::Dataset,
                           buffer::Array{$T, 3},
                           bands::Vector{Cint},
                           xoffset::Integer,
                           yoffset::Integer,
                           xsize::Integer,
                           ysize::Integer,
                           access::GDALRWFlag   = GDAL.GF_Read,
                           pxspace::Integer     = 0,
                           linespace::Integer   = 0,
                           bandspace::Integer   = 0,
                           extraargs = Ptr{GDAL.GDALRasterIOExtraArg}(C_NULL))
            (dataset == C_NULL) && error("Can't read invalid rasterband")
            xbsize, ybsize, zbsize = size(buffer)
            nband = length(bands); @assert nband == zbsize
            result = ccall((:GDALDatasetRasterIOEx,GDAL.libgdal),GDAL.CPLErr,
                           (GDALDataset,GDAL.GDALRWFlag,Cint,Cint,Cint,Cint,
                            Ptr{Cvoid},Cint,Cint,GDAL.GDALDataType,Cint,
                            Ptr{Cint},GDAL.GSpacing,GDAL.GSpacing,GDAL.GSpacing,
                            Ptr{GDAL.GDALRasterIOExtraArg}),dataset.ptr,access,
                            xoffset,yoffset,xsize,ysize,pointer(buffer),xbsize,
                            ybsize,$GT,nband,pointer(bands),pxspace,linespace,
                            bandspace,extraargs)
            @cplerr result "Access in DatasetRasterIO failed."
            buffer
        end

        function rasterio!(rasterband::RasterBand,
                           buffer::Array{$T,2},
                           xoffset::Integer,
                           yoffset::Integer,
                           xsize::Integer,
                           ysize::Integer,
                           access::GDALRWFlag   = GDAL.GF_Read,
                           pxspace::Integer     = 0,
                           linespace::Integer   = 0,
                           extraargs = Ptr{GDAL.GDALRasterIOExtraArg}(C_NULL))
            (rasterband == C_NULL) && error("Can't read invalid rasterband")
            xbsize, ybsize = size(buffer)
            result = ccall((:GDALRasterIOEx,GDAL.libgdal),GDAL.CPLErr,
                           (GDALRasterBand,GDAL.GDALRWFlag,Cint,Cint,Cint,Cint,
                            Ptr{Cvoid},Cint,Cint,GDAL.GDALDataType,GDAL.GSpacing,
                            GDAL.GSpacing,Ptr{GDAL.GDALRasterIOExtraArg}),
                            rasterband.ptr,access,xoffset,yoffset,xsize,ysize,
                            pointer(buffer),xbsize,ybsize,$GT,pxspace,linespace,
                            extraargs)
            @cplerr result "Access in RasterIO failed."
            buffer
        end
    end)
end

"""
Read a block of image data efficiently.

This method accesses a "natural" block from the raster band without resampling,
or data type conversion. For a more generalized, but potentially less efficient
access use RasterIO().

### Parameters
* `xoffset` the horizontal block offset, with zero indicating the left most
            block, 1 the next block and so forth.
* `yoffset` the vertical block offset, with zero indicating the top most block,
            1 the next block and so forth.
* `buffer`  the buffer into which the data will be read. The buffer must be
            large enough to hold GetBlockXSize()*GetBlockYSize() words of type
            GetRasterDataType().
"""
function readblock!(rb::RasterBand, xoffset::Integer, yoffset::Integer, buffer)
    result = GDAL.readblock(rb.ptr, xoffset, yoffset, buffer)
    @cplerr result "Failed to read block at ($xoffset,$yoffset)"
    buffer
end

"""
Write a block of image data efficiently.

This method accesses a "natural" block from the raster band without resampling,
or data type conversion. For a more generalized, but potentially less efficient
access use RasterIO().

### Parameters
* `xoffset` the horizontal block offset, with zero indicating the left most
            block, 1 the next block and so forth.
* `yoffset` the vertical block offset, with zero indicating the left most block,
            1 the next block and so forth.
* `buffer`  the buffer from which the data will be written. The buffer must be
            large enough to hold GetBlockXSize()*GetBlockYSize() words of type
            GetRasterDataType().
"""
function writeblock!(rb::RasterBand, xoffset::Integer, yoffset::Integer, buffer)
    result = GDAL.writeblock(rb.ptr, xoffset, yoffset, buffer)
    @cplerr result "Failed to write block at ($xoffset,$yoffset)"
    buffer
end


text:
#using Optim # numerical estimation routines moved to NormalTransforms

"""
    LogitNormal(μ,σ)

The *logit normal distribution* is the distribution of
of a random variable whose logit has a [`Normal`](@ref) distribution.
Or inversely, when applying the logistic function to a Normal random variable
then the resulting random variable follows a logit normal distribution.

If ``X \\sim \\operatorname{Normal}(\\mu, \\sigma)`` then
``\\operatorname{logistic}(X) \\sim \\operatorname{LogitNormal}(\\mu,\\sigma)``.

The probability density function is

```math
f(x; \\mu, \\sigma) = \\frac{1}{x \\sqrt{2 \\pi \\sigma^2}}
\\exp \\left( - \\frac{(\\text{logit}(x) - \\mu)^2}{2 \\sigma^2} \\right),
\\quad x > 0
```

where the logit-Function is
```math
\\text{logit}(x) = \\ln\\left(\\frac{x}{1-x}\\right)
\\quad 0 < x < 1
```

```julia
LogitNormal()          # Logit-normal distribution with zero logit-mean and unit scale
LogitNormal(mu)        # Logit-normal distribution with logit-mean mu and unit scale
LogitNormal(mu, sig)   # Logit-normal distribution with logit-mean mu and scale sig

params(d)            # Get the parameters, i.e. (mu, sig)
median(d)            # Get the median, i.e. logistic(mu)
```

The following properties have no analytical solution but numerical
approximations. In order to avoid package dependencies for
numerical optimization, they are currently not implemented.

```julia
mean(d)
var(d)
std(d)
mode(d)
```

Similarly, skewness, kurtosis, and entropy are not implemented.

External links

* [Logit normal distribution on Wikipedia](https://en.wikipedia.org/wiki/Logit-normal_distribution)

"""
struct LogitNormal{T<:Real} <: ContinuousUnivariateDistribution
    μ::T
    σ::T
    LogitNormal{T}(μ::T, σ::T) where {T} = new{T}(μ, σ)
end

function LogitNormal(μ::T, σ::T; check_args=true) where {T <: Real}
    check_args && @check_args(LogitNormal, σ > zero(σ))
    return LogitNormal{T}(μ, σ)
end

LogitNormal(μ::Real, σ::Real) = LogitNormal(promote(μ, σ)...)
LogitNormal(μ::Integer, σ::Integer) = LogitNormal(float(μ), float(σ))
LogitNormal(μ::T) where {T} = LogitNormal(μ, one(T))
LogitNormal() = LogitNormal(0.0, 1.0, check_args=false)

# minimum and maximum not defined for logitnormal
# but see https://github.com/JuliaStats/Distributions.jl/pull/457
@distr_support LogitNormal 0.0 1.0
#insupport(d::Union{D,Type{D}},x::Real) where {D<:LogitNormal} = 0.0 < x < 1.0

#### Conversions
convert(::Type{LogitNormal{T}}, μ::S, σ::S) where
  {T <: Real, S <: Real} = LogitNormal(T(μ), T(σ))
convert(::Type{LogitNormal{T}}, d::LogitNormal{S}) where
  {T <: Real, S <: Real} = LogitNormal(T(d.μ), T(d.σ), check_args=false)

#### Parameters

params(d::LogitNormal) = (d.μ, d.σ)
location(d::LogitNormal) = d.μ
scale(d::LogitNormal) = d.σ
@inline partype(d::LogitNormal{T}) where {T<:Real} = T

#### Statistics

# meanlogitx(d::LogitNormal) = d.μ
# varlogitx(d::LogitNormal) = abs2(d.σ)
# stdlogitx(d::LogitNormal) = d.σ

# mean, mode, and variance
# moved to NormalTransforms because
# they depend on the Optim package.

# skewness, kurtosis, entropy
# not implemented

median(d::LogitNormal) = logistic(d.μ)

#### Evalution

#TODO check pd and logpdf
function pdf(d::LogitNormal{T}, x::Real) where T<:Real
    if zero(x) < x < one(x)
        return normpdf(d.μ, d.σ, logit(x)) / (x * (1-x))
    else
        return T(0)
    end
end

function logpdf(d::LogitNormal{T}, x::Real) where T<:Real
    if zero(x) < x < one(x)
        lx = logit(x)
        return normlogpdf(d.μ, d.σ, lx) - log(x) - log1p(-x)
    else
        return -T(Inf)
    end
end

cdf(d::LogitNormal{T}, x::Real) where {T<:Real} =
    x ≤ 0 ? zero(T) : x ≥ 1 ? one(T) : normcdf(d.μ, d.σ, logit(x))
ccdf(d::LogitNormal{T}, x::Real) where {T<:Real} =
    x ≤ 0 ? one(T) : x ≥ 1 ? zero(T) : normccdf(d.μ, d.σ, logit(x))
logcdf(d::LogitNormal{T}, x::Real) where {T<:Real} =
    x ≤ 0 ? -T(Inf) : x ≥ 1 ? zero(T) : normlogcdf(d.μ, d.σ, logit(x))
logccdf(d::LogitNormal{T}, x::Real) where {T<:Real} =
    x ≤ 0 ? zero(T) : x ≥ 1 ? -T(Inf) : normlogccdf(d.μ, d.σ, logit(x))

quantile(d::LogitNormal, q::Real) = logistic(norminvcdf(d.μ, d.σ, q))
cquantile(d::LogitNormal, q::Real) = logistic(norminvccdf(d.μ, d.σ, q))
invlogcdf(d::LogitNormal, lq::Real) = logistic(norminvlogcdf(d.μ, d.σ, lq))
invlogccdf(d::LogitNormal, lq::Real) = logistic(norminvlogccdf(d.μ, d.σ, lq))

function gradlogpdf(d::LogitNormal{T}, x::Real) where T<:Real
    #TODO check
    (μ, σ) = params(d)
    0 < x < 1 ? - ((log(x) - μ) / ((σ^2) + 1) * x * (1-x)) : zero(T)
end
function heslogpdf(d::LogitNormal{T}, x::Real) where T<:Real
    #TODO check
    (μ, σ) = params(d)
    0 < x < 1 ? ((2 * x - 1) * log(x) + (1 - 2 * μ) * x + μ - 1) / (σ^2 + 1) : zero(T)
end

# mgf(d::LogitNormal)
# cf(d::LogitNormal)

#### Sampling

rand(rng::AbstractRNG, d::LogitNormal) = logistic(randn(rng) * d.σ + d.μ)

## Fitting

function fit_mle(::Type{<:LogitNormal}, x::AbstractArray{T}) where T<:Real
    lx = logit.(x)
    μ, σ = mean_and_std(lx)
    LogitNormal(μ, σ)
end

ext: jl


text:
@testset "Disk" begin

D = Disk()
@test D.id == 1
@test D.rmin ≈ 0
@test D.rmax ≈ 1
@test D.inverted == false
@test GRAYTR.area(D) ≈ π
@test GRAYTR.can_intersect(D)

# Create with transformation
D = Disk(GRAYTR.scaling(2.0, 1.0, 1.0))
@test GRAYTR.obj_bounds(D) == GRAYTR.BoundingBox(Point3(-1, -1, 0), Point3(1, 1, 0))
@test GRAYTR.world_bounds(D) == GRAYTR.BoundingBox(Point3(-2, -1, 0), Point3(2, 1, 0))

# Transform after creation
T = GRAYTR.scaling(1.0, 2.0, 1.0)
D = T(Disk())
@test GRAYTR.obj_bounds(D) == GRAYTR.BoundingBox(Point3(-1, -1, 0), Point3(1, 1, 0))
@test GRAYTR.world_bounds(D) == GRAYTR.BoundingBox(Point3(-1, -2, 0), Point3(1, 2, 0))

# intersectP
D = Disk(translation(0, 0, -1))
R = GRAYTR.Ray(Point3(0, 0, 1), Vector3(0, 0, -1))
@test GRAYTR.intersectP(R, D)

R = GRAYTR.Ray(Point3(0, 0, 1), Vector3(0, 0, -1), 10, Inf, 1)
@test !GRAYTR.intersectP(R, D)

R = GRAYTR.Ray(Point3(0, 2, 1), Vector3(0, 0, -1))
@test !GRAYTR.intersectP(R, D)

R = GRAYTR.Ray(Point3(0.9*sqrt(2), 0.9*sqrt(2), 1), Vector3(0, 0, -1))
@test !GRAYTR.intersectP(R, D)

R = GRAYTR.Ray(Point3(0, 0, -1), Vector3(0, 0, 1))
@test GRAYTR.intersectP(R, D)

R = GRAYTR.Ray(Point3(0, 0, 0), Vector3(0, 1, 0))
@test !GRAYTR.intersectP(R, D)

# intersect

R = GRAYTR.Ray(Point3(0, 0, 1), Vector3(0, 0, -1))
dg, t, e = GRAYTR.shape_intersect(R, D)
@test dg.p ≈ Point3(0, 0, -1)
@test dg.n ≈ Normal3(0, 0, 1)
@test t ≈ 2.0

R = GRAYTR.Ray(Point3(0.9*sqrt(2), 0.9*sqrt(2), 1), Vector3(0, 0, -1))
dg, t, e = GRAYTR.shape_intersect(R, D)
@test dg == nothing
@test isnan(t)
@test isnan(e)

R = GRAYTR.Ray(Point3(1, 0, 0), Vector3(-1, 0, 0))
dg, t, e = GRAYTR.shape_intersect(R, D)
@test dg == nothing
@test isnan(t)
@test isnan(e)

R = GRAYTR.Ray(Point3(0, 0, 1), Vector3(0, 0, -1), 10, Inf, 1)
dg, t, e = GRAYTR.shape_intersect(R, D)
@test dg == nothing
@test isnan(t)
@test isnan(e)

end

max_stars_repo_name: dronir/GRAYTR

max_issues_repo_name: dronir/GRAYTR

max_forks_repo_name: dronir/GRAYTR


text:
"""
    function length_conversion(value, from_type, to_type)

A function that converts a value from a measurement unit to another one

Accepted units are: millimeter(s), centimeter(s), meter(s), kilometer(s),
inch(es), feet, foot, yard(s), mile(s). Abbreviations are also supported.

# Examples/Tests (optional but recommended)
```julia
length_conversion(10, "METERS", "cm")
length_conversion(12345, "yards", "FEET")
```

Because the algorithm converts values to meters, and from meters to the final type,
some imperial system units may have errors:
length_conversion(1, "yards", "FEET") returns 3.000000096, intead of 3.

# Contributed by:- [Fernanda Kawasaki](https://github.com/fernandakawasaki)
"""

# Lookup table that returns conversion of 1 unit of type to meters
METER_CONVERSION = Dict{String,Float64}(
    "mm" => 0.001,
    "cm" => 0.01,
    "m" => 1,
    "km" => 1000,
    "in" => 0.0254,
    "ft" => 0.3048,
    "yd" => 0.9144,
    "mi" => 1609.34,
)

# Lookup table that returns the conversion of 1 meter to type
TYPE_CONVERSION = Dict{String,Float64}(
    "mm" => 1000,
    "cm" => 100,
    "m" => 1,
    "km" => 0.001,
    "in" => 39.3701,
    "ft" => 3.28084,
    "yd" => 1.09361,
    "mi" => 0.000621371,
)

NAME_CONVERSION = Dict{String,String}(
    "millimeter" => "mm",
    "millimeters" => "mm",
    "centimeter" => "cm",
    "centimeters" => "cm",
    "meter" => "m",
    "meters" => "m",
    "kilometer" => "km",
    "kilometers" => "km",
    "inch" => "in",
    "inches" => "in",
    "feet" => "ft",
    "foot" => "ft",
    "yard" => "yd",
    "yards" => "yd",
    "mile" => "mi",
    "miles" => "mi",
)

function normalize_type(type)
    l_type = lowercase(type)
    if !haskey(METER_CONVERSION, l_type)
        if !haskey(NAME_CONVERSION, l_type)
            throw(
                error(
                    "Invalid 'type' value: $(type)\n",
                    "Supported values are: $(keys(NAME_CONVERSION))\n",
                    "Supported abbreviations are: $(keys(METER_CONVERSION))\n",
                ),
            )
        end
        return NAME_CONVERSION[l_type]
    end
    return l_type
end

function length_conversion(value, from_type, to_type)
    from_type_norm = normalize_type(from_type)
    to_type_norm = normalize_type(to_type)
    value_in_meters = value * METER_CONVERSION[from_type_norm]
    return value_in_meters * TYPE_CONVERSION[to_type_norm]
end

max_issues_repo_name: AugustoCL/Julia

max_forks_repo_name: AugustoCL/Julia


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.
module Multilens

using LinearAlgebra
import Unitful
using Unitful:uconvert,ustrip
using Unitful.DefaultSymbols:mm,μm,nm #Unitful.DefaultSymbols exports a variable T which can cause major confusion with type declarations since T is a common parametric type symbol. Import only what is needed.
using Colors
using StaticArrays
import DataFrames
import ..Repeat
import SpecialFunctions
import Plots
import ...OpticSim
using ...OpticSim:plane_from_points,centroid,pointonplane,focallength,ParaxialLens
import ...OpticSim.Geometry
using ...OpticSim.Repeat:region,HexBasis1,HexBasis3,LatticeBasis,LatticeCluster,clustersize,ClusterWithProperties,cluster_coordinates_from_tile_coordinates, AbstractLatticeCluster,elementbasis,euclideandiameter
using ...OpticSim.Data

include("HexClusters.jl")
include("HexTilings.jl")
include("Analysis.jl")
include("DisplayGeneration.jl")
include("LensletAssignment.jl")
include("Example.jl")

end # module
export Multilens


text:
# minimal working example for fitting a sem

using LinearAlgebra, Optim, Zygote, ReverseDiff, AutoGrad, Yota,
    FiniteDiff, BenchmarkTools

pars = [1.0, 2]
A = rand(2,2); A = A*A'
pre = copy(A)
prod = copy(A)
a = cholesky(A)

function floss(pars, A, a, pre, prod)
    A[1] = pars[2]
    A[2] = pars[1]
    A[3] = pars[1]
    A[4] = pars[2]
    a = cholesky!(A; check = false)
    if !isposdef(a) return Inf end
    pre = LinearAlgebra.inv(a)
    #mul!(prod, pre, A)
    F = logdet(a) + tr(pre)
end

function floss2(pars, A, a, pre, prod)
    A = [pars[2] pars[1]
        pars[1] pars[2]]
    B = copy(A)
    a = cholesky!(A)
    pre = LinearAlgebra.inv(a)
    prod = pre*B
    F = logdet(a) + tr(prod)
end

floss(pars, A, a, pre, prod)
floss2(pars, A, a, pre, prod)

@benchmark floss(pars, A, a, pre, prod)

# reverse diff
# Zygote; ReverseDiff; AutoGrad; Yota

# FiniteDiff
fgrad =
    FiniteDiff.finite_difference_jacobian(
        x -> floss2(x, A, a, pre, prod), pars)

# Zygote
zgrad =
    Zygote.gradient(x -> floss2(x, A, a, pre, prod), pars)

x = Param([1,2,3])		# user declares parameters
x => P([1,2,3])			# they are wrapped in a struct
value(x) => [1,2,3]		# we can get the original value
sum(abs2,x) => 14		# they act like regular values outside of differentiation
y = @diff sum(abs2,x)	        # if you want the gradients
y => T(14)			# you get another struct
value(y) => 14			# which represents the same value
grad(y,x) => [2,4,6]

a_pars = Param(pars)
agrad = @diff floss2(a_pars, A, a, pre, prod)

#Yota
ygrad =
    Yota.grad(x -> floss2(x, A, a, pre, prod), pars)

# ReverseDiff
# pre-record a GradientTape for `f` using inputs of shape 100x100 with Float64 elements
const floss_tape =
    ReverseDiff.GradientTape(x -> floss(x, A, a, pre, prod), (pars))

const floss_tape22 =
    ReverseDiff.GradientTape(x -> floss2(x, A, a, pre, prod), (pars))
# compile `f_tape` into a more optimized representation

const compiled_f_tape = ReverseDiff.compile(floss_tape)
const compiled_f_tape22 = ReverseDiff.compile(floss_tape22)

results = similar(pars)
ReverseDiff.gradient!(results, compiled_f_tape, pars)
ReverseDiff.gradient!(results, compiled_f_tape22, pars)

### => Yota and AutoGrad give errors

#########################
function floss_r(pars, A, a, pre, prod)
    A = [pars[2] pars[1]
        pars[1] pars[2]]
    B = copy(A)
    a = cholesky!(A)
    pre = LinearAlgebra.inv(a)
    mul!(prod, pre, B)
    F = logdet(a) + tr(prod)
end

function floss_z(pars, A, a, pre, prod)
    A = [pars[2] pars[1]
        pars[1] pars[2]]
    B = copy(A)
    a = cholesky(A)
    pre = LinearAlgebra.inv(a)
    prod = pre*B
    F = logdet(a) + tr(prod)
end

const floss_tape22222 =
    ReverseDiff.GradientTape(x -> floss_r(x, A, a, pre, prod), (pars))
# compile `f_tape` into a more optimized representation
const compiled_f_tape22222 = ReverseDiff.compile(floss_tape22222)

results = similar(pars)

ReverseDiff.gradient!(results, compiled_f_tape22222, pars)

zgrad =
    Zygote.gradient(x -> floss_z(x, A, a, pre, prod), pars)

#########################
function floss_r(pars, A, a, pre, prod)
    A = [pars[2] pars[1]
        pars[1] pars[2]]
    #B = copy(A)
    a = cholesky!(A; check = false)
    if !isposdef(a) return Inf end
    pre = LinearAlgebra.inv(a)
    #mul!(prod, pre, B)
    F = logdet(a) + tr(pre)
end

@btime optimize(x -> floss_r(x, A, a, pre, prod), pars, BFGS())

fgrad =
    FiniteDiff.finite_difference_jacobian(
        x -> floss_r(x, A, a, pre, prod), pars)

const floss_t =
    ReverseDiff.GradientTape(x -> floss_r(x, A, a, pre, prod), (pars))
# compile `f_tape` into a more optimized representation
const compiled_f_t= ReverseDiff.compile(floss_t)

results = similar(pars)

g(G, x) = ReverseDiff.gradient!(G, compiled_f_t, x)

@btime optimize(x -> floss(x, A, a, pre, prod), g, pars, BFGS())

@benchmark g(results, pars)
@benchmark floss_r(pars, A, a, pre, prod)

function floss_z(pars, A, a, pre, prod)
    A = [pars[2] pars[1]
        pars[1] pars[2]]
    #B = copy(A)
    if !isposdef(A) return Inf end
    pre = LinearAlgebra.inv(A)
    #prod = pre*B
    F = logdet(A) + tr(pre)
end

floss_z(pars, A, a, pre, prod)

zgrad =
    Zygote.gradient(x -> floss_z(x, A, a, pre, prod), pars)

function g_z(G, x)
    check = Zygote.gradient(x -> floss_z(x, A, a, pre, prod), x)[1]
    isnothing(check) ? G .= 0.0 : G.= check[1]
end

@btime optimize(x -> floss(x, A, a, pre, prod), g_z, pars, BFGS())

### with build function

using ModelingToolkit

@variables v[1:2]
A = [v[2] v[1]
    v[1] v[2]]

b, b_ = eval.(ModelingToolkit.build_function(A, v))

b(pars)
B = zeros(2,2)
b_(B, pars)

function floss_r(pars, A, a, pre, prod)
    A = b(pars)
    #B = copy(A)
    a = cholesky!(Hermitian(A); check = false)
    if !isposdef(a) return Inf end
    ld = logdet(a)
    pre = LinearAlgebra.inv(a)
    #mul!(prod, pre, B)
    F = ld + tr(pre)
end

floss_r(pars, B, a, pre, prod)

const floss_t_5 =
    ReverseDiff.GradientTape(x -> floss_r(x, B, a, pre, prod), (pars))
# compile `f_tape` into a more optimized representation
const compiled_f_t_4 = ReverseDiff.compile(floss_t_4)

results = similar(pars)
ReverseDiff.gradient!(results, compiled_f_t_4, pars)

#### Test generated function vs matrix multiplication

#matrix

function fmat(x)
    S =[x[1]  0     0     0     0     0     0     0     0     0     0     0     0     0
    0     x[2]  0     0     0     0     0     0     0     0     0     0     0     0
    0     0     x[3]  0     0     0     0     0     0     0     0     0     0     0
    0     0     0     x[4]  0     0     0     x[15] 0     0     0     0     0     0
    0     0     0     0     x[5]  0     x[16] 0     x[17] 0     0     0     0     0
    0     0     0     0     0     x[6]  0     0     0     x[18] 0     0     0     0
    0     0     0     0     x[16] 0     x[7]  0     0     0     x[19] 0     0     0
    0     0     0     x[15] 0     0     0     x[8]  0     0     0     0     0     0
    0     0     0     0     x[17] 0     0     0     x[9]  0     x[20] 0     0     0
    0     0     0     0     0     x[18] 0     0     0     x[10] 0     0     0     0
    0     0     0     0     0     0     x[19] 0     x[20] 0     x[11] 0     0     0
    0     0     0     0     0     0     0     0     0     0     0     x[12] 0     0
    0     0     0     0     0     0     0     0     0     0     0     0     x[13] 0
    0     0     0     0     0     0     0     0     0     0     0     0     0     x[14]]

    F =[1.0 0 0 0 0 0 0 0 0 0 0 0 0 0
        0 1 0 0 0 0 0 0 0 0 0 0 0 0
        0 0 1 0 0 0 0 0 0 0 0 0 0 0
        0 0 0 1 0 0 0 0 0 0 0 0 0 0
        0 0 0 0 1 0 0 0 0 0 0 0 0 0
        0 0 0 0 0 1 0 0 0 0 0 0 0 0
        0 0 0 0 0 0 1 0 0 0 0 0 0 0
        0 0 0 0 0 0 0 1 0 0 0 0 0 0
        0 0 0 0 0 0 0 0 1 0 0 0 0 0
        0 0 0 0 0 0 0 0 0 1 0 0 0 0
        0 0 0 0 0 0 0 0 0 0 1 0 0 0]

    A =[0  0  0  0  0  0  0  0  0  0  0     1     0     0
        0  0  0  0  0  0  0  0  0  0  0     x[21] 0     0
        0  0  0  0  0  0  0  0  0  0  0     x[22] 0     0
        0  0  0  0  0  0  0  0  0  0  0     0     1     0
        0  0  0  0  0  0  0  0  0  0  0     0     x[23] 0
        0  0  0  0  0  0  0  0  0  0  0     0     x[24] 0
        0  0  0  0  0  0  0  0  0  0  0     0     x[25] 0
        0  0  0  0  0  0  0  0  0  0  0     0     0     1
        0  0  0  0  0  0  0  0  0  0  0     0     0     x[26]
        0  0  0  0  0  0  0  0  0  0  0     0     0     x[27]
        0  0  0  0  0  0  0  0  0  0  0     0     0     x[28]
        0  0  0  0  0  0  0  0  0  0  0     0     0     0
        0  0  0  0  0  0  0  0  0  0  0     x[29] 0     0
        0  0  0  0  0  0  0  0  0  0  0     x[30] x[31] 0]

    invia = inv(I-A)
    implied = F*invia*S*invia'*F'
    return logdet(implied)
end

# generated
using ModelingToolkit, SparseArrays, sem

@ModelingToolkit.variables x[1:31]

S =[x[1]  0     0     0     0     0     0     0     0     0     0     0     0     0
    0     x[2]  0     0     0     0     0     0     0     0     0     0     0     0
    0     0     x[3]  0     0     0     0     0     0     0     0     0     0     0
    0     0     0     x[4]  0     0     0     x[15] 0     0     0     0     0     0
    0     0     0     0     x[5]  0     x[16] 0     x[17] 0     0     0     0     0
    0     0     0     0     0     x[6]  0     0     0     x[18] 0     0     0     0
    0     0     0     0     x[16] 0     x[7]  0     0     0     x[19] 0     0     0
    0     0     0     x[15] 0     0     0     x[8]  0     0     0     0     0     0
    0     0     0     0     x[17] 0     0     0     x[9]  0     x[20] 0     0     0
    0     0     0     0     0     x[18] 0     0     0     x[10] 0     0     0     0
    0     0     0     0     0     0     x[19] 0     x[20] 0     x[11] 0     0     0
    0     0     0     0     0     0     0     0     0     0     0     x[12] 0     0
    0     0     0     0     0     0     0     0     0     0     0     0     x[13] 0
    0     0     0     0     0     0     0     0     0     0     0     0     0     x[14]]

F =[1.0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 1 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 1 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 1 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 1 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 1 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 1 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 1 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 1 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 1 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 1 0 0 0]

A =[0  0  0  0  0  0  0  0  0  0  0     1     0     0
    0  0  0  0  0  0  0  0  0  0  0     x[21] 0     0
    0  0  0  0  0  0  0  0  0  0  0     x[22] 0     0
    0  0  0  0  0  0  0  0  0  0  0     0     1     0
    0  0  0  0  0  0  0  0  0  0  0     0     x[23] 0
    0  0  0  0  0  0  0  0  0  0  0     0     x[24] 0
    0  0  0  0  0  0  0  0  0  0  0     0     x[25] 0
    0  0  0  0  0  0  0  0  0  0  0     0     0     1
    0  0  0  0  0  0  0  0  0  0  0     0     0     x[26]
    0  0  0  0  0  0  0  0  0  0  0     0     0     x[27]
    0  0  0  0  0  0  0  0  0  0  0     0     0     x[28]
    0  0  0  0  0  0  0  0  0  0  0     0     0     0
    0  0  0  0  0  0  0  0  0  0  0     x[29] 0     0
    0  0  0  0  0  0  0  0  0  0  0     x[30] x[31] 0]

S = sparse(S)

#F
F = sparse(F)

#A
A = sparse(A)

invia = sem.neumann_series(A)

impcov = F*invia*S*permutedims(invia)*permutedims(F)

impcov = Array(impcov)

impcov = simplify.(impcov)

f_ = eval(ModelingToolkit.build_function(impcov, x)[1])

function fgen(par)
    impcov = f_(par)
    return tr(impcov)
end

### test

using BenchmarkTools

start_val = vcat(
    fill(1.0, 11),
    fill(0.05, 3),
    fill(0.0, 6),
    fill(1.0, 8),
    fill(0, 3)
    )

@btime fmat(start_val)

@btime fgen(start_val)

const mat_t =
    ReverseDiff.GradientTape(x -> fmat(x), (start_val))
# compile `f_tape` into a more optimized representation
const compiled_mat_t= ReverseDiff.compile(mat_t)

results = similar(start_val)

gmat(G, x) = ReverseDiff.gradient!(G, compiled_mat_t, x)

const gen_t3 =
    ReverseDiff.GradientTape(x -> fgen(x), (start_val))
# compile `f_tape` into a more optimized representation
const compiled_gen_t2 = ReverseDiff.compile(gen_t3)

results = similar(start_val)

ggen(G, x) = ReverseDiff.gradient!(G, compiled_gen_t2, x)

@btime ggen(results, start_val)

@btime gmat(results, start_val)

using ForwardDiff

@btime ForwardDiff.gradient(fgen, start_val)

FiniteDiff.finite_difference_gradient(fgen, start_val)

B = rand(100,100)

@btime inv(B)

#193.399 μs

Threads.nthreads()

@sync for i in 1:length(imply)
    @spawn for j in 1:length(imply[i] -> length(observed))
        F[k] = model[k](parameter)
    end
    @spawn for j in 1:length(imply[i] -> length(observed))
    end
    @spawn for j in 1:length(imply[i] -> length(observed))
    end
    @spawn for j in 1:length(imply[i] -> length(observed))
    end
end

rand(10,10)


text:
# This file was generated, do not modify it. # hide
tree = machine(DecisionTreeRegressor(), X, y)
e = evaluate!(tree, resampling=Holdout(fraction_train=0.8),
              measure=[rms, rmslp1])
e |> pprint # use PrettyPrinting


text:
export coulomb

using OffsetArrays

#using StaticArrays # Surprisingly, didn't speed up at all
"""
    coulomb(aexpn,ax,ay,az,aI,aJ,aK,
        bexpn,bx,by,bz,bI,bJ,bK,
        cexpn,cx,cy,cz,cI,cJ,cK,
        dexpn,dx,dy,dz,dI,dJ,dK)

Compute the coulomb repulsion between four primitive Gaussian basis
functions defined by their exponents, xyz coordinates, and IJK powers.
"""
function coulomb(aexpn,ax,ay,az,aI,aJ,aK,
    bexpn,bx,by,bz,bI,bJ,bK,
    cexpn,cx,cy,cz,cI,cJ,cK,
    dexpn,dx,dy,dz,dI,dJ,dK)

    rab2 = dist2(ax-bx,ay-by,az-bz)
    rcd2 = dist2(cx-dx,cy-dy,cz-dz)

    px,py,pz = gaussian_product_center(aexpn,ax,ay,az,bexpn,bx,by,bz)
    qx,qy,qz = gaussian_product_center(cexpn,cx,cy,cz,dexpn,dx,dy,dz)
    rpq2 = dist2(px-qx,py-qy,pz-qz)

    g1 = aexpn+bexpn
    g2 = cexpn+dexpn
    delta = 0.25*(1/g1+1/g2)

    Bx = Barray(aI,bI,cI,dI,px,ax,bx,qx,cx,dx,g1,g2,delta)
    By = Barray(aJ,bJ,cJ,dJ,py,ay,by,qy,cy,dy,g1,g2,delta)
    Bz = Barray(aK,bK,cK,dK,pz,az,bz,qz,cz,dz,g1,g2,delta)

    x = 0.25*rpq2/delta

    s = 0
    for I in 0:(aI+bI+cI+dI)
        for J in 0:(aJ+bJ+cJ+dJ)
            for K in 0:(aK+bK+cK+dK)
                s += Bx[I]*By[J]*Bz[K]*Fgamma(I+J+K,x)
            end
        end
    end
    return 2*pi^(2.5)/(g1*g2*sqrt(g1+g2))*exp(-aexpn*bexpn*rab2/g1)*exp(-cexpn*dexpn*rcd2/g2)*s
end

"""
    coulomb(a::PGBF,b::PGBF,c::PGBF,d::PGBF)

Compute the coulomb repulsion between four primitive Gaussian basis
functions a,b,c,d.
"""
function coulomb(a::PGBF,b::PGBF,c::PGBF,d::PGBF)
    return a.norm*b.norm*c.norm*d.norm*coulomb(a.expn,a.xyz...,a.I,a.J,a.K,
        b.expn,b.xyz...,b.I,b.J,b.K,
        c.expn,c.xyz...,c.I,c.J,c.K,
        d.expn,d.xyz...,d.I,d.J,d.K)
end

fB(i,l1,l2,p,a,b,r,g) = binomial_prefactor(i,l1,l2,p-a,p-b)*B0(i,r,g)
B0(i,r,g) = fact_ratio2(i,r)*(4g)^(r-i)
#fact_ratio2(a,b) = factorial(a)/factorial(b)/factorial(a-2b)
fact_ratio2(a,b) = prod((b+1):a)/factorial(a-2b) # faster

function Bterm(i1,i2,r1,r2,u,l1,l2,l3,l4,
        Px,Ax,Bx,Qx,Cx,Dx,
        gamma1,gamma2,delta)
    # THO eq. 2.22
    return fB(i1,l1,l2,Px,Ax,Bx,r1,gamma1)*
        (-1)^i2*fB(i2,l3,l4,Qx,Cx,Dx,r2,gamma2)*
        (-1)^u*fact_ratio2(i1+i2-2*(r1+r2),u)*
        (Qx-Px)^(i1+i2-2*(r1+r2)-2*u)/delta^(i1+i2-2*(r1+r2)-u)
end

function Barray(l1,l2,l3,l4,p,a,b,q,c,d,g1,g2,delta)
    Imax = l1+l2+l3+l4+1
    B = OffsetArray(zeros(Float64,Imax),0:(Imax-1))
    for i1 in 0:(l1+l2)
        for i2 in 0:(l3+l4)
            for r1 in 0:div(i1,2)
                for r2 in 0:div(i2,2)
                    for u in 0:(div(i1+i2,2)-r1-r2)
                        I = i1+i2-2*(r1+r2)-u
                        B[I] += Bterm(i1,i2,r1,r2,u,l1,l2,l3,l4,
                                      p,a,b,q,c,d,g1,g2,delta)
                    end
                end
            end
        end
    end
    return B
end

"""
    coulomb(a::CGBF,b::CGBF,c::CGBF,d::CGBF)

Compute the coulomb repulsion between four contracted Gaussian basis
functions a,b,c,d.
"""
coulomb(a::CGBF,b::CGBF,c::CGBF,d::CGBF) = contract(coulomb,a,b,c,d)

function all_twoe_ints(bfs,ERI=coulomb)
    n = length(bfs)
    totlen = div(n*(n+1)*(n*n+n+2),8)
    ints2e = zeros(Float64,totlen)
    Threads.@threads for (i,j,k,l) in collect(iiterator(n))
        ints2e[iindex(i,j,k,l)] = ERI(bfs[i],bfs[j],bfs[k],bfs[l])
    end
    return ints2e
end

function make2JmK(D::Array{Float64,2},Ints::Array{Float64,1})
    n = size(D,1)
    G = Array{Float64}(n,n)
    D1 = reshape(D,n*n)
    temp = Array{Float64}(n*n)
    for (i,j) in pairs(n)
        kl = 1
        for (k,l) in rpairs(n)
            temp[kl] = 2*Ints[iindex(i,j,k,l)]-Ints[iindex(i,k,j,l)]
            kl += 1
        end
        G[i,j] = G[j,i] = dot(D1,temp)
    end
    return G
end

dmat(U::Array{Float64,2},nocc::Int64) = U[:,1:nocc]*U[:,1:nocc]'


text:
include("polymodule.jl")
include("variable.jl")
include("constraint.jl")
include("functions.jl")

include("Mock/mock_tests.jl")


text:
if isdefined(@__MODULE__, :DEBUG)
    MQuerySymbol(x...) = Symbol(join(["MQuery", x...], "_"))
else
    MQuerySymbol(x...) = Symbol(join(["MQuery", x...], "."))
end

ARG = MQuerySymbol("ARG") # we just need limited mangled symbols here.
TYPE_ROOT = MQuerySymbol("TYPE_ROOT")

IN_TYPES = MQuerySymbol("IN", "TYPES")
IN_FIELDS = MQuerySymbol("IN", "FIELDS")
IN_SOURCE = MQuerySymbol("IN", "SOURCE")

OUT_TYPES = MQuerySymbol("OUT", "TYPES")
OUT_FIELDS = MQuerySymbol("MQuery", "OUT.FIELDS")
OUT_SOURCE = MQuerySymbol("MQuery", "OUT.SOURCE")
RECORD = MQuerySymbol("RECORD")

N = MQuerySymbol("N")
GROUPS = MQuerySymbol("GROUPS")
GROUP_KEY = MQuerySymbol("GROUP_KEY")

FN = MQuerySymbol("FN")
FN_RETURN_TYPES = MQuerySymbol("FN", "RETURN_TYPES")
FN_OUT_FIELDS = MQuerySymbol("FN", "OUT_FIELDS")

AGG = MQuerySymbol("AGG")
AGG_TYPES = MQuerySymbol("AGG", "TYPES")

_gen_sym_count = 0
function gen_sym()
    global _gen_sym_count
    let sym = MQuerySymbol("TMP", _gen_sym_count)
        _gen_sym_count = _gen_sym_count  + 1
        sym
    end
end

function gen_sym(a :: Union{Symbol, Int, String})
    global _gen_sym_count
    MQuerySymbol("Symbol", a)
end


text:
function deprecatesolver(solver)
	error("""
	NEOS$(solver)Solver(server; kwargs...) has been removed from NEOS.jl

	You should use
		NEOSSolver(server; solver=:$(solver), kwargs...)
	instead. See the README for more details.
	""")
end
NEOSCPLEXSolver(server; kwargs...)    = deprecatesolver("CPLEX")
NEOSMOSEKSolver(server; kwargs...)    = deprecatesolver("MOSEK")
NEOSXpressSolver(server; kwargs...)   = deprecatesolver("Xpress")
NEOSSYMPHONYSolver(server; kwargs...) = deprecatesolver("SYMPHONY")

# NEOSSolver{:Solver, :Format}
mutable struct NEOSSolver{S,F} <: AbstractMathProgSolver
	server::NEOSServer
	requires_email::Bool
	solves_sos::Bool
	provides_duals::Bool
	template::String
	params::Dict{String,Any}
	gzipmodel::Bool
	print_results::Bool
	print_level::Int
	result_file::String
end

function NEOSSolver(server::NEOSServer=NEOSServer();
		solver::Symbol=:CPLEX, format=:MPS,
		email::String="",  gzipmodel::Bool=true,
		print_results::Bool=false, result_file::String="",
		print_level::Int = 0,
		kwargs...
	)
	if !haskey(TEMPLATES, (solver, format))
		error("Solver $(solver) with format $(format) not supported.")
	end

	if email != ""
		addemail!(server, email)
	end
	params=Dict{String,Any}()
	for (key, value) in kwargs
		params[string(key)] = value
	end
	(template, requires_email, solves_sos, provides_duals) = TEMPLATES[(solver, format)]
	return NEOSSolver{solver, format}(server, requires_email, solves_sos,
		provides_duals, template, params, gzipmodel, print_results, print_level, result_file)
end

function addparameter!(s::NEOSSolver, param::String, value)
	s.params[param] = value
end
addparameter!(s::NEOSSolver, param::Symbol, value) = addparameter!(s, string(param), value)

addemail!(s::NEOSSolver, email::String) = addemail!(s.server, email)

# TEMPLATE[(solver, format)] = (template, requireemail, solvesos, provideduals)
const TEMPLATES = Dict{Tuple{Symbol, Symbol}, Tuple{String, Bool, Bool, Bool}}(
    (:CPLEX, :MPS) => ("""
        <document>
        <client>NEOS.jl</client>
        <category>milp</category>
        <solver>CPLEX</solver>
        <inputMethod>MPS</inputMethod>
        <email></email>
        <MPS></MPS>
        <options></options>
        <post></post>
        </document>
        """, true, true, true),
    (:Xpress, :MPS) => ("""
    	<document>
    	<client>NEOS.jl</client>
    	<category>lp</category>
    	<solver>FICO-Xpress</solver>
    	<inputMethod>MPS</inputMethod>
    	<email></email>
    	<MPS></MPS>
    	<algorithm></algorithm>
    	<par></par>
    	</document>
    	""", true, true, false),
    (:MOSEK, :MPS) => ("""
    	<document>
    	<client>NEOS.jl</client>
    	<category>milp</category>
    	<solver>MOSEK</solver>
    	<inputMethod>MPS</inputMethod>
    	<email></email>
    	<MPS></MPS>
    	<param></param>
    	<wantsol><![CDATA[yes]]></wantsol>
    	<wantint><![CDATA[yes]]></wantint>
    	</document>
    	""", false, false, true),
    (:SYMPHONY, :MPS) => ("""
    	<document>
    	<client>NEOS.jl</client>
    	<category>milp</category>
    	<solver>SYMPHONY</solver>
    	<inputMethod>MPS</inputMethod>
    	<email></email>
    	<MPS></MPS>
    	<options></options>
    	</document>
    	""", false, false, false),
    (:CPLEX, :NL) => ("""
    	<document>
    	<client>NEOS.jl</client>
    	<category>milp</category>
    	<solver>CPLEX</solver>
    	<inputMethod>NL</inputMethod>
    	<email></email>
    	<model></model>
    	<options></options>
    	</document>
    	""", true, false, false)
)


text:
#===============================================================================
# Some practice using various julia interfaces
#
# Author: Carl Cortright
# Date: 12/31/2016
#
===============================================================================#
# A simple exampe of an interable squares object with finite length
immutable squares
  count::Int
end

Base.start(::squares) = 1
Base.next(S::squares, state) = (state*state, state+1)
Base.done(S::squares, state) = state > S.count
Base.eltype(::Type{squares}) = Int
Base.length(S::squares) = S.count

# We can now iterate over a squares object in a for loop
s = squares(100)
for i in s
  println(i)
end

println(25 in squares(10))

# We can also define the getindex() function
function Base.getindex(S::squares, i::Int)
  1 <= i <= S.count || throw(BoundsError(S, i))
  return i*i
end

s = squares(1000000)
println(s[982031])
# If we want to support the end syntax
Base.endof(S::squares) = length(s)
endval = s[end]
println("s[end] = $endval")

# We can also define as an abstract array to make things faster
immutable SquaresVector <: AbstractArray{Int, 1}
  count::Int
end
# Note {Int, 1} defines eltype and ndims
Base.length(S::SquaresVector) = (S.count,)
Base.linearindexing{T <: SquaresVector}(::Type{T}) = Base.LinearFast()
Base.getindex(S::SquaresVector, i::Int) = i*i

max_issues_repo_name: carlcortright/Playground


text:
abstract type AbstractSaturated <: Function end

"""
    Saturated(f, lb, ub)

Function wrapper that saturates the output at lower bound `lb` and
upper bound `ub`.
"""
Base.@kwdef @concrete struct Saturated <: AbstractSaturated
    f
    lb
    ub
end

(f::Saturated)(args...; kwargs...) = clamp(f.f(args...; kwargs...), f.lb, f.ub)

Base.@kwdef @concrete struct SoftSaturated <: AbstractSaturated
    f
    lb
    ub
    halfwidth
end

(f::SoftSaturated)(args...; kwargs...) = soft_clamp(f.f(args...; kwargs...), f.lb, f.ub; f.halfwidth)

"""
    WithControls(sys_fun, controls::NamedTuple)

ODE function wrapper that applies control inputs through keyword arguments.

## Example
```julia
function f(x, p, t; u)
    A, B = p.A, p.B
    return A*x + B*u
end

ctrl_f = WithControls(f; u = (x,p,t) -> -p.kp*x)

p = (A=2.5, B=1.5, kp=2.5)  # Parameters
x0 = 1.0                    # Initial Condition
tspan = (0.0, 10.0)         # Time span

prob = ODEProblem(ctrl_f, x0, tspan, p)
sol = solve(prob)
```
"""
@concrete struct WithControls <: AbstractSaturated
    sys_fun
    controls <: NamedTuple
end
WithControls(sys_fun; kwargs...) = WithControls(sys_fun, NamedTuple(kwargs))

(sys::WithControls)(dx, x, p, t) = sys.sys_fun(dx, x, p, t; map(f->maybe_apply(f,x,p,t), sys.controls)...)
(sys::WithControls)(x, p, t) = sys.sys_fun(x, p, t; map(f->maybe_apply(f,x,p,t), sys.controls)...)

DifferentialEquations.isinplace(f::WithControls, args...) = isinplace(f.sys_fun, args...)

"""
    maybe_apply(f, args...)

Maybe apply a function, maybe don't. I don't know. I'm not your boss.
"""
maybe_apply(f::Function, args...) = f(args...)
maybe_apply(val, args...) = val

function saturation_callback(sat_fun::AbstractSaturated)
    @unpack f, lb, ub = sat_fun

    function condition(vars, t, integrator)

    end
end

# Needed to make the symbolic stuff work
ComponentArrays.ComponentArray(x::Num) = x

reachable_constraint(f, prob::ODEProblem) = reachable_constraint(f, prob.u0, prob.p, prob.tspan[1])
# reachable_constraint(f, prob::RA.InitialValueProblem) = reachable_constraint(f, prob.x0.center.u, prob.x0.center.p, prob.x0.center.t)
function reachable_constraint(f::Function, u0::ComponentArray, p::ComponentArray, t)
    x0 = make_symbolic_x0(u0, p, t)
    exprs = f(x0.u, x0.p, x0.t)
    return HalfSpace(exprs, x0)
end
function reachable_constraint(fs::AbstractArray, u0::ComponentArray, p::ComponentArray, t)
    x0 = make_symbolic_x0(u0, p, t)
    exprs = map(f->f(x0.u, x0.p, x0.t), fs)
    return HPolyhedron(exprs, x0)
end

function make_symbolic_x0(u0, p, t)
    @variables U[1:length(u0)], P[1:length(p)], T
    u_ca = ComponentArray(collect(U), getaxes(u0))
    p_ca = ComponentArray(collect(P), getaxes(p))
    return ComponentArray(u=u_ca, p=p_ca, t=T)
end

function make_reach_function(f, unflatten_u, unflatten_p, unflatten_t; inplace=true)
    return if inplace
        function (dx, x, _, t)
            U = unflatten_u(x)
            dU = unflatten_u(dx)
            P = unflatten_p(x)
            dP = unflatten_p(dx)
            T = unflatten_t(x)

            f(dU, U, P, T) # Update state
            dP .= zero(x[1])    # Params have no dynamics
            dx[end] = one(x[1]) # Update time variable
            return nothing
        end
    else
        function (x, _, t)
            U = unflatten_u(x)
            P = unflatten_p(x)
            T = unflatten_t(x)

            dU = f(U, P, T)
            return [dU; zeros(typeof(x[1]), length(P)); one(x[1])]
        end
    end
end

function reach_prepare(u0, p, t)
    nu = length(u0)
    np = length(p)
    u0, unflatten_u = _reach_prepare(u0)
    p, unflatten_p = _reach_prepare(p, nu)
    t, unflatten_t = _reach_prepare(t, nu+np)
    x0 = prod(u0) × prod(p) × t |> concretize
    # center = ComponentArray(u=unflatten_u(x0.center), p=unflatten_p(x0.center), t=unflatten_t(x0.center))
    # radius = ComponentArray(u=unflatten_u(x0.radius), p=unflatten_p(x0.radius), t=unflatten_t(x0.radius))
    # x0 = Hyperrectangle(center, radius)
    return x0, unflatten_u, unflatten_p, unflatten_t
end

function reach_prepare(u0, p)
    nu = length(u0)
    u0, unflatten_u = _reach_prepare(u0)
    p, unflatten_p = _reach_prepare(p, nu)
    x0 = prod(u0) × prod(p) |> concretize
    # center = ComponentArray(u=unflatten_u(x0.center), p=unflatten_p(x0.center))
    # radius = ComponentArray(u=unflatten_u(x0.radius), p=unflatten_p(x0.radius))
    # x0 = Hyperrectangle(center, radius)
    return x0, unflatten_u, unflatten_p
end

function _reach_prepare(p, offset=0)
    p = to_reach_interval(p)
    i = offset + 1
    unflatten = x -> x[i]
    return p, unflatten
end
function _reach_prepare(p::AbstractArray, offset=0)
    np = length(p)
    p = recursive_convert(to_reach_interval, p)
    i = (1:np) .+ offset
    unflatten = x -> view(x, i)
    return p, unflatten
end
function _reach_prepare(p::ComponentArray, offset=0)
    ax_p = getaxes(p)
    np = length(p)
    p = recursive_convert(to_reach_interval, p)
    i = (1:np) .+ offset
    unflatten = x -> ComponentArray(view(x, i), ax_p)
    return p, unflatten
end


text:
using Parameters

mutable struct Root{T}
    name::String
    self_dim::Int64
    children::Vector{String}
    q::Vector{T}
    q̇::Vector{T}
    q̈::Vector{T}
    f::Vector{T}
    M::Matrix{T}
end

function Root(self_dim::Int64)
    name = "root"
    children = String[]
    q = zeros(Float64, self_dim)
    q̇ = zero(q)
    q̈ = zero(q)
    f = zero(q)
    M = zeros(Float64, self_dim, self_dim)
    Root(name, children, self_dim, q, q̇, q̈, f, M)
end

"""ノード"""
mutable struct Node{T}
    name::String
    parent::String
    self_dim::Int64
    parent_dim::Int64
    children::Vector{String}
    x::Vector{T}
    ẋ::Vector{T}
    J::Matrix{T}
    J̇::Matrix{T}
    f::Vector{T}
    M::Matrix{T}
end

function Node(name::String, parent::String, self_dim::Int64, parent_dim::Int64)
    children = String[]
    x = zeros(Float64, self_dim)
    ẋ = zero(x)
    J = zeros(Float64, self_dim, parent_dim)
    J̇ = zero(J)
    f = zero(x)
    M = zeros(Float64, self_dim, parent_dim)
    Node(name, parent, self_dim, parent_dim, children, x, ẋ, J, J̇, f, M)
end

mutable struct Leaf{T}
    name::String
    self_dim::Int64
    parent_dim::Int64
    parent::String
    x::Vector{T}
    ẋ::Vector{T}
    J::Matrix{T}
    J̇::Matrix{T}
    f::Vector{T}
    M::Matrix{T}
end

function Leaf(name::String, parent::String, self_dim::Int64, parent_dim::Int64)
    x = zeros(Float64, self_dim)
    ẋ = zero(x)
    J = zeros(Float64, self_dim, parent_dim)
    J̇ = zero(J)
    f = zero(x)
    M = zeros(Float64, self_dim, parent_dim)
    Leaf(name, parent, self_dim, parent_dim, x, ẋ, J, J̇, f, M)
end

function pushforward!(root::Root{T}, rmp_tree::Dict{String, Any}) where {T}
    for child in root.children
        update!(root.x, root.ẋ, rmp_tree[child].x, rmp_tree[child].ẋ, rmp_tree[child].J, rmp_tree[child].J̇)
        pushforward!(rmp_tree[child], rmp_tree)
    end
end

function pushforward!(node::Node{T}, rmp_tree::Dict{String, Any}) where {T}
    for child in node.children
        pushforward!(rmp_tree[child], rmp_tree)
    end
end

function pushforward!(leaf::Leaf{T}, rmp_tree::Dict{String, Any}) where {T}
    println("this is leaf")
end


text:
using Weave
using Test

# Make a dummy codehunk with figure
chunk = Weave.CodeChunk("plot(x)", 1, 1, "", Dict())
options = merge(Weave.rcParams[:chunk_defaults], chunk.options)
merge!(chunk.options, options)
chunk.figures = ["figs/figures_plot1.png"]

@test Weave.formatfigures(chunk, Weave.md2tex) == "\\includegraphics{figs/figures_plot1.png}\n"
@test Weave.formatfigures(chunk, Weave.tex) == "\\includegraphics{figs/figures_plot1.png}\n"
@test Weave.formatfigures(chunk, Weave.texminted) == "\\includegraphics{figs/figures_plot1.png}\n"
@test Weave.formatfigures(chunk, Weave.pandoc) == "![](figs/figures_plot1.png)\\ \n\n"
@test Weave.formatfigures(chunk, Weave.github) == "![](figs/figures_plot1.png)\n"
@test Weave.formatfigures(chunk, Weave.hugo) == "{{
}}"
@test Weave.formatfigures(chunk, Weave.multimarkdown) == "![][figs/figures_plot1.png]\n\n[figs/figures_plot1.png]: figs/figures_plot1.png \n"
@test Weave.formatfigures(chunk, Weave.md2html) == "\n"

chunk.options[:out_width] = "100%"
@test Weave.formatfigures(chunk, Weave.adoc) == "image::figs/figures_plot1.png[width=100%]\n"
@test Weave.formatfigures(chunk, Weave.rst) == ".. image:: figs/figures_plot1.png\n   :width: 100%\n\n"

chunk.options[:fig_cap] = "Nice plot"
@test Weave.formatfigures(chunk, Weave.tex) == "\\begin{figure}[!h]\n\\center\n\\includegraphics[width=1.0\\linewidth]{figs/figures_plot1.png}\n\\caption{Nice plot}\n\\end{figure}\n"
@test Weave.formatfigures(chunk, Weave.pandoc) == "![Nice plot](figs/figures_plot1.png){width=100%}\n"
@test Weave.formatfigures(chunk, Weave.md2tex) == "\\begin{figure}[!h]\n\\center\n\\includegraphics[width=1.0\\linewidth]{figs/figures_plot1.png}\n\\caption{Nice plot}\n\\end{figure}\n"
@test Weave.formatfigures(chunk, Weave.md2html) == "
\n\n
Nice plot
\n
\n"
@test Weave.formatfigures(chunk, Weave.rst) == ".. figure:: figs/figures_plot1.png\n   :width: 100%\n\n   Nice plot\n\n"
@test Weave.formatfigures(chunk, Weave.multimarkdown) == "![Nice plot][figs/figures_plot1.png]\n\n[figs/figures_plot1.png]: figs/figures_plot1.png width=100%\n"
@test Weave.formatfigures(chunk, Weave.adoc) == "image::figs/figures_plot1.png[width=100%,title=\"Nice plot\"]"

chunk.options[:label] = "somefig"
@test Weave.formatfigures(chunk, Weave.pandoc) == "![Nice plot](figs/figures_plot1.png){width=100% #fig:somefig}\n"
@test Weave.formatfigures(chunk, Weave.tex) == "\\begin{figure}[!h]\n\\center\n\\includegraphics[width=1.0\\linewidth]{figs/figures_plot1.png}\n\\caption{Nice plot}\n\\label{fig:somefig}\n\\end{figure}\n"
@test Weave.formatfigures(chunk, Weave.tex) == Weave.formatfigures(chunk, Weave.md2tex)

chunk.options[:label] = nothing
chunk.options[:fig_cap] = nothing
chunk.options[:fig_env] = "center"
chunk.options[:fig_pos] = ""
@test Weave.formatfigures(chunk, Weave.tex) == "\\begin{center}\n\\includegraphics[width=1.0\\linewidth]{figs/figures_plot1.png}\n\\end{center}\n"

chunk.options[:out_width] = "50%"
chunk.options[:out_height] = "75 %"
@test Weave.formatfigures(chunk, Weave.tex) == "\\begin{center}\n\\includegraphics[width=0.5\\linewidth,height=0.75\\paperheight]{figs/figures_plot1.png}\n\\end{center}\n"

chunk.options[:out_width] = "A%"
chunk.options[:out_height] = "0.5\\textwidth"
@test Weave.formatfigures(chunk, Weave.tex) == "\\begin{center}\n\\includegraphics[width=A%,height=0.5\\textwidth]{figs/figures_plot1.png}\n\\end{center}\n"


text:
module Management

using CoreExtended

abstract type IManagerReference end

list = SortedDict{Symbol, IManagerReference}(Forward)

""" TODO """
type Manager{T} <: IManagerReference
	list::SortedDict{Symbol, T}
	selected::Union{Void,T}

	create::Function
	remove::Function
	init::Function
	reset::Function
	link::Function
	unlink::Function

	Manager{T}() where T = new(
		SortedDict{Symbol, T}(Forward),
		nothing,
		(id::Symbol) -> T(id),
		(obj::T) -> nothing,
		(obj::T) -> nothing,
		() -> nothing,
		(obj::T) -> nothing,
		() -> nothing
	)
end

""" TODO """
function presetManager(T::DataType)
	m=createManager(T)

	eval(T.name.module,Expr(:toplevel,:(
		m=Management.list[Symbol($(T.name))];
		const Typ = Management.getType(m);

		create(k) = Management.create(m,k);
		getSelected() = Management.getSelected(m);
		setSelected(o) = Management.setSelected(m,o);

		get(k) =  Management.get(m,k);
		set(k,o) = Management.set(m,k,o);
		del(k) = Management.del(m,k);

		reset() = Management.reset(m);

		isInvalid(o) = Management.isInvalid(m,o);
		isLinked(o) = Management.isLinked(m,o);

		link() = Management.link(m);
		unlink() = Management.unlink(m);

		select(o) = Management.select(m,o);
	)))

	m
end

""" TODO """
function createManager(T::DataType)
	k = Symbol(T.name)
	if !haskey(list,k)
		e=Manager{T}()
		list[k]=e
	else
		e=list[k]
	end
	e
end

""" TODO """
function create{T}(m::Manager{T}, k::Symbol)
	if !haskey(m.list,k)
		e=m.create(k)
		m.list[k]=e
		setSelected(m,e)
		m.init(e)
	else
		e=m.list[k]
		setSelected(m,e)
	end
	e
end

""" TODO """
getSelected{T}(m::Manager{T}) = m.selected

""" TODO """
setSelected{T}(m::Manager{T}, obj::Union{Void,T}) = (m.selected = obj)

""" TODO """
getType{T}(m::Manager{T}) = Union{Void,T}

""" TODO """
get{T}(m::Manager{T}, k::Symbol) = m.list[k]

""" TODO """
set{T}(m::Manager{T}, k::Symbol, obj::Union{Void,T}) = (m.list[k] = obj)

""" TODO """
function del{T}(m::Manager{T}, k::Symbol)
	obj=get(m,k)
	selected=getSelected(m)
	if selected == obj unlink() end
	m.remove(obj)
	set(m,k,nothing)
end

""" TODO """
function reset{T}(m::Manager{T})
	unlink(m)
	m.reset()
end;

""" TODO """
isInvalid{T}(m::Manager{T}, obj::Union{Void,T}) = obj == nothing

""" TODO """
isLinked{T}(m::Manager{T}, obj::Union{Void,T}) =	!isInvalid(m,obj) && obj == getSelected(m)

""" TODO """
unlink{T}(m::Manager{T}) = (setSelected(m,nothing); m.unlink())

""" TODO """
function link{T}(m::Manager{T})
	obj=getSelected(m)
	if !isInvalid(m,obj) m.link(obj)
	else unlink()
	end
end

""" TODO """
function select{T}(m::Manager{T}, obj::Union{Void,T})
	other=getSelected(m)
	if other == obj return end
	setSelected(m,obj)
	link(m)
end

end #Management

max_stars_repo_name: Gilga/JGE

max_issues_repo_name: Gilga/JGE

max_forks_repo_name: Gilga/JGE


text:
#csrlsvlu
for (fname, elty, relty) in ((:cusolverSpScsrlsvluHost, :Float32, :Float32),
                             (:cusolverSpDcsrlsvluHost, :Float64, :Float64),
                             (:cusolverSpCcsrlsvluHost, :Complex64, :Float32),
                             (:cusolverSpZcsrlsvluHost, :Complex128, Float64))
    @eval begin
        function csrlsvlu!(A::SparseMatrixCSC{$elty},
                           b::Vector{$elty},
                           x::Vector{$elty},
                           tol::$relty,
                           reorder::Cint,
                           inda::SparseChar)
            cuinda = cusparseindex(inda)
            n = size(A,1)
            if size(A,2) != n
                throw(DimensionMismatch("LU factorization is only possible for square matrices!"))
            end
            if size(A,2) != length(b)
                throw(DimensionMismatch("second dimension of A, $(size(A,2)), must match the length of b, $(length(b))"))
            end
            if length(x) != length(b)
                throw(DimensionMismatch("length of x, $(length(x)), must match the length of b, $(length(b))"))
            end
            Mat = transpose(A)
            cudesca = cusparseMatDescr_t(CUSPARSE_MATRIX_TYPE_GENERAL, CUSPARSE_FILL_MODE_LOWER, CUSPARSE_DIAG_TYPE_NON_UNIT, cuinda)
            singularity = zeros(Cint,1)
            statuscheck(ccall(($(string(fname)),libcusolver), cusolverStatus_t,
                              (cusolverSpHandle_t, Cint, Cint,
                               Ptr{cusparseMatDescr_t}, Ptr{$elty}, Ptr{Cint},
                               Ptr{Cint}, Ptr{$elty}, $relty, Cint, Ptr{$elty},
                               Ptr{Cint}),
                              cusolverSphandle[1], n, length(A.nzval), &cudesca,
                              Mat.nzval, convert(Vector{Cint},Mat.colptr),
                              convert(Vector{Cint},Mat.rowval), b, tol, reorder,
                              x, singularity))
            if singularity[1] != -1
                throw(Base.LinAlg.SingularException(singularity[1]))
            end
            x
        end
    end
end

#csrlsvqr
for (fname, elty, relty) in ((:cusolverSpScsrlsvqr, :Float32, :Float32),
                             (:cusolverSpDcsrlsvqr, :Float64, :Float64),
                             (:cusolverSpCcsrlsvqr, :Complex64, :Float32),
                             (:cusolverSpZcsrlsvqr, :Complex128, Float64))
    @eval begin
        function csrlsvqr!(A::CudaSparseMatrixCSR{$elty},
                           b::CudaVector{$elty},
                           x::CudaVector{$elty},
                           tol::$relty,
                           reorder::Cint,
                           inda::SparseChar)
            cuinda = cusparseindex(inda)
            n = size(A,1)
            if size(A,2) != n
                throw(DimensionMismatch("QR factorization is only possible for square matrices!"))
            end
            if size(A,2) != length(b)
                throw(DimensionMismatch("second dimension of A, $(size(A,2)), must match the length of b, $(length(b))"))
            end
            if length(x) != length(b)
                throw(DimensionMismatch("length of x, $(length(x)), must match the length of b, $(length(b))"))
            end
            cudesca = cusparseMatDescr_t(CUSPARSE_MATRIX_TYPE_GENERAL, CUSPARSE_FILL_MODE_LOWER, CUSPARSE_DIAG_TYPE_NON_UNIT, cuinda)
            singularity = Array{Cint}(1)
            statuscheck(ccall(($(string(fname)),libcusolver), cusolverStatus_t,
                              (cusolverSpHandle_t, Cint, Cint,
                               Ptr{cusparseMatDescr_t}, Ptr{$elty}, Ptr{Cint},
                               Ptr{Cint}, Ptr{$elty}, $relty, Cint, Ptr{$elty},
                               Ptr{Cint}),
                              cusolverSphandle[1], n, A.nnz, &cudesca, A.nzVal,
                              A.rowPtr, A.colVal, b, tol, reorder, x, singularity))
            if singularity[1] != -1
                throw(Base.LinAlg.SingularException(singularity[1]))
            end
            x
        end
    end
end

#csrlsvchol
for (fname, elty, relty) in ((:cusolverSpScsrlsvchol, :Float32, :Float32),
                             (:cusolverSpDcsrlsvchol, :Float64, :Float64),
                             (:cusolverSpCcsrlsvchol, :Complex64, :Float32),
                             (:cusolverSpZcsrlsvchol, :Complex128, Float64))
    @eval begin
        function csrlsvchol!(A::CudaSparseMatrixCSR{$elty},
                             b::CudaVector{$elty},
                             x::CudaVector{$elty},
                             tol::$relty,
                             reorder::Cint,
                             inda::SparseChar)
            cuinda = cusparseindex(inda)
            n      = size(A,1)
            if size(A,2) != n
                throw(DimensionMismatch("Cholesky factorization is only possible for square matrices!"))
            end
            if size(A,2) != length(b)
                throw(DimensionMismatch("second dimension of A, $(size(A,2)), must match the length of b, $(length(b))"))
            end
            if length(x) != length(b)
                throw(DimensionMismatch("length of x, $(length(x)), must match the length of b, $(length(b))"))
            end

            cudesca     = cusparseMatDescr_t(CUSPARSE_MATRIX_TYPE_GENERAL, CUSPARSE_FILL_MODE_LOWER, CUSPARSE_DIAG_TYPE_NON_UNIT, cuinda)
            singularity = zeros(Cint,1)
            statuscheck(ccall(($(string(fname)),libcusolver), cusolverStatus_t,
                              (cusolverSpHandle_t, Cint, Cint,
                               Ptr{cusparseMatDescr_t}, Ptr{$elty}, Ptr{Cint},
                               Ptr{Cint}, Ptr{$elty}, $relty, Cint, Ptr{$elty},
                               Ptr{Cint}),
                              cusolverSphandle[1], n, A.nnz, &cudesca, A.nzVal,
                              A.rowPtr, A.colVal, b, tol, reorder, x, singularity))
            if singularity[1] != -1
                throw(Base.LinAlg.SingularException(singularity[1]))
            end
            x
        end
    end
end

#csrlsqvqr
for (fname, elty, relty) in ((:cusolverSpScsrlsqvqrHost, :Float32, :Float32),
                             (:cusolverSpDcsrlsqvqrHost, :Float64, :Float64),
                             (:cusolverSpCcsrlsqvqrHost, :Complex64, :Float32),
                             (:cusolverSpZcsrlsqvqrHost, :Complex128, Float64))
    @eval begin
        function csrlsqvqr!(A::SparseMatrixCSC{$elty},
                            b::Vector{$elty},
                            x::Vector{$elty},
                            tol::$relty,
                            inda::SparseChar)
            cuinda = cusparseindex(inda)
            m,n    = size(A)
            if m < n
                throw(DimensionMismatch("csrlsqvqr only works when the first dimension of A, $m, is greater than or equal to the second dimension of A, $n"))
            end
            if size(A,2) != length(b)
                throw(DimensionMismatch("second dimension of A, $(size(A,2)), must match the length of b, $(length(b))"))
            end
            if length(x) != length(b)
                throw(DimensionMismatch("length of x, $(length(x)), must match the length of b, $(length(b))"))
            end
            cudesca  = cusparseMatDescr_t(CUSPARSE_MATRIX_TYPE_GENERAL, CUSPARSE_FILL_MODE_LOWER, CUSPARSE_DIAG_TYPE_NON_UNIT, cuinda)
            p        = zeros(Cint,n)
            min_norm = zeros($relty,1)
            rankA    = zeros(Cint,1)
            Mat      = transpose(A)
            statuscheck(ccall(($(string(fname)),libcusolver), cusolverStatus_t,
                              (cusolverSpHandle_t, Cint, Cint, Cint,
                               Ptr{cusparseMatDescr_t}, Ptr{$elty}, Ptr{Cint},
                               Ptr{Cint}, Ptr{$elty}, $relty, Ptr{Cint},
                               Ptr{$elty}, Ptr{Cint}, Ptr{$relty}),
                              cusolverSphandle[1], m, n, length(A.nzval),
                              &cudesca, Mat.nzval, convert(Vector{Cint},Mat.colptr),
                              convert(Vector{Cint},Mat.rowval), b,
                              tol, rankA, x, p, min_norm))
            x, rankA[1], p, min_norm[1]
        end
    end
end

#csreigvsi
for (fname, elty, relty) in ((:cusolverSpScsreigvsi, :Float32, :Float32),
                             (:cusolverSpDcsreigvsi, :Float64, :Float64),
                             (:cusolverSpCcsreigvsi, :Complex64, :Float32),
                             (:cusolverSpZcsreigvsi, :Complex128, Float64))
    @eval begin
        function csreigvsi(A::CudaSparseMatrixCSR{$elty},
                           μ_0::$elty,
                           x_0::CudaVector{$elty},
                           tol::$relty,
                           maxite::Cint,
                           inda::SparseChar)
            cuinda = cusparseindex(inda)
            m,n    = size(A)
            if m != n
                throw(DimensionMismatch("A must be square!"))
            end
            if n != length(x_0)
                throw(DimensionMismatch("second dimension of A, $(size(A,2)), must match the length of x_0, $(length(x_0))"))
            end
            cudesca = cusparseMatDescr_t(CUSPARSE_MATRIX_TYPE_GENERAL, CUSPARSE_FILL_MODE_LOWER, CUSPARSE_DIAG_TYPE_NON_UNIT, cuinda)
            x       = copy(x_0)
            μ       = CudaArray(zeros($elty,1))
            statuscheck(ccall(($(string(fname)),libcusolver), cusolverStatus_t,
                              (cusolverSpHandle_t, Cint, Cint,
                               Ptr{cusparseMatDescr_t}, Ptr{$elty}, Ptr{Cint},
                               Ptr{Cint}, $elty, Ptr{$elty}, Cint,
                               $relty, Ptr{$elty}, Ptr{$elty}),
                              cusolverSphandle[1], n, A.nnz, &cudesca, A.nzVal,
                              A.rowPtr, A.colVal, μ_0, x_0, maxite, tol, μ, x))
            to_host(μ)[1], x
        end
    end
end

#csreigs
for (fname, elty, relty) in ((:cusolverSpScsreigsHost, :Complex64, :Float32),
                             (:cusolverSpDcsreigsHost, :Complex128, :Float64),
                             (:cusolverSpCcsreigsHost, :Complex64, :Complex64),
                             (:cusolverSpZcsreigsHost, :Complex128, :Complex128))
    @eval begin
        function csreigs(A::SparseMatrixCSC{$relty},
                         lbc::$elty,
                         ruc::$elty,
                         inda::SparseChar)
            cuinda = cusparseindex(inda)
            m,n    = size(A)
            if m != n
                throw(DimensionMismatch("A must be square!"))
            end
            cudesca = cusparseMatDescr_t(CUSPARSE_MATRIX_TYPE_GENERAL, CUSPARSE_FILL_MODE_LOWER, CUSPARSE_DIAG_TYPE_NON_UNIT, cuinda)
            numeigs = zeros(Cint,1)
            Mat     = A.'
            statuscheck(ccall(($(string(fname)),libcusolver), cusolverStatus_t,
                              (cusolverSpHandle_t, Cint, Cint,
                               Ptr{cusparseMatDescr_t}, Ptr{$relty}, Ptr{Cint},
                               Ptr{Cint}, $elty, $elty, Ptr{Cint}),
                              cusolverSphandle[1], n, length(A.nzval), &cudesca,
                              Mat.nzval, convert(Vector{Cint},Mat.colptr),
                              convert(Vector{Cint},Mat.rowval), lbc, ruc, numeigs))
            numeigs[1]
        end
    end
end


text:
@testset "Matrix multiplication" begin

    # test default settings
    @test get_multiplication_mode() == Dict(:multiplication => :fast)

    A = [2..4 -2..1; -1..2 2..4]
    set_multiplication_mode(:slow)
    @test A * A == [0..18 -16..8; -8..16 0..18]

    set_multiplication_mode(:rank1)
    @test A * A == [0..18 -16..8; -8..16 0..18]

    set_multiplication_mode(:fast)
    @test A * A == [-2..19.5 -16..10; -10..16 -2..19.5]
    @test A * mid.(A) == [5..12.5 -8..2; -2..8 5..12.5]
    @test mid.(A) * A == [5..12.5 -8..2; -2..8 5..12.5]
end


text:
"""
    datasetpath(name)

Return the folder that dataset `name` is stored.

If it hasn't been downloaded yet, you will be asked if you want to
download it. See [`Datasets.DATASETS`](#) for a list of available datasets.
"""
function datasetpath(name)
    datadeppath = @datadep_str "fastai-$name"
    return Path(joinpath(datadeppath, name))
end

function loadfolderdata(
        dir;
        pattern="**",
        splitfn = nothing,
        filterfn = nothing,
        loadfn = nothing)
    data = FileDataset(dir, pattern)
    if filterfn !== nothing
        data = filterobs(filterfn, data)
    end
    if splitfn !== nothing
        data = groupobs(splitfn, data)
    end
    if loadfn !== nothing
        if splitfn === nothing
            data = mapobs(loadfn, data)
        else
            data = Dict(zip(keys(data), map(d -> mapobs(loadfn, d), values(data))))
        end
    end
    return data
end

parentname(f) = f |> pathparent |> pathname
grandparentname(f) = f |> pathparent |> pathparent |> pathname
matches(re::Regex) = f -> matches(re, f)
matches(re::Regex, f) = !isnothing(match(re, f))
const RE_IMAGEFILE = r".*\.(gif|jpe?g|tiff?|png|webp|bmp)$"i
isimagefile(f) = matches(RE_IMAGEFILE, f)

function getclassessegmentation(dir::AbstractPath)
    classes = readlines(open(joinpath(dir, "codes.txt")))
    return classes
end
getclassessegmentation(name::String) = getclassessegmentation(datasetpath(name))

#=
"""
    loadtaskdata(dir, ImageSegmentation; [split = false])

Load a data container for `ImageSegmentation` with observations
`(input = image, target = mask)`.

If `split` is `true`, returns a tuple of the data containers split by
the name of the grandparent folder.

"""
function loadtaskdata(
        dir,
        ::Type{FastAI.ImageSegmentation};
        split=false,
        kwargs...)
    imagedata = mapobs(loadfile, filterobs(isimagefile, FileDataset(joinpath(dir, "images"))))
    maskdata = mapobs(maskfromimage ∘ loadfile, filterobs(isimagefile, FileDataset(joinpath(dir, "labels"))))
    return mapobs((input = obs -> obs[1], target = obs -> obs[2]), (imagedata, maskdata))
end
=#

maskfromimage(a::AbstractArray{<:Gray{T}}, classes) where T = maskfromimage(reinterpret(T, a), classes)
maskfromimage(a::AbstractArray{<:Normed{T}}, classes) where T = maskfromimage(reinterpret(T, a), classes)
function maskfromimage(a::AbstractArray{I}, classes) where {I<:Integer}
    a .+= one(I)
    return IndirectArray(a, classes)
end

"""
    loadmask(file, classes)

Load a segmentation mask from an image file. Returns an efficiently stored
array of type `eltype(classes)`.

"""
loadmask(file, classes) = maskfromimage(loadfile(file), classes)


text:
function [Z,dZ]=Topo(x,zeta,Param)
switch Param.TopoS
  case 'AgnesiCart'
    h=Param.hC/(((x-Param.x0C)/Param.aC)^2+1);
  otherwise
    h=0;
end
Z=zeta+(Param.H-zeta)*h/Param.H;
dZ=1-h/Param.H;
end


text:
# Parameter transformations

## TransformVariables transformations for θ <--> η
function transformation(problem::ProblemLike,eq::SdeRoots)
  rr,ri = splitroots(eq.r)
  ni = length(ri)
  nr = length(rr)+ni
  ir = problem.parint
  tr = as(Vector,as(Real,-ir.last,-ir.first),nr)
  return if ni==0
    as((r_re=tr,))
  else
    ii = problem.freqint
    as((r_re=tr,r_im=as(Vector,as(Real,-ii.last,-ii.first),ni)))
  end
end

function transformation(problem::ProblemLike,::Sde)
  i = problem.parint
  return as((a=as(Vector,as(Real,i.first,i.last),problem.p),))
end

transformation(problem::ProblemLike,proc::SdeProcess) = transformation(problem,proc.eq)

transformation(problem::ProblemProb,args...) = transformation(problem.like,args...)

## Equation to parameter tuples: proc <--> θ

"Estimation parameters `θ` for `TransformVariables`."
paramstuple(eq::Sde) = (a=eq.a,)

function paramstuple(eq::SdeRoots)
  rr,ri = splitroots(eq.r)
  r_re = [real(ri);real(rr)]
  θ = isempty(ri) ? (r_re=r_re,) : (r_re=r_re,r_im=imag(ri))
  return θ
end

paramstuple(::Problem,proc::SdeProcess) = paramstuple(proc.eq)

sde(::Problem{E},args...;kwargs...) where {E} = sdetype(E)(args...;kwargs...)

function SdeProcess(problem::ProblemLike,θ)
  _,σϵ = logpdf_ml_σϵ(problem,θ)
  proc = SdeProcess(sde(problem,θ,σϵ))
  return proc
end

SdeProcess(problem::ProblemProb,θ) = SdeProcess(problem.like,θ)

## Transform processes

# Some of the functionality `transform_ss` can probably more elegantly be handled by
# implementing `convert` and/or constructors combined with a function that transforms the state.

function transform_ss(D::Type,proc::SdeProcess{T,U}) where {T,U<:AbstractSde{T}}
  if !(U <: D)  error("transform_ss from $U to $D not defined")  end
  return proc
end

"
$SIGNATURES

Transform process `proc` to the state space of `D`.
"
function transform_ss(::Type{SdeRoots},proc::SdeProcess{<:Any,<:Sde})
  eq′ = SdeRoots(proc.eq)
  e = SdeRootsEngine(eq′)
  @unpack μ,Σ = proc.state
  μ′ = e.U\μ
  Σ′ = e.U\Σ/e.U'
  state′ = State(μ′,Σ′)
  proc′ = SdeProcess(eq′,state′)
  return proc′
end

function transform_ss(::Type{Sde},proc::SdeProcess{<:Any,<:SdeRoots},tol::Float64=def_tol)
  eq′ = Sde(proc.eq,tol)
  e = SdeRootsEngine(proc.eq)
  @unpack μ,Σ = proc.state
  μ′ = real(e.U*μ)
  Σ′ = real(e.U*Σ*e.U')
  state′ = State(μ′,Σ′)
  return SdeProcess(eq′,state′)
end

function transform_ss(::Type{VectorSde},proc::SdeProcess{<:Any,<:Sde})
  eq = convert(VectorSde,proc.eq)
  return SdeProcess(eq,proc.state)
end

function transform_ss(::Type{Sde},proc::SdeProcess{<:Any,<:VectorSde})
  a = -proc.eq.A[end,:]
  σϵ = √proc.eq.Q[end,end]
  eq′ = Sde(a,σϵ)
  return SdeProcess(eq′,proc.state)
end

# Custom adjoint for transform to guarantee real-valued adjoint
using TransformVariables:VectorTransform,TransformTuple,@argcheck,transform_with,NOLOGJAC

Zygote.@adjoint TransformVariables.transform(t::TransformTuple{<:Union{NamedTuple{(:r_re, :r_im)},NamedTuple{(:r_re,)}}},η) = begin
  θ,Bc = Zygote.pullback(t,η) do t,η
    # Code duplicate for `transform(t,η)`. We cannot use transform directly because the adjoint would call itself.
    # Is there a better way?
    @argcheck dimension(t) == length(η)
    first(transform_with(NOLOGJAC, t, η, firstindex(η)))
  end
  B = function(θ̄)
    t̄,η̄c = Bc(θ̄)
    η̄ = real.(η̄c)
    return t̄,η̄
  end
  return θ,B
end


text:
abstract type Output end

# mergeoutput{O<:Output}(::Vector{O}; kwargs...) = nothing
Base.close(::Output) = nothing

include("filtered.jl")
include("split.jl")

max_stars_repo_name: pwl/FASTQDemultiplexer

max_issues_repo_name: pwl/FASTQDemultiplexer

max_forks_repo_name: pwl/FASTQDemultiplexer


text:
using Flux
using Flux:LSTMCell, RNNCell
using Distributions

using Convex

# lstm
mutable struct lstm
    nu
    nv
    ny
    cell
    output
end

function lstm(nu, nv, ny)
    cell = LSTMCell(nu, nv)
    output = Dense(nv, ny)
    return lstm(nu, nv, ny, cell, output)
end

function (model::lstm)(x0, u)
    h0 = x0[1:model.nv, :]
    c0 = x0[model.nv + 1:2 * model.nv, :]

    recurrent = Flux.Recur(model.cell, (h0, c0))
    states = recurrent.(u)

    return model.output.(states), vcat(recurrent.state...)
end

function init_state(model::lstm, batches)
    return zeros(2 * model.nv, batches)
end
Flux.trainable(model::lstm) = (Flux.params(model.cell)..., Flux.params(model.output)...)

# rnn
# lstm
mutable struct rnn
    nu
    nv
    ny
    cell
    output
end

function rnn(nu, nv, ny)
    cell = RNNCell(nu, nv)
    output = Dense(nv, ny)
    return rnn(nu, nv, ny, cell, output)
end

function (model::rnn)(x0, u)
    recurrent = Flux.Recur(model.cell, x0)
    states = recurrent.(u)
    return model.output.(states), recurrent.state
end
function init_state(model::rnn, batches)
    return zeros(model.nv, batches)
end
Flux.trainable(model::rnn) = (Flux.params(model.cell)..., Flux.params(model.output)...)

# Constant memory RNN

mutable struct mRNNCell
    W
    B
    b
end
function mRNNCell(nu, nv)
    W = randn(nv, nv) / sqrt(nv + nv)
    B = randn(nv, nu) / sqrt(nv + nu)
    b =  randn(nv) / sqrt(nv)
    return mRNNCell(W, B, b)
end

function (M::mRNNCell)(z, u)
    xn = M.W * z + M.B * u .+ M.b
    J = [Diagonal(xn[:, i] .>= 0.0) for i in 1:batches]
    α = norm.(J[i] * M.W for i in 1:batches)[:, :]'
    zn = Flux.relu.(xn) ./ α
    return res, res
end
Flux.trainable(M::mRNNCell) = (M.W, M.B, M.b)

mutable struct mrnn
    nu
    nv
    ny
    cell
    output
end

function mrnn(nu, nv, ny)
    cell = mRNNCell(nu, nv)
    output = Dense(nv, ny)
    return mrnn(nu, nv, ny, cell, output)
end

function (model::mrnn)(x0, u)
    recurrent = Flux.Recur(model.cell, x0)
    states = recurrent.(u)
    return model.output.(states), recurrent.state
end

function init_state(model::mrnn, batches)
    return zeros(model.nv, batches)
end
Flux.trainable(model::mrnn) = (Flux.params(model.cell)..., Flux.params(model.output)...)

# LTI Models
mutable struct LTI
    A
    B
    C
    D
end

function simulate_states(model::LTI, x0, u)
    x = zeros(length(x0), length(u))
    x[:, 1] = x0
    for tt in 1:length(u)-1
        x[:, tt+1] = model.A*x[:, tt] + model.B*u[tt]
    end
    return x
end

function (model::LTI)(x0, u)
    states = simulate_states(model, x0, u)
    output(xt, ut) = model.C*xt + model.D*ut
    return output.(states, u)
end

mutable struct explicit_robust_rnn
    A
    B1
    B2
    C1
    D12
    bv
    C2
    D21
    D22
    by
end

mutable struct implicit_robust_rnn
    E
    B1
    B2
    C1
    D12
    bv
    C2
    D21
    D22
    by
    Λ
    P
end

function project_stable(exp_model::explicit_robust_rnn)
    𝔸 = exp_model.A
    𝔹1 = exp_model.B1
    𝔹2 = exp_model.B2
    ℂ1 = exp_model.C1
    𝔻12 = exp_model.D12
    ℂ2 = exp_model.C2
    𝔻21 = exp_model.D21
    𝔻22 = exp_model.D22

    Λ = Diagonal(Variable(nv))
    P = Variable(nx, nx)
    E = Variable(nx, nx)

    F = Variable(nx, nx)
    B1 = Variable(nx, nv)
    B2 = Variable(nx, nu)

    C1 = Variable(nv, nx)
    D12 = Variable(nv, nu)

    P = Variable(nx, nx)
    Λ = Diagonal(Variable(nv))

    # Correctness
    error1 = E * [𝔸 𝔹1 𝔹2] - [F B1 B2]
    error2 = Λ * [ℂ1 𝔻12] - [C1 D12]
    J = norm(error1) + norm(error2)

    # Stability constraint
    M = [E+E'-P -C1' F'
         -C1 2Λ B1'
         F B1 P]

    constraints = [M⪰1E-4]
    solver = () -> SCS.Optimizer(verbose=1)
    problem = minimize(γ, Constraints...)
    Convex.solve!(problem, solver)
end

function make_robust_rnn_plot()

end

max_stars_repo_name: imanchester/REN

max_issues_repo_name: imanchester/REN

max_forks_repo_name: imanchester/REN


text:
#using SparseArrays

mutable struct woodbury_identity
    U::AbstractMatrix
    V::AbstractMatrix
    invA_U::AbstractMatrix
    factored_capacitance_matrix
    factored_A::abstract_linear_system_solver

    # (A + U*V')^-1
    function woodbury_identity(factored_A::abstract_linear_system_solver, U::AbstractArray, V::AbstractArray)
        invA_U = false
        capacitance_matrix = false
        try
          invA_U = ls_solve(factored_A, U)
          n = size(U,1);
          k = size(U,2);

          capacitance_matrix = (eye(k) + V * invA_U);
          factored_capacitance_matrix = lufact(capacitance_matrix);

          return new(U, V, invA_U, factored_capacitance_matrix, factored_A);
        catch e
            println("ERROR woodbury_identity")
            #@show full(factored_A._SparseMatrix)
            #@show capacitance_matrix
            #@show invA_U'
            throw(e)
        end
    end
end

function evaluate(W::woodbury_identity, sol::Array{Float64,1})
    return W.factored_A._SparseMatrix * sol + W.U * (W.V * sol)
end

function rel_error(W::woodbury_identity, rhs::Array{Float64,1}, sol::Array{Float64,1})
    return LinearAlgebra.norm(evaluate(W, sol) - rhs,1)/LinearAlgebra.norm(rhs,1);
end

function ls_solve(W::woodbury_identity, rhs::Array{Float64,1})
  try
    rhs_err = rhs
    sol = zeros(length(rhs))
    for i = 1:1
        # iterative refinement
        my_temp = ls_solve(W.factored_A, rhs_err)
        Δd = (my_temp - W.invA_U * (W.factored_capacitance_matrix \ (W.V * my_temp)))
        sol = sol + Δd[:]
        rhs_err = rhs - evaluate(W,sol)
        #@show i, LinearAlgebra.norm(rhs_err,1)
    end

    # catch errors
    tol = 1e-8;
    err = rel_error(W, rhs, sol) # why is this rel error?
    if err > tol
        my_warn("numerical instability using Woodbury, using direct factorization instead of woodbury.")
        sol = ls_solve_direct(W, rhs)
        err = rel_error(W, rhs, sol)

        if (err > tol)
            error("numerical stability using direct factorization, computation skipped")
        end
    end

    return sol
  catch e
      println("ERROR ls_solve")
      throw(e)
  end
end

function ls_solve_direct(W::woodbury_identity, rhs::Array{Float64,1}) # if there is numerical errors
    mat = W.factored_A._SparseMatrix + SparseArrays.sparse(W.U) * SparseArrays.sparse(W.V)
    return lufact(mat) \ rhs;
end

function numerical_error(A, sol, rhs)
    return LinearAlgebra.norm(A * sol - rhs,1);
end

function e_(i::Int64,n::Int64)
    @assert(i <= n)
    vec = zeros(n);
    vec[i] = 1.0;
    return vec;
end


text:
using Test, ARFIMA, Statistics

N = 1000
σ = 2.3

x = arfima(N, σ, nothing) # white noise
@test abs(std(x) - σ) < 1
@test length(x) == N

x = arfima(N, σ, nothing, SVector(0.2, 0.2))
@test length(x) == N
@test all(!isinf, x)
x = arfima(N, σ, nothing, SVector(0.2,0.2), SVector(0.1))
@test length(x) == N
@test all(!isinf, x)
x = arfima(N, σ, 1, SVector(0.02,0.02), SVector(0.01))
@test length(x) == N
@test all(!isinf, x)
x = arfima(N, σ, 0.25, SVector(0.02,0.02), SVector(0.01))
@test length(x) == N
@test all(!isinf, x)


text:
mutable struct SimulatedCamera <: Camera
    trigger_source::Channel{UInt64}
    image_source::Channel
    isopen::Bool
    isrunning::Bool
    current_id::Int
    SimulatedCamera(image_source::Channel) = new(Channel{UInt64}(Inf), image_source, false, false, 0)
    SimulatedCamera(trigger_source::Channel{UInt64}, image_source::Channel) = new(trigger_source, image_source, false, false, 0)
    function SimulatedCamera(trigger_period::Real, image_source::Channel)
        function produce_triggers!(trigger_source::Channel{UInt64})
            while true
                sleep(trigger_period)
                put!(trigger_source, time_ns())
            end
        end
        trigger_source = Channel(produce_triggers!; ctype = UInt64)
        new(trigger_source, image_source, false, false, 0)
    end
end

isopen(camera::SimulatedCamera) = camera.isopen
open!(camera::SimulatedCamera) = camera.isopen = true
close!(camera::SimulatedCamera) = camera.isopen = false

isrunning(camera::SimulatedCamera) = camera.isrunning
start!(camera::SimulatedCamera) = camera.isrunning = true

function stop!(camera::SimulatedCamera)
    close(camera.image_source)
    close(camera.trigger_source)
    camera.isrunning = false
end

function take!(camera::SimulatedCamera)
    @debug "Taking image from $(camera.image_source)"
    image = take!(camera.image_source)
    id = camera.current_id += 1
    @debug "Taking trigger from $(camera.trigger_source)"
    timestamp = take!(camera.trigger_source)
    return AcquiredImage(image, id, timestamp)
end

function trigger!(camera::SimulatedCamera)
    put!(camera.trigger_source, time_ns())
end

function Base.show(io::IO, camera::SimulatedCamera)
    write(io, "$(nameof(SimulatedCamera))(trigger_source:$(camera.trigger_source), image_source:$(camera.image_source), current_id:$(camera.current_id))")
end


text:
# debugging and plotting utilities for classification

mutable struct TuningParameters
  Mfair::Int # Cluster distribution approximation accuracy, more is better, computation becomes prohibative beyond ~1000
  EMiters::Int # expectation maximization iterations to refine Cluster distribution density estimate and sample classification, diminishing returns beyond ~30
  TuningParameters(;Mfair=200, EMiters=30) = new(Mfair, EMiters)
end

mutable struct ClassificationSystem
  numCategories::Int
  categoryLabels::Dict{Int, String}
  colors::Dict{Int, String}
  expertBelief::Dict{Int, BallTreeDensity}
  temporalBelief::Dict{Int, BallTreeDensity}
  currentBelief::Dict{Int, BallTreeDensity}
  assignment::Vector{Int}
end

mutable struct SampleData
  samples::Array{Float64,2}
  measurementDims::Int
  numSamples::Int
end

mutable struct DataGroundTruth
  assignment::Array{Int,1}
  clustersizes::Dict{Int,Int}
  clusters::Dict{Int,Array{Float64,2}}
  plt_lyr_cluster::Dict{Int, Any}
  plt_lyr_cluster_nocolor::Dict{Int, Any}
  DataGroundTruth() = new(zeros(Int, 0),  Dict{Int,Int}(),  Dict{Int,Array{Float64,2}}(),  Dict{Int,Any}(),  Dict{Int,Any}() )
  DataGroundTruth(x...) = new(x[1],x[2],x[3],x[4],x[5])
end

# Common struct to store debug information during running of the algorithm
mutable struct DebugResults
  ASSIGNED::Array{Array{Int,1},1}
  PL_MEAS::Array{Any,1}
  ACCUR_C::Dict{Int,Array{Float64,1}}
  REL_ACCUR_C::Dict{Int,Array{Float64,1}}
  INDV_MISASSIGN_C::Dict{Int,Array{Float64,1}}
  DebugResults() = new(Array{Array{Int,1},1}(),
                      [],
                      Dict{Int,Array{Float64,1}}(),
                      Dict{Int,Array{Float64,1}}(),
                      Dict{Int,Array{Float64,1}}() )
end
function defaultDebugResults()
  dbg = DebugResults()
  dbg.ACCUR_C[1] = Float64[]
  dbg.ACCUR_C[2] = Float64[]

  dbg.REL_ACCUR_C[1] = Float64[]
  dbg.REL_ACCUR_C[2] = Float64[]

  dbg.INDV_MISASSIGN_C[1] = Float64[]
  dbg.INDV_MISASSIGN_C[2] = Float64[]

  return dbg
end

mutable struct ClassificationStats
  POPFRAC::Array{Float64,2}
  SEQKLDIVERG::Dict{Int, Vector{Float64}}
  ESTBELIEF::Dict{Int, Vector{BallTreeDensity}}
  ClassificationStats() = new(zeros(0,0),
                              Dict{Int, Vector{Float64}}(),
                              Dict{Int, Vector{BallTreeDensity}}()   )
end
function defaultClassificationStats(params::TuningParameters, cs::ClassificationSystem)
  stats = ClassificationStats()
  for clbl in cs.categoryLabels
    stats.SEQKLDIVERG[clbl[1]] = Float64[]
    stats.ESTBELIEF[clbl[1]] = BallTreeDensity[]
  end
  stats.POPFRAC = zeros(cs.numCategories,params.EMiters) #Array{Float64,2}()
  stats
end

# attempt at balancing populations of different size (non-textbook, there is no textbook)
# using Dirichlet distribution as conjugate prior of categorical distribution
# a is fractions of current population classification estimates, sum(a) = 1, a.>= 0
# b is current likelihood estimate of being classified label length(b), and should ==length(a); sum(b) = 1, b.>=0
# yes, function needs a better name
function sdc2(a,b)
  aa = a.+0.6*abs(maximum(a)-minimum(a))#1#./maximum(a)
  p = rand(Dirichlet(aa))
  pp = p.*b
  pp /= sum(pp)
  rand(Categorical(pp))
end

function packDebugResults!(dbg::DebugResults, assigned, GT)

  push!(dbg.ACCUR_C[1], 100*abs(sum(assigned .== 1)-GT.clustersizes[1])/(GT.clustersizes[1]+0.0))
  push!(dbg.ACCUR_C[2], 100*abs(sum(assigned .== 2)-GT.clustersizes[2])/(GT.clustersizes[2]+0.0))

  push!(dbg.REL_ACCUR_C[1], 100*abs(sum(assigned .== 1)-GT.clustersizes[1])/(GT.clustersizes[1]+GT.clustersizes[2]))
  push!(dbg.REL_ACCUR_C[2], 100*abs(sum(assigned .== 2)-GT.clustersizes[2])/(GT.clustersizes[1]+GT.clustersizes[2]))

  push!(dbg.INDV_MISASSIGN_C[1], 100*abs(sum(assigned .== 1)-GT.clustersizes[1])/(GT.clustersizes[1]+0.0) )
  push!(dbg.INDV_MISASSIGN_C[2], 100*abs(sum(assigned .== 2)-GT.clustersizes[2])/(GT.clustersizes[2]+0.0) )

  push!(dbg.ASSIGNED, assigned)

  nothing
end

function EMClassificationRun!(
      params::TuningParameters,
      cs::ClassificationSystem,
      data::SampleData;
      debug::NothingUnion{DebugResults}=nothing,
      groundtruth::NothingUnion{DataGroundTruth}=nothing )

  # allocate some common variables
  stats = defaultClassificationStats(params, cs)
  idxs = collect(1:data.numSamples)
  weights = zeros(cs.numCategories, data.numSamples)

  ### [1.0] Classification algorithm part 1 of 2 -- the measurement, generalized expectation maximization (EM) procedure
  # [1.1] loop over cluster assignment sample drawing process
  for m in 1:params.EMiters # EM iterations per measurement cycle, this is a tuning parameter -> more iterations should assymptote to true assignment distribution
    # [1.2] numerically approximate a Bayesian prior -- (interpret as Dirichlet sampling process with some implicit concentration parameter)

    # find the joint distribution using expert prior, EM measurement iteration (HMM prediction not shown in this example yet)
    if m>1  for clbl in cs.categoryLabels
      cs.currentBelief[clbl[1]] = cs.currentBelief[clbl[1]] * cs.expertBelief[clbl[1]]
    end end

    # evaluate the point likelihood of being in each cluster
    for i in 1:cs.numCategories
      weights[i,:] = evaluateDualTree(cs.currentBelief[i], data.samples)
    end

    # [1.3] assign cluster as sampling from categorical distribution
    popFrac = [sum(cs.assignment .== 1); sum(cs.assignment .== 2)]/(data.numSamples+0.0)
    stats.POPFRAC[:,m] = popFrac'
    for i in 1:data.numSamples
      p = vec(weights[:,i])
      p = p/sum(p) # renormalize to valid probability distribution
      sel = false ?   rand( Categorical( p ) )   :   round(Int,sdc2(popFrac, p))
      cs.assignment[i] =  sel
    end

    # build new intermediate cluster density estimates
    # but subsample to fair chance and practically computable sets first
    clusters = Dict{Int, Array{Float64,2}}()
    subSampledClusters = Dict{Int, Array{Float64,2}}()
    for clbl in cs.categoryLabels
      msk = cs.assignment .== clbl[1]
      @show size(msk), size(data.samples)
      clusters[clbl[1]] = data.samples[:, msk ]
      len = size(clusters[clbl[1]],2)

      useidx = StatsBase.sample(idxs[1:size(clusters[clbl[1]],2)], ProbabilityWeights(1.0/len*ones(len)), params.Mfair, replace=false)
      subSampledClusters[clbl[1]] = clusters[clbl[1]][:,useidx]

      newBel = kde!(subSampledClusters[clbl[1]])
      diverg = kld(newBel, cs.currentBelief[clbl[1]])[1]
      cs.currentBelief[clbl[1]] = newBel
      push!(stats.SEQKLDIVERG[clbl[1]], diverg)
      push!(stats.ESTBELIEF[clbl[1]], deepcopy(cs.currentBelief[clbl[1]]))
    end

    # measure classification accuracy
    if debug != nothing && groundtruth != nothing
      packDebugResults!(debug, cs.assignment, groundtruth)
    end
  end
  stats
end

function dispersePoints(X::Array{Float64,2}, dispersion)
  r, c = size(X)
  Y = zeros(r,c)
  for i in 1:c
    Y[:,i] = X[:,i] + rand(dispersion)
  end
  return Y
end

#