text:
"""
Delete named dataset.

The driver will attempt to delete the named dataset in a driver specific
fashion. Full featured drivers will delete all associated files, database
objects, or whatever is appropriate. The default behaviour when no driver
specific behaviour is provided is to attempt to delete the passed name as a
single file.

It is unwise to have open dataset handles on this dataset when it is deleted.

### Parameters
* `pszFilename`     name of dataset to delete.

### Returns
`CE_None` on success, or `CE_Failure` if the operation fails.
"""
_deletedataset(dataset::GDALDriverH, pszFilename::Ptr{UInt8}) =
    GDALDeleteDataset(dataset, pszFilename)::CPLErr

function deletedataset(dataset::GDALDriverH, filename::String)
    result = _deletedataset(dataset, pointer(filename))
    (result == CE_Failure) && error("Failed to delete dataset")
    @assert result == CE_None
end

"""
Rename a dataset.

This may including moving the dataset to a new directory or filesystem.

It is unwise to have open dataset handles on this dataset when it is being
renamed.

### Parameters
* `pszNewName`  new name for the dataset.
* `pszOldName`  old name for the dataset.

### Returns
`CE_None` on success, or `CE_Failure` if the operation fails.
"""
_renamedataset(dataset::GDALDriverH,
               pszNewName::Ptr{UInt8},
               pszOldName::Ptr{UInt8}) =
    GDALRenameDataset(dataset, pszNewName, pszOldName)::CPLErr

function renamedataset(dataset::GDALDriverH,
                       newname::String,
                       oldname::String)
    result = _renamedataset(dataset, pointer(newname), pointer(oldname))
    (result == CE_Failure) && error("Failed to rename $oldname to $newname")
    @assert result == CE_None
end

"""
Copy the files associated with a dataset.

### Parameters
* `pszNewName`  new name for the dataset.
* `pszOldName`  old name for the dataset.

### Returns
`CE_None` on success, or `CE_Failure` if the operation fails.
"""
_copydatasetfiles(dataset::GDALDriverH,
                  pszNewName::Ptr{UInt8},
                  pszOldName::Ptr{UInt8}) =
    GDALCopyDatasetFiles(dataset, pszNewName, pszOldName)::CPLErr

function copydatasetfiles(dataset::GDALDriverH,
                          newname::String,
                          oldname::String)
    result = _copydatasetfiles(dataset, pointer(newname), pointer(oldname))
    (result == CE_Failure) && error("Failed to copy $oldname to $newname")
    @assert result == CE_None
end

"""
Validate the list of creation options that are handled by a driver.

This is a helper method primarily used by `Create()` and `CreateCopy()` to
validate that the passed in list of creation options is compatible with the
`GDAL_DMD_CREATIONOPTIONLIST` metadata item defined by some drivers.

See also: `GDALGetDriverCreationOptionList()`

If the `GDAL_DMD_CREATIONOPTIONLIST` metadata item is not defined, this
function will return `TRUE`. Otherwise it will check that the keys and values
in the list of creation options are compatible with the capabilities declared
by the `GDAL_DMD_CREATIONOPTIONLIST` metadata item. In case of incompatibility
a (non fatal) warning will be emited and `FALSE` will be returned.

### Parameters
* `hDriver`     the handle of the driver with whom the lists of creation option
must be validated
* `options`     the list of creation options. An array of strings, whose last
element is a `NULL` pointer

### Returns
`TRUE` if the list of creation options is compatible with the `Create()` and
`CreateCopy()` method of the driver, `FALSE` otherwise.
"""
_validatecreationoptions(hDriver::GDALDriverH, options::Ptr{Ptr{UInt8}}) =
    GDALValidateCreationOptions(hDriver, options)::Cint

function validatecreationoptions(hDriver::GDALDriverH,
                                 options::Vector{String})
    poptions = Ptr{Ptr{UInt8}}(pointer(options))
    Bool(_validatecreationoptions(hDriver, poptions))
end

"Fetch raster width in pixels."
_getrasterxsize(dataset::GDALDatasetH) = GDALGetRasterXSize(dataset)::Cint

"Fetch raster height in pixels."
_getrasterysize(dataset::GDALDatasetH) = GDALGetRasterYSize(dataset)::Cint

"Fetch the number of raster bands on this dataset."
_getrastercount(dataset::GDALDatasetH) = GDALGetRasterCount(dataset)::Cint

"""
Fetch a band object for a dataset from its index (from 1 to `GetRasterCount()`)
"""
_getrasterband(dataset::GDALDatasetH, i::Integer) =
  GDALGetRasterBand(dataset, i)::GDALRasterBandH

getrasterband(dataset::GDALDatasetH, i::Int) =
  _getrasterband(dataset, i)

"""
Add a band to a dataset.

This method will add a new band to the dataset if the underlying format
supports this action. Most formats do not.

Note that the new `GDALRasterBand` is not returned. It may be fetched after
successful completion of the method through `GetRasterBand(GetRasterCount())`
as the newest band will always be the last band.

### Parameters
* `eType`           the data type of the pixels in the new band.
* `papszOptions`    a list of NAME=VALUE option strings. The supported options
are format specific. `NULL` may be passed by default.

### Returns
`CE_None` on success or `CE_Failure` on failure.
"""
_addband(hDS::GDALDatasetH,
         eType::GDALDataType,
         papszOptions::Ptr{Ptr{UInt8}}) =
    GDALAddBand(hDS, eType, papszOptions)::CPLErr

function addband(hDS::GDALDatasetH,
                 eType::GDALDataType,
                 options::Vector{String})
    poptions = Ptr{Ptr{UInt8}}(pointer(options))
    result = _addband(hDS, eType, poptions)
    (result == CE_Failure) && error("Failed to add band")
end

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

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

The `nPixelSpace`, `nLineSpace` and `nBandSpace` parameters allow reading into
or writing from various organization of buffers.

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

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

### Returns
`CE_Failure` if the access fails, otherwise `CE_None`.
"""
_datasetrasterio(hDS::GDALDatasetH,
                 eRWFlag::GDALRWFlag,
                 nDSXOff::Integer,
                 nDSYOff::Integer,
                 nDSXSize::Integer,
                 nDSYSize::Integer,
                 pData::Ptr{Void},
                 nBXSize::Integer,
                 nBYSize::Integer,
                 eBDataType::GDALDataType,
                 nBandCount::Integer,
                 panBandCount::Ptr{Cint},
                 nPixelSpace::Integer = 0,
                 nLineSpace::Integer = 0,
                 nBandSpace::Integer = 0) =
    GDALDatasetRasterIO(hDS, eRWFlag, nDSXOff, nDSYOff, nDSXSize, nDSYSize,
                        pData, nBXSize, nBYSize, eBDataType, nBandCount,
                        panBandCount, nPixelSpace, nLineSpace,
                        nBandSpace)::CPLErr

function rasterio!{T <: Real}(hDS::GDALDatasetH,
                              buffer::Array{T, 3},
                              bands::Vector{Cint},
                              width::Integer,
                              height::Integer,
                              xoffset::Integer = 0,
                              yoffset::Integer = 0,
                              access::GDALRWFlag = GF_Read,
                              nPixelSpace::Integer = 0,
                              nLineSpace::Integer = 0,
                              nBandSpace::Integer = 0)
    xsize, ysize, zsize = size(buffer)
    nband = length(bands)
    @assert nband == zsize
    result = _datasetrasterio(hDS, access, xoffset, yoffset, width, height,
                              Ptr{Void}(pointer(buffer)), xsize, ysize,
                              _gdaltype(T), nband, pointer(bands),
                              nPixelSpace, nLineSpace, nBandSpace)
    (result == CE_Failure) && error("Failed to access raster bands")
    buffer
end

function rasterio!{T <: Real}(hDS::GDALDatasetH,
                              buffer::Array{T, 3},
                              bands::Vector{Cint},
                              access::GDALRWFlag = GF_Read,
                              nPixelSpace::Integer = 0,
                              nLineSpace::Integer = 0,
                              nBandSpace::Integer = 0)
    rasterio!(hDS, buffer, bands, _getrasterxsize(hDS), _getrasterysize(hDS),
              0, 0, access, nPixelSpace, nLineSpace, nBandSpace)
end

function rasterio!{T <: Real, U <: Integer}(hDS::GDALDatasetH,
                                            buffer::Array{T, 3},
                                            bands::Vector{Cint},
                                            rows::UnitRange{U},
                                            cols::UnitRange{U},
                                            access::GDALRWFlag = GF_Read,
                                            nPixelSpace::Integer = 0,
                                            nLineSpace::Integer = 0,
                                            nBandSpace::Integer = 0)
    width = cols[end] - cols[1] + 1
    width < 0 && error("invalid window width")
    height = rows[end] - rows[1] + 1
    height < 0 && error("invalid window height")
    rasterio!(hDS, buffer, bands, width, height, cols[1], rows[1], access,
              nPixelSpace, nLineSpace, nBandSpace)
end

"Fetch a format specific internally meaningful handle."
_getinternalhandle(dataset::GDALDatasetH, request::Ptr{UInt8}) =
    GDALGetInternalHandle(dataset, request)::Ptr{Void}

getinternalhandle(dataset::GDALDatasetH, request::String) =
    _getinternalhandle(dataset, Ptr{UInt8}(pointer(request)))

"""
Add one to dataset reference count. The reference is one after instantiation.

### Returns
the post-increment reference count.
"""
_referencedataset(dataset::GDALDatasetH) = GDALReferenceDataset(dataset)::Cint

"""
Subtract one from dataset reference count. The reference is one after
instantiation. Generally when the reference count has dropped to zero the
dataset may be safely deleted (closed).

### Returns
the post-decrement reference count.
"""
_dereferencedataset(dataset::GDALDatasetH) =
    GDALDereferenceDataset(dataset)::Cint

"""
Build raster overview(s).

If the operation is unsupported for the indicated dataset, then CE_Failure is
returned, and CPLGetLastErrorNo() will return CPLE_NotSupported.

### Parameters
* `pszResampling`   one of "NEAREST", "GAUSS", "CUBIC", "AVERAGE", "MODE",
"AVERAGE_MAGPHASE" or "NONE" controlling the downsampling method applied.
* `nOverviews`      number of overviews to build.
* `panOverviewList` the list of overview decimation factors to build.
* `nListBands`      number of bands to build overviews for in `panBandList`.
Build for all bands if this is 0.
* `panBandList`     list of band numbers.
* `pfnProgress`     a function to call to report progress, or `NULL`.
* `pProgressData`   application data to pass to the progress function.

### Returns
`CE_None` on success or `CE_Failure` if the operation doesn't work.

For example, to build overview level 2, 4 and 8 on all bands the following
call could be made:

```C
   int       anOverviewList[3] = { 2, 4, 8 };
   poDataset->BuildOverviews( "NEAREST", 3, anOverviewList, 0, NULL,
                              GDALDummyProgress, NULL );
```
"""
_buildoverviews(arg1::GDALDatasetH,
                pszResampling::Ptr{UInt8},
                nOverviews::Integer,
                panOverviewList::Ptr{Cint},
                nListBands::Integer,
                panBandList::Ptr{Cint},
                pfnProgress::GDALProgressFunc,
                pProgressData::Ptr{Void}) =
    GDALBuildOverviews(pszResampling, nOverviews, panOverviewList, nListBands,
                       panBandList, pfnProgress, pProgressData)::CPLErr

function buildoverviews(hDS::GDALDatasetH,
                        resampling::String,
                        overviewlist::Vector{Cint},
                        bandList::Vector{Cint})
    result = _buildoverviews(hDS, pointer(resampling),
                             length(overviewlist),
                             pointer(overviewlist),
                             length(bandlist),
                             pointer(bandlist),
                             C_NULL, C_NULL)
    (result == CE_Failure) && error("Failed to build overviews")
end

"""
Fetch all open GDAL dataset handles.

NOTE: This method is not thread safe. The returned list may change at any time
and it should not be freed.

### Parameters
* `pnCount`     integer into which to place the count of dataset pointers being
returned.

### Returns
a pointer to an array of dataset handles.
"""
_getopendatasets(hDS::Ptr{Ptr{GDALDatasetH}}, pnCount::Ptr{Cint}) =
    GDALGetOpenDatasets(hdS, pnCount)

"Return access flag."
_getaccess(dataset::GDALDatasetH) = GDALGetAccess(dataset)::Cint

"""
Flush all write cached data to disk.

Any raster (or other GDAL) data written via GDAL calls, but buffered
internally will be written to disk.

