text:
"""
Fetch the "natural" block size of this band.

GDAL contains a concept of the natural block size of rasters so that
applications can organized data access efficiently for some file formats.
The natural block size is the block size that is most efficient for accessing
the format. For many formats this is simple a whole scanline in which case
`*pnXSize` is set to `GetXSize()`, and *pnYSize is set to 1.

However, for tiled images this will typically be the tile size.

Note that the X and Y block sizes don't have to divide the image size evenly,
meaning that right and bottom edge blocks may be incomplete. See `ReadBlock()`
for an example of code dealing with these issues.
"""
function getblocksize(rasterband::RasterBand)
    xy = Array(Cint, 2); x = pointer(xy); y = x + sizeof(Cint)
    GDAL.getblocksize(rasterband, x, y)
    xy
end

"Fetch the pixel data type for this band."
getdatatype(rasterband::RasterBand) =
    _jltype[GDAL.getrasterdatatype(rasterband)]

"Fetch the width in pixels of this band."
width(rasterband::RasterBand) = GDAL.getrasterbandxsize(rasterband)

"Fetch the height in pixels of this band."
height(rasterband::RasterBand) = GDAL.getrasterbandysize(rasterband)

"Find out if we have update permission for this band."
getaccess(rasterband::RasterBand) = GDAL.getrasteraccess(rasterband)

"Fetch a band object for a dataset from its index"
fetchband(dataset::Dataset, i::Integer) =
    GDAL.getrasterband(dataset, i)

"""Fetch the band number (1+) within its dataset, or 0 if unknown.

This method may return a value of 0 to indicate overviews, or free-standing
`GDALRasterBand` objects without a relationship to a dataset.
"""
indexof(rasterband::RasterBand) = GDAL.getbandnumber(rasterband)

"""
Fetch the handle to its dataset handle, or `NULL` if this cannot be determined.

Note that some `GDALRasterBands` are not considered to be a part of a dataset,
such as overviews or other "freestanding" bands.
"""
getdataset(rasterband::RasterBand) =
    GDAL.getbanddataset(rasterband)

"""
Return raster unit type.

Return a name for the units of this raster's values. For instance, it might be
"m" for an elevation model in meters, or "ft" for feet.
"""
getunittype(rasterband::RasterBand) = GDAL.getrasterunittype(rasterband)

"Set unit type."
setunittype(rasterband::RasterBand, unitstring::AbstractString) =
    GDAL.setrasterunittype(rasterband, unitstring)

"""
Fetch the raster value offset.

This (in combination with `GetScale()`) is used to transform raw pixel values
into the units returned by `GetUnits()`. For e.g. this might be used to store
elevations in `GUInt16` bands with a precision of 0.1, starting from -100.

    Units value = (raw value * scale) + offset

For file formats that don't know this intrinsically, a value of 0 is returned.
"""
getoffset(rasterband::RasterBand) =
    GDAL.getrasteroffset(rasterband, C_NULL)

"Set scaling offset."
setoffset(rasterband::RasterBand, offset::Cdouble) =
    GDAL.setrasteroffset(band, offset)

"""
Fetch the raster value scale.

This value (in combination with the `GetOffset()` value) is used to transform
raw pixel values into the units returned by `GetUnits()`. For example this
might be used to store elevations in GUInt16 bands with a precision of 0.1,
and starting from -100.

    Units value = (raw value * scale) + offset

For file formats that don't know this intrinsically a value of one is returned.
"""
getscale(rasterband::RasterBand) = GDAL.getrasterscale(rasterband, C_NULL)

"Set scaling ratio."
setscale(rasterband::RasterBand, scale::Cdouble) =
    GDAL.setrasterscale(rasterband, scale)

"""
Fetch the no data value for this band.

If there is no out of data value, an out of range value will generally be
returned. The no data value for a band is generally a special marker value
used to mark pixels that are not valid data. Such pixels should generally
not be displayed, nor contribute to analysis operations.

### Parameters
* `pbSuccess`   pointer to a boolean to use to indicate if a value is actually
associated with this layer. May be `NULL` (default).

### Returns
the nodata value for this band.
"""
getnodatavalue(rasterband::RasterBand) =
    GDAL.getrasternodatavalue(rasterband, C_NULL)

"Set the no data value for this band."
function setnodatavalue(rasterband::RasterBand, value::Real)
    result = GDAL.setrasternodatavalue(rasterband, value)
    (result == GDAL.CE_Failure) && error("Could not set nodatavalue")
end

# """
#     GDALDeleteRasterNoDataValue(GDALRasterBandH hBand) -> CPLErr
# Remove the no data value for this band.
# """
# function deleterasternodatavalue{T <: GDALRasterBandH}(arg1::Ptr{T})
#     ccall((:GDALDeleteRasterNoDataValue,libgdal),CPLErr,(Ptr{GDALRasterBandH},),arg1)
# end

# """
#     GDALGetRasterCategoryNames(GDALRasterBandH hBand) -> char **
# Fetch the list of category names for this raster.
# """
# function getrastercategorynames{T <: GDALRasterBandH}(arg1::Ptr{T})
#     bytestring(unsafe_load(ccall((:GDALGetRasterCategoryNames,libgdal),Ptr{Cstring},(Ptr{GDALRasterBandH},),arg1)))
# end

# """
#     GDALSetRasterCategoryNames(GDALRasterBandH hBand,
#                                char ** papszNames) -> CPLErr
# Set the category names for this band.
# """
# function setrastercategorynames{T <: GDALRasterBandH}(arg1::Ptr{T},arg2)
#     ccall((:GDALSetRasterCategoryNames,libgdal),CPLErr,(Ptr{GDALRasterBandH},Ptr{Cstring}),arg1,arg2)
# end

# """
#     GDALGetRasterMinimum(GDALRasterBandH hBand,
#                          int * pbSuccess) -> double
# Fetch the minimum value for this band.
# """
# function getrasterminimum{T <: GDALRasterBandH}(arg1::Ptr{T},pbSuccess)
#     ccall((:GDALGetRasterMinimum,libgdal),Cdouble,(Ptr{GDALRasterBandH},Ptr{Cint}),arg1,pbSuccess)
# end

# """
#     GDALGetRasterMaximum(GDALRasterBandH hBand,
#                          int * pbSuccess) -> double
# Fetch the maximum value for this band.
# """
# function getrastermaximum{T <: GDALRasterBandH}(arg1::Ptr{T},pbSuccess)
#     ccall((:GDALGetRasterMaximum,libgdal),Cdouble,(Ptr{GDALRasterBandH},Ptr{Cint}),arg1,pbSuccess)
# end

# """
#     GDALGetRasterStatistics(GDALRasterBandH hBand,
#                             int bApproxOK,
#                             int bForce,
#                             double * pdfMin,
#                             double * pdfMax,
#                             double * pdfMean,
#                             double * pdfStdDev) -> CPLErr
# Fetch image statistics.
# """
# function getrasterstatistics{T <: GDALRasterBandH}(arg1::Ptr{T},bApproxOK::Integer,bForce::Integer,pdfMin,pdfMax,pdfMean,pdfStdDev)
#     ccall((:GDALGetRasterStatistics,libgdal),CPLErr,(Ptr{GDALRasterBandH},Cint,Cint,Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble}),arg1,bApproxOK,bForce,pdfMin,pdfMax,pdfMean,pdfStdDev)
# end

# """
#     GDALComputeRasterStatistics(GDALRasterBandH hBand,
#                                 int bApproxOK,
#                                 double * pdfMin,
#                                 double * pdfMax,
#                                 double * pdfMean,
#                                 double * pdfStdDev,
#                                 GDALProgressFunc pfnProgress,
#                                 void * pProgressData) -> CPLErr
# Compute image statistics.
# """
# function computerasterstatistics{T <: GDALRasterBandH}(arg1::Ptr{T},bApproxOK::Integer,pdfMin,pdfMax,pdfMean,pdfStdDev,pfnProgress::Ptr{GDALProgressFunc},pProgressData)
#     ccall((:GDALComputeRasterStatistics,libgdal),CPLErr,(Ptr{GDALRasterBandH},Cint,Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{GDALProgressFunc},Ptr{Void}),arg1,bApproxOK,pdfMin,pdfMax,pdfMean,pdfStdDev,pfnProgress,pProgressData)
# end

# """
#     GDALSetRasterStatistics(GDALRasterBandH hBand,
#                             double dfMin,
#                             double dfMax,
#                             double dfMean,
#                             double dfStdDev) -> CPLErr
# Set statistics on band.
# """
# function setrasterstatistics{T <: GDALRasterBandH}(hBand::Ptr{T},dfMin::Real,dfMax::Real,dfMean::Real,dfStdDev::Real)
#     ccall((:GDALSetRasterStatistics,libgdal),CPLErr,(Ptr{GDALRasterBandH},Cdouble,Cdouble,Cdouble,Cdouble),hBand,dfMin,dfMax,dfMean,dfStdDev)
# end

# """
#     GDALComputeRasterMinMax(GDALRasterBandH hBand,
#                             int bApproxOK,
#                             double adfMinMax) -> void
# Compute the min/max values for a band.
# """
# function computerasterminmax{T <: GDALRasterBandH}(hBand::Ptr{T},bApproxOK::Integer,adfMinMax::Array_2_Cdouble)
#     ccall((:GDALComputeRasterMinMax,libgdal),Void,(Ptr{GDALRasterBandH},Cint,Array_2_Cdouble),hBand,bApproxOK,adfMinMax)
# end

# """
#     GDALGetRasterHistogram(GDALRasterBandH hBand,
#                            double dfMin,
#                            double dfMax,
#                            int nBuckets,
#                            int * panHistogram,
#                            int bIncludeOutOfRange,
#                            int bApproxOK,
#                            GDALProgressFunc pfnProgress,
#                            void * pProgressData) -> CPLErr
# Compute raster histogram.
# """
# function getrasterhistogram{T <: GDALRasterBandH}(hBand::Ptr{T},dfMin::Real,dfMax::Real,nBuckets::Integer,panHistogram,bIncludeOutOfRange::Integer,bApproxOK::Integer,pfnProgress::Ptr{GDALProgressFunc},pProgressData)
#     ccall((:GDALGetRasterHistogram,libgdal),CPLErr,(Ptr{GDALRasterBandH},Cdouble,Cdouble,Cint,Ptr{Cint},Cint,Cint,Ptr{GDALProgressFunc},Ptr{Void}),hBand,dfMin,dfMax,nBuckets,panHistogram,bIncludeOutOfRange,bApproxOK,pfnProgress,pProgressData)
# end

# """
#     GDALGetRasterHistogramEx(GDALRasterBandH hBand,
#                              double dfMin,
#                              double dfMax,
#                              int nBuckets,
#                              GUIntBig * panHistogram,
#                              int bIncludeOutOfRange,
#                              int bApproxOK,
#                              GDALProgressFunc pfnProgress,
#                              void * pProgressData) -> CPLErr
# Compute raster histogram.
# """
# function getrasterhistogramex{T <: GDALRasterBandH}(hBand::Ptr{T},dfMin::Real,dfMax::Real,nBuckets::Integer,panHistogram,bIncludeOutOfRange::Integer,bApproxOK::Integer,pfnProgress::Ptr{GDALProgressFunc},pProgressData)
#     ccall((:GDALGetRasterHistogramEx,libgdal),CPLErr,(Ptr{GDALRasterBandH},Cdouble,Cdouble,Cint,Ptr{GUIntBig},Cint,Cint,Ptr{GDALProgressFunc},Ptr{Void}),hBand,dfMin,dfMax,nBuckets,panHistogram,bIncludeOutOfRange,bApproxOK,pfnProgress,pProgressData)
# end

# """
#     GDALGetDefaultHistogram(GDALRasterBandH hBand,
#                             double * pdfMin,
#                             double * pdfMax,
#                             int * pnBuckets,
#                             int ** ppanHistogram,
#                             int bForce,
#                             GDALProgressFunc pfnProgress,
#                             void * pProgressData) -> CPLErr
# Fetch default raster histogram.
# """
# function getdefaulthistogram{T <: GDALRasterBandH}(hBand::Ptr{T},pdfMin,pdfMax,pnBuckets,ppanHistogram,bForce::Integer,pfnProgress::Ptr{GDALProgressFunc},pProgressData)
#     ccall((:GDALGetDefaultHistogram,libgdal),CPLErr,(Ptr{GDALRasterBandH},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cint},Ptr{Ptr{Cint}},Cint,Ptr{GDALProgressFunc},Ptr{Void}),hBand,pdfMin,pdfMax,pnBuckets,ppanHistogram,bForce,pfnProgress,pProgressData)
# end

# """
#     GDALGetDefaultHistogramEx(GDALRasterBandH hBand,
#                               double * pdfMin,
#                               double * pdfMax,
#                               int * pnBuckets,
#                               GUIntBig ** ppanHistogram,
#                               int bForce,
#                               GDALProgressFunc pfnProgress,
#                               void * pProgressData) -> CPLErr
# Fetch default raster histogram.
# """
# function getdefaulthistogramex{T <: GDALRasterBandH}(hBand::Ptr{T},pdfMin,pdfMax,pnBuckets,ppanHistogram,bForce::Integer,pfnProgress::Ptr{GDALProgressFunc},pProgressData)
#     ccall((:GDALGetDefaultHistogramEx,libgdal),CPLErr,(Ptr{GDALRasterBandH},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cint},Ptr{Ptr{GUIntBig}},Cint,Ptr{GDALProgressFunc},Ptr{Void}),hBand,pdfMin,pdfMax,pnBuckets,ppanHistogram,bForce,pfnProgress,pProgressData)
# end