The default implementation of this method just calls the `FlushCache()` method
on each of the raster bands and the `SyncToDisk()` method on each of the
layers. Conceptionally, calling `FlushCache()` on a dataset should include any
work that might be accomplished by calling `SyncToDisk()` on layers in that
dataset.

Using this method does not prevent use from calling `GDALClose()` to properly
close a dataset and ensure that important data not addressed by `FlushCache()`
is written in the file.
"""
_flushcache(dataset::GDALDatasetH) = GDALFlushCache(dataset)

"""
Adds a mask band to the dataset.

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 one band.
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

### Parameters
* `nFlags`  ignored. `GMF_PER_DATASET` will be assumed.

### Returns
`CE_None` on success or `CE_Failure` on an error.
"""
_createdatasetmaskband(hDS::GDALDatasetH,
                       nFlags::Integer = GMF_PER_DATASET) =
    GDALCreateDatasetMaskBand(hDS, nFlags)::CPLErr

function createdatasetmaskband(hDS::GDALDatasetH)
    result = _createdatasetmaskband(hDS)
    (result == CE_Failure) && error("Failed to create dataset mask band")
end

"""
Copy all dataset raster data.

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

This function is primarily intended to support implementation of driver
specific `CreateCopy()` functions. It implements efficient copying, in
particular "chunking" the copy in substantial blocks and, if appropriate,
performing the transfer in a pixel interleaved fashion.

Currently the only papszOptions value supported are : "INTERLEAVE=PIXEL" to
force pixel interleaved operation and "COMPRESSED=YES" to force alignment on
target dataset block sizes to achieve best compression. More options may be
supported in the future.

### Parameters
* `hSrcDS`          the source dataset
* `hDstDS`          the destination dataset
* `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.
"""
_datasetcopywholeraster(hSrcDS::GDALDatasetH,
                        hDstDS::GDALDatasetH,
                        papszOptions::Ptr{Ptr{UInt8}},
                        pfnProgress::GDALProgressFunc,
                        pProgressData::Ptr{Void}) =
    GDALDatasetCopyWholeRaster(hSrcDS, hDstDS, papszOptions, pfnProgress,
                               pProgressData)::CPLErr

function datasetcopywholeraster(sourceDS::GDALDatasetH,
                                destDS::GDALDatasetH,
                                options::Vector{String})
    result = _datasetcopywholeraster(sourceDS, destDS,
                                     Ptr{Ptr{UInt8}}(pointer(options)))
    (result == CE_Failure) && error("Failed to copy whole raster")
end


text:
function findtrackgraph(Atrack)
    #TRACKGRAPH forms graph of locations from Track Associative array.

        # Find 1 hop and >1 hop tracks.
        AtrackHop = sum(Atrack,1)
        Hop1 = getcol(AtrackHop == 1)
        Hop2 = getcol(AtrackHop > 1)

        # Get track list.  Naturally comes out sorted by p.
        t1,p1,x1 = find(Atrack[:,Hop1])
        AtrackGraph1 = Assoc(x1,x1,1,(+))

        t2,p2,x2 = find(Atrack[:,Hop2])

        # Create matrices and shifted matrices.
        p22 = circshift(p2,-1)
        x22 = circshift(x2,-1)

        # Find where p21 and p22 are the same.
        test = p2 .== p22
        x21 = x2[test]
        x22 = x22[test]

        AtrackGraph2 = Assoc(x21,x22,1)

        AtrackGraph = AtrackGraph1 + AtrackGraph2

    end

ext: jl


text:
using CoreNLPClient

res1 = CoreNLP("http://localhost:9000", "I am a Julia programmer. I love Julia programming.")

res2 = CoreNLP("localhost:9000", "I am a Julia programmer. I love Julia programming.")

res3 = CoreNLP("localhost", "I am a Julia programmer. I love Julia programming.")

res4 = CoreNLP("I am a Julia programmer. I love Julia programming.")

annotations = getNLPAnnotations("I am a Julia programmer. I love Julia programming.")


text:
struct Data end
struct Grid end
struct StaggeredField end
struct Variable end

data_kind(::Type{<:PlainVariableBlockHeader}) = Data()
data_kind(::Type{<:PointVariableBlockHeader}) = Data()
data_kind(::Type{<:PlainMeshBlockHeader}) = Grid()
data_kind(::Type{<:PointMeshBlockHeader}) = Grid()

discretization_type(::Type{<:PlainVariableBlockHeader}) = StaggeredField()
discretization_type(::Type{<:PlainMeshBlockHeader}) = StaggeredField()
discretization_type(::Type{<:PointVariableBlockHeader}) = Variable()
discretization_type(::Type{<:PointMeshBlockHeader}) = Variable()


text:
export ScoreQuadratic

"""
    ScoreQuadratic <: ScoringRule

 Calculate Quadratic Scores for a set of outcomes and predictions

## References
* https://en.wikipedia.org/wiki/Scoring_rule
*

## Parameters
* `normalise::false`: if true, normalise scores so that for binary choice
     + S( [0.5,0.5], 1) = 0.0
     + S( [0.0,1.0], 1) = 1.0

```
"""
struct ScoreQuadratic <: ScoringRule
    normalise::Bool
end
ScoreQuadratic(; normalise::Bool=false) = ScoreQuadratic(normalise)

function score_direction(srule::ScoreQuadratic)
    return 1
end

function scoring_function(srule::ScoreQuadratic, predicted_probabilities::Array{Float64,1}, outcome::Int)
    if abs( sum(predicted_probabilities) - 1.0 ) > 1.0e-2
        error("probabilities must sum to 1")
    end
    result = 2*predicted_probabilities[outcome] - sum( predicted_probabilities.^2 )
    if srule.normalise
        return 2*result - 1
    else
        return result
    end
end


text:
 	centuries = ("I", "II", "III", "IV", "V", "VI", "VII", "VIII", "IX", "X", "XI", "XII", "XIII", "XIV", "XV", "XVI","XVII", "XVIII", "XIX", "XX", "XXI")

 	year_top =2100

 	unix_time()= rand(0: Dates.value((Dates.DateTime(Dates.now())) - (Dates.DateTime(1970,1,1))))

 	function date_time()
 		fec = Dates.DateTime(rand(1970:year_top),rand(01:12),1,rand(01:23),rand(00:59),rand(00:59)) + Dates.Day(rand(01:31))
 		return Dates.format(fec, "Y-m-d HH:MM:SS")
 	end

 	date_time_ad()= Dates.format(iso8601(), "Y-m-d HH:MM:SS")

 	iso8601()=
 	Dates.DateTime(rand(1970:year_top),rand(01:12),1,rand(01:23),rand(00:59),rand(00:59)) + Dates.Day(rand(01:31))

 	function date(pattern="Y-m-d")
 		x = Dates.DateTime(rand(1970:year_top),rand(01:12),1,rand(01:23),rand(00:59),rand(00:59)) + Dates.Day(rand(01:31))

 		return Dates.format(x, pattern)

 	end

 	function time(pattern="H:M:S")
 		x= Dates.DateTime(1,1,1,rand(01:23),rand(00:59),rand(00:59))

 		return Dates.format(x, pattern)
 	end

 	function date_time_between(start_date="-30y", end_date="now")

 		try
 			start_date = Dates.DateTime(Dates.now()) + Dates.Year(parse(Int,start_date[1:end-1]))
 		catch
 			start_date = Dates.DateTime(Dates.now())
 		end

 		try
 		end_date = Dates.DateTime(Dates.now()) + Dates.Year(parse(Int,end_date[1:end-1]))
	 	catch
	 	end_date = Dates.DateTime(Dates.now())
	 	end

	 dr=(start_date:Day(1):end_date)[rand(1:end)]
	 return Dates.format(dr, "Y-m-d HH:MM:SS")
	end

function random_datetime(bf_now,af_now,time_start,now_time,time_finish)
    if bf_now && af_now
		return(time_start:Day(1):time_finish)[rand(1:end)]
	elseif (bf_now==false) && af_now
		return(now_time:Day(1):time_finish)[rand(1:end)]
	else
		return(time_start:Day(1):now_time)[rand(1:end)]
	end
	return  now_time
end

function date_time_this_century(;before_now=true, after_now=false)
	r = Dates.DateTime(Dates.now())
	this_century_start = Dates.DateTime((Int(Dates.year(r)) - (Dates.year(r) % 100)), 1, 1)
	this_century_finish = (Dates.DateTime(Dates.year(this_century_start) + 100, 1, 1))-Dates.Day(1)

	return Dates.format(random_datetime(before_now,after_now,this_century_start,r,this_century_finish), "Y-m-d HH:MM:SS")
end

function  date_time_this_decade(;before_now=true, after_now=false)
	r = Dates.DateTime(Dates.now())
	this_decade_start = Dates.DateTime((Int(Dates.year(r)) - (Dates.year(r) % 10)), 1, 1)
	this_decade_finish = (Dates.DateTime(Dates.year(this_decade_start) + 10, 1, 1))-Dates.Day(1)

	return Dates.format(random_datetime(before_now,after_now,this_decade_start,r,this_decade_finish), "Y-m-d HH:MM:SS")
end

function date_time_this_year(;before_now=true, after_now=false)
	r = Dates.DateTime(Dates.now())
	this_year_start = Dates.DateTime(Dates.year(r), 1, 1)
	this_year_finish = this_year_start+Dates.Year(1)-Dates.Day(1)

	return Dates.format(random_datetime(before_now,after_now,this_year_start,r,this_year_finish), "Y-m-d HH:MM:SS")
end

function date_time_this_month(;before_now=true, after_now=false)
	r = Dates.DateTime(Dates.now())
	this_month_start = r-Dates.Day(r)+Dates.Day(1)
	this_month_finish = this_month_start+Dates.Month(1)-Dates.Day(1)

	return Dates.format(random_datetime(before_now,after_now,this_month_start,r,this_month_finish), "Y-m-d HH:MM:SS")
end

ap = ("AM","PM")
am_pm()= ap[rand(1:2)]

day_of_month()= rand(1:31)

day_of_week()= Dates.dayname(Date(2015,07,rand(12:18)))

months()=rand(1:12)

month_name()= Dates.monthname(Date(2015,rand(1:12),15))

year()=rand(1920:year_top)

century()= centuries[rand(1:length(centuries))]

function timezone()
	executor(data["en"]["faker"]["address"]["time_zone"])
end


text:
using Documenter, TrajTools

makedocs(
    sitename="TrajTools Documentation",
    # prettyURLs disabled to make local browsing easier -> should be enabled for online hosting
    format = Documenter.HTML(prettyurls = false),
    pages = [
        "Home" => "index.md",
        "Reading trajectories" => "trajecio.md",
        ]
    )


text:
function openhtmldoc(filepath::AbstractString)
    NOSHOW && return
    is_linux() && run(`xdg-open $filepath`)
    is_apple() && run(`open $filepath`)
    is_windows() && run(`cmd /c start $filepath`)
end

_basic(p::Plot) = "Plot(\"$(p.title)\" with $(length(p.datacolumns)) datacolumns)"

function Base.show(io::IO, ::MIME"text/html", p::Plot)
    print(io, renderplot(p, true))
    print(io, "
", _basic(p), "
")
end

Base.show(io::IO, ::MIME"text/plain", p::Plot) = print(io, _basic(p))

# If in the notebook, activate a notebook display
if !isdefined(Main, :IJulia)
    mutable struct BokehDisplay <: Display; end
    pushdisplay(BokehDisplay())

    function Base.display(d::BokehDisplay, p::Plot)
        display("text/plain", p)
        AUTOOPEN && showplot(p)
    end
end

function showplot(p::Plot, filename::NullString=nothing)
    genplot(p, filename)
    openhtmldoc(filename == nothing ? p.filename : filename)
end

showplot() = showplot(CURPLOT)
showplot(filename::AbstractString) = showplot(CURPLOT, filename)

function setupnotebook()
    jspath, csspath = _bokehjs_paths(!DEBUG)
    jscss = _render_jscss(jspath, csspath, true)
    display("text/html", jscss)
    display("text/html", "
BokehJS successfully loaded.
")
end


text:
"""
  BEDQuery
   A query object holds the coordinates
"""

type BEDQuery
  qid::String
  seq_id::String
  start_pos::Integer
  stop_pos::Integer
end

"""
   make_query_path
"""
function make_query_str(queries::Array{BEDQuery} )
    query_array = String[]
    for query in queries
        str="$(query.seq_id):$(query.start_pos)-$(query.stop_pos)"
        push!(query_array, str)
    end
    return query_array
end

"""
    bedquery2list
"""
function bedquery2list(Q::BEDQuery)
    return (Q.qid,Q.seq_id,Q.start_pos,Q.stop_pos)
end

"""
    initiate

    Get query seq_id,start_pos,stop_pos and new vector
"""

function initiate(query::BEDQuery)
    query_seq_id    = query.seq_id
    query_start_pos = query.start_pos
    query_stop_pos  = query.stop_pos
    query_length    = query_stop_pos - query_start_pos
    V=fill(float32(0.0),query_length)
    return (query_seq_id,query_start_pos,query_stop_pos,query_length,V)
end

"""
    Query
"""

function query(filepath::String,query::Array{BEDQuery})
    query_idx=1
    query_str_array = make_query_str(query)
    cmd = `tabix -p bed $filepath $query_str_array`
    D = Dict{String,Array{Float32}}()
    (query_seq_id,query_start_pos,
     query_stop_pos,query_length,V) = initiate(query[query_idx])
    for line in eachline(cmd)
        (seq_id,start,stop,score)=split( chomp(line), '\t' )
        start_pos = int32(start)
        stop_pos  = int32(stop)

        idx_start_pos = ( start_pos - query_start_pos ) + 1
        idx_stop_pos  = ( stop_pos  - query_start_pos ) + 1

        try
           for i in (idx_start_pos:idx_stop_pos )
               if i < 1
                  continue
               end
               if i > query_length
                  continue
               end
               V[i]= float32(score)
               println("$i ($idx_start_pos) $seq_id, $start,$stop,$score")
           end
        catch e
           println("$seq_id, $start,$stop,$score")
           error(e)
        end

        # If current stop position is greater than the current query stop position
        # then the next interval must be from the next query
        if ( ( stop_pos >= query[query_idx].stop_pos) | (seq_id != query_seq_id) )
            D[ query_str_array[query_idx] ]=copy(V)
            if query_idx < length(query)
                query_idx += 1
                (query_seq_id,query_start_pos,
                 query_stop_pos,query_length,V) = initiate( query[query_idx] )
            end

        end
    end
    D[ query_str_array[query_idx] ]=copy(V)
    return D
end

"""
   bgzip_and_index

   bgzips - then uses tabix -p bed to create an index
"""
function bgzip_and_index(filepath; outputpath=filepath*".gz")
   cmd_bgzip = `bgzip -c $filepath` |> outputpath
   readall(cmd_bgzip)
   cmd_tabix = `tabix -f -p bed $outputpath`
   readall(cmd_tabix)
end

#println("CHANGE AT $seq_id, $start,$stop,$score")
# println("NEXT INTERVAL $seq_id, $query_start_pos,$query_stop_pos,$score")
#println("IS $start_pos GREATER THAN $(query[query_idx].stop_pos) ?")


text:
using AnyMODResults
using Test
using Suppressor
using CSV
using DataFrames
using JSON

function compare_table_values(loaded_Result, test_result)

    error_flag = false

    for x in eachindex(loaded_Result[:,1]), y in names(loaded_Result)[2:end]
        target = loaded_Result[x,y]
        result = test_result[x,y]
        if ismissing(target)
            ismissing(result) || (error_flag = true)
        elseif target isa AbstractString
            String(target) == String(result) || (error_flag = true)
        elseif target isa Number && result isa Number
            target ≈ result || (error_flag = true)
        else
            error_flag = true
        end

        if error_flag
            @warn "Test failed" target=target actual_result=result
            return false
        end
    end

    return true
end

include("testset1.jl")
include("testset2.jl")
include("testset3.jl")


text:
"""
$(SIGNATURES)

Computes standardized matrix on GPU. Theoretically this function
works with CPU as well. But it is not as fast as the
`get_standardized_matrix` function, which utilizes multi threads on CPU.

# Arguments
- `m` : a matrix to be standardized

# Output:
returns the standardized matrix of `m`

"""
function get_standardized_matrix_gpu(m::AbstractArray{<:Real,2})
    return (m .- mean(m, dims=1)) ./ std(m, corrected=false, dims=1)
end

"""
$(SIGNATURES)
Computes standardized matrix on CPU utilizing multi-threads.

# Arguments
- `m` : a matrix to be standardized

# Output:
returns the standardized matrix of `m`

"""
function get_standardized_matrix(mat::AbstractArray{<:Real,2})
    m = copy(mat)
    Threads.@threads for col in 1:size(m)[2]
        summ = 0.0f0
        rows = size(m)[1]
        for row in 1:rows
            summ += m[row, col]
        end
        mean = summ/rows
        sums = 0.0f0
        for row in 1:rows
            sums += (m[row,col] - mean)^2
        end
        std = sqrt(sums/rows)
        for row in 1:rows
            m[row,col] = (m[row,col]-mean)/std
        end
    end
    return m
end

"""
$(SIGNATURES)

# Arguments
- `x` : input matrix

# Output:
returns px

"""
function calculate_px(x::AbstractArray{<:Real,2})
    XtX = transpose(x)*x
    result = x*inv(XtX)*transpose(x)
    # display(result)
    return result
end


text:
using YAActL

function fib(D::Dict{Int,BigInt}, n::Int)
    get!(D, n) do
        n == 0 && return big(0)
        n == 1 && return big(1)
        return fib(D, n-1) + fib(D, n-2)
    end
end

myfib = Actor(fib, Dict{Int,BigInt}())

call!(myfib, 1000)


text:
include("utils.jl")

using Statistics
using Dates

include("data_graph.jl")
include("model.jl")

import CUDA
CUDA.has_cutensor()
import IterTools

include("train.jl")

using Profile
using BenchmarkTools: @btime
import CSV
using DataFrames: DataFrame, Not, All

function readdir_nohidden(dir)
    files = readdir(dir)
    return filter((x)->!startswith(x, "."), files)
end

function load_most_recent(model_dir)
    if !isdir(model_dir) return nothing, 1 end
    # get the most recent one
    files = readdir_nohidden(model_dir)
    if length(files) == 0 return nothing, 1 end

    # step-1000.bson
    steps = map(files) do fname
        step_str = match(r"step-(\d+).bson", fname).captures[1]
        parse(Int, step_str)
    end
    max_step = maximum(steps)
    most_recent = files[argmax(steps)]
    @info "Loading $most_recent .."
    @load joinpath(model_dir, most_recent) model
    return model, max_step
end

function keep_most_recent!(model_dir)
    @info "cleaning $model_dir .."
    # if there's more than one bson file, keep the most recent, and delete the rest
    if !isdir(model_dir) return nothing, 1 end
    # get the most recent one
    files = readdir_nohidden(model_dir)
    if length(files) == 0 return nothing, 1 end

    # step-1000.bson
    steps = map(files) do fname
        step_str = match(r"step-(\d+).bson", fname).captures[1]
        parse(Int, step_str)
    end
    max_step = maximum(steps)
    most_recent = files[argmax(steps)]
    @info "will keep only $most_recent"
    # delete all other files
    for f in files
        if f != most_recent
            @info "Removing $f .."
            rm(joinpath(model_dir, f))
        end
    end
end

"""These functions takes care of data creation. It is fine as long as the data
is created only once, before first use.