# """
#     GDALSetDefaultHistogram(GDALRasterBandH hBand,
#                             double dfMin,
#                             double dfMax,
#                             int nBuckets,
#                             int * panHistogram) -> CPLErr
# Set default histogram.
# """
# function setdefaulthistogram{T <: GDALRasterBandH}(hBand::Ptr{T},dfMin::Real,dfMax::Real,nBuckets::Integer,panHistogram)
#     ccall((:GDALSetDefaultHistogram,libgdal),CPLErr,(Ptr{GDALRasterBandH},Cdouble,Cdouble,Cint,Ptr{Cint}),hBand,dfMin,dfMax,nBuckets,panHistogram)
# end

# """
#     GDALSetDefaultHistogramEx(GDALRasterBandH hBand,
#                               double dfMin,
#                               double dfMax,
#                               int nBuckets,
#                               GUIntBig * panHistogram) -> CPLErr
# Set default histogram.
# """
# function setdefaulthistogramex{T <: GDALRasterBandH}(hBand::Ptr{T},dfMin::Real,dfMax::Real,nBuckets::Integer,panHistogram)
#     ccall((:GDALSetDefaultHistogramEx,libgdal),CPLErr,(Ptr{GDALRasterBandH},Cdouble,Cdouble,Cint,Ptr{GUIntBig}),hBand,dfMin,dfMax,nBuckets,panHistogram)
# end

# """
#     GDALGetRandomRasterSample(GDALRasterBandH,
#                               int,
#                               float *) -> int
# """
# function getrandomrastersample{T <: GDALRasterBandH}(arg1::Ptr{T},arg2::Integer,arg3)
#     ccall((:GDALGetRandomRasterSample,libgdal),Cint,(Ptr{GDALRasterBandH},Cint,Ptr{Cfloat}),arg1,arg2,arg3)
# end

# """
#     GDALComputeBandStats(GDALRasterBandH hBand,
#                          int nSampleStep,
#                          double * pdfMean,
#                          double * pdfStdDev,
#                          GDALProgressFunc pfnProgress,
#                          void * pProgressData) -> CPLErr
# """
# function computebandstats{T <: GDALRasterBandH}(hBand::Ptr{T},nSampleStep::Integer,pdfMean,pdfStdDev,pfnProgress::Ptr{GDALProgressFunc},pProgressData)
#     ccall((:GDALComputeBandStats,libgdal),CPLErr,(Ptr{GDALRasterBandH},Cint,Ptr{Cdouble},Ptr{Cdouble},Ptr{GDALProgressFunc},Ptr{Void}),hBand,nSampleStep,pdfMean,pdfStdDev,pfnProgress,pProgressData)
# end

# """
#     GDALOverviewMagnitudeCorrection(GDALRasterBandH hBaseBand,
#                                     int nOverviewCount,
#                                     GDALRasterBandH * pahOverviews,
#                                     GDALProgressFunc pfnProgress,
#                                     void * pProgressData) -> CPLErr
# """
# function overviewmagnitudecorrection{T <: GDALRasterBandH}(hBaseBand::Ptr{T},nOverviewCount::Integer,pahOverviews,pfnProgress::Ptr{GDALProgressFunc},pProgressData)
#     ccall((:GDALOverviewMagnitudeCorrection,libgdal),CPLErr,(Ptr{GDALRasterBandH},Cint,Ptr{GDALRasterBandH},Ptr{GDALProgressFunc},Ptr{Void}),hBaseBand,nOverviewCount,pahOverviews,pfnProgress,pProgressData)
# end

# """
#     GDALGetDefaultRAT(GDALRasterBandH hBand) -> GDALRasterAttributeTableH
# Fetch default Raster Attribute Table.
# """
# function getdefaultrat{T <: GDALRasterBandH}(hBand::Ptr{T})
#     checknull(ccall((:GDALGetDefaultRAT,libgdal),Ptr{GDALRasterAttributeTableH},(Ptr{GDALRasterBandH},),hBand))
# end

# """
#     GDALSetDefaultRAT(GDALRasterBandH hBand,
#                       GDALRasterAttributeTableH hRAT) -> CPLErr
# Set default Raster Attribute Table.
# """
# function setdefaultrat{T <: GDALRasterBandH}(arg1::Ptr{T},arg2::Ptr{GDALRasterAttributeTableH})
#     ccall((:GDALSetDefaultRAT,libgdal),CPLErr,(Ptr{GDALRasterBandH},Ptr{GDALRasterAttributeTableH}),arg1,arg2)
# end

# """
#     GDALAddDerivedBandPixelFunc(const char * pszName,
#                                 GDALDerivedPixelFunc pfnPixelFunc) -> CPLErr
# This adds a pixel function to the global list of available pixel functions for derived bands.
# ### Parameters
# * **pszFuncName**: Name used to access pixel function
# * **pfnNewFunction**: Pixel function associated with name. An existing pixel function registered with the same name will be replaced with the new one.
# ### Returns
# CE_None, invalid (NULL) parameters are currently ignored.
# """
# function addderivedbandpixelfunc(pszName,pfnPixelFunc::Ptr{GDALDerivedPixelFunc})
#     ccall((:GDALAddDerivedBandPixelFunc,libgdal),CPLErr,(Cstring,Ptr{GDALDerivedPixelFunc}),pszName,pfnPixelFunc)
# end

"""
Copy all raster band raster data.

This function copies the complete raster contents of one band to another
similarly configured band. The source and destination bands must have the same
width and height. The bands do not have to have the same data type.

It implements efficient copying, in particular "chunking" the copy in
substantial blocks.

Currently the only `papszOptions` value supported is : "COMPRESSED=YES" to
force alignment on target dataset block sizes to achieve best compression.
More options may be supported in the future.

### Parameters
* `hSrcBand`        the source band
* `hDstBand`        the destination band
* `papszOptions`    transfer hints in "StringList" Name=Value format.
* `pfnProgress`     progress reporting function.
* `pProgressData`   callback data for progress function.

### Returns
`CE_None` on success, or `CE_Failure` on failure.
"""
copywholeraster(source::RasterBand, dest::RasterBand, options) =
    GDAL.rasterbandcopywholeraster(source, dest, pointer(options))

copywholeraster(source::RasterBand, dest::RasterBand) =
    GDAL.rasterbandcopywholeraster(source, dest, C_NULL,C_NULL,C_NULL)

"""
Return the number of overview layers available.

### Returns
overview count, zero if none.
"""
noverview(rasterband::RasterBand) = GDAL.getoverviewcount(rasterband)

"Fetch overview raster band object."
fetchoverview(rasterband::RasterBand, i::Integer) =
    GDAL.getoverview(rasterband, i-1)

"""
Fetch best overview satisfying `nsamples` number of samples.

Returns the most reduced overview of the given band that still satisfies the
desired number of samples `nsamples`. This function can be used with zero as the
number of desired samples to fetch the most reduced overview. The same band as was
passed in will be returned if it has not overviews, or if none of the overviews
have enough samples.
"""
getsampleoverview(rasterband::RasterBand, nsamples::Integer) =
    GDAL.getrastersampleoverview(rasterband, nsamples)

# """
#     GDALGetRasterSampleOverviewEx(GDALRasterBandH hBand,
#                                   GUIntBig nDesiredSamples) -> GDALRasterBandH
# Fetch best sampling overview.
# """
# function getrastersampleoverviewex{T <: GDALRasterBandH}(arg1::Ptr{T},arg2::GUIntBig)
#     checknull(ccall((:GDALGetRasterSampleOverviewEx,libgdal),Ptr{GDALRasterBandH},(Ptr{GDALRasterBandH},GUIntBig),arg1,arg2))
# end

"Color Interpretation value for band"
getcolorinterp(rasterband::RasterBand) =
    GDAL.getrastercolorinterpretation(rasterband)

"Set color interpretation of a band."
setcolorinterp(rasterband::RasterBand, color::GDAL.GDALColorInterp) =
    GDAL.setrastercolorinterpretation(rasterband, color)

"""
Fetch the color table associated with band.

If there is no associated color table, the return result is `NULL`. The
returned color table remains owned by the `GDALRasterBand`, and can't be
depended on for long, nor should it ever be modified by the caller.
"""
getcolortable(rasterband::RasterBand) =
    GDAL.getrastercolortable(rasterband)

"""
Set the raster color table.

The driver will make a copy of all desired data in the colortable. It remains
owned by the caller after the call.

### Parameters
* `poCT`    the color table to apply. This may be NULL to clear the color table
(where supported).

### Returns
`CE_None` on success, or `CE_Failure` on failure. If the action is unsupported
by the driver, a value of `CE_Failure` is returned, but no error is issued.
"""
setcolortable(rasterband::RasterBand, colortable::Ptr{GDAL.GDALColorTableH}) =
    GDAL.setrastercolortable(rasterband, colortable)

clearcolortable(rasterband::RasterBand) =
    GDAL.setrastercolortable(rasterband, C_NULL)

"""
Generate downsampled overviews.

This function will generate one or more overview images from a base image using
the requested downsampling algorithm. Its primary use is for generating
overviews via BuildOverviews(), but it can also be used to generate downsampled
images in one file from another outside the overview architecture.

The output bands need to exist in advance.

This function will honour properly `NODATA_VALUES` tuples (special dataset
metadata) so that only a given RGB triplet (in case of a RGB image) will be
considered as the nodata value and not each value of the triplet independantly
per band.

### Parameters
* `hSrcBand`        the source (base level) band.
* `nOverviewCount`  the number of downsampled bands being generated.
* `pahOvrBands`     the list of downsampled bands to be generated.
* `pszResampling`   Resampling algorithm (eg. "AVERAGE").
* `pfnProgress`     progress report function.
* `pProgressData`   progress function callback data.

### Returns
`CE_None` on success or `CE_Failure` on failure.
"""
function regenerateoverviews(rasterband::RasterBand,
                             noverviews::Integer,
                             overviewbands::Vector{Ptr{GDAL.GDALRasterBandH}},
                             resampling::AbstractString = "NEAREST")
    result = GDAL.regenerateoverviews(rasterband, noverviews, overviewbands,
                                      resampling,
                                      ProgressFunc(C_NULL),
                                      C_NULL)
    (result == GDAL.CE_Failure) && error("Failed to regenerate overviews")
end

# "Advise driver of upcoming read requests."
# _rasteradviseread(hRB::GDALRasterBandH,
#                   nDSXOff::Integer,
#                   nDSYOff::Integer,
#                   nDSXSize::Integer,
#                   nDSYSize::Integer,
#                   nBXSize::Integer,
#                   nBYSize::Integer,
#                   eBDataType::GDALDataType,
#                   papszOptions::Ptr{Ptr{UInt8}}) =
#     GDALRasterAdviseRead(hRB, nDSXOff, nDSYOff, nDSXSize, nDSYSize, nBXSize,
#                          nBYSize, eBDataType, papszOptions)::CPLErr

"""
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().

See the `GetLockedBlockRef()` method for a way of accessing internally cached
block oriented data without an extra copy into an application buffer.

### Parameters
* `nXBlockOff`  the horizontal block offset, with zero indicating the left most
block, 1 the next block and so forth.
* `nYBlockOff`  the vertical block offset, with zero indicating the left most
block, 1 the next block and so forth.
* `pImage`      the buffer into which the data will be read. The buffer must be
large enough to hold `GetBlockXSize()*GetBlockYSize()` words of type
`GetRasterDataType()`.

### Returns
`CE_None` on success or `CE_Failure` on an error.
"""
function readblock!(band::RasterBand, nxoff::Integer, nyoff::Integer,
                    buffer::Array)
    result = GDAL.readblock(band, nxoff, nyoff, pointer(buffer))
    (result != GDAL.CE_None) && error("Failed to read block")
    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()`.

See `ReadBlock()` for an example of block oriented data access.

### Parameters
* `nXBlockOff`  the horizontal block offset, with zero indicating the left most
block, 1 the next block and so forth.
* `nYBlockOff`  the vertical block offset, with zero indicating the left most
block, 1 the next block and so forth.
* `pImage`      the buffer from which the data will be written. The buffer must
be large enough to hold `GetBlockXSize()*GetBlockYSize()` words of type
`GetRasterDataType()`.

### Returns
`CE_None` on success or `CE_Failure` on an error.
"""
function writeblock!(band::RasterBand, nxoff::Integer, nyoff::Integer,
                     image::Array)
    result = GDAL.writeblock(band, nxoff, nyoff, pointer(image))
    (result != GDAL.CE_None) && error("Failed to write block")
    image
end

"""
Check for arbitrary overviews.

This returns `TRUE` if the underlying datastore can compute arbitrary overviews
efficiently, e.g. with OGDI over a network. Datastores with arbitrary overviews
don't generally have any fixed overviews, but the `RasterIO()` method can be
used in downsampling mode to get overview data efficiently.

### Returns
`TRUE` if arbitrary overviews available (efficiently), otherwise `FALSE`.
"""
hasarbitraryoverviews(band::RasterBand) =
    Bool(GDAL.hasarbitraryoverviews(band))

"""
Fetch the list of category names for this raster.

The return list is a "StringList" in the sense of the CPL functions. That is a
NULL terminated array of strings. Raster values without associated names will
have an empty string in the returned list. The first entry in the list is for
raster values of zero, and so on.

The returned stringlist should not be altered or freed by the application.
It may change on the next GDAL call, so please copy it if it is needed for any
period of time.

### Returns
list of names, or `NULL` if none.
"""
getcategorynames(band::RasterBand) =
    loadstringlist(GDAL.C.GDALGetRasterCategoryNames(band))

"Set the category names for this band."
function setrastercategorynames(band::RasterBand, names::Vector{ASCIIString})
    result = GDAL.setrastercategorynames(band, Ptr{Ptr{UInt8}}(pointer(names)))
    (result != GDAL.CE_None) && error("Failed to set category names for this band")
end

"""
Flush raster data cache.

This call will recover memory used to cache data blocks for this raster band,
and ensure that new requests are referred to the underlying driver.

### Returns
`CE_None` on success.
"""
function flushcache(band::RasterBand)
    result = GDAL.flushrastercache(band)
    (result != GDAL.CE_None) && error("Failed to flush raster data cache")
end

"""
Fill this band with a constant value.

GDAL makes no guarantees about what values pixels in newly created files are
set to, so this method can be used to clear a band to a specified "default"
value. The fill value is passed in as a double but this will be converted to
the underlying type before writing to the file. An optional second argument
allows the imaginary component of a complex constant value to be specified.

### Parameters
* `dfRealValue`         Real component of fill value
* `dfImaginaryValue`    Imaginary component of fill value, defaults to zero

### Returns
`CE_Failure` if the write fails, otherwise `CE_None`
"""
function fillraster(band::RasterBand,
                    realvalue::Cdouble,
                    imagvalue::Cdouble = 0)
    result = GDAL.fillraster(band, realvalue, imagvalue)
    (result != GDAL.CE_None) && error("Failed to fill raster band")
end

"""
Return the mask band associated with the band. (Since: GDAL 1.6.0)

The `GDALRasterBand` class includes a default implementation of `GetMaskBand()`
that returns one of four default implementations:

- If a corresponding .msk file exists it will be used for the mask band.
- If the dataset has a `NODATA_VALUES` metadata item, an instance of the new
GDALNoDataValuesMaskBand class will be returned. `GetMaskFlags()` will return
`GMF_NODATA | GMF_PER_DATASET`.
- If the band has a nodata value set, an instance of the new
`GDALNodataMaskRasterBand` class will be returned. `GetMaskFlags()` will
return `GMF_NODATA`.
- If there is no nodata value, but the dataset has an alpha band that seems to
apply to this band (specific rules yet to be determined) and that is of type
`GDT_Byte` then that alpha band will be returned, and the flags
`GMF_PER_DATASET` and `GMF_ALPHA` will be returned in the flags.
- If neither of the above apply, an instance of the new
`GDALAllValidRasterBand` class will be returned that has 255 values for all
pixels. The null flags will return `GMF_ALL_VALID`.

Note that the `GetMaskBand()` should always return a `GDALRasterBand` mask,
even if it is only an all 255 mask with the flags indicating `GMF_ALL_VALID`.

See also: http://trac.osgeo.org/gdal/wiki/rfc15_nodatabitmask

### Returns
a valid mask band.
"""
getmaskband(band::RasterBand) = GDAL.getmaskband(band)

"""
Return the status flags of the mask band associated with the band.

The GetMaskFlags() method returns an bitwise OR-ed set of status flags with the
following available definitions that may be extended in the future:

* `GMF_ALL_VALID` (`0x01`):    There are no invalid pixels, all mask values
will be 255. When used this will normally be the only flag set.
* `GMF_PER_DATASET` (`0x02`):  The mask band is shared between all bands on
the dataset.
- `GMF_ALPHA` (`0x04`):        The mask band is actually an alpha band and may
have values other than 0 and 255.
- `GMF_NODATA` (`0x08`):       Indicates the mask is actually being generated
from nodata values. (mutually exclusive of `GMF_ALPHA`)

The `GDALRasterBand` class includes a default implementation of `GetMaskBand()`
that returns one of four default implementations:

- If a corresponding .msk file exists it will be used for the mask band.
- If the dataset has a `NODATA_VALUES` metadata item, an instance of the new
`GDALNoDataValuesMaskBand` class will be returned. `GetMaskFlags()` will return
`GMF_NODATA | GMF_PER_DATASET`.
- If the band has a nodata value set, an instance of the new
`GDALNodataMaskRasterBand` class will be returned. `GetMaskFlags()` will return
`GMF_NODATA`.
- If there is no nodata value, but the dataset has an alpha band that seems to
apply to this band (specific rules yet to be determined) and that is of type
`GDT_Byte` then that alpha band will be returned, and the flags `GMF_PER_DATASET`
and `GMF_ALPHA` will be returned in the flags.
- If neither of the above apply, an instance of the new `GDALAllValidRasterBand`
class will be returned that has 255 values for all pixels. The null flags will
return `GMF_ALL_VALID`.

See also: http://trac.osgeo.org/gdal/wiki/rfc15_nodatabitmask

### Returns
a valid mask band.
"""
getmaskflags(band::RasterBand) = GDAL.getmaskflags(band)

"""
Adds a mask band to the current band.

The default implementation of the `CreateMaskBand()` method is implemented
based on similar rules to the `.ovr` handling implemented using the
`GDALDefaultOverviews` object. A `TIFF` file with the extension `.msk` will be
created with the same basename as the original file, and it will have as many
bands as the original image (or just one for `GMF_PER_DATASET`). The mask
images will be deflate compressed tiled images with the same block size as the
original image if possible.

If you got a mask band with a previous call to `GetMaskBand()`, it might be
invalidated by `CreateMaskBand()`. So you have to call `GetMaskBand()` again.

See also: http://trac.osgeo.org/gdal/wiki/rfc15_nodatabitmask

### Returns
`CE_None` on success or `CE_Failure` on an error.
"""
function createmaskband(band::RasterBand, nFlags::Integer)
    result = GDAL.createmaskband(band, nFlags)
    (result != GDAL.CE_None) && error("Failed to create mask band")
end


text:
"""
`@with(d, x)`

Replace all symbols in expression `x` with the respective column in `d`. In this context,
`_` refers to the whole table `d`. To use actual symbols, escape them with `^`, as in `^(:a)`.
Use `cols(c)` to refer to column `c` where `c` is a variable that evaluates to a symbol. `c` must be available in
the scope where the macro is called.

## Examples

```jldoctest with
julia> t = table((a = [1,2,3], b = ["x","y","z"]));

julia> @with t mean(:a)
2.0

julia> @with t mean(:a)*length(_)
6.0

julia> @with t join(:b)
"xyz"

julia> @with t @show ^(:a) != :a
:a != getfield(JuliaDBMeta.columns(t), :a) = true
true

julia> c = :a
:a

julia> @with t cols(c)
3-element Array{Int64,1}:
 1
 2
 3
```

Note that you can use this syntax to modify columns in place as well:

```jldoctest with
julia> @with t :b .= :b .* string.(:a)
3-element Array{String,1}:
 "x1"
 "y2"
 "z3"