"""
function spec2ds(spec::DataSpec; merge=false)
    # FIXME merge is not used, only to match the signature
    batch_size = spec.bsize
    ds, test_ds = load_sup_ds(spec, batch_size)
    # DEBUG move on GPU on demand to reduce memory usage
    # ds, test_ds = (ds, test_ds) .|> CuDataSetIterator
    return ds, test_ds
end

function spec2ds(specs::Array{DataSpec, N} where N; merge=false)
    if merge
        # separate into groups with same d
        gs = IterTools.groupby((x)->x.d, specs)
        dses = map(gs) do g
            dses = map(g) do spec
                ds, test_ds = load_sup_ds(spec, spec.bsize)
            end
            ds = [ds[1] for ds in dses]
            test_ds = [ds[2] for ds in dses]
            # do the merge
            ds = merge_dses(ds)
            test_ds = merge_dses(test_ds)
            ds, test_ds
        end
    else
        # FIXME this does not seem to be used in training?
        dses = map(specs) do spec
            ds, test_ds = load_sup_ds(spec)
            # ds, test_ds = (ds, test_ds) .|> CuDataSetIterator
        end
    end
    return [ds[1] for ds in dses], [ds[2] for ds in dses]
end

function exp_train(spec, model_fn;
                   prefix,
                   train_steps=1e5, test_throttle=10,
                   merge=false)
    train_steps=Int(train_steps)

    # setting expID
    expID = "$prefix-$(dataspec_to_id(spec))"

    model_dir = joinpath("saved_models", expID)

    # FIXME load model if already trained
    most_recent_model, from_step = load_most_recent(model_dir)
#     keep_most_recent!(model_dir)
    if !isnothing(most_recent_model)
        @info "using trained model, starting at step $from_step"
        # FIXME it does not seem to be smooth at the resume point
        model = most_recent_model |> gpu
    else
        model = model_fn() |> gpu
        from_step = 1
    end

    if from_step >= train_steps
        @info "already trained"
        return expID
    end

    ds, test_ds = spec2ds(spec, merge=merge)
    @show ds
    @show test_ds
    x, y = next_batch!(test_ds) |> gpu

    @show size(x)
    @show size(y)

    @info "warming up model with x .."
    model(x)
    @info "warming up gradient .."
    @show size(x)
    gradient(()->sum(model(x)), Flux.params(model))
    @info "actual training .."

    # TODO lr decay
    opt = ADAM(1e-4)

    # when continual training, new files are created
    # CAUTION there will be some overlapping
    logger = TBLogger("tensorboard_logs/train-$expID",
                      tb_append, min_level=Logging.Info)
    test_logger = TBLogger("tensorboard_logs/test-$expID",
                           tb_append, min_level=Logging.Info)

    print_cb = create_print_cb(logger=logger)
    # do not log training loss
#     print_cb = (a,b)->nothing
    test_cb = create_test_cb(model, test_ds, "test_ds", logger=test_logger)
    save_cb = create_save_cb("saved_models/$expID", model)

    @info "training .."

    sup_train!(model, opt, ds,
               # print tensorboard logs every 1 second
               print_cb=Flux.throttle(print_cb, 10),
               # test throttle = 10
               test_cb=Flux.throttle(test_cb, test_throttle),
               # save every 60 seconds TODO use a ultra large number
               #
               # I'll save all the model during training, if that's too large,
               # I'll consider:
               #
               # 1. delete old ones, i.e. rotating
               # 2. use a larger interval
        # TODO the first save should be 60 later
               save_cb=Flux.throttle(save_cb, 600),
               from_step=from_step,
               train_steps=train_steps)

    # final save and test
    test_cb(train_steps)
    save_cb(train_steps)

    return expID
end

# this should be a map, from (expID, testID) to values
# - expID is the model name, manually specified
# - testID is the DataSpec's string representation
# Also, I should save it as a text file, and support external editing
# I should also record the time of the test
#
# so, better just use a csv file
#
# expID, testID, date, prec, recall, shd
_results = nothing
function load_results!()
    global _results
    if isfile("result.csv")
        df = CSV.read("result.csv", DataFrame)
        df.date
        _results = Dict((df.expID .=> df.testID)
                        .=>
                        zip(df.prec, df.recall, df.shd,
                            df.time, df.date))
    else
        _results = Dict()
    end
end

function test()
    df[in.(df.test_d, Ref([200, 300, 400])), All(:model, :train_d, :test_d,  :prec, :recall, :shd, :time)]
end

# However, although the date is string, it is saved without quotes. Then, the
# read CSV.read() is too smart to that it converts that into DateTime. However,
# it cannot convert DateTime to string when it saves it.
#
# FIXME this is causing ambiguous issues
Base.convert(::Type{String}, date::DateTime) = "$date"

function result2csv(res, fname)
    df = DataFrame(expID=String[], testID=String[],
                   prec=Float64[], recall=Float64[],
                   shd=Float64[], time=Float64[],
                   date=String[])
    for (key,value) in res
        push!(df, (key[1], key[2], value[1], value[2], value[3], value[4], value[5]))
    end
    # save df
    CSV.write(fname, df)
end

function save_results!()
    global _results
    result2csv(_results, "result.csv")
end

function test()
    result = Dict(("hello"=>"world")=>(0.2, 0.3, 3),
                  ("1"=>"2")=>(0.2, 0.3, 3))'
    result2csv(result, "result.csv")
end

function exp_test(expID, spec, testID=nothing; use_raw=false)
    global _results
    if isnothing(_results)
        load_results!()
    end

    if isnothing(testID)
        testID = dataspec_to_id(spec)

        if use_raw
            testID = "raw-$testID"
        end
    end

    if !haskey(_results, expID=>testID)
        model_dir = joinpath("saved_models", expID)
        model, _ = load_most_recent(model_dir)
        if isnothing(model)
            error("Cannot load model in $model_dir")
        end

        model = gpu(model)
        if use_raw
            test_ds = load_sup_ds(spec, 16, use_raw=true);
            metrics, t = sup_test_raw(model, test_ds, nbatch=16)
        else
            # FIXME the merge has to be true, otherwise batchsize is messed up, e.g. 100 for d=80
            ds, test_ds = spec2ds(spec, merge=true)
            @show ds
            @show test_ds
            # DEBUG TODO not using all data for testing
            metrics, t = sup_test(model, test_ds, nbatch=4)
        end

        # add this to result
        _results[expID=>testID] = (metrics.prec, metrics.recall,
                                   metrics.shd,
                                   t,
                                   now())

        # finally, save the results
        @info "Saving results .."
        save_results!()
    end
end

function test_profile()
    Profile.clear()
    Profile.init(n = 10^7, delay = 0.01)

    open("/tmp/prof.txt", "w") do s
        # I want to ignore all low-cost operations (e.g. <10)
        Profile.print(IOContext(s, :displaysize => (24, 200)), mincount=10)
    end

    # other visualization methods
    #
    # Profile.print()
    # import ProfileView
    # ProfileView.view()
end

function inf_one(model, x)
    # assume x is (dim, dim), no batch
    reshape(model(reshape(x, size(x)..., 1)), size(x,1), size(x,1))
end

max_stars_repo_name: lihebi/DAG-EQ

max_issues_repo_name: lihebi/DAG-EQ

max_forks_repo_name: lihebi/DAG-EQ


text:
# Задача.
# ДАНО: робот - в югозападном углу поля, на которм имеются
#       прямолинейные перегородки, не соприкасющиеся друг с другом и с внешней рамкой
# РЕЗУЛЬТАТ: функция должна вернуть число имеющихся на поле перегородок
# (в упрощенном варианте - число горизонтальных перегородок)

mutable struct  CountbordersRobot <: AbstractBorderRobot
    robot::Robot
    state::Bool
    count::Int
    border_side::HorizonSide
    CountbordersRobot(r, border_side=Nord) = new(r,0,0, border_side)
end

get(robot::CountbordersRobot) = robot.robot

function try_move!(robot::CountbordersRobot,side)
    try_move!(get(robot), side)
    if isborder(robot, robot.border_side)
        if robot.state == 0
            robot.state = 1
            robot.count += 1
        end
    else
        robot.state = 0
    end
end

get_num_borders(robot::CountbordersRobot) = robot.count
#--------------------------------------------

function count_borders!(robot::Robot)
    # робот - в юго-западном углу

    robot = CountbordersRobot(robot)
    snake(robot)
    return get_num_borders(robot) # - число горизонтальных перегородок
end

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

#ЗАМЕЧАНИЕ.
#Если определить функцию snake! так, чтобы она возвращала ссылку на робота нового типа:
function snake!(robot::AbstractRobot, (next_row_side::HorizonSide, move_side::HorizonSide)=(Nord, Ost)) # - это обобщенная функция
    # Робот - в (inverse(next_row_side), inverse(move_side)) - углу поля
    movements!(robot, move_side)
    while !isborder(robot, next_row_side)
        move!(robot,next_row_side)
        move_side = inverse(move_side)
        movements!(robot,move_side)
    end
    return robot
end

#то можно было бы определить функцию count_borders! совсем кратко:
count_borders!(robot::Robot) = get_num_borders(snake!(CountbordersRobot(robot)))


text:
# calculate the stoichiometry coefficients
function stoich_coeff(species::Vector{T}, reaction::DiffEqBiological.ReactionStruct, e_order::Vector{Int}) where T <: AlgebraSet
    coeff = 1
    for ind in eachindex(e_order)
        coeff *= DiffEqBiological.get_stoch_diff(reaction, Symbol(species[ind]))^e_order[ind]
    end
    coeff
end

# calculate the sign coefficients
function sign_coeff(a_ords::Vector{Int}, b_ords::Vector{Int})
    *((-1).^(a_ords .- b_ords)...)
end

# calculate the exponentiated mean terms
function α_coeff(raw_means::Vector{T}, n_orders::Vector{Int}, k_orders::Vector{Int}) where T <: AlgebraSet
    *(raw_means.^(n_orders .- k_orders)...)
end

# calculate the bionomial coefficients
function binomial_coeff(a_ords::Vector{Int}, b_ords::Vector{Int})
    *(binomial.(a_ords, b_ords)...)
end

# calculate the 1 / n! coefficients
function inv_factorial_calc(orders::Vector{Int})
    1 / *(factorial.(orders)...)
end

# to rewrite the mass rate calculation of DiffEqBiological, so that it can be differentiated later on.
function mass_rate_MA(reaction::DiffEqBiological.ReactionStruct; algebra::Algebra=sympy())
    #println(SymPy.Sym(string(reaction.rate_org)))
    rate = algebra.Sym(string(reaction.rate_org))
    if reaction.is_pure_mass_action
        for sub in reaction.substrates
             for i in 1:sub.stoichiometry
                rate *= (algebra.Sym(sub.reactant) + 1 - i) / i
            end
        end
    end
    rate
end

# convert raw moments to central moments
function raw_to_cen(momentstruct::MomentStruct, n_order::Vector{Int})
    eq = zero(eltype(momentstruct.species))
    for ord in Iterators.product([0:n_max for n_max in n_order]...)
        __n =collect(ord)
        eq += binomial_coeff(n_order, __n) * α_coeff(momentstruct.species, n_order, __n) * momentstruct.cen_moments[__n]
    end
    eq
end

# convert central moments to raw moments (currently unused)
function cen_to_raw(momentstruct::MomentStruct, n_order::Vector{Int})
	eq = zero(eltype(momentstruct.species))
	for ord in Iterators.product([0:n_max for n_max in n_order]...)
		__n =collect(ord)
		if sum(__n) == 1
			eq += binomial_coeff(n_order, __n) * sign_coeff(n_order, __n) * α_coeff(momentstruct.species, n_order, __n) * momentstruct.species[findall(isone,__n)[1]]
		else
			eq += binomial_coeff(n_order, __n) * sign_coeff(n_order, __n) * α_coeff(momentstruct.species, n_order, __n) * momentstruct.raw_moments[__n]
		end
	end
	eq
end

# compute taylor expansion
function taylor_expand(expr::AlgebraSet, momentstruct::MomentStruct; algebra::Algebra=sympy())
    eq = zero(algebra)
    for orders in keys(momentstruct.cen_moments)
        diff_args = reduce(append!, repeat([momentstruct.species[i]], orders[i]) for i in eachindex(orders))
        factorial_coeff = inv_factorial_calc(orders)
        eq += foldl(algebra.diff, (expr, diff_args...))*momentstruct.cen_moments[orders]*factorial_coeff
    end
    eq
end

# generate ode expressions and function from symbolic arrays
function get_ode(network::DiffEqBase.AbstractReactionNetwork, moments::MomentStruct, dμ::Vector{T}, dM::Vector{T}) where T <: AlgebraSet
    filtered_cen_moments = filter(p->sum(p.first)>1, moments.cen_moments)
    central_dict = OrderedDict{Symbol, Int}(Symbol(pair[2])=>i+moments.num_species for (i,pair) in enumerate(filtered_cen_moments))
    ode_expr = Vector{Expr}(undef, moments.num_species + length(central_dict))
    for i = 1:moments.num_species
        ex = dμ[i]
        ode_expr[i] = :(internal_var___du[$i] = $(Meta.parse(string(ex))))
    end
    for i = (moments.num_species+1):(moments.num_species+length(central_dict))
        ex = dM[i-moments.num_species]
        ode_expr[i] = :(internal_var___du[$i] = $(Meta.parse(string(ex))))
    end
    combined_dict = merge(speciesmap(network), central_dict)
    ode_func = DiffEqBiological.make_func(ode_expr, combined_dict, paramsmap(network))
    ode_expr, ODEFunction(eval(ode_func))
end


text:
export Material, Lambertian, Metal, Dielectric, scatter

function random_unit_vector(::Type{T}) where T
    # normalize(randn(Vec3{T}))
    u = T(2) * rand(T) - T(1)
    ϕ = T(2) * T(π) * rand(T)
    v = sqrt(T(1) - u^2)
    Vec3(cos(ϕ) * v, sin(ϕ) * v, u)
end

random_in_unit_sphere(::Type{T}) where {T} = cbrt(rand(T)) * random_unit_vector(T)

abstract type Material end

struct Lambertian{T} <: Material
    albedo::RGB{T}
end

function scatter(r_in::Ray{T}, rec::HitRecord{T, Lambertian{T}}) where T
    scatter_direction = rec.normal + random_unit_vector(T)
    # Catch degenerate scatter direction
    if iszero(scatter_direction)
        scatter_direction = rec.normal
    end
    return (attenuation=rec.material.albedo, scattered=Ray(rec.p, scatter_direction))
end

struct Metal{T} <: Material
    albedo::RGB{T}
    fuzz::T
end

reflect(v, n) = v - 2*(v⋅n)*n

function scatter(r_in::Ray{T}, rec::HitRecord{T, Metal{T}}) where T
    reflected = reflect(normalize(r_in.dir), rec.normal)
    scattered = Ray(rec.p, reflected + rec.material.fuzz * random_in_unit_sphere(T))
    if scattered.dir ⋅ rec.normal > 0
        return (;attenuation=rec.material.albedo, scattered)
    else
        return nothing
    end
end

struct Dielectric{T} <: Material
    ir::T  # Index of Refraction
end

function refract(uv::Vec3{T}, n, etai_over_etat) where T
    cos_theta = min(-uv ⋅ n, T(1))
    r_out_perp =  etai_over_etat * (uv + cos_theta*n)
    r_out_parallel = -sqrt(abs(T(1) - norm_sqr(r_out_perp))) * n
    return r_out_perp + r_out_parallel
end

function scatter(r_in::Ray{T}, rec::HitRecord{T, Dielectric{T}}) where T
    attenuation = RGB(T(1))
    refraction_ratio = rec.front_face ? (T(1)/rec.material.ir) : rec.material.ir
    unit_direction = normalize(r_in.dir)
    refracted = refract(unit_direction, rec.normal, refraction_ratio)
    scattered = Ray(rec.p, refracted)
    return (;attenuation, scattered)
end


text:
@testset "API keys" begin
  println(testname("API key tests"))
  print("  [  ] Make api key ")
  result = make_api_key(cl)
  @test result["ok"] == true
  println("\r  [OK] Make api key")

  api_key = result["key"]

  print("  [  ] Set permissions ")
  result = set_permissions(db; key=api_key, roles=["_reader", "_writer"])
  @test result["ok"] == true
  println("\r  [OK] Set permissions")

  print("  [  ] Set permissions with no permissions should fail ")
  @test_throws ErrorException set_permissions(db)
  println("\r  [OK] Set permissions with no permissions should fail")

  print("  [  ] Set permissions with no roles should fail ")
  @test_throws ErrorException set_permissions(db; key=api_key)
  println("\r  [OK] Set permissions with no roles should fail")

  print("  [  ] View permissions ")
  permissions = get_permissions(db)
  @test haskey(permissions, "cloudant") && haskey(permissions["cloudant"], api_key)
  println("\r  [OK] View permissions")

  print("  [  ] Modify existing permissions ")
  result = set_permissions(db, permissions; key=api_key, roles=["_reader"])
  @test result["ok"] == true
  println("\r  [OK] Modify existing permissions")

  print("  [  ] Delete API key ")
  @test delete_api_key(db, api_key)
  println("\r  [OK] Delete API key")

   print("  [  ] Delete non-existing API key ")
  @test delete_api_key(db, "NOSUCHKEY") == false
  println("\r  [OK] Delete non-existing API key")
end

max_stars_repo_name: xpqz/couchzilla

max_issues_repo_name: xpqz/couchzilla

max_forks_repo_name: xpqz/couchzilla


text:
include("langevin.jl");

using NLopt;
using Cubature;
using LinearAlgebra;

 _ALGO_DICT_ = Dict("PRAXIS" => :LN_PRAXIS, "SBPLX" => :LN_SBPLX,
                         "ISRES" => :GN_ISRES, "Newton" => :NEWTON,
                         "SLSQP" => :LD_SLSQP, "LBFGS" => :LD_LBFGS,
                         "AUGLAG" => :AUGLAG, "Anneal" => :ANNEAL);

function algostr_to_symbol(algostr::AbstractString)
  if algostr in keys(_ALGO_DICT_)
    return _ALGO_DICT_[algostr];
  else
    error("Do not understand the algorithm: \"" * algostr * "\"");
  end
  return :AJKLLKJALDSKJFAS
end

function parse_algostr(algostr::AbstractString)
  if startswith(algostr, "try")
    steps = map(strip, split(algostr[4:end], "then"));
    return map(parse_algostr, steps);
  end
  return map(algostr_to_symbol, map(strip, split(algostr, ",")));
end

function parse_max_iters(max_iters_str::AbstractString)
  if startswith(max_iters_str, "try")
    steps = map(strip, split(max_iters_str[4:end], "then"));
    return map(parse_max_iters, steps);
  end
  return map(x -> parse(Int, x), map(strip, split(max_iters_str, ",")));
end

function initial_guess(γx::Real, γz::Real)
  γ = sqrt(γx^2 + γz^2);
  return [(γ > 0) ? N*langevin_inv(γ)*csch(langevin_inv(γ))/(4*pi) : 0;
          (γz > 0) ? γz / γ * langevin_inv(γ) : 0;
          (γx > 0) ? γx / γ * langevin_inv(γ) : 0];
end

function initial_guess(γ::Real)
  return [γ > 0 ? N*langevin_inv(γ)*csch(langevin_inv(γ))/(4*pi) : 0;
          γ > 0 ? langevin_inv(γ) : 0;
          0];
end

function initial_guess(ω::Real, γs::Tuple{Real,Real})
  sqrtω = sqrt(abs(ω));
  erfsω = ω > 0 ? erf(sqrtω) : erfi(sqrtω);
  expω = exp(ω);
  sqrtpi = sqrt(pi);
  γx, γz = γs;
  return [
          N * sqrtω / (2*pi*sqrtpi*erfsω);
          (2*expω*sqrtpi*γz*ω*erfω) / (-2*sqrtω + expω*sqrtpi*erfsω);
          (4*expω*sqrtpi*γx*ω*erfω) / (2*sqrtω + expω*sqrtpi*(-1+2*ω)*erfsω)
         ];
end

function initial_guess(ωs::Tuple{Real,Real}, γ::Real)
  ωx, ωz = ωs;
  ω = ωx + ωz;
  sqrtωx = sqrt(abs(ωx));
  sqrtωz = sqrt(abs(ωz));
  sqrtω = sqrt(abs(ω));
  erfsω = ω > 0 ? erf(sqrtω) : erfi(sqrtω);
  expω = exp(ω);
  sqrtpi = sqrt(pi);
  return [
          N * sqrtω / (2*pi*sqrtpi*erfsω);
          (2*expω*sqrtpi*sqrtωz*sqrtω*erfsω) / (-2*sqrtω + expω*sqrtpi*erfsω);
          (4*expω*sqrtpi*sqrtωx*sqrtω*erfsω) / (2*sqrtω + expω*sqrtpi*(-1+2*ω)*erfsω)
         ];
end

 _RHS_ECOORD_ = Float64[N; γz*N; γx*N];

function residuals(Is::Vector)
  global _RHS_;
  map(i -> Is[i][1] - _RHS_[i], 1:3)
end

ω(θs::Vector) = ω(θs[1], θs[2]);

ρ(ϕ::Real, θ::Real, c::Real, λ::Real, α::Real) = c * exp(-ω(ϕ, θ) + λ*cos(θ)
                                                         + α*cos(ϕ)*sin(θ));
ρ(θs::Vector, xs::Vector) = ρ(θs[1], θs[2], xs[1], xs[2], xs[3]);
u(ϕ::Real, θ::Real) = kT * ω(ϕ, θ) + _u0_;
u(θs::Vector) = kT * ω(θs[1], θs[2]) + _u0_;

integrand_functions(xs::Vector) = (
  [
    (θs::Vector) -> ρ(θs, xs)*sin(θs[2]),
    (θs::Vector) -> ρ(θs, xs)*sin(θs[2])*cos(θs[2]),
    (θs::Vector) -> ρ(θs, xs)*sin(θs[2])*sin(θs[2])*cos(θs[1])
  ]
);

 _LOWER_BOUNDS_ = [0.0, 0.0];
 _UPPER_BOUNDS_ = [2*pi, pi];

function integrate_sys(xs::Vector)
  return map(integrand -> pcubature(integrand, _LOWER_BOUNDS_, _UPPER_BOUNDS_;
                                    reltol=RELTOL_INT, maxevals=MAXEVALS_INT),
             integrand_functions(xs));
end

function gradc(xs::Vector)
  global RELTOL_INT;
  global MAXEVALS_INT;

  grad = zeros(3);
  grad[1], = pcubature(θs -> ρ(θs, xs)/xs[1]*sin(θs[2]), _LOWER_BOUNDS_, _UPPER_BOUNDS_; reltol=RELTOL_INT, maxevals=MAXEVALS_INT);
  grad[2], = pcubature(θs -> ρ(θs, xs)/xs[1]*sin(θs[2])*cos(θs[2]), _LOWER_BOUNDS_, _UPPER_BOUNDS_; reltol=RELTOL_INT, maxevals=MAXEVALS_INT);
  grad[3], = pcubature(θs -> ρ(θs, xs)/xs[1]*sin(θs[2])*sin(θs[2])*cos(θs[1]), _LOWER_BOUNDS_, _UPPER_BOUNDS_; reltol=RELTOL_INT, maxevals=MAXEVALS_INT);

  return grad;
end

function gradλ(xs::Vector)
  global RELTOL_INT;
  global MAXEVALS_INT;

  grad = zeros(3);
  grad[1], = pcubature(θs -> ρ(θs, xs)*sin(θs[2])*cos(θs[2]), _LOWER_BOUNDS_, _UPPER_BOUNDS_; reltol=RELTOL_INT, maxevals=MAXEVALS_INT);
  grad[2], = pcubature(θs -> ρ(θs, xs)*sin(θs[2])*cos(θs[2])*cos(θs[2]), _LOWER_BOUNDS_, _UPPER_BOUNDS_; reltol=RELTOL_INT, maxevals=MAXEVALS_INT);
  grad[3], = pcubature(θs -> ρ(θs, xs)*sin(θs[2])*sin(θs[2])*cos(θs[1])*cos(θs[2]), _LOWER_BOUNDS_, _UPPER_BOUNDS_; reltol=RELTOL_INT, maxevals=MAXEVALS_INT);

  return grad;
end

function gradα(xs::Vector)
  global RELTOL_INT;
  global MAXEVALS_INT;

  grad = zeros(3);
  grad[1], = pcubature(θs -> ρ(θs, xs)*sin(θs[2])*sin(θs[2])*cos(θs[1]), _LOWER_BOUNDS_, _UPPER_BOUNDS_; reltol=RELTOL_INT, maxevals=MAXEVALS_INT);
  grad[2], = pcubature(θs -> ρ(θs, xs)*sin(θs[2])*cos(θs[2])*sin(θs[2])*cos(θs[1]), _LOWER_BOUNDS_, _UPPER_BOUNDS_; reltol=RELTOL_INT, maxevals=MAXEVALS_INT);
  grad[3], = pcubature(θs -> ρ(θs, xs)*sin(θs[2])*(sin(θs[2])*cos(θs[1]))^2, _LOWER_BOUNDS_, _UPPER_BOUNDS_; reltol=RELTOL_INT, maxevals=MAXEVALS_INT);

  return grad;
end

function hessian(xs::Vector)
  h = zeros(3, 3);

  h[:, 1] = gradc(xs);
  h[:, 2] = gradλ(xs);
  h[:, 3] = gradα(xs);

  return h;
end

function newton_raphson_iteration(xs::Vector, rs::Vector)
  h = hessian(xs);
  return xs - inv(h) * rs;
end

function newton_raphson(xs::Vector; reltol::Real=1e-10, abstol::Real=1e-8, max_iters::Int=100000)
  initial_Is = integrate_sys(xs);
  current_residuals = residuals(initial_Is);
  residual_norm = norm(current_residuals);
  init_residual_norm = residual_norm;

  if residual_norm < abstol
    return (xs, current_residuals, true, 0);
  end

  global print_every_percent_compl;
  iters_per_msg = ceil(Int, print_every_percent_compl * max_iters);

  for iter = 1:max_iters

    lin_alg_error_flag = false;
    try
      xs = newton_raphson_iteration(xs, current_residuals);
    catch
      lin_alg_error_flag = true;
    end

    # if Hessian is singular, or Newton iteration fails, perturb guess
    if lin_alg_error_flag
      @warn "hessian was singular; perturbing guess"
      xs += rand(3);
    end

    Is = integrate_sys(xs);
    current_residuals = residuals(Is);
    residual_norm = norm(current_residuals);

    relnorm = residual_norm / init_residual_norm;

    if residual_norm < abstol || relnorm < reltol
      return (xs, current_residuals, true, iter);
    end

    if iter % iters_per_msg == 0
      println("iter: $iter / $max_iters; Rnorm = $residual_norm");
    end

  end

  return (xs, current_residuals, false, max_iters);
end

function simulated_annealing(xs::Vector; anneal_factor::Real=1e-3,
                             init_T::Real=10.0, final_T::Real=1e-7,
                             relax_steps::Int=100,
                             init_delta::Real=0.5,
                             delta_adapt_factor::Real=1.2,
                             adapt_acc_thrshld::Real=0.7,
                             reltol::Real=1e-10, abstol::Real=1e-8,
                             max_iters::Int=1000000)

  initial_Is = integrate_sys(xs);
  current_residuals = residuals(initial_Is);
  residual_norm = norm(current_residuals);
  init_residual_norm = residual_norm;
  Is = initial_Is;

  if residual_norm < abstol
    return (xs, current_residuals, true, 0);
  end

  acc = 0;
  delta = init_delta;
  T = init_T;

  for iter = 1:max_iters

    # check to see if relaxation should occur
    if iter % relax_steps == 0
      T *= (1.0 - anneal_factor);
      if T < final_T
        break;
      end

      # adapt step size to current acceptance ratio
      if acc / relax_steps > adapt_acc_thrshld
        delta *= delta_adapt_factor;
      elseif acc / relax_steps < (1 - adapt_acc_thrshld)
        delta /= delta_adapt_factor;
      end

      # reset acceptance count
      acc = 0;
    end

    # trial move
    choice = rand(1:3);
    dx = (rand(Float64)*2 - 1) * delta;
    new_xs = copy(xs);
    new_xs[choice] += dx;

    new_Is = integrate_sys(new_xs);
    new_residuals = residuals(new_Is);

    # metropolis acceptance criteria
    if rand(Float64) < exp(-(dot(new_residuals, new_residuals) -
                             dot(current_residuals, current_residuals)) / T)

      # accept move
      acc += 1;
      copy!(xs, new_xs);
      current_residuals = new_residuals;
      Is = copy(new_Is);

      # check for convergence
      residual_norm = norm(current_residuals);
      relnorm = residual_norm / init_residual_norm;
      if residual_norm < abstol || relnorm < reltol
        return (xs, current_residuals, true, iter);
      end

    end

  end

  return (xs, current_residuals, false, iter);

end

_GRADS_ = [gradc, gradλ, gradα];

NLopt_raint_functions = map(i -> begin
  (xs::Vector, grad::Vector) -> begin
    if length(grad) > 0
      grad = gradc(xs);
    end

    return (pcubature(integrand_functions(xs)[i], _LOWER_BOUNDS_, _UPPER_BOUNDS_;
                      reltol=RELTOL_INT, maxevals=MAXEVALS_INT)[1] - _RHS_[i]);
  end
end, 1:3);

function NLopt_optimization_function!(xs::Vector, grad::Vector)

  Is = integrate_sys(xs);
  current_residuals = residuals(Is);

  if length(grad) > 0
    h = hessian(xs);
    grad = vec(transpose(current_residuals) * h);
  end

  return 0.5 * dot(current_residuals, current_residuals);
end

 _CONV_SMBLS_ = [:SUCCESS; :STOPVAL_REACHED; :FTOL_REACHED; :XTOL_REACHED];

function run_solver(algo::Symbol;
                    x0s::Vector = DEFAULT_INIT_GUESS, reltol::Real=1e-10,
                    abstol::Real=1e-8, max_iters::Int=10000)
  solver = init_solver(algo, x0s, reltol, abstol, max_iters);
  return solve(solver);
end

function run_solver(algos::Vector{Symbol}, maxs_iters::Vector{Int};
                    x0s::Vector = DEFAULT_INIT_GUESS, reltol::Real=1e-10,
                    abstol::Real=1e-8)

  if length(algos) != length(maxs_iters)
    error("Length of algorithm and maximum iteration vectors must be equal");
  end

  total_iterations = 0;
  xs, rs, converged = x0s, [NaN; NaN; NaN], false;

  for (algo, max_iters) in zip(algos, maxs_iters)
    (xs, rs, converged, iterations) = run_solver(algo; x0s=xs,
                                                        reltol=reltol,
                                                        abstol=abstol,
                                                        max_iters=max_iters);
    total_iterations += iterations;
  end

  return (xs, rs, converged, total_iterations);
end

function run_solver(algo_steps::Vector{Vector{Symbol}},
                    max_iters_by_step::Vector{Vector{Int}};
                    x0s::Vector = DEFAULT_INIT_GUESS, reltol::Real=1e-10,
                    abstol::Real=1e-6)

  if length(algo_steps) != length(max_iters_by_step)
    error("Length of algorithm and maximum iteration vectors must be equal");
  end

  total_iterations = 0;
  xs, rs, converged = x0s, [NaN; NaN; NaN], false;

  for (algo_step, max_iters_step) in zip(algo_steps, max_iters_by_step)
    (xs, rs, converged, iterations) = run_solver(algo_step, max_iters_step;
                                                        x0s=x0s,
                                                        reltol=reltol,
                                                        abstol=abstol);
    total_iterations += iterations;
    if converged
      return (xs, rs, converged, total_iterations);
    end
  end

  return (xs, rs, converged, total_iterations);
end

 _NLOPT_SOLVERS_ = [:LN_PRAXIS; :LN_SBPLX; :GN_ISRES; :LD_SLSQP;
                         :LD_LBFGS; :AUGLAG];
 _NLOPT_REQ_BOUNDS_ = [:GN_ISRES];
 _NLOPT_REQ_RAINTS_ = [:GN_ISRES; :LD_SLSQP; :AUGLAG];
 _NLOPT_DEFAULT_LOWER_BOUNDS_ = [-N; -100.0; -100.0];
 _NLOPT_DEFAULT_UPPER_BOUNDS_ = [100.0*N; 100.0; 100.0];

function init_solver(algo::Symbol, x0s::Vector, reltol::Real, abstol::Real,
                     max_iters::Int)

  global _NUMEVALS_;
  global _ITERS_PER_PRINT_;
  global print_every_percent_compl;

  _NUMEVALS_ = 0;
  _ITERS_PER_PRINT_ = ceil(Int, print_every_percent_compl * max_iters);

  if algo in _NLOPT_SOLVERS_
    opt = Opt(algo, 3);
    min_objective!(opt, NLopt_optimization_function!);
    if algo in _NLOPT_REQ_BOUNDS_
      lower_bounds!(_NLOPT_DEFAULT_LOWER_BOUNDS_);
      upper_bounds!(_NLOPT_DEFAULT_UPPER_BOUNDS_);
    end
    if algo in _NLOPT_REQ_RAINTS_
      for i=1:3
        equality_raint!(opt, NLopt_raint_functions[i]);
      end
    end
    if algo == :AUGLAG
      local_opt = init_solver(:LD_LBFGS, x0s, reltol, abstol, max_iters);
      local_optimizer!(opt, local_opt[1]);
    end
    ftol_rel!(opt, reltol);
    ftol_abs!(opt, abstol);
    maxeval!(opt, max_iters);
    return (opt, x0s);
  elseif algo == :NEWTON
    return ((xs) -> newton_raphson(xs; reltol=reltol, abstol=abstol, max_iters=max_iters),
            x0s);
  elseif algo == :ANNEAL
    return ((xs) -> simulated_annealing(xs; anneal_factor=ANNEAL_FACTOR,
                             init_T=ANNEAL_INIT_T, final_T=ANNEAL_FINAL_T,
                             relax_steps=ANNEAL_RELAX_STEPS,
                             init_delta=ANNEAL_INIT_DELTA,
                             delta_adapt_factor=ANNEAL_ADAPT_FACTOR,
                             adapt_acc_thrshld=ANNEAL_ACC_THRSHLD,
                             reltol=reltol, abstol=abstol,
                             max_iters=max_iters),
            x0s);
  else
    error("Do not understand algorithm symbol: ", algo);
  end
end

function solve(opt_pair::Tuple{NLopt.Opt, Vector})
  (minf, minx, ret) = optimize(opt_pair[1], opt_pair[2]);
  if ret == :XTOL_REACHED; @warn "Convergence given by x-tolerance"; end
  return (minx, sqrt(2 * minf), ret in _CONV_SMBLS_, _NUMEVALS_);
end

solve(newton_pair::Tuple{Function, Vector}) = newton_pair[1](newton_pair[2]);

max_stars_repo_name: grasingerm/polymer-stats

max_issues_repo_name: grasingerm/polymer-stats

max_forks_repo_name: grasingerm/polymer-stats


text:
@load_symbol _active_info :inform_active_info _libinform
@load_symbol _active_info_ensemble :inform_active_info_ensemble _libinform

activeinfo(xs :: Vector{Any}, k; b=2) = convert(Float64,NaN)

function activeinfo(xs :: Vector{Int}, k; b=max(2, maximum(xs)+1))
    if isempty(xs)
        convert(Float64, NaN)
    else
        ccall(_active_info, Float64, (Ptr{Int}, Int, Int, Int),
            xs, length(xs), b, k)
    end
end

function activeinfo(xs :: Array{Int,2}, k; b=max(2, maximum(xs)+1))
    ccall(_active_info_ensemble, Float64, (Ptr{Int}, Int, Int, Int, Int),
        xs, size(xs,2), size(xs,1), b, k)
end


text:
module IBKR

using PyCall
using Unitful
using UnitfulAssets
using Dates

using ..Emporos: CashQuantity, CashType, Position, Account
using ..Emporos: Asset, Stock, fullsymbol
using ..Emporos: @pywrap

const TWS_DEFAULT_PORT = 7497

const ib_insync = PyNULL()

include("types.jl")
include("util.jl")
include("data.jl")

function __init__()
    copy!(ib_insync, pyimport_conda("ib_insync", "ib-insync", "conda-forge"))

    # TODO
    # PyCall.pytype_mapping(ib_insync.Stock, Stock)
end

end


text:
function coupleproximal!(system::CVSystem,n::Int64)
    if (system.heart.lv.P[n+1]/mmHgToPa < system.branches.P[1][1,n+1] &&
        system.heart.av.zeta[n+1] == 0)
        AV = "closed";
    elseif (system.heart.av.zeta[n+1] >= 0 && system.heart.lv.P[n+1]/mmHgToPa >
        system.branches.P[1][1,n+1])
        AV = "opening";
    elseif system.heart.lv.P[n+1]/mmHgToPa < system.branches.P[1][1,n+1]
        AV = "closing";
    end

    if AV == "closed"
        # update left-running invariant by extrapolation
        ret1 = Float64[0];
        CVModule.rootinvariant!(ret1,system.branches.Q[1][1,n+1],system.branches.A[1][1,n+1],
            system.branches.Q[1][2,n+1],system.branches.A[1][2,n+1],system.branches.c0[1][end],
            system.branches.beta[1][end],system.solverparams.rho,system.solverparams.mu,
            system.solverparams.diffusioncoeff,system.branches.k[1],system.solverparams.h);
        system.branches.W2root = ret1[1];
        # update right-running invariant
        system.branches.W1root = -system.branches.W2root;
        # update proximal A, Q w/ invariants
        system.branches.A[1][1,n+2] = ((2*system.solverparams.rho/
            system.branches.beta[1][end])^2*(0.125*(system.branches.W1root-
            system.branches.W2root) + system.branches.c0[1][end])^4);
        system.branches.Q[1][1,n+2] = (system.branches.A[1][1,n+2]*0.5*(
            system.branches.W1root + system.branches.W2root));
        system.heart.lv.V[n+2] = (system.heart.lv.V[n+1]+
            system.solverparams.h*(system.heart.la.Q[n+1]-
            system.branches.Q[1][1,n+2]));
    else
        # ventricular elastance at next time step
        CVModule.elastancefn!(system,n+1);
        # update W2 by extrapolation
        ret1 = Float64[0];
        CVModule.rootinvariant!(ret1,system.branches.Q[1][1,n+1],system.branches.A[1][1,n+1],
            system.branches.Q[1][2,n+1],system.branches.A[1][2,n+1],system.branches.c0[1][end],
            system.branches.beta[1][end],system.solverparams.rho,system.solverparams.mu,
            system.solverparams.diffusioncoeff,system.branches.k[1],system.solverparams.h);
        system.branches.W2root = ret1[1];
        # non-dimensionalizing parameters
        Vs = 100*cm3Tom3;
        vs = system.branches.c0[1][end];
        ts = system.heart.av.leff/vs;
        if system.heart.av.zeta[n+1] == 0
            zs = 1e-5;
        else
            zs = system.heart.av.zeta[n+1];
        end
        # tolerance on dζ/dt (needs to be lower for higher HR)
        tolz = 1e-7 + (0.8 - system.heart.activation.th[end])*2.5e-5;
        # Newton iterations
        yout = zeros(5);
        iters = Int64[0];
        CVModule.newtonav!(yout,iters,Vs,vs,zs,ts,tolz,AV,system.heart.av.Kvo,system.heart.av.Kvc,
            system.heart.lv.P[n+1],system.branches.P[1][1,n+1],system.branches.Q[1][1,n+1],
            system.heart.lv.V[n+1],system.heart.av.zeta[n+1],system.branches.W1root,
            system.branches.W2root,system.heart.lv.E[n+2],system.heart.lv.V0,
            system.branches.beta[1][end],system.solverparams.rho,system.branches.c0[1][end],
            system.branches.A0[1][end],system.heart.av.Aann,system.heart.av.leff,
            system.heart.av.Ks,system.heart.activation.th[end],system.solverparams.h,CVModule.fav,CVModule.Jav,
            system.solverparams.maxiter,system.solverparams.maxval,
            system.solverparams.epsJ,system.solverparams.epsN)
        system.solverparams.totaliter += iters[1];
        system.branches.A[1][1,n+2] = yout[1];
        system.branches.Q[1][1,n+2] = yout[2];
        system.heart.lv.V[n+2] = yout[3];
        system.heart.av.zeta[n+2] = yout[4];
        system.branches.W1root = yout[5];
    end
end

max_stars_repo_name: dcanuto/CardioModeling

max_issues_repo_name: dcanuto/CardioModeling

max_forks_repo_name: dcanuto/CardioModeling


text:
module Function
using  ...Analyzer.Function
using  ...PatternStructure.Trees
using  ...PatternStructure.Checks
using  ...PatternStructure.SlurpTypes
using  ...PatternStructure.Special
using  ...Matching.Environment
using  ...Helper
export matcher

function matcher(pattern, mod=current_module())
  pattern_tree = analyze(pattern, mod).child

  return match!(ex, variables=Variables()) =
    matchtree(pattern_tree, ex, newstate(variables))
end

struct MatchState
  inslurp     :: Bool
  variables   :: Variables
  slurpranges :: SlurpRanges
end
newstate(vars)    = MatchState(false, vars, SlurpRanges())
enterslurp(state) = MatchState(true, state.variables, SlurpRanges())

struct SlurpInfo
  node        :: PatternNode
  parent      :: PatternTree
  children    :: Vector
  position    :: Int
  starting_ai :: Int
  parentstate :: MatchState
  state       :: MatchState
end

#----------------------------------------------------------------------------
# matchtree: Tries to match the expression with a pattern tree.
#----------------------------------------------------------------------------

matchtree(leaf::PatternLeaf, ex, st) = true

function matchtree(gate::PatternGate, ex, st)
  match_check(gate.check, ex, st) &&
    matchtree(gate.child, ex, st)
end

function matchtree(node::PatternNode, ex, st)
  exprhead(ex) == nodehead(node) || return false

  args = node.step(ex)
  match_children(node, args, 1, 1, st)
end

#----------------------------------------------------------------------------
# match_check: matches a pattern check with an expression.
#----------------------------------------------------------------------------

function match_check(bd::Binding, ex, st)
  st.inslurp ? true :  match_variable!(st.variables, bd.name, ex)
end

match_check(check::EqualityCheck{T}, ex::T, st) where T = check.value == ex
match_check(check::TypeCheck{T},     ex::T, st) where T = true
match_check(check::PredicateCheck,   ex,    st) = check.predicate(ex)

match_check(check, ex, st) = false

#----------------------------------------------------------------------------
# match_children: matches a vector of patterns with a vector of expressions.
#----------------------------------------------------------------------------

function match_children(node, args, ci, ai, st)
  matches = true

  for cj in ci:endof(node.children)
    matches  || return false
    matches, ai = match_child(node, args, cj, ai, false, st)
  end
  return matches && ai > length(args)
end

function match_child(node, args, ci, ai, slurpchild, st)
  is_slurp(node.children[ci]) ?
    (match_slurp(node, args, ci, ai, slurpchild, st)) :
    (match_nonslurp_child(node, args, ci, ai, st), ai+1)
  #end
end

function match_nonslurp_child(node, args, ci, ai, st)
  ai <= length(args) &&
  matchtree(node.children[ci], args[ai], st)
end

#----------------------------------------------------------------------------
# slurps: Patterns of the form (P...). Can be matched greedily or lazily.
#----------------------------------------------------------------------------

function match_slurp(node, args, ci, ai, slurpchild, st)
  slurpchild && match_directly_nested_slurp(node)

  local matches, mi

  slurp  = node.children[ci]
  mstate = st.inslurp ? st : enterslurp(st)

  if ai <= length(args)
     s_info = SlurpInfo(slurp, node, slurp.children, ci, ai, st, mstate)
     matches, mi = match_slurp_impl(slurp.head, s_info, args)
  else
    matches = true
    mi = ai
  end

  if matches
    add_slurp_range!(st.slurpranges, slurp, ai:mi-1)

    if !st.inslurp
      extract!(mstate.variables, slurp, args, ai, mi-1, mstate.slurpranges)
    end
  end

  matches, mi
end

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

function match_after_slurp(slurp, args, ai)
  match_children(slurp.parent, args, slurp.position+1, ai, slurp.parentstate)
end

function match_as_many_as_possible(slurp, args)
  mi       = slurp.starting_ai # match index
  args_len = length(args)
  matches  = true

  while matches && mi <= args_len
    oldmi = mi
    for ci in eachindex(slurp.children)
      matches || break
      matches, mi = match_child(slurp.node, args, ci, mi, true, slurp.state)
    end
    !matches && begin mi = oldmi end
  end
  return mi
end

function put_back_until_match(slurp, args, mi)
  match_size = length(slurp.children)
  matches    = false

  while slurp.starting_ai < mi && !matches
    matches =  match_after_slurp(slurp, args, mi)
    matches || begin mi -= match_size end
  end
  return mi, matches
end

function match_slurp_impl(::GenericGreedySlurp, slurp, args)
  mi = match_as_many_as_possible(slurp, args)
  mi, matches = put_back_until_match(slurp, args, mi)

  return (mi == slurp.starting_ai || matches), mi
end

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

function match_slurp_impl(::GenericLazySlurp, slurp, args)
  match_size  = length(slurp.children)
  args_length = length(args)

  mi = slurp.starting_ai
  while !(match_after_slurp(slurp, args, mi))
    for ci in eachindex(slurp.children)
      matches, mi = match_child(slurp.node, args, ci, mi, true, slurp.state)
      matches || return false, 0
    end
  end
  return true, mi
end

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

function match_slurp_impl(h::SimpleLastSlurp, slurp, args)
  return length(args)>=h.post, length(args)-h.post+1
end

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

function slurp_until(found, slurp, args, si)
  ai = slurp.starting_ai
  ci = slurp.position
  st = slurp.parentstate

  found      = filter(x->x>=ai+si-1, found)
  matches(f) = match_children(slurp.parent, args, ci+1, f[1]-si+1, st)

  while !isempty(found) &&
        !matches(found)
    found = found[2:end]
  end
  return !isempty(found), isempty(found) ? 0 : found[1]-si+1
end

function match_slurp_impl(h::SimpleGreedySlurpUntil, slurp, args)
  slurp_until(reverse(findin(args, h.until)), slurp, args, h.index)
end

function match_slurp_impl(h::SimpleLazySlurpUntil, slurp, args)
  slurp_until(findin(args, h.until), slurp, args, h.index)
end

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

match_nested_slurp(node) =
  error("Nesting slurps in another slurp's arguments isn't supported. (found in $(node))")

#----------------------------------------------------------------------------
# Utility functions
#----------------------------------------------------------------------------

const Iterable = Union{Vector, Tuple}

exprhead(ex::Expr)      = ex.head in [:kw, :(=)] ? :assign : ex.head
exprhead(ex::QuoteNode) = :quote
exprhead(ex::Iterable)  = :iterable
exprhead(ex::Any)       = :notexpr

#----------------------------------------------------------------------------
end


text:
"""
    AdaptiveParticleSwarm(n_particles = 3,
                          c1 = 2.0,
                          c2 = 2.0,
                          prob_shift = 0.25,
                          rng = Random.GLOBAL_RNG)

Instantiate an adaptive particle swarm optimization tuning strategy. A swarm is
initiated by sampling hyperparameters with their customizable priors, and new
models are generated by referencing each member's and the swarm's best models so
far.

### Supported Ranges and Discrete Hyperparameter Handling

See [`ParticleSwarm`](@ref) for more information about supported ranges and how
discrete hyperparameters are handled.

### Algorithm

Hyperparameter ranges are sampled and concatenated into position vectors for
each swarm particle. Velocity is initiated to be zeros, and in each iteration,
every particle's position is updated to approach its personal best and the
swarm's best models so far with the equations:

\$vₖ₊₁ = w⋅vₖ + c₁⋅rand()⋅(pbest - xₖ) + c₂⋅rand()⋅(gbest - xₖ)\$

\$xₖ₊₁ = xₖ + vₖ₊₁\$

Coefficients `w`, `c1`, `c2` are adaptively adjusted at each iteration by
determining the evolutionary phase of the swarm. We calculate the evolutionary
factor by comparing the mean distance from each particle to other members of the
swarm. This factor is then used to classify whether the swarm is in exploration,
exploitation, convergence, or jumping out phase and calibrate the tuning
hyperparameters accordingly. For more information, refer to "Adaptive Particle
Swarm Optimiztion" by Zhan, Zhang, Li, and Chung. Note that we omit the elitist
learning strategy in the paper.