julia> t
Table with 3 rows, 2 columns:
a  b
───────
1  "x1"
2  "y2"
3  "z3"
```
"""
macro with(args...)
    esc(with_helper(args...))
end

function with_helper(args...)
    d = gensym()
    func, _ = extract_anonymous_function(last(args), replace_column)
    Expr(:call, :(JuliaDBMeta._pipe), func, replace_keywords(args[1:end-1])...)
end

replace_column(d, x) = Expr(:call, :getfield, :(JuliaDBMeta.columns($d)), x)
replace_column(d) = d

ext: jl


text:
"""
A full drainage basin, with a series of rivers and of reservoirs.

TODO: how to handle topology?
TODO: unused for now!
"""
immutable DrainageBasin
  name::AbstractString
  variant::AbstractString

  reservoirs::Array{Reservoir, 1}

  function DrainageBasin(name::AbstractString, reservoirs::Array{Reservoir, 1})
    if length(reservoirs) < 1
      error("List of reservoirs does not even contain a reservoir.")
    end
    obj = new(name, "vanilla", reservoirs)

    # Check consistency of scenarios.
    countScenariosWithConsistency(obj)
    countTimeStepsWithConsistency(obj)

    return obj
  end

  function DrainageBasin(name::AbstractString, reservoir::Reservoir)
    obj = new(name, "vanilla", [reservoir])

    # Check consistency of scenarios.
    countScenariosWithConsistency(obj)
    countTimeStepsWithConsistency(obj)

    return obj
  end
end

### Generic getters.

getName(b::DrainageBasin) = b.name
getVariant(b::DrainageBasin) = b.variant
getReservoirs(b::DrainageBasin) = b.reservoirs

# Counting things.

"""
Counts the number of scenarios within this basin.
"""
countScenarios(b::DrainageBasin) = countScenarios(getReservoirs(b)[1])

"""
Counts the number of scenarios within this basin and ensures the whole data structure is consistent.
"""
function countScenariosWithConsistency(b::DrainageBasin)
  cnt = -1

  for reservoir in getReservoirs(b)
    if cnt < 0
      cnt = countScenariosWithConsistency(reservoir)
    elseif cnt != countScenariosWithConsistency(reservoir)
      error("Inconsistent data in basin " * getName(b) * ": different number of scenarios within the basin (while checking reservoir " * getName(reservoir) * ").")
    end
  end

  if cnt < 0
    error("Inconsistent data in basin " * getName(b) * ": no time series found in the rivers.")
  end

  return cnt
end

"""
Counts the number of time steps within this basin.
"""
countTimeSteps(b::DrainageBasin) = countTimeSteps(getReservoirs(b)[1])

"""
Counts the number of time steps within this basin and ensures the whole data structure is consistent.
"""
function countTimeStepsWithConsistency(b::DrainageBasin)
  cnt = -1

  for reservoir in getReservoirs(b)
    if cnt < 0
      cnt = countTimeStepsWithConsistency(reservoir)
    elseif cnt != countTimeStepsWithConsistency(reservoir)
      error("Inconsistent data in basin " * getName(b) * ": different number of time steps within the basin (while checking reservoir " * getName(reservoir) * ").")
    end
  end

  if cnt < 0
    error("Inconsistent data in basin " * getName(b) * ": no time series found in the rivers.")
  end

  return cnt
end

"""
Retrieves the periodicity of this basin data.
"""
getPeriod(b::DrainageBasin) = getPeriod(getReservoirs(b)[1])

"""
Retrieves the periodicity of this basin data and ensures the whole data structure is consistent.
"""
function getPeriodWithConsistency(b::DrainageBasin)
  period = :Unknown

  for reservoir in getReservoirs(b)
    if period == :Unknown
      period = getPeriod(reservoir)
    elseif period != getPeriod(reservoir)
      error("Inconsistent data in basin " * getName(b) * ": different period within the basin (while checking reservoir " * getName(reservoir) * ").")
    end
  end

  if period == :Unknown
    error("Inconsistent data in basin " * getName(b) * ": no time series found in the rivers.")
  end

  return period::Period
end


text:
"maybe of type `T` or nothing"
const Maybe{T} = Union{Nothing, T}

"""
    PartialDefault{F}

Type for non-constant default value, it depends on the value of another field that has default value.
"""
struct PartialDefault{F}
    lambda::F
    vars::Vector{Symbol}
    expr::Expr
end

(f::PartialDefault)(x) = f.lambda(x)

function Base.show(io::IO, x::PartialDefault)
    print(io, JLFunction(;head=:->, args=x.vars, body=x.expr))
end

"""
    create(::Type{T}; kwargs...) where T

Create an instance of option type `T` from `kwargs`. Similar
to the default keyword argument constructor, but one can use this to create
custom keyword argument constructor with extra custom keywords.
"""
function create(::Type{T}; kwargs...) where T
    error("$T is not an option type")
end


text:
begin # Create (Deterministic) Limit Order Generator
    MyOrderSubTypes = (Int64,Float32,Int64,Int64) # define types for Order Size, Price, Order IDs, Account IDs
    MyOrderType = Order{MyOrderSubTypes...}
    MyLOBType = OrderBook{MyOrderSubTypes...}
    using Base.Iterators: cycle, take, zip, flatten
    orderid_iter = Base.Iterators.countfrom(1)
    sign_iter = cycle([1,-1,-1,1,1,-1])
    side_iter = ( s>0 ? SELL_ORDER : BUY_ORDER for s in sign_iter )
    spread_iter = cycle([3 2 3 2 2 2 3 2 3 4 2 2 3 2 3 2 3 3 2 2 3 2 5 2 2 2 2 2 4 2 3 6 5 6 3 2 3 5 4]*1e-2)
    price_iter = ( Float32(100.0 + sgn*δ) for (δ,sgn) in zip(spread_iter,sign_iter) )
    size_iter = cycle([2,5,3,4,10,15,1,6,13,11,4,1,5])
    # zip them all together
    lmt_order_info_iter = zip(orderid_iter,price_iter,size_iter,side_iter)
end

begin # Create (Deterministic) Market Order Generator
    mkt_size_iter = cycle([10,20,30,15,25,5,7])
    mkt_side_iter = cycle([SELL_ORDER,BUY_ORDER,BUY_ORDER,SELL_ORDER,BUY_ORDER,SELL_ORDER])
    mkt_order_info_iter = zip(mkt_size_iter,mkt_side_iter)
end

@testset "Submit and Cancel 1" begin # Add and delete all orders, verify book is empty, verify account tracking
    ob = MyLOBType() #Initialize empty book
    order_info_lst = take(lmt_order_info_iter,50000)
    # Add a bunch of orders
    for (orderid, price, size, side) in order_info_lst
        submit_limit_order!(ob,orderid,side,price,size,10101)
    end
    @test length(ob.acct_map[10101]) == 50000 # Check account order tracking
    # Cancel them all
    for (orderid, price, size, side) in order_info_lst
        cancel_order!(ob,orderid,side,price)
    end
    # Check emptiness
    @test isempty(ob.bid_orders)
    @test isempty(ob.ask_orders)
    @test isempty(ob.acct_map[10101])
end

@testset "MO Liquidity Wipe" begin # Wipe out book completely, try MOs on empty book
    ob = MyLOBType() #Initialize empty book
    # Add a bunch of orders
    for (orderid, price, size, side) in Base.Iterators.take( lmt_order_info_iter, 50 )
        submit_limit_order!(ob,orderid,BUY_ORDER,price,size)
    end
    mo_matches, mo_ltt = submit_market_order!(ob,SELL_ORDER,100000)

    # Tests
    @test length( mo_matches ) == 50
    @test mo_ltt > 0
    @test isempty(submit_market_order!(ob,SELL_ORDER,10000)[1] )
    @test isempty(ob.bid_orders)
    @test isempty(ob.ask_orders)
end

@testset "Order match exact - bid" begin # Test correctness in order matching system / Stat calculation (:BID)
    ob = MyLOBType() #Initialize empty book
    # record order book info before
    order_lst_tmp = Base.Iterators.take( Base.Iterators.filter( x-> x[4]===BUY_ORDER, lmt_order_info_iter), 7 ) |> collect

    # Add a bunch of orders
    for (orderid, price, size, side) in order_lst_tmp
        submit_limit_order!(ob,orderid,side,price,size)
    end

    orders_before = Iterators.flatten(q.queue for (k,q) in ob.bid_orders.book) |> collect

    # record information from before
    expected_bid_volm_before = sum( x[3] for x in order_lst_tmp )
    expected_bid_n_orders_before =  length(order_lst_tmp)

    # execute MO
    mo_matches, mo_ltt = submit_market_order!(ob,SELL_ORDER,30)
    mo_match_sizes = [o.size for o in mo_matches]

    # record what is expected to be seen
    expected_bid_volm_after = expected_bid_volm_before - 30
    expected_bid_n_orders_after = expected_bid_n_orders_before - 5
    expected_best_bid_after = Float32(99.97)

    # record what is expected of MO result
    expected_mo_match_size = [5,15,6,1,2,1]

    # Compute realized values
    book_info_after = book_depth_info(ob,1000)
    realized_bid_volm_after  = sum(book_info_after[:BID][:volume])
    realized_bid_n_orders_after = sum(book_info_after[:BID][:orders])
    realized_best_bid_after = first(book_info_after[:BID][:price])

    # Check all expected vs realized values
    @test realized_bid_volm_after == expected_bid_volm_after
    @test realized_bid_n_orders_after == expected_bid_n_orders_after
    @test realized_best_bid_after == expected_best_bid_after
    @test mo_match_sizes == expected_mo_match_size
    @test mo_ltt == 0
end

@testset "Test MO, LO insert, LO cancel outputs" begin
    ob = MyLOBType() #Initialize empty book
    order_info_lst = take(lmt_order_info_iter,500)
    # Add a bunch of orders
    for (orderid, price, size, side) in order_info_lst
        submit_limit_order!(ob,orderid,side,price,size)
    end

    # Test that inserting LO returns correctly
    lmt_info = (10_000, BUY_ORDER, 99.97f0, 3,)
    lmt_obj, _, _ = submit_limit_order!(ob,lmt_info...)
    @test lmt_info[1:4] == (lmt_obj.orderid,lmt_obj.side,lmt_obj.price,lmt_obj.size)

    # Add something on top
    submit_limit_order!(ob, 10_001, BUY_ORDER, 100.02f0, 900,)

    # Test that cancelling present order returns correctly
    lmt_obj_cancel = cancel_order!(ob,lmt_obj)
    @test lmt_obj_cancel == lmt_obj

    # Test that missing order returns correctly
    lmt_obj_cancel_2 = cancel_order!(ob,lmt_obj_cancel)
    @test isnothing(lmt_obj_cancel_2)

    # Test that complete MO returns correctly
    mo_match_list, mo_ltt = submit_market_order!(ob,SELL_ORDER,100)
    @test typeof(mo_match_list) <: Vector{<:Order}
    @test mo_ltt == 0

    mo_match_list, mo_ltt = submit_market_order!(ob,BUY_ORDER,797+13)
    @test length(mo_match_list) == 128
    @test mo_ltt == 13

    mo_match_list, mo_ltt = submit_market_order!(ob,BUY_ORDER,13)
    @test isempty(mo_match_list)
    @test mo_ltt == 13

end

@testset "Test Account Tracking" begin
    ob = MyLOBType() #Initialize empty book

    # Add a bunch of orders
    for (orderid, price, size, side) in take(lmt_order_info_iter,100)
        submit_limit_order!(ob,orderid,side,price,size)
    end

    # Add order with an account ID
    acct_id = 1313
    order_id0 = 10001
    my_acct_orders = MyOrderType[]
    push!(my_acct_orders,submit_limit_order!(ob,order_id0,SELL_ORDER,100.03f0,50,acct_id)[1])
    push!(my_acct_orders,submit_limit_order!(ob,order_id0+1,BUY_ORDER,99.98f0,20,acct_id)[1])
    push!(my_acct_orders,submit_limit_order!(ob,order_id0+2,BUY_ORDER,99.97f0,30,acct_id)[1])

    # Throw some more nameless orders on top
    for (orderid, price, size, side) in take(lmt_order_info_iter,20)
        submit_limit_order!(ob,orderid,side,price,size)
    end

    # Get account list from book
    book_acct_list = collect(get_acct(ob,acct_id))
    @test (order_id0 .+ collect(0:2)) == [first(x) for x in book_acct_list] # Test correct ids
    @test my_acct_orders == [last(x) for x in book_acct_list] # Test correct orders
    @test isnothing(get_acct(ob,0))

    # Delete some orders and maintain checks
    to_canc = popat!(my_acct_orders,2)
    canc_order = cancel_order!(ob,to_canc)
    @test to_canc == canc_order
    book_acct_list = collect(get_acct(ob,acct_id))
    @test to_canc ∉ book_acct_list

end

# using BenchmarkTools
# # Add a bunch of orders

# ob = MyLOBType() #Initialize empty book

# order_info_lst = take(lmt_order_info_iter,Int64(10_000)) |> collect
# for (orderid, price, size, side) in order_info_lst
#     submit_limit_order!(ob,orderid,side,price,size)
# end

# (orderid, price, size, side), _ =  Iterators.peel(lmt_order_info_iter)
# @benchmark submit_limit_order!($ob,$orderid,$side,$price,$size,)
# @benchmark (submit_market_order!($ob,BUY_ORDER,1000);)

# @code_typed submit_limit_order!(ob,orderid,side,price,size)

# ob = MyLOBType() # initialize order book

# # fill book with random limit orders
# randspread() = ceil(-0.03*log(rand()),digits=2)
# for i=1:1000
#     submit_limit_order!(ob,2i,BUY_ORDER,99.0-randspread(),rand(1:25))
#     submit_limit_order!(ob,3i,SELL_ORDER,99.0+randspread(),rand(1:25))
# end

# @benchmark submit_limit_order!(ob,$2,$(rand([BUY_ORDER,SELL_ORDER])),$(99.0+rand([1,-1])*randspread()),$(rand(1:25)))


text:
export simloglik_posterior!, SharedPosterior

@GenGlobal g_SharedPosterior

struct SharedPosterior
    drill::SharedMatrix{Float64}
    royalty::SharedMatrix{Float64}
    produce::SharedMatrix{Float64}
    function SharedPosterior(d,r,p)
        size(d) == size(r) == size(p) ||
            throw(DimensionMismatch())
        return new(d,r,p)
    end
end

@getFieldFunction SharedPosterior drill royalty produce

function fill!(s::SharedPosterior, x)
    fill!(drill(s), x)
    fill!(royalty(s), x)
    fill!(produce(s), x)
    return nothing
end

function SharedPosterior(M,N; kwargs...)
    d = SharedMatrix{Float64}(M,N; kwargs...)
    r = SharedMatrix{Float64}(M,N; kwargs...)
    p = SharedMatrix{Float64}(M,N; kwargs...)
    s = SharedPosterior(d,r,p)
    fill!(s, 0)
    return s
end

SharedPosterior(d::DataFull, M; kwargs...) = SharedPosterior(M, num_i(d); kwargs...)

function view(s::SharedPosterior, i)
    d = view(drill(s), :, i)
    r = view(royalty(s), :, i)
    p = view(produce(s), :, i)
    return (d,r,p)
end

# SML for unit i
function simloglik_posterior!(grptup::NTuple{3,ObservationGroup}, thetas,
    sim::SimulationDrawsVector, posteriors::NTuple{3,AbstractVector};
    kwargs...
) where {N}

    qm = _qm(sim)

    for (k,(theta, grp, post)) in enumerate(zip(thetas, grptup, posteriors))
        fill!(qm, 0)
        simloglik!(theta, grp, theta, sim, false; kwargs...)
        post .= qm
    end
    return nothing
end

@noinline function simloglik_posterior!(i, theta, sim, posteriors, data::DataFull; kwargs...)

    grptup = getindex.(data, i)
    simi = view(sim, i)
    ptup = view(posteriors, i)

    θρ = theta_ρ(data, theta)
    update!(simi, θρ)

    thetasvw = split_thetas(data, theta)

    simloglik_posterior!(grptup, thetasvw, simi, ptup; kwargs...)
    return nothing
end

function simloglik_posterior!(i, theta; kwargs...)
    data = get_g_BaseDataSetofSets()
    sim = get_g_SimulationDrawsMatrix()# ::SimulationDraws{Float64,2,Matrix{Float64}}
    posteriors = get_g_SharedPosterior()
    simloglik_posterior!(i, theta, sim, posteriors::SharedPosterior, data; kwargs...)
    return nothing
end


text:
export PrimitiveMeasure

abstract type PrimitiveMeasure <: AbstractMeasure end

basemeasure(μ::PrimitiveMeasure) = μ

logdensity(μ::PrimitiveMeasure, x) = 0.0

logdensity(μ::M, ν::M, x) where {M<:PrimitiveMeasure} = 0.0


text:
using .Threads

const N = nthreads()
const t1 = zeros(Int, N)
const s1 = zeros(Int, N)

@threads for i in 1:N
    t1[i] = threadid()
end

@sync for i in 1:N
    Threads.@spawn s1[i] = threadid()
end

for j in 1:5
    @threads for i in 1:N
        t1[i] = threadid()
    end
    println(t1, " - ", Set(t1))
end

for j in 1:5
    @sync for i in 1:N
        Threads.@spawn s1[i] = threadid()
    end
    println(s1, " - ", Set(s1))
    s1[:] .= 0
end


text:
using MyExample
using DelimitedFiles
using LinearAlgebra: Eigen, eigen

pre = "Data/0.5_(15, 15, 15)_(1, 1, 1)"

Y = SparseTensor(pre*"_data.txt")
K1 = readdlm(pre*"_k1.txt")
K2 = readdlm(pre*"_k2.txt")
K3 = readdlm(pre*"_k3.txt")
kernels = [K1,K2,K3]
kernelseigen = eigen.(kernels)

model = NSKRRegressor(size(Y), [10^(-2), 10^(-2), 10^(-2)], "zeros")
fit(model, kernelseigen, Y)
predict(model, kernels)

predict_LOO(model, kernelseigen, Y, (0,))
predict_LOO(model, kernelseigen, Y, (1,2,))

# LOO without indication same result as just predict (which uses the hatmatrixes,
# bit like fit_predict, but without setting the A_ coefficients (faster)).
predict(model, kernelseigen, Y) == predict_LOO(model, kernelseigen, Y, ())

max_stars_repo_name: padwulf/MyExample

max_issues_repo_name: padwulf/MyExample

max_forks_repo_name: padwulf/MyExample


text:
## BEFORE RUNNING THIS SCRIPT, MAKE SURE ros is in the environment

if parsed_args["imshow"]
  using ImageView
  using Gtk.ShortNames
end

using TimeZones
using LinearAlgebra
using ImageMagick
using FileIO
using ImageFeatures, ImageDraw, Images, CoordinateTransformations, Rotations
using DataStructures
using ColorTypes, FixedPointNumbers
using FreeTypeAbstraction # for drawTagID!
using AprilTags
using BSON
using JSON2
using Dates
using DataInterpolations

# using CuArrays
using Flux
using RoME
using Caesar

# source /opt/ros/melodic/setup.bash

## Project path variables

projdir = @__DIR__

## work with multiple processes

using Distributed

# populate remote machines on cluster later (was need for Python workaround and left)
prcs115 = Int[]

# ## Process on this machine but not remote machines
# using Caesar

## Prepare python version

using Pkg
@everywhere using Pkg

# ENV["PYTHON"] = "/usr/bin/python3.6"
@everywhere begin
  ENV["PYTHON"] = "/usr/bin/python"
  Pkg.build("PyCall")
end

using PyCall
@everywhere using PyCall

# for prc in prcs115
#   @info "Load rospy at worker $prc"
#   fut = Distributed.@spawnat prc begin
#     # import subprocess
#     # rc = subprocess.call("/opt/ros/melodic/setup.bash",shell=True)
#     # """
#     py"""
#     import sys
#     sys.path.insert(0, "/opt/ros/melodic/lib/python2.7/dist-packages")
#     """
#     py"""
#     import rospy
#     """
#   end
#   @show fetch(fut)
# end

## prepare ROS integration

# this will trigger inclusion of RobotOS functionality in Caesar too.
using RobotOS

# standard types
# @rosimport sensor_msgs.msg: CompressedImage
# @rosimport geometry_msgs.msg: PoseStamped
#
# rostypegen()

## start common processes on local machine

if nprocs() < parsed_args["localprocs"]
  addprocs(parsed_args["localprocs"])
end

## Start processes on remote machines

if 0 < parsed_args["remoteprocs"]
  machines = []
  allmc = split(ENV[parsed_args["remoteserver"]], ';')
  for mc in allmc
    push!(machines, (mc,parsed_args["remoteprocs"]))
  end
  prcs115 = addprocs(machines)
end

## Load solver libraries everywhere

using Flux
using RoME
using DataInterpolations
using FileIO, ImageMagick
using BSON
@everywhere using BSON
@everywhere using Flux
@everywhere using RoME
@everywhere using IncrementalInference, DistributedFactorGraphs, TransformUtils, DataInterpolations, Rotations
@everywhere using TimeZones
@everywhere using FileIO, ImageMagick

## Constant parameters

# using JSON2
const ImgType = Array{ColorTypes.RGB{FixedPointNumbers.Normed{UInt8,8}},2}

# Consume rosbag with car data
runfile = parsed_args["folder_name"] # "labrun7"
bagfile = joinpath(ENV["HOME"],"data/racecar/$runfile.bag")

leftimgtopic = "/zed/left/image_rect_color"
# rightimgtopic = "/zed/right/image_rect_color"
zedodomtopic = "/zed/odom"
joysticktopic = "/joy"

# from bagfile -- TODO -- make this simpler
fx = 341.4563903808594
fy = 341.4563903808594
cx = 329.19091796875
cy = 196.3658447265625

K = [-fx 0  cx;
      0 fy cy]

# guessing
# 758016/3 = 252672
# (c, a / c) = (672, 376.0)

WEIRDOFFSET = Dict{Symbol, Int}() #:cmdVal => -11345)
detector = AprilTagDetector()

## load visualization if requested

if parsed_args["vis2d"]

using Cairo
using Gadfly
using RoMEPlotting
Gadfly.set_default_plot_size(35cm, 20cm)

end

#


text:
module ArrayInterface

using Requires
using LinearAlgebra
using SparseArrays

using Base: @propagate_inbounds, tail, OneTo, LogicalIndex, Slice

Base.@pure __parameterless_type(T) = Base.typename(T).wrapper
parameterless_type(x) = parameterless_type(typeof(x))
parameterless_type(x::Type) = __parameterless_type(x)

"""
    parent_type(x)

Returns the parent array that `x` wraps.
"""
parent_type(x) = parent_type(typeof(x))
parent_type(::Type{<:SubArray{T,N,P}}) where {T,N,P} = P
parent_type(::Type{<:Base.ReshapedArray{T,N,P}}) where {T,N,P} = P
parent_type(::Type{Adjoint{T,S}}) where {T,S} = S
parent_type(::Type{Transpose{T,S}}) where {T,S} = S
parent_type(::Type{Symmetric{T,S}}) where {T,S} = S
parent_type(::Type{<:LinearAlgebra.AbstractTriangular{T,S}}) where {T,S} = S
parent_type(::Type{<:PermutedDimsArray{T,N,I1,I2,A}}) where {T,N,I1,I2,A} = A
parent_type(::Type{Slice{T}}) where {T} = T
parent_type(::Type{T}) where {T} = T

"""
    known_length(::Type{T})

If `length` of an instance of type `T` is known at compile time, return it.
Otherwise, return `nothing`.
"""
known_length(x) = known_length(typeof(x))
known_length(::Type{<:NamedTuple{L}}) where {L} = length(L)
known_length(::Type{T}) where {T<:Slice} = known_length(parent_type(T))
known_length(::Type{<:Tuple{Vararg{Any,N}}}) where {N} = N
known_length(::Type{<:Number}) = 1
function known_length(::Type{T}) where {T}
    if parent_type(T) <: T
        return nothing
    else
        return known_length(parent_type(T))
    end
end
@inline function known_length(::Type{<:SubArray{T,N,P,I}}) where {T,N,P,I}
    return _known_length(ntuple(i -> known_length(I.parameters[i]), Val(N)))
end
_known_length(x::Tuple{Vararg{<:Union{Int,Nothing}}}) = nothing
_known_length(x::Tuple{Vararg{Int}}) = prod(x)

"""
    can_change_size(::Type{T}) -> Bool

Returns `true` if the Base.size of `T` can change, in which case operations
such as `pop!` and `popfirst!` are available for collections of type `T`.
"""
can_change_size(x) = can_change_size(typeof(x))
can_change_size(::Type{T}) where {T} = false
can_change_size(::Type{<:Vector}) = true
can_change_size(::Type{<:AbstractDict}) = true
can_change_size(::Type{<:Base.ImmutableDict}) = false

function ismutable end

"""
    ismutable(x::DataType)

Query whether a type is mutable or not, see
https://github.com/JuliaDiffEq/RecursiveArrayTools.jl/issues/19.
"""
ismutable(x) = ismutable(typeof(x))

function ismutable(::Type{T}) where {T<:AbstractArray}
    if parent_type(T) <: T
        return true
    else
        return ismutable(parent_type(T))
    end
end
ismutable(::Type{<:AbstractRange}) = false
ismutable(::Type{<:AbstractDict}) = true
ismutable(::Type{<:Base.ImmutableDict}) = false
function ismutable(::Type{T}) where {T}
    if parent_type(T) <: T
        return T.mutable
    else
        return ismutable(parent_type(T))
    end
end

# Piracy
function Base.setindex(x::AbstractArray,v,i...)
  _x = copy(x)
  _x[i...] = v
  _x
end

function Base.setindex(x::AbstractVector,v,i::Int)
  n = length(x)
  x .* (i .!== 1:n) .+ v .* (i .== 1:n)
end

function Base.setindex(x::AbstractMatrix,v,i::Int,j::Int)
  n,m = Base.size(x)
  x .* (i .!== 1:n) .* (j  .!== i:m)' .+ v .* (i .== 1:n) .* (j  .== i:m)'
end

"""
    can_setindex(x::DataType)

Query whether a type can use `setindex!`.
"""
can_setindex(x) = true
can_setindex(x::AbstractArray) = can_setindex(typeof(x))
can_setindex(::Type{<:AbstractRange}) = false

"""
    aos_to_soa(x)

Converts an array of structs formulation to a struct of array.
"""
aos_to_soa(x) = x

"""
    fast_scalar_indexing(x)

Query whether an array type has fast scalar indexing.
"""
fast_scalar_indexing(x) = true
fast_scalar_indexing(x::AbstractArray) = fast_scalar_indexing(typeof(x))
fast_scalar_indexing(::Type{<:LinearAlgebra.AbstractQ}) = false
fast_scalar_indexing(::Type{<:LinearAlgebra.LQPackedQ}) = false

"""
    allowed_getindex(x,i...)

A scalar `getindex` which is always allowed.
"""
allowed_getindex(x,i...) = x[i...]

"""
    allowed_setindex!(x,v,i...)

A scalar `setindex!` which is always allowed.
"""
allowed_setindex!(x,v,i...) = Base.setindex!(x,v,i...)

"""
    isstructured(x::DataType)

Query whether a type is a representation of a structured matrix.
"""
isstructured(x) = false
isstructured(x::AbstractArray) = isstructured(typeof(x))
isstructured(::Symmetric) = true
isstructured(::Hermitian) = true
isstructured(::UpperTriangular) = true
isstructured(::LowerTriangular) = true
isstructured(::Tridiagonal) = true
isstructured(::SymTridiagonal) = true
isstructured(::Bidiagonal) = true
isstructured(::Diagonal) = true

"""
    has_sparsestruct(x::AbstractArray)

Determine whether `findstructralnz` accepts the parameter `x`.
"""
has_sparsestruct(x) = false
has_sparsestruct(x::AbstractArray) = has_sparsestruct(typeof(x))
has_sparsestruct(x::Type{<:AbstractArray}) = false
has_sparsestruct(x::Type{<:SparseMatrixCSC}) = true
has_sparsestruct(x::Type{<:Diagonal}) = true
has_sparsestruct(x::Type{<:Bidiagonal}) = true
has_sparsestruct(x::Type{<:Tridiagonal}) = true
has_sparsestruct(x::Type{<:SymTridiagonal}) = true

"""
    issingular(A::AbstractMatrix)

Determine whether a given abstract matrix is singular.
"""
issingular(A::AbstractMatrix) = issingular(Matrix(A))
issingular(A::AbstractSparseMatrix) = !issuccess(lu(A, check=false))
issingular(A::Matrix) = !issuccess(lu(A, check=false))
issingular(A::UniformScaling) = A.λ == 0
issingular(A::Diagonal) = any(iszero,A.diag)
issingular(B::Bidiagonal) = any(iszero, A.dv)
issingular(S::SymTridiagonal) = diaganyzero(iszero, ldlt(S).data)
issingular(T::Tridiagonal) = !issuccess(lu(A, check=false))
issingular(A::Union{Hermitian,Symmetric}) = diaganyzero(bunchkaufman(A, check=false).LD)
issingular(A::Union{LowerTriangular,UpperTriangular}) = diaganyzero(A.data)
issingular(A::Union{UnitLowerTriangular,UnitUpperTriangular}) = false
issingular(A::Union{Adjoint,Transpose}) = issingular(parent(A))
diaganyzero(A) = any(iszero, view(A, diagind(A)))

"""
    findstructralnz(x::AbstractArray)

Return: (I,J) #indexable objects
Find sparsity pattern of special matrices, the same as the first two elements of findnz(::SparseMatrixCSC).
"""
function findstructralnz(x::Diagonal)
  n = Base.size(x,1)
  (1:n,1:n)
end

abstract type MatrixIndex end

struct BidiagonalIndex <: MatrixIndex
  count::Int
  isup::Bool
end

struct TridiagonalIndex <: MatrixIndex
  count::Int#count==nsize+nsize-1+nsize-1
  nsize::Int
  isrow::Bool
end

struct BandedMatrixIndex <: MatrixIndex
  count::Int
  rowsize::Int
  colsize::Int
  bandinds::Array{Int,1}
  bandsizes::Array{Int,1}
  isrow::Bool
end

function _bandsize(bandind,rowsize,colsize)
  -(rowsize-1) <= bandind <= colsize-1 || throw(ErrorException("Invalid Bandind"))
  if (bandind*(colsize-rowsize)>0) & (abs(bandind)<=abs(colsize-rowsize))
    return min(rowsize,colsize)
  elseif bandind*(colsize-rowsize)<=0
    return min(rowsize,colsize)-abs(bandind)
  else
    return min(rowsize,colsize)-abs(bandind)+abs(colsize-rowsize)
  end
end

function BandedMatrixIndex(rowsize,colsize,lowerbandwidth,upperbandwidth,isrow)
  upperbandwidth>-lowerbandwidth || throw(ErrorException("Invalid Bandwidths"))
  bandinds=upperbandwidth:-1:-lowerbandwidth
  bandsizes=[_bandsize(band,rowsize,colsize) for band in bandinds]
  BandedMatrixIndex(sum(bandsizes),rowsize,colsize,bandinds,bandsizes,isrow)
end

struct BlockBandedMatrixIndex <: MatrixIndex
  count::Int
  refinds::Array{Int,1}
  refcoords::Array{Int,1}#storing col or row inds at ref points
  isrow::Bool
end

function BlockBandedMatrixIndex(nrowblock,ncolblock,rowsizes,colsizes,l,u)
  blockrowind=BandedMatrixIndex(nrowblock,ncolblock,l,u,true)
  blockcolind=BandedMatrixIndex(nrowblock,ncolblock,l,u,false)
  sortedinds=sort([(blockrowind[i],blockcolind[i]) for i in 1:length(blockrowind)],by=x->x[1])
  sort!(sortedinds,by=x->x[2],alg=InsertionSort)#stable sort keeps the second index in order
  refinds=Array{Int,1}()
  refrowcoords=Array{Int,1}()
  refcolcoords=Array{Int,1}()
  rowheights=cumsum(pushfirst!(copy(rowsizes),1))
  blockheight=0
  blockrow=1
  blockcol=1
  currenti=1
  lastrowind=sortedinds[1][1]-1
  lastcolind=sortedinds[1][2]
  for ind in sortedinds
    rowind,colind=ind
    if colind==lastcolind
      if rowind>lastrowind
        blockheight+=rowsizes[rowind]
      end
    else
      for j in blockcol:blockcol+colsizes[lastcolind]-1
        push!(refinds,currenti)
        push!(refrowcoords,blockrow)
        push!(refcolcoords,j)
        currenti+=blockheight
      end
      blockcol+=colsizes[lastcolind]
      blockrow=rowheights[rowind]
      blockheight=rowsizes[rowind]
    end
    lastcolind=colind
    lastrowind=rowind
  end
  for j in blockcol:blockcol+colsizes[lastcolind]-1
    push!(refinds,currenti)
    push!(refrowcoords,blockrow)
    push!(refcolcoords,j)
    currenti+=blockheight
  end
  push!(refinds,currenti)#guard
  push!(refrowcoords,-1)
  push!(refcolcoords,-1)
  rowindobj=BlockBandedMatrixIndex(currenti-1,refinds,refrowcoords,true)
  colindobj=BlockBandedMatrixIndex(currenti-1,refinds,refcolcoords,false)
  rowindobj,colindobj
end

struct BandedBlockBandedMatrixIndex <: MatrixIndex
  count::Int
  refinds::Array{Int,1}
  refcoords::Array{Int,1}#storing col or row inds at ref points
  reflocalinds::Array{BandedMatrixIndex,1}
  isrow::Bool
end

function BandedBlockBandedMatrixIndex(nrowblock,ncolblock,rowsizes,colsizes,l,u,lambda,mu)
  blockrowind=BandedMatrixIndex(nrowblock,ncolblock,l,u,true)
  blockcolind=BandedMatrixIndex(nrowblock,ncolblock,l,u,false)
  sortedinds=sort([(blockrowind[i],blockcolind[i]) for i in 1:length(blockrowind)],by=x->x[1])
  sort!(sortedinds,by=x->x[2],alg=InsertionSort)#stable sort keeps the second index in order
  rowheights=cumsum(pushfirst!(copy(rowsizes),1))
  colwidths=cumsum(pushfirst!(copy(colsizes),1))
  currenti=1
  refinds=Array{Int,1}()
  refrowcoords=Array{Int,1}()
  refcolcoords=Array{Int,1}()
  reflocalrowinds=Array{BandedMatrixIndex,1}()
  reflocalcolinds=Array{BandedMatrixIndex,1}()
  for ind in sortedinds
    rowind,colind=ind
    localrowind=BandedMatrixIndex(rowsizes[rowind],colsizes[colind],lambda,mu,true)
    localcolind=BandedMatrixIndex(rowsizes[rowind],colsizes[colind],lambda,mu,false)
    push!(refinds,currenti)
    push!(refrowcoords,rowheights[rowind])
    push!(refcolcoords,colwidths[colind])
    push!(reflocalrowinds,localrowind)
    push!(reflocalcolinds,localcolind)
    currenti+=localrowind.count
  end
  push!(refinds,currenti)
  push!(refrowcoords,-1)
  push!(refcolcoords,-1)
  rowindobj=BandedBlockBandedMatrixIndex(currenti-1,refinds,refrowcoords,reflocalrowinds,true)
  colindobj=BandedBlockBandedMatrixIndex(currenti-1,refinds,refcolcoords,reflocalcolinds,false)
  rowindobj,colindobj
end

Base.firstindex(ind::MatrixIndex)=1
Base.lastindex(ind::MatrixIndex)=ind.count
Base.length(ind::MatrixIndex)=ind.count
function Base.getindex(ind::BidiagonalIndex,i::Int)
  1 <= i <= ind.count || throw(BoundsError(ind, i))
  if ind.isup
    ii=i+1
  else
    ii=i+1+1
  end
  convert(Int,floor(ii/2))
end

function Base.getindex(ind::TridiagonalIndex,i::Int)
  1 <= i <= ind.count || throw(BoundsError(ind, i))
  offsetu= ind.isrow ? 0 : 1
  offsetl= ind.isrow ? 1 : 0
  if 1 <= i <= ind.nsize
    return i
  elseif ind.nsize < i <= ind.nsize+ind.nsize-1
    return i-ind.nsize+offsetu
  else
    return i-(ind.nsize+ind.nsize-1)+offsetl
  end
end

function Base.getindex(ind::BandedMatrixIndex,i::Int)
  1 <= i <= ind.count || throw(BoundsError(ind, i))
  _i=i
  p=1
  while _i-ind.bandsizes[p]>0
    _i-=ind.bandsizes[p]
    p+=1
  end
  bandind=ind.bandinds[p]
  startfromone=!xor(ind.isrow,(bandind>0))
  if startfromone
    return _i
  else
    return _i+abs(bandind)
  end
end

function Base.getindex(ind::BlockBandedMatrixIndex,i::Int)
  1 <= i <= ind.count || throw(BoundsError(ind, i))
  p=1
  while i-ind.refinds[p]>=0
    p+=1
  end
  p-=1
  _i=i-ind.refinds[p]
  if ind.isrow
    return ind.refcoords[p]+_i
  else
    return ind.refcoords[p]
  end
end

function Base.getindex(ind::BandedBlockBandedMatrixIndex,i::Int)
  1 <= i <= ind.count || throw(BoundsError(ind, i))
  p=1
  while i-ind.refinds[p]>=0
    p+=1
  end
  p-=1
  _i=i-ind.refinds[p]+1
  ind.reflocalinds[p][_i]+ind.refcoords[p]-1
end

function findstructralnz(x::Bidiagonal)
  n= Base.size(x,1)
  isup= x.uplo=='U' ? true : false
  rowind=BidiagonalIndex(n+n-1,isup)
  colind=BidiagonalIndex(n+n-1,!isup)
  (rowind,colind)
end

function findstructralnz(x::Union{Tridiagonal,SymTridiagonal})
  n= Base.size(x,1)
  rowind=TridiagonalIndex(n+n-1+n-1,n,true)
  colind=TridiagonalIndex(n+n-1+n-1,n,false)
  (rowind,colind)
end

function findstructralnz(x::SparseMatrixCSC)
  rowind,colind,_=findnz(x)
  (rowind,colind)
end

abstract type ColoringAlgorithm end

"""
    fast_matrix_colors(A)

Query whether a matrix has a fast algorithm for getting the structural
colors of the matrix.
"""
fast_matrix_colors(A) = false
fast_matrix_colors(A::AbstractArray) = fast_matrix_colors(typeof(A))
fast_matrix_colors(A::Type{<:Union{Diagonal,Bidiagonal,Tridiagonal,SymTridiagonal}}) = true

"""
    matrix_colors(A::Union{Array,UpperTriangular,LowerTriangular})

The color vector for dense matrix and triangular matrix is simply
`[1,2,3,..., Base.size(A,2)]`.
"""
function matrix_colors(A::Union{Array,UpperTriangular,LowerTriangular})
    eachindex(1:Base.size(A,2)) # Vector Base.size matches number of rows
end

function _cycle(repetend,len)
    repeat(repetend,div(len,length(repetend))+1)[1:len]
end

function matrix_colors(A::Diagonal)
    fill(1, Base.size(A,2))
end

function matrix_colors(A::Bidiagonal)
    _cycle(1:2, Base.size(A,2))
end

function matrix_colors(A::Union{Tridiagonal,SymTridiagonal})
    _cycle(1:3, Base.size(A,2))
end

"""
  lu_instance(A) -> lu_factorization_instance

Returns an instance of the LU factorization object with the correct type
cheaply.
"""
function lu_instance(A::Matrix{T}) where T
  noUnitT = typeof(zero(T))
  luT = LinearAlgebra.lutype(noUnitT)
  ipiv = Vector{LinearAlgebra.BlasInt}(undef, 0)
  info = zero(LinearAlgebra.BlasInt)
  return LU{luT}(similar(A, 0, 0), ipiv, info)
end

"""
  lu_instance(a::Number) -> a

Returns the number.
"""
lu_instance(a::Number) = a

"""
    lu_instance(a::Any) -> lu(a, check=false)

Returns the number.
"""
lu_instance(a::Any) = lu(a, check=false)

"""
safevec(v)

It is a form of `vec` which is safe for all values in vector spaces, i.e., if it
is already a vector, like an AbstractVector or Number, it will return said
AbstractVector or Number.
"""
safevec(v) = vec(v)
safevec(v::Number) = v
safevec(v::AbstractVector) = v

"""
zeromatrix(u::AbstractVector)

Creates the zero'd matrix version of `u`. Note that this is unique because
`similar(u,length(u),length(u))` returns a mutable type, so it is not type-matching,
while `fill(zero(eltype(u)),length(u),length(u))` doesn't match the array type,
i.e., you'll get a CPU array from a GPU array. The generic fallback is
`u .* u' .* false`, which works on a surprising number of types, but can be broken
with weird (recursive) broadcast overloads. For higher-order tensors, this
returns the matrix linear operator type which acts on the `vec` of the array.
"""
function zeromatrix(u)
  x = safevec(u)
  x .* x' .* false
end

"""
restructure(x,y)

Restructures the object `y` into a shape of `x`, keeping its values intact. For
simple objects like an `Array`, this simply amounts to a reshape. However, for
more complex objects such as an `ArrayPartition`, not all of the structural
information is adequately contained in the type for standard tools to work. In
these cases, `restructure` gives a way to convert for example an `Array` into
a matching `ArrayPartition`.
"""
function restructure(x,y)
  out = similar(x,eltype(y))
  vec(out) .= vec(y)
end

function restructure(x::Array,y)
  reshape(convert(Array,y), Base.size(x)...)
end

abstract type AbstractDevice end
abstract type AbstractCPU <: AbstractDevice end
struct CPUPointer <: AbstractCPU end
struct CheckParent end
struct CPUIndex <: AbstractCPU end
struct GPU <: AbstractDevice end
"""
device(::Type{T})

Indicates the most efficient way to access elements from the collection in low-level code.
For `GPUArrays`, will return `ArrayInterface.GPU()`.
For `AbstractArray` supporting a `pointer` method, returns `ArrayInterface.CPUPointer()`.
For other `AbstractArray`s and `Tuple`s, returns `ArrayInterface.CPUIndex()`.
Otherwise, returns `nothing`.
"""
device(A) = device(typeof(A))
device(::Type) = nothing
device(::Type{<:Tuple}) = CPUIndex()
# Relies on overloading for GPUArrays that have subtyped `StridedArray`.
device(::Type{<:StridedArray}) = CPUPointer()
function device(::Type{T}) where {T <: AbstractArray}
    P = parent_type(T)
    T === P ? CPUIndex() : device(P)
end

"""
defines_strides(::Type{T}) -> Bool