New models are then generated for evaluation by mutating the fields of a deep
copy of `model`. If the corresponding range has a specified `scale` function,
then the transformation is applied before the hyperparameter is returned. If
`scale` is a symbol (eg, `:log`), it is ignored.
"""
mutable struct AdaptiveParticleSwarm{R<:AbstractRNG} <: AbstractParticleSwarm
    n_particles::Int
    c1::Float64
    c2::Float64
    prob_shift::Float64
    rng::R
end

# Constructor

function AdaptiveParticleSwarm(;
    n_particles=3,
    c1=2.0,
    c2=2.0,
    prob_shift=0.25,
    rng::R=Random.GLOBAL_RNG
) where {R}
    swarm = AdaptiveParticleSwarm{R}(n_particles, c1, c2, prob_shift, rng)
    message = MLJTuning.clean!(swarm)
    isempty(message) || @warn message
    return swarm
end

# Validate tuning hyperparameters

function MLJTuning.clean!(tuning::AdaptiveParticleSwarm)
    warning = ""
    if tuning.n_particles < 3
        warning *= "AdaptiveParticleSwarm requires at least 3 particles. " *
                   "Resetting n_particles=3. "
        tuning.n_particles = 3
    end
    c1, c2 = tuning.c1, tuning.c2
    if !(1.5 ≤ c1 ≤ 2.5) || !(1.5 ≤ c2 ≤ 2.5) || (c1 + c2 > 4)
        c1, c2 = _clamp_coefficients(c1, c2)
        warning *= "AdaptiveParticleSwarm requires 1.5 ≤ c1 ≤ 2.5, 1.5 ≤ c2 ≤ 2.5, and " *
                   "c1 + c2 ≤ 4. Resetting coefficients c1=$(c1), c2=$(c2). "
        tuning.c1 = c1
        tuning.c2 = c2
    end
    if !(0 ≤ tuning.prob_shift < 1)
        warning *= "AdaptiveParticleSwarm requires 0 ≤ prob_shift < 1. " *
                   "Resetting prob_shift=0.25. "
        tuning.prob_shift = 0.25
    end
    return warning
end

# Helper function to clamp swarm coefficients in the interval [1.5, 2.5] with a sum of less
# than or equal to 4

function _clamp_coefficients(c1, c2)
    c1 = min(max(c1, 1.5), 2.5)
    c2 = min(max(c2, 1.5), 2.5)
    scale = 4. / (c1 + c2)
    if scale < 1
        c1 *= scale
        c2 *= scale
    end
    return c1, c2
end

# Initial state

function MLJTuning.setup(tuning::AdaptiveParticleSwarm, model, ranges, n, verbosity)
    # state, evolutionary phase, swarm coefficients
    return (initialize(ranges, tuning), nothing, tuning.c1, tuning.c2)
end

# New models

function MLJTuning.models(
    tuning::AdaptiveParticleSwarm,
    model,
    history,
    (state, phase, c1, c2),
    n_remaining,
    verbosity
)
    n_particles = tuning.n_particles
    if !isnothing(history)
        sig = MLJTuning.signature(history[1].measure[1])
        measurements = similar(state.pbest)
        map(history[end-n_particles+1:end]) do h
            measurements[h.metadata] = sig * h.measurement[1]
        end
        pbest!(state, measurements, tuning)
        gbest!(state)
        f = _evolutionary_factor(state.X, argmin(state.pbest))
        phase = _evolutionary_phase(f, phase)
        w, c1, c2 = _adapt_parameters(tuning.rng, c1, c2, f, phase)
        move!(tuning.rng, state, w, c1, c2)
    end
    retrieve!(state, tuning)
    fields = getproperty.(state.ranges, :field)
    new_models = map(1:n_particles) do i
        clone = deepcopy(model)
        for (field, param) in zip(fields, getindex.(state.parameters, i))
            recursive_setproperty!(clone, field, param)
        end
        (clone, i)
    end
    return new_models, (state, phase, c1, c2)
end

# Helper function to calculate the evolutionary factor and phase

function _evolutionary_factor(X, gbest_i)
    n_particles = size(X, 1)
    dists = zeros(n_particles, n_particles)
    for i in 1:n_particles
        for j in i+1:n_particles
            dists[j, i] = dists[i, j] = norm(X[i, :] - X[j, :])
        end
    end
    mean_dists = sum(dists, dims=2) / (n_particles - 1)
    min_dist, max_dist = extrema(mean_dists)
    gbest_dist = mean_dists[gbest_i]
    f = (gbest_dist - min_dist) / max(max_dist - min_dist, sqrt(eps()))
    return f
end

function _evolutionary_phase(f, phase)
    # Classify evolutionary phase
    μs = [μ₁(f), μ₂(f), μ₃(f), μ₄(f)]
    if phase === nothing # first iteration
        phase = argmax(μs)
    else
        next_phase = mod1(phase + 1, 4)
        # switch to next phase if possible
        if μs[next_phase] > 0
            phase = next_phase
        # stay in current phase is possible, else pick the most likely phase
        elseif μs[phase] == 0
            phase = argmax(μs)
        end
    end
    return phase
end

# Helper functions to calculate probabilities of the four evolutionary states

μ₁(f) = f ≤ 0.4 ? 0.0         :
        f ≤ 0.6 ? 5 * f - 2   :
        f ≤ 0.7 ? 1           :
        f ≤ 0.8 ? -10 * f + 8 :
        0.0

μ₂(f) = f ≤ 0.2 ? 0.0        :
        f ≤ 0.3 ? 10 * f - 2 :
        f ≤ 0.4 ? 1.0        :
        f ≤ 0.6 ? -5 * f + 3 :
        0.0

μ₃(f) = f ≤ 0.1 ? 1.0          :
        f ≤ 0.3 ? -5 * f + 1.5 :
        0.0

μ₄(f) = f ≤ 0.7 ? 0.0         :
        f ≤ 0.9 ? 5 * f - 3.5 :
        1.0

# Adaptive control of swarm's parameters

function _adapt_parameters(rng, c1, c2, f, phase)
    w = 1.0 / (1.0 + 1.5*exp(-2.6 * f)) # update inertia
    δ = rand(rng) * 0.05 + 0.05 # coefficient acceleration
    if phase === 1 # exploration
        c1 += δ
        c2 -= δ
    elseif phase === 2 # exploitation
        δ *= 0.5
        c1 += δ
        c2 -= δ
    elseif phase === 3 # convergence
        δ *= 0.5
        c1 += δ
        c2 += δ
    else # jumping out
        c1 -= δ
        c2 += δ
    end
    c1, c2 = _clamp_coefficients(c1, c2)
    return w, c1, c2
end

function MLJTuning.tuning_report(
    tuning::AdaptiveParticleSwarm,
    history,
    (state, phase, c1, c2)
)
    fields = getproperty.(state.ranges, :field)
    scales = MLJBase.scale.(state.ranges)
    return (; plotting = MLJTuning.plotting_report(fields, scales, history))
end


text:
using AtomicLevels
using UnicodeFun
using Combinatorics: combinations
using WignerSymbols
using Test

@testset "Terms" begin
    @testset "Construction" begin
        @test T"1S" == Term(0, 0, p"even")
        @test T"1S" == Term(0, 0, 1)
        @test T"1Se" == Term(0, 0, p"even")
        @test T"1So" == Term(0, 0, p"odd")
        @test T"1So" == Term(0, 0, -1)
        @test T"2So" == Term(0, 1//2, p"odd")
        @test T"4P" == Term(1, 3//2, p"even")
        @test T"3D" == Term(2, 1, p"even")
        @test T"3Do" == Term(2, 1, p"odd")
        @test T"1[54]" == Term(54, 0, p"even")
        @test T"1[3/2]" == Term(3//2, 0, p"even")
        @test T"2[3/2]o" == Term(3//2, 1//2, p"odd")
        @test T"2Z" == Term(20, 1//2, p"even")

        for T in [T"1S", T"1Se", T"1So", T"2So", T"4P", T"3D", T"3Do",
                  T"1[54]", T"1[3/2]", T"2[3/2]o", T"2Z"]
            @test parse(Term, string(T)) == T
        end

        @test_throws DomainError Term(HalfInteger(-1//2), HalfInteger(1//2), p"even")
        @test_throws DomainError Term(3//2, -1//2, p"odd")
        @test_throws DomainError Term(-2, 1//2, 1)

        @test_throws ArgumentError parse(Term, "1[4/3]")
        @test_throws ArgumentError parse(Term, "1[43/]")
        @test_throws ArgumentError parse(Term, "1[/43/]")
        @test_throws ArgumentError parse(Term, "1[/43]")
        @test_throws ArgumentError parse(Term, "P")
        @test_throws ArgumentError parse(Term, "asdf")
    end

    @testset "Properties" begin
        @test multiplicity(T"1S") == 1
        @test multiplicity(T"1So") == 1
        @test multiplicity(T"2S") == 2
        @test multiplicity(T"2So") == 2
        @test multiplicity(T"3S") == 3
        @test multiplicity(T"3So") == 3

        @test weight(T"1S") == 1
        @test weight(T"1P") == 3
        @test weight(T"2S") == 2
        @test weight(T"2P") == 6

        @test T"1S" > T"1So"
        @test T"1So" < T"1S"
        @test T"1S" < T"2S"
        @test T"1S" < T"1P"
        @test T"1S" < T"3S"
        @test T"1P" < T"3S"
    end

    @testset "Pretty printing" begin
        map([T"1S" => "¹S",
             T"2So" => "²Sᵒ",
             T"4[3/2]" => "⁴[3/2]"]) do (t,s)
            @test "$(t)" == s
        end
    end

    @testset "Orbital terms" begin
        function test_single_orbital_terms(orb::Orbital, occ::Int, ts::Vector{Term})
            cts = sort(terms(orb, occ))
            ts = sort(ts)
            ccts = copy(cts)
            for t in cts
                if t in ts
                    @test count(e -> e == t, ccts) == count(e -> e == t, ts)
                    ccts = filter!(e -> e != t, ccts)
                    ts = filter!(e -> e != t, ts)
                end
            end
            if length(ts) != 0
                println("fail, terms: ", join(string.(cts), ", "))
                println("missing: ", join(string.(ts), ", "))
                println("should not be there: ", join(string.(ccts), ", "))
                println("==========================================================")
            end
            @test length(ts) == 0
        end

        function test_single_orbital_terms(orb::Orbital, occs::Tuple{Int,Int}, ts::Vector{Term})
            map(occs) do occ
                test_single_orbital_terms(orb, occ, ts)
            end
        end

        function get_orbital(o::AbstractString)
            n = something(findfirst(isequal(o[1]), AtomicLevels.spectroscopic), 0)
            ℓ = n - 1
            orb = Orbital(n, ℓ)
            g = 2(2ℓ + 1)
            occ = length(o)>1 ? parse(Int, o[2:end]) : 1
            if g != occ && g/2 != occ
                occ = (occ, g-occ)
            end
            orb, occ
        end

        function test_orbital_terms(o_ts::Pair{<:ST,<:ST}) where {ST<:AbstractString}
            orb,occ = get_orbital(o_ts[1])
            m = match(r"([0-9]+)([A-Z])", o_ts[2])
            L = something(findfirst(isequal(lowercase(m[2])[1]), AtomicLevels.spectroscopic), 0)-1
            S = (parse(Int, m[1])-1)//2
            test_single_orbital_terms(orb, occ, [Term(L, S, parity(orb)^occ[1])])
        end

        function test_orbital_terms(o_ts::Pair{<:ST,<:Vector{ST}}) where {ST<:AbstractString}
            orb,occ = get_orbital(o_ts[1])

            p = parity(orb)^occ[1]
            ts = o_ts[2]

            p1 = r"([0-9]+)\(((?:(?:[A-Z]|\[[0-9]+\])[0-9]*)+)\)"
            p2 = r"([0-9]+)([A-Z])"
            ts = map(ts) do t
                t = replace(t, " " => "")
                m = match(p1, t)
                if m != nothing
                    S = (parse(Int, m[1]) - 1)//2
                    map(eachmatch(r"([A-Z]|\[[0-9]+\])([0-9]*)", m[2])) do mm
                        term = parse(Term, "$(m[1])$(mm[1])")
                        nterms = isempty(mm[2]) ? 1 : parse(Int, mm[2])
                        [Term(term.L, term.S, p) for j in 1:nterms]
                    end
                else
                    m = match(p2, t)
                    term = parse(Term, t)
                    Term(term.L, term.S, p)
                end
            end
            test_single_orbital_terms(orb, occ, vcat(vcat(ts...)...))
        end

        # Table taken from Cowan, p. 110
        # Numbers following term symbols indicate the amount of times
        # different terms with the same (L,S) occur.
        test_orbital_terms("s" => "2S")
        for o in ["s2", "p6", "d10", "f14"]
            test_orbital_terms(o => "1S")
        end
        test_orbital_terms("p" => "2P")
        test_orbital_terms("p2" => ["1(SD)", "3P"])
        test_orbital_terms("p3" => ["2(PD)", "4S"])
        test_orbital_terms("d" => "2D")
        test_orbital_terms("d2" => ["1(SDG)", "3(PF)"])
        test_orbital_terms("d3" => ["2(PD2FGH)", "4(PF)"])
        test_orbital_terms("d4" => ["1(S2D2FG2I)", "3(P2DF2GH)", "5D"])
        test_orbital_terms("d5" => ["2(SPD3F2G2HI)", "4(PDFG)", "6S"])
        test_orbital_terms("f" => "2F")
        test_orbital_terms("f2" => ["1(SDGI)", "3(PFH)"])
        test_orbital_terms("f3" => ["2(PD2F2G2H2IKL)", "4(SDFGI)"])
        test_orbital_terms("f4" => ["1(S2D4FG4H2I3KL2N)", "3(P3D2F4G3H4I2K2LM)", "5(SDFGI)"])
        test_orbital_terms("f5" => ["2(P4D5F7G6H7I5K5L3M2NO)", "4(SP2D3F4G4H3I3K2LM)", "6(PFH)"])
        test_orbital_terms("f6" => ["1(S4PD6F4G8H4I7K3L4M2N2Q)", "3(P6D5F9G7H9I6K6L3M3NO)", "5(SPD3F2G3H2I2KL)", "7F"])
        test_orbital_terms("f7" => ["2(S2P5D7F10G10H9I9K7L5M4N2OQ)", "4(S2P2D6F5G7H5I5K3L3MN)", "6(PDFGHI)", "8S"])
        # # Data below are from Xu2006
        test_orbital_terms("g9" => [
            "2(S8 P19 D35 F40 G52 H54 I56 K53 L53 M44 N40 O32 Q26 R19 T15 U9 V7 W4 X2 YZ)",
            "4(S6 P16 D24 F34 G38 H40 I42 K39 L35 M32 N26 O20 Q16 R11 T7 U5 V3 WX)",
            "6(S3P3D9F8G12H10I12K9L9M6N6O3Q3RT)",
            "8(PDFGHIKL)", "10(S)"
        ])
        test_orbital_terms("h11" => [
            "2(S36 P107 D173 F233 G283 H325 I353 K370 L376 M371 N357 O335 Q307 R275 T241 U207 V173 W142 X114 Y88 Z68 [21]50 [22]36 [23]25 [24]17 [25]11 [26]7 [27]4 [28]2 [29] [30])",
            "4(S37 P89 D157 F199 G253 H277 I309 K313 L323 M308 N300 O271 Q251 R216 T190 U155 V131 W101 X81 Y59 Z45 [21]30 [22]22 [23]13 [24]9 [25]5 [26]3 [27] [28])",
            "6(S12 P35 D55 F76 G90 H101 I109 K111 L109 M105 N97 O87 Q77 R65 T53 U43 V33 W24 X18 Y12 Z8 [21]5 [22]3 [23] [24])",
            "8(S4 P4 D12 F11 G17 H15 I19 K16 L18 M14 N14 O10 Q10 R6 T6 U3 V3WX)",
            "10(PDFGHIKLMN)", "12S"
        ])

        @test_throws DomainError terms(o"2p", 7)

        @testset "Reference implementation" begin
            # This is an independent implementation that calculates all the terms (L, S) terms
            # of a given ℓ^w orbital is LS coupling. It works by considering the distribution
            # of the M_L and M_S quantum numbers of the many-particle basis states of the ℓ^w
            # orbital. It is much less performant than the implementation of terms() in AtomicLevels.
            function ls_terms_reference(ℓ::Integer, w::Integer, parity::Parity)
                ℓ >= 0 || throw(DomainError("ℓ must be non-negative"))
                w >= 1 || throw(DomainError("w must be positive"))
                # First, let's build a histogram of (M_L, M_S) values of all the product
                # basis states.
                lsbasis = [(ml, ms) for ms = HalfInteger(-1//2):HalfInteger(1//2), ml = -ℓ:ℓ]
                @assert length(lsbasis) == 2*(2ℓ+1)
                Lmax, Smax = w * ℓ, w * HalfInteger(1//2)
                NL, NS = convert(Int, 2*Lmax + 1), convert(Int, 2*Smax + 1)
                hist =  zeros(Int, NL, NS)
                for c in combinations(1:length(lsbasis), w)
                    ML, MS = 0, 0
                    for i in c
                        ML += lsbasis[i][1]
                        MS += lsbasis[i][2]
                    end
                    i = convert(Int, ML + Lmax) + 1
                    j = convert(Int, MS + Smax) + 1
                    hist[i, j] += 1
                end
                # Find the valid (L, S) terms by removing maximal rectangular bases from the
                # 2D histogram. The width and breadth of the rectangles determines the (L,S)
                # term that generates the corresponding (M_L, M_S) values.
                terms = Term[]
                Lmid, Smid = div(NL, 2) + (isodd(NL) ? 1 : 0), div(NS, 2) + (isodd(NS) ? 1 : 0)
                for i = 1:Lmid
                    while !all(hist[i,:] .== 0)
                        is = i:(NL-i+1)
                        for j = 1:Smid
                            js = j:(NS-j+1)
                            any(hist[is, js] .== 0) && continue
                            L = convert(HalfInteger, Lmax - i + 1)
                            S = convert(HalfInteger, Smax - j + 1)
                            push!(terms, Term(L, S, parity))
                            hist[is, js] .-= 1
                            break
                        end
                    end
                    @assert all(hist[NL-i+1, :] .== 0) # make sure everything is symmetric
                end
                return terms
            end

            # We can't go too high in ℓ, since ls_terms_reference becomes quite slow.
            for ℓ = 0:5, w = 1:(4ℓ+2)
                o = Orbital(:X, ℓ)
                p = parity(o)^w
                reference_terms = ls_terms_reference(ℓ, w, p)
                @test sort(terms(o, w)) == sort(reference_terms)
            end
        end
    end

    @testset "Count terms" begin
        @test count_terms(o"1s", 1, T"2S") == 1
        @test count_terms(o"1s", 2, T"1S") == 1
        @test count_terms(o"3d", 1, T"2D") == 1
        @test count_terms(o"3d", 3, T"2D") == 2
        @test count_terms(o"3d", 5, T"2D") == 3
    end
end