Is strides(::T) defined?
"""
defines_strides(::Type) = false
defines_strides(x) = defines_strides(typeof(x))
defines_strides(::Type{<:StridedArray}) = true
defines_strides(::Type{A}) where {A <: Union{<:Transpose,<:Adjoint,<:SubArray,<:PermutedDimsArray}} = defines_strides(parent_type(A))

"""
can_avx(f)

Returns `true` if the function `f` is guaranteed to be compatible with
`LoopVectorization.@avx` for supported element and array types. While a return
value of `false` does not indicate the function isn't supported, this allows a
library to conservatively apply `@avx` only when it is known to be safe to do so.

```julia
function mymap!(f, y, args...)
    if can_avx(f)
        @avx @. y = f(args...)
    else
        @. y = f(args...)
    end
end
```
"""
can_avx(::Any) = false

"""
    insert(collection, index, item)

Returns a new instance of `collection` with `item` inserted into at the given `index`.
"""
Base.@propagate_inbounds function insert(collection, index, item)
  @boundscheck checkbounds(collection, index)
  ret = similar(collection, length(collection) + 1)
  @inbounds for i in firstindex(ret):(index - 1)
      ret[i] = collection[i]
  end
  @inbounds ret[index] = item
  @inbounds for i in (index + 1):lastindex(ret)
      ret[i] = collection[i - 1]
  end
  return ret
end

function insert(x::Tuple{Vararg{Any,N}}, index::Integer, item) where {N}
  @boundscheck if !checkindex(Bool, StaticInt{1}():StaticInt{N}(), index)
    throw(BoundsError(x, index))
  end
  return unsafe_insert(x, Int(index), item)
end

@inline function unsafe_insert(x::Tuple, i::Int, item)
  if i === 1
    return (item, x...)
  else
    return (first(x), unsafe_insert(tail(x), i - 1, item)...)
  end
end

"""
    deleteat(collection, index)

Returns a new instance of `collection` with the item at the given `index` removed.
"""
@propagate_inbounds function deleteat(collection::AbstractVector, index)
  @boundscheck if !checkindex(Bool, eachindex(collection), index)
    throw(BoundsError(collection, index))
  end
  return unsafe_deleteat(collection, index)
end
@propagate_inbounds function deleteat(collection::Tuple{Vararg{Any,N}}, index) where {N}
  @boundscheck if !checkindex(Bool, StaticInt{1}():StaticInt{N}(), index)
    throw(BoundsError(collection, index))
  end
  return unsafe_deleteat(collection, index)
end

function unsafe_deleteat(src::AbstractVector, index::Integer)
  dst = similar(src, length(src) - 1)
  @inbounds for i in indices(dst)
    if i < index
      dst[i] = src[i]
    else
      dst[i] = src[i + 1]
    end
  end
  return dst
end

@inline function unsafe_deleteat(src::AbstractVector, inds::AbstractVector)
  dst = similar(src, length(src) - length(inds))
  dst_index = firstindex(dst)
  @inbounds for src_index in indices(src)
    if !in(src_index, inds)
      dst[dst_index] = src[src_index]
      dst_index += one(dst_index)
    end
  end
  return dst
end

@inline function unsafe_deleteat(src::Tuple, inds::AbstractVector)
  dst = Vector{eltype(src)}(undef, length(src) - length(inds))
  dst_index = firstindex(dst)
  @inbounds for src_index in OneTo(length(src))
    if !in(src_index, inds)
      dst[dst_index] = src[src_index]
      dst_index += one(dst_index)
    end
  end
  return Tuple(dst)
end

@inline unsafe_deleteat(x::Tuple{T}, i::Integer) where {T} = ()
@inline unsafe_deleteat(x::Tuple{T1,T2}, i::Integer) where {T1,T2} = isone(i) ? (x[2],) : (x[1],)
@inline function unsafe_deleteat(x::Tuple, i::Integer)
  if i === one(i)
    return tail(x)
  elseif i == length(x)
    return Base.front(x)
  else
    return (first(x), unsafe_deleteat(tail(x), i - one(i))...)
  end
end

function __init__()

  @require SuiteSparse="4607b0f0-06f3-5cda-b6b1-a6196a1729e9" begin
    lu_instance(jac_prototype::SparseMatrixCSC) = SuiteSparse.UMFPACK.UmfpackLU(Ptr{Cvoid}(), Ptr{Cvoid}(), 1, 1,
                                                                                jac_prototype.colptr[1:1],
                                                                                jac_prototype.rowval[1:1],
                                                                                jac_prototype.nzval[1:1],
                                                                                0)
  end

  @require StaticArrays="90137ffa-7385-5640-81b9-e52037218182" begin
    ismutable(::Type{<:StaticArrays.StaticArray}) = false
    can_setindex(::Type{<:StaticArrays.StaticArray}) = false
    ismutable(::Type{<:StaticArrays.MArray}) = true
    ismutable(::Type{<:StaticArrays.SizedArray}) = true

    function lu_instance(_A::StaticArrays.StaticMatrix{N,N}) where {N}
      A = StaticArrays.SArray(_A)
      L = LowerTriangular(A)
      U = UpperTriangular(A)
      p = StaticArrays.SVector{N,Int}(1:N)
      return StaticArrays.LU(L, U, p)
    end

    function restructure(x::StaticArrays.SArray,y::StaticArrays.SArray)
      reshape(y,StaticArrays.Size(x))
    end

    function restructure(x::StaticArrays.SArray{S},y) where S
      StaticArrays.SArray{S}(y)
    end

    known_first(::Type{<:StaticArrays.SOneTo}) = 1
    known_last(::Type{StaticArrays.SOneTo{N}}) where {N} = N
    known_length(::Type{StaticArrays.SOneTo{N}}) where {N} = N

    device(::Type{<:StaticArrays.MArray}) = CPUPointer()
    contiguous_axis(::Type{<:StaticArrays.StaticArray}) = Contiguous{1}()
    contiguous_batch_size(::Type{<:StaticArrays.StaticArray}) = ContiguousBatch{0}()
    stride_rank(::Type{T}) where {N, T <: StaticArrays.StaticArray{<:Any,<:Any,N}} = StrideRank{ntuple(identity, Val{N}())}()
    dense_dims(::Type{<:StaticArrays.StaticArray{S,T,N}}) where {S,T,N} = DenseDims{ntuple(_ -> true, Val(N))}()
    defines_strides(::Type{<:StaticArrays.MArray}) = true
    @generated function size(A::StaticArrays.StaticArray{S}) where {S}
        t = Expr(:tuple); Sp = S.parameters
        for n in 1:length(Sp)
            push!(t.args, Expr(:call, Expr(:curly, :StaticInt, Sp[n])))
        end
        t
    end
    @generated function strides(A::StaticArrays.StaticArray{S}) where {S}
        t = Expr(:tuple, Expr(:call, Expr(:curly, :StaticInt, 1))); Sp = S.parameters; x = 1
        for n in 1:length(Sp)-1
            push!(t.args, Expr(:call, Expr(:curly, :StaticInt, (x *= Sp[n]))))
        end
        t
    end
    @require Adapt="79e6a3ab-5dfb-504d-930d-738a2a938a0e" begin
      function Adapt.adapt_storage(::Type{<:StaticArrays.SArray{S}},xs::Array) where S
          StaticArrays.SArray{S}(xs)
      end
    end
  end

  @require LabelledArrays="2ee39098-c373-598a-b85f-a56591580800" begin
    ismutable(::Type{<:LabelledArrays.LArray{T,N,Syms}}) where {T,N,Syms} = ismutable(T)
    can_setindex(::Type{<:LabelledArrays.SLArray}) = false
  end

  @require Tracker="9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c" begin
    ismutable(::Type{<:Tracker.TrackedArray}) = false
    can_setindex(::Type{<:Tracker.TrackedArray}) = false
    fast_scalar_indexing(::Type{<:Tracker.TrackedArray}) = false
    aos_to_soa(x::AbstractArray{<:Tracker.TrackedReal,N}) where N = Tracker.collect(x)
  end

    @require CuArrays="3a865a2d-5b23-5a0f-bc46-62713ec82fae" begin
    @require Adapt="79e6a3ab-5dfb-504d-930d-738a2a938a0e" begin
      include("cuarrays.jl")
    end
    @require DiffEqBase="2b5f629d-d688-5b77-993f-72d75c75574e" begin
      # actually do QR
      lu_instance(A::CuArrays.CuMatrix{T}) where T = CuArrays.CUSOLVER.CuQR(similar(A, 0, 0), similar(A, 0))
    end
  end

  @require CUDA="052768ef-5323-5732-b1bb-66c8b64840ba" begin
    @require Adapt="79e6a3ab-5dfb-504d-930d-738a2a938a0e" begin
      include("cuarrays2.jl")
    end
    @require DiffEqBase="2b5f629d-d688-5b77-993f-72d75c75574e" begin
      # actually do QR
      lu_instance(A::CUDA.CuMatrix{T}) where T = CUDA.CUSOLVER.CuQR(similar(A, 0, 0), similar(A, 0))
    end
  end

  @require BandedMatrices="aae01518-5342-5314-be14-df237901396f" begin
    function findstructralnz(x::BandedMatrices.BandedMatrix)
      l,u=BandedMatrices.bandwidths(x)
      rowsize,colsize= Base.size(x)
      rowind=BandedMatrixIndex(rowsize,colsize,l,u,true)
      colind=BandedMatrixIndex(rowsize,colsize,l,u,false)
      (rowind,colind)
    end

    has_sparsestruct(::Type{<:BandedMatrices.BandedMatrix}) = true
    is_structured(::Type{<:BandedMatrices.BandedMatrix}) = true
    fast_matrix_colors(::Type{<:BandedMatrices.BandedMatrix}) = true

    function matrix_colors(A::BandedMatrices.BandedMatrix)
        l,u=BandedMatrices.bandwidths(A)
        width=u+l+1
        _cycle(1:width, Base.size(A,2))
    end

  end

  @require BlockBandedMatrices="ffab5731-97b5-5995-9138-79e8c1846df0" begin
    @require BlockArrays="8e7c35d0-a365-5155-bbbb-fb81a777f24e" begin
      function findstructralnz(x::BlockBandedMatrices.BlockBandedMatrix)
        l,u=BlockBandedMatrices.blockbandwidths(x)
        nrowblock=BlockBandedMatrices.blocksize(x,1)
        ncolblock=BlockBandedMatrices.blocksize(x,2)
        rowsizes=BlockArrays.blocklengths(axes(x,1))
        colsizes=BlockArrays.blocklengths(axes(x,2))
        BlockBandedMatrixIndex(nrowblock,ncolblock,rowsizes,colsizes,l,u)
      end

      function findstructralnz(x::BlockBandedMatrices.BandedBlockBandedMatrix)
        l,u=BlockBandedMatrices.blockbandwidths(x)
        lambda,mu=BlockBandedMatrices.subblockbandwidths(x)
        nrowblock=BlockBandedMatrices.blocksize(x,1)
        ncolblock=BlockBandedMatrices.blocksize(x,2)
        rowsizes=BlockArrays.blocklengths(axes(x,1))
        colsizes=BlockArrays.blocklengths(axes(x,2))
        BandedBlockBandedMatrixIndex(nrowblock,ncolblock,rowsizes,colsizes,l,u,lambda,mu)
      end

      has_sparsestruct(::Type{<:BlockBandedMatrices.BlockBandedMatrix}) = true
      has_sparsestruct(::Type{<:BlockBandedMatrices.BandedBlockBandedMatrix}) = true
      is_structured(::Type{<:BlockBandedMatrices.BlockBandedMatrix}) = true
      is_structured(::Type{<:BlockBandedMatrices.BandedBlockBandedMatrix}) = true
      fast_matrix_colors(::Type{<:BlockBandedMatrices.BlockBandedMatrix}) = true
      fast_matrix_colors(::Type{<:BlockBandedMatrices.BandedBlockBandedMatrix}) = true

      function matrix_colors(A::BlockBandedMatrices.BlockBandedMatrix)
          l,u=BlockBandedMatrices.blockbandwidths(A)
          blockwidth=l+u+1
          nblock=BlockBandedMatrices.blocksize(A,2)
          cols=BlockArrays.blocklengths(axes(A,2))
          blockcolors=_cycle(1:blockwidth,nblock)
          #the reserved number of colors of a block is the maximum length of columns of blocks with the same block color
          ncolors=[maximum(cols[i:blockwidth:nblock]) for i in 1:blockwidth]
          endinds=cumsum(ncolors)
          startinds=[endinds[i]-ncolors[i]+1 for i in 1:blockwidth]
          colors=[(startinds[blockcolors[i]]:endinds[blockcolors[i]])[1:cols[i]] for i in 1:nblock]
          vcat(colors...)
      end

      function matrix_colors(A::BlockBandedMatrices.BandedBlockBandedMatrix)
          l,u=BlockBandedMatrices.blockbandwidths(A)
          lambda,mu=BlockBandedMatrices.subblockbandwidths(A)
          blockwidth=l+u+1
          subblockwidth=lambda+mu+1
          nblock=BlockBandedMatrices.blocksize(A,2)
          cols=BlockArrays.blocklengths(axes(A,2))
          blockcolors=_cycle(1:blockwidth,nblock)
          #the reserved number of colors of a block is the min of subblockwidth and the largest length of columns of blocks with the same block color
          ncolors=[min(subblockwidth,maximum(cols[i:blockwidth:nblock])) for i in 1:min(blockwidth,nblock)]
          endinds=cumsum(ncolors)
          startinds=[endinds[i]-ncolors[i]+1 for i in 1:min(blockwidth,nblock)]
          colors=[_cycle(startinds[blockcolors[i]]:endinds[blockcolors[i]],cols[i]) for i in 1:nblock]
          vcat(colors...)
      end
    end
  end
  @require OffsetArrays="6fe1bfb0-de20-5000-8ca7-80f57d26f881" begin
      size(A::OffsetArrays.OffsetArray) = size(parent(A))
      strides(A::OffsetArrays.OffsetArray) = strides(parent(A))
      # offsets(A::OffsetArrays.OffsetArray) = map(+, A.offsets, offsets(parent(A)))
      device(::OffsetArrays.OffsetArray) = CheckParent()
      contiguous_axis(A::OffsetArrays.OffsetArray) = contiguous_axis(parent(A))
      contiguous_batch_size(A::OffsetArrays.OffsetArray) = contiguous_batch_size(parent(A))
      stride_rank(A::OffsetArrays.OffsetArray) = stride_rank(parent(A))
  end
end

include("static.jl")
include("ranges.jl")
include("indexing.jl")
include("stridelayout.jl")

end


text:
module Thermo

__precompile__()

export EquilibriumData, Component, Mixture
export SRK, PR, RK, VdW, Equation
export PRESSURE, TEMPERATURE, LIQUID, VAPOR
export CriticalPoint, EquilibriumPoint
export solve, dump_json, fit

# include("equations.jl")
include("equilibrium.jl")

end


text:
using DocumentDistances

sdd = SinkhornDocumentDistance()

# From Kusner's WMD paper they used this illustrative example of two
# quite similar sentences (documents) which are similar semantically
# but not syntactically.
d1 = ["obama", "speaks", "media", "illinois"]
d2 = ["president", "greets", "press", "chicago"]
d12 = evaluate(sdd, d1, d2)

# To compare to something let's take a different document which we expect
# to have a larger distance to the ones above than they have to each other.
d3 = ["lawyer", "tanks", "car", "africa"]

d13 = evaluate(sdd, d1, d3)
d23 = evaluate(sdd, d2, d3)

@assert d12 < d13
@assert d12 < d23


text:
function makeshow(CV, CVstr, fields)
    exs = [:(showcompact(io, $(fields[d]))) for d = 1:length(fields)]
    exc = [d < length(fields) ? (:(print(io, ','))) : (:(print(io, ')'))) for d = 1:length(fields)]
    exboth = hcat(exs, exc)'
    ex = Expr(:block, exboth...)
    eval(quote
        function show{T,f}(io::IO, c::$CV{FixedPointNumbers.UfixedBase{T,f}})
            print(io, "$($CVstr){Ufixed", f, "}(")
            $ex
        end
    end)
end

for (CV, CVstr, fields) in ((RGB,  "RGB",  (:(c.r),:(c.g),:(c.b))),
                            (XYZ,  "XYZ",  (:(c.x),:(c.y),:(c.z))),
                            (RGBA, "RGBA", (:(c.c.r),:(c.c.g),:(c.c.b),:(c.alpha))))
    makeshow(CV, CVstr, fields)
end

# Displaying color swatches
# -------------------------

function writemime(io::IO, ::MIME"image/svg+xml", c::ColorValue)
    write(io,
        """
        <?xml version"1.0" encoding="UTF-8"?>
        <!DOCTYPE svg PUBLIC "-//W3C//DTD SVG 1.1//EN"
         "http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd">
        <svg xmlns="http://www.w3.org/2000/svg" version="1.1"
             width="25mm" height="25mm" viewBox="0 0 1 1">
             <rect width="1" height="1"
                   fill="#$(hex(c))" stroke="none"/>
        </svg>
        """)
end

function writemime{T <: ColorValue}(io::IO, ::MIME"image/svg+xml",
                                    cs::AbstractVecOrMat{T})
    m,n = ndims(cs) == 2 ? size(cs) : (1,length(cs))

    xsize,xpad = n > 50 ? (250/n,0) : n > 18 ? (5.,1) : n > 12 ? (10.,1) : n > 1 ? (15.,1) : (25.,0)
    ysize,ypad = m > 28 ? (150/m,0) : m > 14 ? (5.,1) : m > 9 ? (10.,1) : m > 1 ? (15.,1) : (25.,0)

    write(io,
        """
        <?xml version="1.0" encoding="UTF-8"?>
        <!DOCTYPE svg PUBLIC "-//W3C//DTD SVG 1.1//EN"
         "http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd">
        <svg xmlns="http://www.w3.org/2000/svg" version="1.1"
             width="$(n*xsize)mm" height="$(m*ysize)mm"
             shape-rendering="crispEdges">
        """)

    for i in 1:m, j in 1:n
        c = ndims(cs) == 2 ? cs[i, j] : cs[j]
        write(io,
            """
            <rect x="$((j-1)*xsize)mm" y="$((i-1)*ysize)mm"
                  width="$(xsize - xpad)mm" height="$(ysize - ypad)mm"
                  fill="#$(hex(c))" stroke="none" />
            """)
    end

    write(io, "</svg>")
end


text:
# utility functions

hertz(value::Real) = value * Hz
hertz{F <: Real}(value::SIUnits.SIQuantity{F,0,0,-1,0,0,0,0,0,0}) = value

tofloat{T}(array::AbstractArray{T}) = map(Float32, array)
tofloat(array::AbstractArray{Float32}) = array

fft_frequencies(samplerate::Real, nfft::Int) = collect(linspace(0f0, samplerate / 2f0, (nfft >> 1) + 1))

"""returns the number of frames when the signal is partitioned into overlapping frames"""
nframes(length::Int, framesize::Int, hopsize::Int) = div(length - framesize, hopsize) + 1

"""Provides a view of a signal partitioned into overlapping frames, similar to MATLAB's buffer"""
immutable FrameView{T<:AbstractVector} <: AbstractVector{Vector}
    original::T
    framesize::Int
    hopsize::Int
    nframes::Int

    function FrameView(original, framesize, hopsize)
        # n = noverlap is a problem - the algorithm will not terminate.
        #@boundscheck((0 < hopsize <= framesize), error("hopsize must be between zero and framesize"))
        new(original, framesize, hopsize, nframes(length(original), framesize, hopsize))
    end
end

function Base.getindex(x::FrameView, i::Int)
    #@boundscheck((i >= 1 && i <= x.k), Base.throw_boundserror(x, i))
    slice(x.original, (i-1) * x.hopsize + (1:x.framesize))
end
Base.start(x::FrameView) = 1
Base.next(x::FrameView, i::Int) = (x[i], i+1)
Base.done(x::FrameView, i::Int) = i > x.nframes
Base.size(x::FrameView) = (x.nframes,)

"""returns the DCT filters"""
function dct(nfilters::Int, ninput::Int)
    basis = Array(Float32, nfilters, ninput)
    samples = (1f0:2f0:2ninput) * π / 2ninput
    for i = 1:nfilters
        basis[i, :] = cos(i * samples)
    end

    basis *= sqrt(2f0/ninput)
    basis
end

"""returns the number of times that the value changes its sign"""
function zero_crossings{T<:Real}(array::AbstractVector{T}, size::Int = length(array), offset::Int = 0)
    result = 0
    previous = 0
    for i in offset + (1:size)
        number = array[i]
        sgn = number == 0 ? 0 : number > 0 ? 1 : -1
        if sgn != previous
            result += 1
        end
        previous = sgn
    end
    result
end


text:
function PresetTimeCallback(tstops,user_affect!;
                            initialize = SciMLBase.INITIALIZE_DEFAULT,
                            filter_tstops = true, kwargs...)

    condition = function (u, t, integrator)
      t in tstops
    end

    # Call f, update tnext, and make sure we stop at the new tnext
    affect! = function (integrator)
        user_affect!(integrator)
        nothing
    end

    # Initialization: first call to `f` should be *before* any time steps have been taken:
    initialize_preset = function (c, u, t, integrator)
        initialize(c, u, t, integrator)
        if filter_tstops
            tdir = integrator.tdir
            _tstops = tstops[@.((tdir*tstops > tdir*integrator.sol.prob.tspan[1]) * (tdir*tstops < tdir*integrator.sol.prob.tspan[2]))]
            add_tstop!.((integrator,), _tstops)
        else
            add_tstop!.((integrator,), tstops)
        end
        if t ∈ tstops
            user_affect!(integrator)
        end
    end
    DiscreteCallback(condition, affect!; initialize = initialize_preset, kwargs...)
end

export PresetTimeCallback


text:
# CoulombModels.jl

"""
Applies Coulomb's model to compute the friction force at the defined state.
"""
function friction(mdl::CoulombModel, nrm::Number, vel::Number)
    if vel == 0.0
        F = mdl.coefficient * nrm
    else
        F = mdl.coefficient * nrm * sign(vel)
    end
    return F
end

# ------------------------------------------------------------------------------
# Required routines to support friction model fitting.

function model_from_array(mdl::CoulombModel, x::Array{T}) where T <: Number
    CoulombModel(x[1])
end

function model_to_array(x::CoulombModel)
    [x.coefficient]
end


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl

given a sequence of observations X₁,...,Xᵢ with each Xᵢ in Rᵈ,
find a density function f maximizing the log likelihood of the observations
(equivalent to maximizing geomean of f.(X))
maximize    (prod_i f(Xᵢ))^(1/n)
subject to  ∫f = 1
            f ≥ 0
=#

import DelimitedFiles

struct DensityEstNative{T <: Real} <: ExampleInstanceNative{T}
    dataset_name::Symbol
    X::Matrix{T}
    deg::Int
    use_wsos::Bool # use WSOS cone formulation, else PSD formulation
    hypogeomean_obj::Bool # use geomean objective, else sum of logs objective
    use_hypogeomean::Bool # use hypogeomean cone, else 3-dim entropy formulation
end

function DensityEstNative{T}(
    dataset_name::Symbol,
    deg::Int,
    use_wsos::Bool,
    hypogeomean_obj::Bool,
    use_hypogeomean::Bool,
    ) where {T <: Real}
    X = DelimitedFiles.readdlm(joinpath(@__DIR__, "data", "$dataset_name.txt"))
    X = convert(Matrix{T}, X)
    return DensityEstNative{T}(dataset_name, X, deg, use_wsos,
        hypogeomean_obj, use_hypogeomean)
end

function DensityEstNative{T}(num_obs::Int, n::Int, args...) where {T <: Real}
    X = randn(T, num_obs, n)
    return DensityEstNative{T}(:Random, X, args...)
end

function build(inst::DensityEstNative{T}) where {T <: Real}
    X = inst.X
    (num_obs, n) = size(X)
    domain = PolyUtils.BoxDomain{T}(-ones(T, n), ones(T, n)) # domain is unit box

    # rescale X to be in unit box
    minX = minimum(X, dims = 1)
    maxX = maximum(X, dims = 1)
    X .-= (minX + maxX) / 2
    X ./= (maxX - minX) / 2

    # setup interpolation
    halfdeg = div(inst.deg + 1, 2)
    (U, _, Ps, V, w) = PolyUtils.interpolate(domain, halfdeg,
        calc_V = true, get_quadr = true)
    F = qr!(Array(V'), ColumnNorm())
    V_X = PolyUtils.make_chebyshev_vandermonde(X, 2halfdeg)
    X_pts_polys = (F \ V_X')'

    cones = Cones.Cone{T}[]

    num_psd_vars = 0
    if inst.use_wsos
        # U variables
        h_poly = zeros(T, U)
        b_poly = T[]
        push!(cones, Cones.WSOSInterpNonnegative{T, T}(U, Ps))
    else
        # U polynomial coefficient variables plus PSD variables
        # there are length(Ps) new PSD variables:
        # store them scaled, lower triangle, row-wise
        psd_var_list = Matrix{T}[]
        nonneg_cone_size = 0
        for i in eachindex(Ps)
            L = size(Ps[i], 2)
            dim = Cones.svec_length(L)
            if dim == 1
                nonneg_cone_size += 1
            else
                if nonneg_cone_size > 0
                    push!(cones, Cones.Nonnegative{T}(nonneg_cone_size))
                end
                push!(cones, Cones.PosSemidefTri{T, T}(dim))
            end
            num_psd_vars += dim
            push!(psd_var_list, zeros(T, U, dim))
            idx = 1
            # we will work with PSD vars with scaled off-diagonals
            # relevant columns (not rows) in A need to be scaled by sqrt(2) also
            for k in 1:L
                for l in 1:(k - 1)
                    psd_var_list[i][:, idx] = Ps[i][:, k] .*
                        Ps[i][:, l] * sqrt(T(2))
                    idx += 1
                end
                psd_var_list[i][:, idx] = Ps[i][:, k] .* Ps[i][:, k]
                idx += 1
            end
        end
        if nonneg_cone_size > 0
            push!(cones, Cones.Nonnegative{T}(nonneg_cone_size))
        end
        A_psd = hcat(psd_var_list...)
        b_poly = zeros(T, U)
        h_poly = zeros(T, num_psd_vars)
    end

    if inst.hypogeomean_obj
        num_hypo_vars = 1
        if inst.use_hypogeomean
            G_likl = [
                -one(T) zeros(T, 1, U)
                zeros(T, num_obs) -X_pts_polys
                ]
            h_likl = zeros(T, 1 + num_obs)
            push!(cones, Cones.HypoGeoMean{T}(1 + num_obs))
            A_ext = zeros(T, 0, num_obs)
        else
            num_ext_geom_vars = 1 + num_obs
            ext_offset = 2 + U + num_psd_vars
            h_likl = zeros(T, 3 * num_obs + 2)
            # order: hypograph u, U f(obs) vars, psd vars, geomean ext vars (y, z)
            G_likl = zeros(T, 3 * num_obs + 2, ext_offset + num_obs)
            # y >= u, e'z >= 0
            G_likl[1, 1] = 1
            G_likl[1, ext_offset] = -1
            G_likl[2, (end - num_obs + 1):end] .= -1
            push!(cones, Cones.Nonnegative{T}(2))
            # f(x) <= y * log(z / y)
            row_offset = 3
            for i in 1:num_obs
                G_likl[row_offset, ext_offset + i] = 1
                G_likl[row_offset + 1, 2:(1 + U)] = -X_pts_polys[i, :]
                G_likl[row_offset + 2, ext_offset] = -1
                row_offset += 3
                push!(cones, Cones.EpiRelEntropy{T}(3))
            end
        end
    else
        num_hypo_vars = num_obs
        h_likl = zeros(T, 3 * num_obs)
        G_likl = zeros(T, 3 * num_obs, num_obs + U)
        offset = 1
        for i in 1:num_obs
            G_likl[offset + 2, (num_obs + 1):(num_obs + U)] = -X_pts_polys[i, :]
            h_likl[offset + 1] = 1
            G_likl[offset, i] = -1
            offset += 3
            push!(cones, Cones.HypoPerLog{T}(3))
        end
    end

    # extended formulation variables for hypogeomean are added after psd ones, so
    # psd vars already accounted for in hypogeomean_obj && !use_hypogeomean path
    if !inst.hypogeomean_obj || inst.use_hypogeomean
        G_likl = hcat(G_likl, zeros(T, size(G_likl, 1), num_psd_vars))
        num_ext_geom_vars = 0
    end
    c = zeros(T, num_hypo_vars + U + num_psd_vars + num_ext_geom_vars)
    @views c[1:num_hypo_vars] .= -1
    h = vcat(h_poly, h_likl)
    b = vcat(b_poly, one(T))

    if inst.use_wsos
        A = zeros(T, 1, num_hypo_vars + U + num_ext_geom_vars)
        A[1, num_hypo_vars .+ (1:U)] = w
        G = zeros(T, U + size(G_likl, 1), size(G_likl, 2))
        G[1:U, num_hypo_vars .+ (1:U)] = Diagonal(-I, U)
        G[(U + 1):end, :] = G_likl
    else
        A = zeros(T, U + 1, num_hypo_vars + U + num_psd_vars + num_ext_geom_vars)
        A[1:U, num_hypo_vars .+ (1:U)] = Diagonal(-I, U)
        A[1:U, (num_hypo_vars + U) .+ (1:num_psd_vars)] = A_psd
        A[U + 1, num_hypo_vars .+ (1:U)] = w
        G = zeros(T, num_psd_vars + size(G_likl, 1), size(G_likl, 2))
        G[1:num_psd_vars, (num_hypo_vars + U) .+ (1:num_psd_vars)] =
            Diagonal(-I, num_psd_vars)
        G[(num_psd_vars + 1):end, :] = G_likl
    end

    model = Models.Model{T}(c, A, b, G, h, cones)
    return model
end


text:
using SMLMAssociationAnalysis_NCB

# PARAMETERS

outputpath = "path/to/your/output/results.jld2"

percentilerankthreshold = 0.1
distancethreshold = 200

# END PARAMETERS

results = FileIO.load(outputpath)["results"]

# generate localization plots (not print quality)
for result ∈ results
    p = localizationsplot(results)
    savefig(p, "$(result.cell) localizations.png")
end

# generate molecule plots (not print quality)
for result ∈ results
    p = moleculesplot(results)
    savefig(p, "$(result.cell) molecules.png")
end

# generate distance-probability plots
for result ∈ results
    p = distanceprobabilityplot(results)
    savefig(p, "$(result.cell) distance-probability.png")
end

# calculate cell summary data (median distance, fraction bound)
percentilerankslessthan10replicate = [x.percentileranks .< percentilerankthreshold for x ∈ results]
bound = [(x.distances .< distancethreshold) .& (x.percentileranks .< percentilerankthreshold) for x ∈ replicateresults]
fractionbound = count.(bound) ./ length.(percentilerankslessthan10replicate)
mediandistances = map(x -> x.mediandistance, replicateresults)

localization1counts = [map(x -> x.group.localizations |> length, y.channels[1].molecules) |> sum for y ∈ replicateresults]
localization2counts = [map(x -> x.group.localizations |> length, y.channels[2].molecules) |> sum for y ∈ replicateresults]
molecule1counts = map(x -> x.channels[1].molecules |> length, replicateresults)
molecule2counts = map(x -> x.channels[2].molecules |> length, replicateresults)