text:
#=--------------------------------------------------------------------

frequencyfilt - Functions for constructing image filters in the
                frequency domain.

Copyright (c) Peter Kovesi
peterkovesi.com

MIT License:

Permission is hereby granted, free of charge, to any person obtaining
a copy of this software and associated documentation files (the
"Software"), to deal in the Software without restriction, including
without limitation the rights to use, copy, modify, merge, publish,
distribute, sublicense, and/or sell copies of the Software, and to
permit persons to whom the Software is furnished to do so, subject to
the following conditions:

The above copyright notice and this permission notice shall be
included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

PK August 2015   Original porting from MATLAB to Julia
   October 2017  Updates for 0.6
   October 2018  Julia 0.7/1.0
---------------------------------------------------------------------=#

export filtergrid, filtergrids, gridangles
export cosineangularfilter, gaussianangularfilter
export lowpassfilter, highpassfilter, bandpassfilter, highboostfilter
export loggabor, monogenicfilters, packedmonogenicfilters
export perfft2
export geoseries
# export homomorphic

#--------------------------------------------------------------------
# filtergrids
"""
Generate grids for constructing frequency domain filters.
```
Usage:  (f, fx, fy) = filtergrids(rows, cols)
        (f, fx, fy) = filtergrids((rows, cols))

Arguments:  rows, cols - Size of image/filter

Returns:             f - Grid of size (rows, cols) containing frequency
                         values from 0 to 0.5,  where f =
                         sqrt(fx^2 + fy^2). The grid is quadrant
                         shifted so that 0 frequency is at f[1,1]

                fx, fy - Grids containing normalised frequency values
                         ranging from -0.5 to 0.5 in x and y directions
                         respectively. fx and fy are quadrant shifted.
```

See also: [`filtergrid`](@ref)  where you are only needing radius

"""
function filtergrids(rows::Integer, cols::Integer)

    # Set up X and Y spatial frequency matrices, fx and fy, with ranges
    # normalised to +/- 0.5 The following code adjusts things appropriately for
    # odd and even values of rows and columns so that the 0 frequency point is
    # placed appropriately.
    if isodd(cols)
        fxrange = (-(cols-1)/2:(cols-1)/2)/cols
    else
        fxrange = (-cols/2:(cols/2-1))/cols
    end

    if isodd(rows)
        fyrange = (-(rows-1)/2:(rows-1)/2)/rows
    else
        fyrange = (-rows/2:(rows/2-1))/rows
    end

    fx = [c for r = fyrange, c = fxrange]
    fy = [r for r = fyrange, c = fxrange]

    # Quadrant shift so that filters are constructed with 0 frequency at
    # the corners
    fx = ifftshift(fx)
    fy = ifftshift(fy)

    # Construct spatial frequency values in terms of normalised radius from
    # centre.
    f = sqrt.(fx.^2 .+ fy.^2)
    return f, fx, fy
end

# Tuple version
function filtergrids(sze::Tuple{Integer,Integer})
    return filtergrids(sze[1], sze[2])
end

#--------------------------------------------------------------------
# filtergrid
"""
Generate grid for constructing frequency domain filters.
```
Usage:  f = filtergrid(rows, cols)
        f = filtergrid((rows, cols))

Arguments:  rows, cols - Size of image/filter

Returns:             f - Grid of size (rows, cols) containing normalised
                         frequency values from 0 to 0.5.  Grid is quadrant
                         shifted so that 0 frequency is at f[1,1]

```
Used by [`phasecongmono`](@ref), [`phasecong3`](@ref), etc etc

See also: [`filtergrids`](@ref)   if you also want normalized frequency grids in
          the x and y directions as well.
"""
function filtergrid(rows::Integer, cols::Integer)

    # Set up X and Y spatial frequency ranges normalised to +/- 0.5
    # The following code adjusts things appropriately for odd and even
    # values of rows and columns so that the 0 frequency point is
    # placed appropriately.
    if isodd(cols)
        fxrange = (-(cols-1)/2:(cols-1)/2)/cols
    else
        fxrange = (-cols/2:(cols/2-1))/cols
    end

    if isodd(rows)
        fyrange = (-(rows-1)/2:(rows-1)/2)/rows
    else
        fyrange = (-rows/2:(rows/2-1))/rows
    end

    # Construct spatial frequency values in terms of normalised radius from
    # centre.
    f = [sqrt(fx^2 + fy^2) for fy in fyrange, fx in fxrange]
    return ifftshift(f)
end

# Tuple version
function filtergrid(sze::Tuple{Integer,Integer})
    return filtergrid(sze[1], sze[2])
end

#--------------------------------------------------------------------
# monogenicfilters
"""
Generate monogenic filter grids.

```
Usage: (H1, H2, f) = monogenicfilters(rows, cols)
       (H1, H2, f) = monogenicfilters((rows, cols))

Arguments:  rows,cols - Size of filters to generate

Returns: H1, H2 - The two monogenic filters.
              f - Frequency grid corresponding to the filters.

where:
       H1 = i*fx./f
       H2 = i*fy./f

```

Note that H1, H2, and f and quadrant shifted to that the 0 frequency
value is at coordinate [1,1].

See also: [`packedmonogenicfilters`](@ref)
"""
function monogenicfilters(rows::Integer, cols::Integer)

    (f, fx, fy) = filtergrids(rows, cols)
    f[1,1] = 1  # Set DC value to 1 to avoid divide by zero

    H1 = im.*fx./f
    H2 = im.*fy./f

    H1[1,1] = 0  # Restore 0 DC value
    H2[1,1] = 0
    f[1,1] = 0
    return H1, H2, f
end

# Tuple version
function monogenicfilters(sze::Tuple{Integer,Integer})
    return monogenicfilters(sze[1], sze[2])
end

#--------------------------------------------------------------------
# packedmonogenicfilters
"""
Monogenic filter where both filters are packed in the one Complex grid.

```
Usage: (H, f) = packedmonogenicfilters(rows, cols)
       (H, f) = packedmonogenicfilters((rows, cols))

Arguments:  rows,cols - Size of filters to generate

Returns:      H - The two monogenic filters packed into the
                  one Complex64 grid.
              f - Frequency grid corresponding to the filter.
```

The two monogenic filters are defined as

```
       H1 = i*fx./f
       H2 = i*fy./f
```

However the two filters can be packed together as a complex valued
matrix, one in the real part and one in the imaginary part.  Do this
by multiplying H2 by i and then adding it to H1.  When the convolution
is performed via the fft the real part of the result will correspond
to the convolution with H1 and the imaginary part with H2.  This
allows the two convolutions to be done as one in the frequency domain,
saving time and memory.

Note that H and f and quadrant shifted to that the 0 frequency
value is at coordinate [1,1].

See also: [`monogenicfilters`](@ref)
"""
function packedmonogenicfilters(rows::Integer, cols::Integer)
    (f, fx, fy) = filtergrids(rows, cols)
    f[1,1] = 1  # Set DC value to 1 to avoid divide by zero

    # Pack the two monogenic filters by multiplying H2 by i and then
    # adding it to H1 (note the subtraction because i*i = -1).
    H = (im.*fx .- fy)./f

    H[1,1] = 0   # Restore 0 DC value
    f[1,1] = 0
    return H, f
end

# Tuple version
function packedmonogenicfilters(sze::Tuple{Integer,Integer})
    return packedmonogenicfilters(sze[1], sze[2])
end

#--------------------------------------------------------------------
# lowpassfilter
"""
Construct a low-pass Butterworth filter.
```
Usage: f = lowpassfilter(sze, cutoff, n)

where: sze    is a two element tuple specifying the size of filter
              to construct (rows, cols).
       cutoff is the cutoff frequency of the filter 0 - 0.5
       n      is the order of the filter, the higher n is the sharper
              the transition is. (n must be an integer >= 1).
              Note that n is doubled so that it is always an even integer.

                      1
      f =    --------------------
                              2n
              1.0 + (w/cutoff)
```
The frequency origin of the returned filter is at the corners.

See also: [`highpassfilter`](@ref), [`highboostfilter`](@ref), [`bandpassfilter`](@ref)
"""
function lowpassfilter(sze::Tuple{Integer, Integer}, cutoff::Real, n::Integer)

    if cutoff < 0 || cutoff > 0.5
	error("cutoff frequency must be between 0 and 0.5")
    end

    f = filtergrid(sze)
    return 1.0 ./ (1.0 .+ (f ./ cutoff).^(2*n))
end

# Compute the low pass filter value at a specified normalised frequency
function lowpassfilter(f::Real, cutoff::Real, n::Integer)
    return 1.0 / (1.0 + (f / cutoff)^(2*n))
end
#--------------------------------------------------------------------
# bandpassfilter
"""
Construct a band-pass Butterworth filter.
```
Usage: f = bandpassfilter(sze, cutin, cutoff, n)

Arguments:
             sze - A 2 element tuple specifying the size of filter
                   to construct (rows, cols).
   cutin, cutoff - The frequencies defining the band pass 0 - 0.5
               n - The order of the filter, the higher n is the sharper
                   the transition is. (n must be an integer >= 1).
Returns:
               f - Frequency domain filter of size==sze, the frequency
                   origin is at the corners.
```
See also: [`lowpassfilter`](@ref), [`highpassfilter`](@ref), [`highboostfilter`](@ref)
"""
function bandpassfilter(sze::Tuple{Integer, Integer}, cutin::Real, cutoff::Real, n::Integer)

    if cutin < 0 || cutin > 0.5 || cutoff < 0 || cutoff > 0.5
	error("Frequencies must be between 0 and 0.5")
    end

    if n < 1
        error("Order of filter must be greater than 1")
    end

    return lowpassfilter(sze, cutoff, n) - lowpassfilter(sze, cutin, n)
end

#--------------------------------------------------------------------
# highboostfilter
"""
Construct a high-boost Butterworth filter.
```
Usage: f = highboostfilter(sze, cutoff, n, boost)

Arguments:
         sze - A 2 element tuple specifying the size of filter
               to construct (rows, cols).
      cutoff - The cutoff frequency of the filter 0 - 0.5
           n - The order of the filter, the higher n is the sharper
               the transition is. (n must be an integer >= 1).
       boost - The ratio that high frequency values are boosted
               relative to the low frequency values.  If boost is less
               than one then a 'lowboost' filter is generated
Returns:
           f - Frequency domain filter of size==sze, the frequency
               origin is at the corners.
```
See also: [`lowpassfilter`](@ref), [`highpassfilter`](@ref), [`bandpassfilter`](@ref)
"""
function highboostfilter(sze::Tuple{Integer, Integer}, cutoff::Real, n::Integer, boost::Real)

    if cutoff < 0 || cutoff > 0.5
	error("cutoff frequency must be between 0 and 0.5")
    end

    if boost >= 1     # high-boost filter
	f = (1.0 - 1.0/boost)*highpassfilter(sze, cutoff, n) .+ 1.0/boost
    else              # low-boost filter
	f = (1.0 - boost)*lowpassfilter(sze, cutoff, n) .+ boost
    end

    return f
end

#--------------------------------------------------------------------
# highpassfilter
"""
Construct a high-pass Butterworth filter.
```
Usage: f = highpassfilter(sze, cutoff, n)

         sze - A 2 element tuple specifying the size of filter
               to construct (rows, cols).
      cutoff - The cutoff frequency of the filter 0 - 0.5
           n - The order of the filter, the higher n is the sharper
               the transition is. (n must be an integer >= 1).
Returns:
           f - Frequency domain filter of size==sze, the frequency
               origin is at the corners.
```
See also: [`lowpassfilter`](@ref), [`highboostfilter`](@ref), [`bandpassfilter`](@ref)
"""
function highpassfilter(sze::Tuple{Integer, Integer}, cutoff::Real, n::Integer)

    if cutoff < 0 || cutoff > 0.5
	error("cutoff frequency must be between 0 and 0.5")
    end

    return 1.0 .- lowpassfilter(sze, cutoff, n)
end
#--------------------------------------------------------------------
# loggabor
"""
The logarithmic Gabor function in the frequency domain.

```
Usage: v = loggabor(f::Real, fo::Real, sigmaOnf::Real)

Arguments:
            f - Frequency to evaluate the function at.
           fo - Centre frequency of filter.
     sigmaOnf - Ratio of the standard deviation of the Gaussian
                describing the log Gabor filter's transfer function
                in the frequency domain to the filter center frequency.

sigmaOnf = 0.75 gives a filter bandwidth of about 1 octave.
sigmaOnf = 0.55 gives a filter bandwidth of about 2 octaves.

```
"""
function loggabor(f::Real, fo::Real, sigmaOnf::Real)
    if f < eps()
        return 0.0
    else
        return exp((-(log(f/fo))^2) / (2 * log(sigmaOnf)^2))
    end
end

#-------------------------------------------------------------
# gridangles
"""
Generate arrays of filter grid angles.
```
Usage: (sintheta, costheta) = gridangles(freq, fx, fy)

Arguments: freq, fx, fy - The output of filtergrids()

Returns:       sintheta - The sine and cosine of the angles in the filtergrid
               costheta

```
See also [`filtergrids`](@ref)

"""
function gridangles(freq::AbstractArray{T1,2},
                    fx::AbstractArray{T2,2}, fy::AbstractArray{T3,2}) where {T1 <: Real, T2 <: Real, T3 <: Real}

    freq[1,1] = 1         # Avoid divide by 0
    sintheta = fx./freq   # sine and cosine of filter grid angles
    costheta = fy./freq
    freq[1,1] = 0         # Restore 0 DC

    return sintheta, costheta
end

#--------------------------------------------------------------------
# cosineangularfilter
"""
Orientation selective frequency domain filter with cosine windowing function.

```
Usage: filter = cosineangularfilter(angl, wavelen, sintheta, costheta)

Arguments:
               angl - Orientation of the filter (radians)
            wavelen - Wavelength of the angular cosine window function.
 sintheta, costheta - Grids as returned by gridangles()
```
See also: [`gaussianangularfilter`](@ref), [`filtergrids`](@ref)
"""
function cosineangularfilter(angl::Real, wavelen::Real,
                             sintheta::Array{T1,2}, costheta::Array{T2,2}) where {T1 <: Real, T2 <: Real}

    sinangl = sin(angl); cosangl = cos(angl)
    fltr = zeros(size(sintheta))

    # For each point in the filter matrix calculate the angular
    # distance from the specified filter orientation.  To overcome
    # the angular wrap-around problem sine difference and cosine
    # difference values are first computed and then the atan2
    # function is used to determine angular distance.
    for n in eachindex(sintheta)
        ds = sintheta[n] * cosangl - costheta[n] * sinangl  # Difference in sine.
        dc = costheta[n] * cosangl + sintheta[n] * sinangl  # Difference in cosine.
        dtheta = abs(atan(ds, dc))                   # Absolute angular distance.

        # Scale theta so that cosine spread function has the right
        # wavelength and clamp to pi.  dtheta has a wavelength of
        # 2pi. If desired wavelength of cosine window function is
        # wavelen we need to multiply dtheta by 2*pi/wavelen.
        dtheta = min(dtheta*2*pi/wavelen, pi)
        # The spread function is cos(dtheta) between -pi and pi.  We add 1,
        # and then divide by 2 so that the value ranges 0-1
        fltr[n] = (cos(dtheta)+1)/2
    end

    return fltr
end

#--------------------------------------------------------------------
# gaussianangularfilter
"""
Orientation selective frequency domain filter with Gaussian windowing function.

```
Usage: filter = gaussianangularfilter(angl, thetaSigma, sintheta, costheta)

Arguments:
               angl - Orientation of the filter (radians)
         thetasigma - Standard deviation of angular Gaussian window function.
 sintheta, costheta - Grids as returned by gridangles()
```

See also: [`cosineangularfilter`](@ref), [`gridangles`](@ref), [`filtergrids`](@ref)
"""
function gaussianangularfilter(angl::Real, thetaSigma::Real,
                               sintheta::Array{T1,2}, costheta::Array{T2,2}) where {T1 <: Real, T2 <: Real}

    sinangl = sin(angl); cosangl = cos(angl)
    fltr = zeros(size(sintheta))

    # For each point in the filter matrix calculate the angular
    # distance from the specified filter orientation.  To overcome
    # the angular wrap-around problem sine difference and cosine
    # difference values are first computed and then the atan2
    # function is used to determine angular distance.
    for n in eachindex(sintheta)
        ds = sintheta[n] * cosangl - costheta[n] * sinangl  # Difference in sine.
        dc = costheta[n] * cosangl + sintheta[n] * sinangl  # Difference in cosine.
        dtheta = atan(ds, dc)                               # Angular distance.

        fltr[n] = exp((-dtheta.^2) / (2 * thetaSigma^2))
    end
    return fltr
end

#--------------------------------------------------------------------
#=
"""
homomorphic - Performs homomorphic filtering on an image.

Function performs homomorphic filtering on an image. This form of
filtering sharpens features and flattens lighting variantions in an
image.  It usually is very effective on images which have large
variations in lighting, for example when a subject appears against
strong backlighting.
```
Usage: newim =
homomorphic(inimage,boost,CutOff,order,lhistogram_cut,uhistogram_cut, hndl)
homomorphic(inimage,boost,CutOff,order,lhistogram_cut,uhistogram_cut)
homomorphic(inimage,boost,CutOff,order,hndl)
homomorphic(inimage,boost,CutOff,order)

Parameters:  (suggested values are in brackets)
        boost    - The ratio that high frequency values are boosted
                   relative to the low frequency values (2).
        CutOff   - Cutoff frequency of the filter (0 - 0.5)
        order    - Order of the modified Butterworth style filter that
                   is used, this must be an integer > 1 (2)
        lhistogram_cut - Percentage of the lower end of the filtered image's
                         histogram to be truncated, this eliminates extreme
                         values in the image from distorting the final result. (0)
        uhistogram_cut - Percentage of upper end of histogram to truncate. (5)
        hndl           - Optional handle to text box for updating
                         messages to be sent to a GUI interface.
```

If lhistogram_cut and uhistogram_cut are not specified no histogram
truncation will be applied.

Suggested values: newim = homomorphic(im, 2, .25, 2, 0, 5)
"""

# June 1999
# December 2001 cleaned up and modified to work with colour images

function him = homomorphic(im, boost, CutOff, order, varargin)

	if ndims(im) == 2  # Greyscale image
	    him = Ihomomorphic(im, boost, CutOff, order, varargin)

	else               # Assume colour image in RGB format
	    hsv = rgb2hsv(im)   # Convert to HSV and apply homomorphic
				 # filtering to just the intensity component.
            hsv(:,:,3) = Ihomomorphic(hsv(:,:,3), boost, CutOff, order, varargin)
	    him = hsv2rgb(hsv)  # Convert back to RGB
	end

#------------------------------------------------------------------------
# Internal function that does the real work
#------------------------------------------------------------------------

function him = Ihomomorphic(im, boost, CutOff, order, varargin)

    # The possible elements in varargin are:
    # {lhistogram_cut, uhistogram_cut, hndl}

    varargin = varargin{:}

    if nargin == 5
	nopparams  = length(varargin)
    end

    if (nopparams == 3)
	dispStatus = 1
	truncate = 1
	lhistogram_cut = varargin{1}
	uhistogram_cut = varargin{2}
	hndl = varargin{3}
    elseif (nopparams == 2)
	dispStatus = 0
	truncate = 1
	lhistogram_cut = varargin{1}
	uhistogram_cut = varargin{2}
    elseif (nopparams == 1)
	dispStatus = 1
	truncate = 0
	hndl = varargin{1}
    elseif (nopparams == 0)
	dispStatus = 0
	truncate = 0
    else
	disp('Usage: newim = homomorphic(inimage,LowGain,HighGain,CutOff,order,lhistogram_cut,uhistogram_cut)')
	error('or    newim = homomorphic(inimage,LowGain,HighGain,CutOff,order)')
    end

    [rows,cols] = size(im)

    im = normalise(im)                        # Rescale values 0-1 (and cast
					       # to `double' if needed).
    FFTlogIm = fft2(log(im+.01))              # Take FFT of log (with offset
                                               # to avoid log of 0).
    h = highboostfilter([rows cols], CutOff, order, boost)
    him = exp(real(ifft2(FFTlogIm.*h)))       # Apply the filter, invert
					       # fft, and invert the log.

    if truncate

	# Problem:
	# The extreme bright values in the image are exaggerated by the filtering.
	# These (now very) bright values have the overall effect of darkening the
	# whole image when we rescale values to 0-255.
	#
	# Solution:
	# Construct a histogram of the image.  Find the level below which a high
	# percentage of the image lies (say 95#).  Saturate the grey levels in
	# the image to this level.

	if dispStatus
	    set(hndl,'String','Calculating histogram and truncating...')
	    drawnow
	else
	    disp('Calculating histogram and truncating...')
	end

	him = histtruncate(him, lhistogram_cut, uhistogram_cut)

    else
	him = normalise(him)  # No truncation, but fix range 0-1
    end

end
=#

#--------------------------------------------------------------------
#  perfft2
"""
2D Fourier transform of Moisan's periodic image component.
```
Usage: (P, S, p, s) = perfft2(img)

Argument: img - Image to be transformed
Returns:    P - 2D fft of periodic image component
            S - 2D fft of smooth component
            p - Periodic component (spatial domain)
            s - Smooth component (spatial domain)
```
Moisan's "Periodic plus Smooth Image Decomposition" decomposes an image
into two components

        img = p + s

where s is the 'smooth' component with mean 0 and p is the 'periodic'
component which has no sharp discontinuities when one moves cyclically
across the image boundaries.

This decomposition is very useful when one wants to obtain an FFT of
an image with minimal artifacts introduced from the boundary
discontinuities.  The image p gathers most of the image information
but avoids periodization artifacts.

The typical use of this function is to obtain a 'periodic only' fft of an
image

      P = perfft2(img)

Displaying the amplitude spectrum of P will yield a clean spectrum without the
typical vertical-horizontal 'cross' arising from the image boundaries that you
would normally see.

Note if you are using the function to perform filtering in the frequency
domain you may want to retain s (the smooth component in the spatial domain)
and add it back to the filtered result at the end.

The computational cost of obtaining the 'periodic only' FFT involves taking an
additional FFT.
"""
function perfft2(img::Array{T,2}) where T <: Real
    #=
    Reference:
    This code is adapted from Lionel Moisan's Scilab function 'perdecomp.sci'
    "Periodic plus Smooth Image Decomposition" 07/2012 available at

    http://www.mi.parisdescartes.fr/~moisan/p+s

    Paper:
    L. Moisan, "Periodic plus Smooth Image Decomposition", Journal of
    Mathematical Imaging and Vision, vol 39:2, pp. 161-179, 2011.
    =#

    (rows,cols) = size(img)

    # Compute the boundary image which is equal to the image discontinuity
    # values across the boundaries at the edges and is 0 elsewhere
    s = zeros(rows, cols)
    s[1,:]   = img[1,:] - img[end,:]
    s[end,:] = -s[1,:]
    s[:,1]   = s[:,1]   + img[:,1] - img[:,end]
    s[:,end] = s[:,end] - img[:,1] + img[:,end]

    # Generate grid upon which to compute the filter for the boundary image in
    # the frequency domain.  Note that cos() is cyclic hence the grid values can
    # range from 0 .. 2*pi rather than 0 .. pi and then pi .. 0

    # Generate FFT of smooth component
    cxrange = 2*pi*(0:cols-1)/cols
    cyrange = 2*pi*(0:rows-1)/rows
    denom = [2*(2 - cos(cx) - cos(cy)) for cy in cyrange, cx in cxrange]
    S = fft(s)./denom

    # The [1,1] element of the filter will be 0 so S[1,1] may be Inf or NaN
    S[1,1] = 0.0         # Enforce 0 mean

    P = fft(img) .- S    # FFT of periodic component

# ** ? Perhaps have a separate version or a flag to request the
# ** spatial versions of p and s
    s = real(ifft(S))
    p = img .- s

    return P, S, p, s
end

#----------------------------------------------------------------------
# geoseries
"""
Generate geometric series.

Useful for generating geometrically scaled wavelengths for specifying
filter banks.

```
Usage 1: s = geoseries(s1, mult, n)

Arguments:      s1 - The starting value in the series.
              mult - The scaling factor between succesive values.
                 n - The desired number of elements in the series.

Usage 2: s = geoseries((s1, sn), n)

Arguments: (s1, sn) - Tuple specifying the 1st and last values
                      in the the series.
                  n - The desired number of elements in the series.
```

Example:
```
      s = geoseries(0.5, 2, 4)
      s =  [0.5000,    1.0000,    2.0000,    4.0000]
```
 Alternatively obtain the same series using
```
           s = geoseries((0.5, 4), 4)
```
"""
function geoseries(s1::Real, mult::Real, n::Integer)
    @assert n > 0  "Number of elements must be a +ve integer"
    @assert s1 > 0  "Starting value must be > 0"
    return s = s1 * mult.^(0:(n-1))
end

function geoseries(s1sn::Tuple{Real, Real}, n::Int)

    # Compute the multiplier from the desired number of elements.
    # max_val = s1*mult^(n-1)
    s1 = s1sn[1]
    sn = s1sn[2]
    @assert s1 > 0  "Starting value must be > 0"
    mult = exp(log(sn/s1)/(n-1))
    return  geoseries(s1, mult, n)
end

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

max_stars_repo_name: peterkovesi/ImagePhaseCongruency

max_issues_repo_name: peterkovesi/ImagePhaseCongruency

max_forks_repo_name: peterkovesi/ImagePhaseCongruency


text:
export foralledges

"""
    foralledges(f::Function, s::Stinger, v::Int64)

Iterates over all the edges edges of a vertex and applies a function to each
edge. The function should take 3 arguments.
`f(current_edge::StingerEdge, vertexid::Int64, etype::Int64)`
"""
function foralledges(f::Function, s::Stinger, v::Int64)
    ebpool_priv_ptr = storageptr(s) + s[ebpool_start] * (sizeof(UInt8)) + sizeof(UInt64) * 2
    current_eb_ptr = ebpool_priv_ptr + getvertex(s, v).edges * (sizeof(StingerEdgeBlock) + sizeof(StingerEdge)*NUMEDGEBLOCKS)
    while current_eb_ptr != ebpool_priv_ptr
        current_eb = unsafe_load(convert(Ptr{StingerEdgeBlock}, current_eb_ptr))
        for i=1:current_eb.high
            current_edge = unsafe_load(convert(Ptr{StingerEdge}, current_eb_ptr+sizeof(StingerEdgeBlock)), i)
            f(current_edge, current_eb.vertexid, current_eb.etype)
        end
        current_eb_ptr = ebpool_priv_ptr + current_eb.next * (sizeof(StingerEdgeBlock) + sizeof(StingerEdge)*NUMEDGEBLOCKS);
    end
end

ext: jl


text:
# Autogenerated wrapper script for libjlnode_jll for aarch64-linux-musl-cxx11
export jlnode_addon, libjlnode

using libnode_jll
JLLWrappers.@generate_wrapper_header("libjlnode")
JLLWrappers.@declare_file_product(jlnode_addon)
JLLWrappers.@declare_library_product(libjlnode, "libjlnode.so")
function __init__()
    JLLWrappers.@generate_init_header(libnode_jll)
    JLLWrappers.@init_file_product(
        jlnode_addon,
        "lib/jlnode_addon.node",
    )

    JLLWrappers.@init_library_product(
        libjlnode,
        "lib/libjlnode.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using MPI, MPIArrays

MPI.Init()
comm = MPI.COMM_WORLD
rank = MPI.Comm_rank(comm)
n_proc = MPI.Comm_size(comm)

# 1
@assert n_proc % 2 == 0
process_grid = (Int(n_proc/2), 2)
A = CyclicMPIArray(Float64, 7, 5, proc_grids=process_grid)
forlocalpart!(x -> fill!(x, rank), A)
sync(A)
if rank == 0
    show(stdout, "text/plain", A)
    println()

    println("pids")
    show(stdout, "text/plain", pids(A))
    println()

    println("blocksizes")
    show(stdout, "text/plain", blocksizes(A))
    println()
end
free(A)

A = CyclicMPIArray(Float64, 7, 9, blocksizes=(2, 2))
forlocalpart!(x -> fill!(x, rank), A)
sync(A)
if rank == 0
    show(stdout, "text/plain", A)
    println()

    println("pids")
    show(stdout, "text/plain", pids(A))
    println()

    println("blocksizes")
    show(stdout, "text/plain", blocksizes(A))
    println()
end
free(A)

MPI.Finalize()


text:
include("abstract_policy.jl")
include("q_based_policy.jl")
include("v_based_policy.jl")
include("tabular_deterministic_policy.jl")
include("tabular_random_policy.jl")
include("exploring_start_policy.jl")
include("off_policy.jl")
include("reinforce.jl")


text:
struct RadauVectorCache{n_stage,NX,NC}
    seed::NTuple{NC,ForwardDiff.Partials{NC,Float64}}
    store_float::Vector{Float64}
    store_complex::Vector{ComplexF64}
    neg_J::Matrix{Float64}
    x_dual::Vector{ForwardDiff.Dual{Nothing,Float64,NC}}
    xx_dual::Vector{ForwardDiff.Dual{Nothing,Float64,NC}}
    xx_0::Vector{Float64}
    function RadauVectorCache{n_stage_,NX,NC}() where {n_stage_,NX,NC}
        seed_  = ForwardDiff.construct_seeds(ForwardDiff.Partials{NC,Float64})
        store_float_ = Vector{Float64}(undef, NX)
        store_complex_ = Vector{ComplexF64}(undef, NX)
        neg_J_ = Matrix{Float64}(undef, NX, NX)
        x_dual_ = Vector{ForwardDiff.Dual{Nothing,Float64,NC}}(undef, NX)
        xx_dual_ = Vector{ForwardDiff.Dual{Nothing,Float64,NC}}(undef, NX)
        xx_0_ = Vector{Float64}(undef, NX)
        return new(seed_, store_float_, store_complex_, neg_J_, x_dual_, xx_dual_, xx_0_)
    end
end

struct RadauCacheTuple{n_stage,NX}
    Ew_stage::NTuple{n_stage,Vector{ComplexF64}}
    delta_Z_stage::NTuple{n_stage,Vector{ComplexF64}}
    X_stage::NTuple{n_stage,Vector{Float64}}
    F_X_stage::NTuple{n_stage,Vector{Float64}}
    inv_C_stage::NTuple{n_stage,Matrix{ComplexF64}}
    function RadauCacheTuple{n_stage_,NX}() where {n_stage_,NX}
        Ew_stage_ = Tuple(Vector{ComplexF64}(undef, NX) for _ = 1:n_stage_)
        delta_Z_stage_ = Tuple(Vector{ComplexF64}(undef, NX) for _ = 1:n_stage_)
        X_stage_ = Tuple(Vector{Float64}(undef, NX) for _ = 1:n_stage_)
        F_X_stage_ = Tuple(Vector{Float64}(undef, NX) for _ = 1:n_stage_)
        inv_C_stage_ = Tuple(Matrix{ComplexF64}(undef, NX, NX) for _ = 1:n_stage_)
        return new(Ew_stage_, delta_Z_stage_, X_stage_, F_X_stage_, inv_C_stage_)
    end
end

struct RadauTable{n_stage}
    A::Matrix{Float64}
    c::NTuple{n_stage,Float64}
    λ::Vector{ComplexF64}
    T::Matrix{ComplexF64}
    T⁻¹::Matrix{ComplexF64}
    bi::Matrix{Float64}
    b̂::Vector{Float64}
    b̂_0::Float64
    function RadauTable(n_rule::Int64)
        (0 <= n_rule) || error("radau_rule stage number needs to be positive, but is $n_rule")
        (n_rule <= 6) || error("RadauIIA rule $n_rule with $(radau_rule_to_stage(n_rule)) stages of order $(radau_rule_to_order(n_rule)) not implemented")

        A, b, c, λ, T, T⁻¹, bi, b̂, b̂_0 = load_radau_table_from_file(n_rule)

        n_stages = radau_rule_to_stage(n_rule)
        return new{n_stages}(A, c, λ, T, T⁻¹, bi, b̂, b̂_0)
    end
end

mutable struct RadauStep
    h::Float64
    h⁻¹::Float64
    tol_a::Float64
    tol_r::Float64
    tol_newton::Float64
    is_has_prev_step::Bool
    hᵏ⁻¹::Float64
    x_err_norm::Float64
    x_err_normᵏ⁺¹::Float64
    h_max::Float64
    h_min::Float64
    exit_flag::Int64
    function RadauStep(;h=1.0e-4, tol_a=1.0e-4, tol_r=1.0e-4, tol_newton=1.0e-16)
        return new(h, 1 / h, tol_a, tol_r, tol_newton, false, -9999.0, -9999.0, -9999.0, 0.01, 1.0e-8, -9999)
    end
end

mutable struct RadauRule
    s::Int64
    n_increase_cooldown::Int64
    θ::Float64
    θᵏ⁻¹::Float64
    k_iter::Int64
    k_iter_max::Int64
    Ψ_k::Float64
    max_rule::Int64
    function RadauRule(NR::Int64)
        return new(1, 10, -9999.0, -9999.0, -9999, 15, 9999.0, NR)
    end
end

mutable struct RadauDenseOutput{n_stage}
    is_has_X::Bool
    h::Float64
    s::Int64
    X_stage::NTuple{n_stage,Vector{Float64}}
    function RadauDenseOutput{n_stage}(NX) where {n_stage}
        X_stage = Tuple(Vector{Float64}(undef, NX) for _ = 1:n_stage)
        return new(false, -9999.0, -9999, X_stage)
    end
end

struct RadauIntegrator{T_object,NX,NC,NR,NSM}
    table::NTuple{NR,RadauTable}
    step::RadauStep
    rule::RadauRule
    dense::RadauDenseOutput{NSM}
    cv::RadauVectorCache{NSM,NX,NC}
    ct::RadauCacheTuple{NSM,NX}
    de_object::T_object
    function RadauIntegrator{NX,NC,NR,NSM}(tol::Float64, de_object_::T_object_) where {T_object_,NX,NC,NR,NSM}
        table_ = Tuple([RadauTable(k) for k = 1:NR])
        n_stage_max_ = radau_rule_to_stage(NR)
        cv_ = RadauVectorCache{n_stage_max_,NX,NC}()
        ct_ = RadauCacheTuple{n_stage_max_,NX}()
        radau_step = RadauStep(tol_newton=tol)
        radau_rule = RadauRule(NR)
        dense = RadauDenseOutput{n_stage_max_}(NX)
        return new{T_object_,NX,NC,NR,NSM}(table_, radau_step, radau_rule, dense, cv_, ct_, de_object_)
    end
end


text:
# Complete spec detailing how `@doc` *should* work.
module MacroSpecTests

using Base.Test

import Markdown

## test case text –––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––

TESTCASE = readall(joinpath(Pkg.dir("Docile"), "test", "test-case.md"))
TESTCASE_PARSED = Markdown.parse(TESTCASE)

# TODO: maybe add more cases

## test driver ––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––

symname(s::Symbol)   = s
symname(s::Method)   = s.func.code.name
symname(s::DataType) = s.name.name
symname(s::Function) = symname(s.env.defs)
symname(s::Module)   = symbol("$s")

function tester(entries, typ, case; verbose = false)

    if verbose
        xdump(entries)
        @show entries typ case
        @show typeof(entries[1])
    end

    obj, ent = entries[1]

    @test isa(obj, typ)
    @test isdefined(symname(obj))

    @test !isempty(ent.docs.content)
    @test isa(ent.docs.content[1], Markdown.Header{1})
    @test isa(ent.docs.content[end], Markdown.Code)

    if case == :docsmeta
        @test haskey(ent.meta, :key)
        @test length(ent.meta) == 2
    elseif case == :docs
        @test !haskey(ent.meta, :key)
        @test length(ent.meta) == 1
    elseif case == :meta
        @test haskey(ent.meta, :key)
        @test haskey(ent.meta, :file)
        @test ent.meta[:file] == "test-case.md"
        @test length(ent.meta) == 3
    else
        error("undefined case in `tester`")
    end
end

## state init –––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––

using Docile
@docstrings

OBJECT_COUNT = 0

## ––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––
info("testing __METADATA__ object")

@test isdefined(:__METADATA__)

@test isa(__METADATA__, Docile.Documentation)
@test isa(__METADATA__.entries, Dict{Any,Docile.Entry})
@test isa(__METADATA__.modname, Module)

@test isempty(__METADATA__.entries)
@test __METADATA__.modname == current_module()

@test length(__METADATA__.entries) == OBJECT_COUNT

## ––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––
info("testing methods")

@doc """
$(TESTCASE)
""" {
    :key => :value
    } ->
function docs_and_meta(x)
end
@query docs_and_meta(1)
tester((@query docs_and_meta(1);), Method, :docsmeta)

@doc """
$(TESTCASE)
""" ->
function docs_no_meta(x)
end
tester((@query docs_no_meta(1);), Method, :docs)

@doc { :file => "test-case.md", :key => :value } ->
function extern_docs_meta(x)
end
tester((@query extern_docs_meta(1);), Method, :meta)

OBJECT_COUNT += 3
@test length(__METADATA__.entries) == OBJECT_COUNT

## ––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––
info("testing single line methods")

@doc """
$(TESTCASE)
""" {
    :key => :value
    } ->
inline_docs_meta(x) = x
tester((@query inline_docs_meta(1);), Method, :docsmeta)

@doc """
$(TESTCASE)
""" ->
inline_docs_no_meta(x) = x
tester((@query inline_docs_no_meta(1);), Method, :docs)

@doc { :file => "test-case.md", :key => :value } ->
inline_extern_docs_meta(x) = x
tester((@query inline_extern_docs_meta(1);), Method, :meta)

OBJECT_COUNT += 3
@test length(__METADATA__.entries) == OBJECT_COUNT

## ––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––
info("testing generic functions")

@doc """
$(TESTCASE)
""" {
    :key => :value
    } ->
inline_docs_meta
tester(query(inline_docs_meta; all = false), Function, :docsmeta)

@doc """
$(TESTCASE)
""" ->
inline_docs_no_meta
tester(query(inline_docs_no_meta; all = false), Function, :docs)

@doc { :file => "test-case.md", :key => :value } ->
inline_extern_docs_meta
tester(query(inline_extern_docs_meta; all = false), Function, :meta)

OBJECT_COUNT += 3
@test length(__METADATA__.entries) == OBJECT_COUNT

## ––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––
info("testing types")

@doc """
$(TESTCASE)
""" {
    :key => :value
    } ->
type DocsAndMeta
end
tester(query(DocsAndMeta), DataType, :docsmeta)

@doc """
$(TESTCASE)
""" ->
type DocsNoMeta
end
tester(query(DocsNoMeta), DataType, :docs)

@doc { :file => "test-case.md", :key => :value } ->
type ExternDocsMeta
end
tester(query(ExternDocsMeta), DataType, :meta)

OBJECT_COUNT += 3
@test length(__METADATA__.entries) == OBJECT_COUNT

## ––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––
info("testing immutable types")

@doc """
$(TESTCASE)
""" {
    :key => :value
    } ->
immutable DocsAndMetaImm
end
tester(query(DocsAndMetaImm), DataType, :docsmeta)

@doc """
$(TESTCASE)
""" ->
immutable DocsNoMetaImm
end
tester(query(DocsNoMetaImm), DataType, :docs)

@doc { :file => "test-case.md", :key => :value } ->
immutable ExternDocsMetaImm
end
tester(query(ExternDocsMetaImm), DataType, :meta)

OBJECT_COUNT += 3
@test length(__METADATA__.entries) == OBJECT_COUNT

## ––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––
info("testing abstract types")

@doc """
$(TESTCASE)
""" {
    :key => :value
    } ->
abstract DocsAndMetaAbs
tester(query(DocsAndMetaAbs), DataType, :docsmeta)

@doc """
$(TESTCASE)
""" ->
abstract DocsNoMetaAbs
tester(query(DocsNoMetaAbs), DataType, :docs)

@doc { :file => "test-case.md", :key => :value } ->
abstract ExternDocsMetaAbs
tester(query(ExternDocsMetaAbs), DataType, :meta)

OBJECT_COUNT += 3
@test length(__METADATA__.entries) == OBJECT_COUNT

## ––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––
info("testing constants")

@doc """
$(TESTCASE)
""" {
    :key => :value
    } ->
const DocsAndMetaConst = 1
tester((@query DocsAndMetaConst;), Symbol, :docsmeta)

@doc """
$(TESTCASE)
""" ->
const DocsNoMetaConst = 1
tester((@query DocsNoMetaConst;), Symbol, :docs)

@doc { :file => "test-case.md", :key => :value } ->
const ExternDocsMetaConst = 1
tester((@query ExternDocsMetaConst;),Symbol, :meta)

OBJECT_COUNT += 3
@test length(__METADATA__.entries) == OBJECT_COUNT

## ––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––
info("testing globals")

@doc """
$(TESTCASE)
""" {
    :key => :value
    } ->
DocsAndMetaGlobal = 1
tester((@query DocsAndMetaGlobal;), Symbol, :docsmeta)

@doc """
$(TESTCASE)
""" ->
DocsNoMetaGlobal = 1
tester((@query DocsNoMetaGlobal;), Symbol, :docs)

@doc { :file => "test-case.md", :key => :value } ->
ExternDocsMetaGlobal = 1
tester((@query ExternDocsMetaGlobal;), Symbol, :meta)

OBJECT_COUNT += 3
@test length(__METADATA__.entries) == OBJECT_COUNT

## ––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––
info("testing macros")

@doc """
$(TESTCASE)
""" {
    :key => :value
    } ->
macro docs_and_meta(x)
end
tester((@query @docs_and_meta;), Symbol, :docsmeta)

@doc """
$(TESTCASE)
""" ->
macro docs_no_meta(x)
end
tester((@query @docs_no_meta;), Symbol, :docs)

@doc { :file => "test-case.md", :key => :value } ->
macro extern_docs_meta(x)
end
tester((@query @extern_docs_meta;), Symbol, :meta)

OBJECT_COUNT += 3
@test length(__METADATA__.entries) == OBJECT_COUNT

## ––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––
info("testing modules")

@doc """
$(TESTCASE)
""" {
    :key => :value
    } ->
MacroSpecTests
tester(query(MacroSpecTests), Module, :docsmeta)

@doc """
$(TESTCASE)
""" ->
MacroSpecTests
tester(query(MacroSpecTests), Module, :docs)

@doc { :file => "test-case.md", :key => :value } ->
MacroSpecTests
tester(query(MacroSpecTests), Module, :meta)

OBJECT_COUNT += 1
@test length(__METADATA__.entries) == OBJECT_COUNT

## ––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––
info("testing doctest")

doctest(MacroSpecTests)
doctest(MacroSpecTests; verbose = true)

end


text:
# SYNTAX TEST "Packages/Julia/Julia.sublime-syntax"

# For information on how this file is used, see
# https://www.sublimetext.com/docs/3/syntax.html#testing
# Run tests by pressing `ctrl(or cmd)+shift+b`, i.e. run the `build` command

#=
  command block
# ^^^^^^^^^^^^^ comment.block.number-equal-sign
=#
# <- comment.block.number-equal-sign punctuation.definition.comment.number-equal-sign

##
## NUMBERS ====
##
  0b101
# ^^^^^ constant.numeric
  0o7
# ^^^ constant.numeric
  0xa3
# ^^^^ constant.numeric
  1e+123
# ^^^^^^ constant.numeric
  12e123
# ^^^^^^ constant.numeric
  1.32e+123
# ^^^^^^^^^ constant.numeric
  .32e+123
# ^^^^^^^^ constant.numeric
  1.e-123
# ^^^^^^^ constant.numeric
  11
# ^^ constant.numeric
  .11
# ^^^ constant.numeric
  11.
# ^^^ constant.numeric
  11.11
# ^^^^^ constant.numeric
  2.a
# ^^ constant.numeric
#   ^ meta.generic-name
# (issue 37)
  123_4_56_7
# ^^^^^^^^^^ constant.numeric
  0xa_3_f
# ^^^^^^^ constant.numeric
  0b1_0_1
# ^^^^^^^ constant.numeric
  1.3_2e+1_2_3
# ^^^^^^^^^^^^ constant.numeric.julia
# (issue 51)
  e2
# ^^ meta.generic-name
  2e
# ^ constant.numeric
#  ^ meta.generic-name
  2e2
# ^^^ constant.numeric
  e+2
# ^ meta.generic-name
#  ^ keyword.operator
#   ^ constant.numeric
  2+e
# ^ constant.numeric
#  ^ keyword.operator
#   ^ meta.generic-name

##
## CONSTANTS ====
##

  true
# ^^^^ constant.language
  false
# ^^^^^ constant.language
  nothing
# ^^^^^^^ constant.language
  missing
# ^^^^^^^ constant.language
  undef
# ^^^^^ constant.language
  NaN
# ^^^ constant.language
  NaN16
# ^^^^^ constant.language
  NaN32
# ^^^^^ constant.language
  NaN64
# ^^^^^ constant.language
  Inf
# ^^^ constant.language
  Inf16
# ^^^^^ constant.language
  Inf32
# ^^^^^ constant.language
  Inf64
# ^^^^^ constant.language
  ℯ
# ^ constant.language
  pi
# ^^ constant.language
  π
# ^ constant.language
  im
# ^^ constant.language
  ARGS
# ^^^^ constant.language
  C_NULL
# ^^^^^^ constant.language
  ENDIAN_BOM
# ^^^^^^^^^^ constant.language
  ENV
# ^^^ constant.language
  LOAD_PATH
# ^^^^^^^^^ constant.language
  PROGRAM_FILE
# ^^^^^^^^^^^^ constant.language
  STDERR
# ^^^^^^ constant.language
  STDIN
# ^^^^^ constant.language
  STDOUT
# ^^^^^^ constant.language
  VERSION
# ^^^^^^^ constant.language

##
## Variables and Functions ====
##

  bar
# ^^^ meta.generic-name.julia

  foo.bar
# ^^^ meta.generic-name
#    ^ punctuation.accessor.dot
#     ^^^ meta.generic-name

  foo.bar = 1
# ^^^ meta.generic-name.julia
#    ^ punctuation.accessor.dot.julia
#     ^^^ meta.generic-name.julia
#         ^ keyword.operator.assignment.julia

  foo.bar(x)
# ^^^ meta.function-call.julia variable.function.julia meta.generic-name.julia
#    ^ meta.function-call.julia punctuation.accessor.dot.julia
#     ^^^ meta.function-call.julia variable.function.julia meta.generic-name.julia
#         ^ meta.function-call.julia meta.function-call.arguments.julia meta.generic-name.julia
  (foo).bar(x)
#  ^^^ meta.group.julia meta.generic-name.julia
#      ^ punctuation.accessor.dot.julia
#       ^^^ meta.function-call.julia variable.function.julia meta.generic-name.julia

  bar(x = 1, y = 2, kwargs...)
# ^^^ meta.function-call variable.function
#    ^^^^^^^^^^^^^^^^^^^^^^^^^ meta.function-call.julia
#     ^ meta.function-call.arguments.julia variable.parameter.julia
#       ^ meta.function-call.arguments.julia keyword.operator.assignment.julia
#            ^ meta.function-call.arguments.julia variable.parameter.julia
#              ^ meta.function-call.arguments.julia keyword.operator.assignment.julia
#                ^ meta.function-call.arguments.julia constant.numeric.julia
#                   ^^^^^^ meta.function-call.arguments.julia meta.generic-name.julia
#                         ^^^ meta.function-call.arguments.julia keyword.operator.splat.julia
  bar(x == 1)
# ^^^ meta.function-call.julia variable.function.julia meta.generic-name.julia
#     ^ meta.function-call.arguments.julia meta.generic-name.julia
#       ^^ meta.function-call.arguments.julia keyword.operator.julia
#          ^ meta.function-call.arguments.julia constant.numeric.julia

  (bar())
# ^^^^^^^ meta.group.julia
#  ^^^ meta.function-call.julia variable.function.julia
#     ^^ meta.function-call.julia punctuation.section.arguments

  bar.()

  (bar).()
#  ^^^ meta.generic-name.julia
#      ^ keyword.operator.broadcast.julia
#       ^^ meta.function-call.julia punctuation.section.arguments

  (bar)(x, y) = 1
#  ^^^ meta.group.julia meta.generic-name.julia
#       ^^^^ meta.function.inline.julia meta.function.parameters.julia

  bar(x, y) = 2
# ^^^ meta.function.inline.julia entity.name.function.julia meta.generic-name.julia
#     ^ meta.function.inline.julia meta.function.parameters.julia variable.parameter.julia
#        ^ meta.function.inline.julia meta.function.parameters.julia variable.parameter.julia

  bar(a = 1, z = abc) = 1
#     ^ meta.function.inline.julia meta.function.parameters.julia variable.parameter.julia
#                ^^^ meta.function.inline.julia meta.function.parameters.julia meta.generic-name.julia
# ^^^ meta.function.inline.julia entity.name.function.julia
  bar(a = 1) == 1
# ^^^ meta.function-call.julia variable.function.julia meta.generic-name.julia
#            ^^ keyword.operator.julia

  foo.bar(x) = 1
# ^^^^^^^ entity.name.function.julia
#    ^ punctuation.accessor.dot.julia
#         ^ meta.function.inline.julia meta.function.parameters.julia variable.parameter.julia

  function foo(x, y)
#          ^^^ meta.function.julia entity.name.function.julia meta.generic-name.julia
#              ^ meta.function.julia meta.function.parameters.julia variable.parameter.julia
#               ^ meta.function.julia meta.function.parameters.julia punctuation.separator.parameters.julia
#                 ^ meta.function.julia meta.function.parameters.julia variable.parameter.julia
  end

  function bar.foo(x, y)
#          ^^^ meta.function.julia entity.name.function.julia meta.generic-name.julia
#             ^ meta.function.julia entity.name.function.julia punctuation.accessor.dot.julia
#              ^^^ meta.function.julia entity.name.function.julia meta.generic-name.julia
  end

  function (foo)(x, y)
#           ^^^ meta.function.julia meta.group.julia entity.name.function.julia meta.generic-name.julia
#                ^^^^ meta.function.julia meta.function.parameters.julia
  end

(
# <- meta.group
  function foo
#          ^^^ meta.function.julia entity.name.function.julia meta.generic-name.julia
  end
#
# <- meta.group
)

(
# <- meta.group

  function (x, y)
#           ^^^^ meta.function.julia meta.function.parameters.julia
  end
#
# <- meta.group
)

  function (f::Mytype)(x, y)
#           ^ meta.function.julia meta.generic-name.julia
#                      ^^^^ meta.function.julia meta.function.parameters.julia
  end

  (f::Mytype)(x, y)
# ^^^^^^^^^^^^^^^^^ meta.function-call.julia
#             ^^^^ meta.function-call.julia meta.function-call.arguments.julia
  (f::Mytype)(x, y) = 1
# ^^^^^^^^^^^^^^^^^ meta.function.inline.julia
#             ^^^^ meta.function.inline.julia meta.function.parameters.julia

  bar().(x = 1, y = 2)
#      ^ keyword.operator.broadcast.julia
#        ^^^^^^^^^^^^ meta.function-call.julia meta.function-call.arguments.julia

  filter!()
# ^^^^^^^ meta.function-call.julia variable.function.julia support.function.julia
  length([1, 2])
# ^^^^^^ meta.function-call.julia variable.function.julia support.function.julia
  length(x::Mytype) = 1
# ^^^^^^ meta.function.inline.julia entity.name.function.julia
  Base.filter!()
# ^^^^ meta.function-call.julia variable.function.julia support.module.julia
#     ^ meta.function-call.julia punctuation.accessor.dot.julia
#      ^^^^^^^ meta.function-call.julia variable.function.julia support.function.julia
  Base.length([1, 2])
# ^^^^ meta.function-call.julia variable.function.julia support.module.julia
#     ^ meta.function-call.julia punctuation.accessor.dot.julia
#      ^^^^^^ meta.function-call.julia variable.function.julia support.function.julia
  Base.length(x::Mytype) = 1
# ^^^^ meta.function.inline.julia entity.name.function.julia
#     ^ meta.function.inline.julia entity.name.function.julia punctuation.accessor.dot.julia
#      ^^^^^^ meta.function.inline.julia entity.name.function.julia

##
## Lambda functions ====
##

  x -> x^2
# ^ meta.function.lambda.julia meta.function.parameters.julia variable.parameter.julia
#   ^^ meta.function.lambda.julia keyword.operator.arrow.julia
  x::Int -> x^2
# ^ meta.function.lambda.julia meta.function.parameters.julia variable.parameter.julia
#        ^^ meta.function.lambda.julia keyword.operator.arrow.julia
  () -> 3^2
# ^^ meta.function.lambda.julia punctuation.section.parameters
#    ^^ meta.function.lambda.julia keyword.operator.arrow.julia
  (x::Int, y) -> x + y
# ^^^^^^^^^^ meta.function.lambda.julia
#     ^^^ meta.function.lambda.julia meta.function.parameters.julia support.type.julia
#             ^^ meta.function.lambda.julia keyword.operator.arrow.julia

  (::Int) -> 1
# ^^^^^^^^^^ meta.function.lambda.julia

  ( (x=1) -> x )
# ^ meta.group.julia
#  ^ - meta.funciton.lambda.julia

##
## UNICODE WORD BOUDARIES ====
##

# Unicode and numbers in names (issue 18)
  β1 = 5
# ^^ meta.generic-name
  β3(x)
# ^^ variable.function
  β2(x) = x
# ^^ entity.name.function
  ∇1 = 5
# ^^ meta.generic-name
  ∇3(x)
# ^^ variable.function
  ∇2(∇2) = x
# ^^ entity.name.function
#    ^^^ meta.function meta.function.parameters

##
## RANGES ====
##

# (issue 14)
  a:b
# ^ meta.generic-name
#  ^ keyword.operator
#   ^ meta.generic-name
  1.:a
# ^^ constant.numeric
#   ^ keyword.operator
#    ^ meta.generic-name
  a:2.
# ^ meta.generic-name
#  ^ keyword.operator
#   ^^ constant.numeric
  23.:31.
# ^^^ constant.numeric
#    ^ keyword.operator
#     ^^^ constant.numeric
  β:f()
# ^ meta.generic-name
#  ^ keyword.operator
#   ^ variable.function
  f():1
# ^ variable.function
#    ^ keyword.operator
#     ^ constant.numeric

##
## TERNARY OPERATORS ====
##

  a ? b :c
# ^ meta.generic-name
#   ^ keyword.operator
#     ^ meta.generic-name
#       ^ keyword.operator
#        ^ meta.generic-name

  a? b : c
#  ^ invalid.operator.julia

  a ?b : c
#   ^ invalid.operator.julia

##
## Types ====
##

  Array{Int, 2}
# ^^^^^ meta.parametric-type.julia support.type.julia
#      ^ meta.parametric-type.julia punctuation.section.parameter.begin.julia
#       ^^^ meta.parametric-type.julia meta.parametric-type.parameters.julia support.type.julia
  Foo{Bar, 2}
# ^^^ meta.parametric-type.julia meta.generic-name.julia
#    ^ meta.parametric-type.julia punctuation.section.parameter.begin.julia
#     ^^^ meta.parametric-type.julia meta.parametric-type.parameters.julia meta.generic-name.julia
  y::MyArray{T, 2} where T
# ^ meta.generic-name.julia
#  ^^ keyword.operator.colons.julia
#    ^^^^^^^ meta.parametric-type.julia meta.generic-name.julia
#                  ^^^^^ meta.where-clause.julia keyword.control.julia
#                        ^ meta.where-clause.julia meta.generic-name.julia
#                          ^ - meta.where-clause.julia
  MyArray{T, 2} where T >: Int
# ^^^^^^^ meta.parametric-type.julia meta.generic-name.julia
#               ^^^^^ meta.where-clause.julia keyword.control.julia
#                     ^ meta.where-clause.julia meta.generic-name.julia
#                       ^^ meta.where-clause.julia keyword.operator.superset.julia
#                          ^^^ meta.where-clause.julia support.type.julia
  MyArray{T, 1} where T <: Array{S} where S
#                     ^ meta.where-clause.julia meta.generic-name.julia
#                       ^^ meta.where-clause.julia keyword.operator.subset.julia
#                          ^^^^^ meta.where-clause.julia meta.parametric-type.julia support.type.julia
#                               ^ meta.where-clause.julia meta.parametric-type.julia punctuation.section.parameter.begin.julia
#                                ^ meta.where-clause.julia meta.parametric-type.julia meta.parametric-type.parameters.julia meta.generic-name.julia
#                                 ^ meta.where-clause.julia meta.parametric-type.julia meta.parametric-type.parameters.julia punctuation.section.parameter.end.julia
#                                   ^^^^^ meta.where-clause.julia meta.where-clause.julia keyword.control.julia
#                                         ^ meta.where-clause.julia meta.where-clause.julia meta.generic-name.julia
  MyArray{T, S} where {T, S}
#                     ^ meta.where-clause.julia meta.parametric-type.julia punctuation.section.parameter.begin.julia
#                      ^ meta.where-clause.julia meta.parametric-type.julia meta.parametric-type.parameters.julia meta.generic-name.julia
#                         ^ meta.where-clause.julia meta.parametric-type.julia meta.parametric-type.parameters.julia meta.generic-name.julia
  MyArray{T, 1} where Int <: T <: Number
#                     ^^^ meta.where-clause.julia support.type.julia
#                         ^^ meta.where-clause.julia keyword.operator.subset.julia
#                            ^ meta.where-clause.julia meta.generic-name.julia
#                              ^^ meta.where-clause.julia keyword.operator.subset.julia
#                                 ^^^^^^ meta.where-clause.julia support.type.julia

  x::Int = 1
# ^ meta.generic-name.julia
#  ^^ keyword.operator.colons.julia
#    ^^^ support.type.julia
  x::Array{S, 2} where S = [2 2; 1 1]
#    ^^^^^ meta.parametric-type.julia support.type.julia
#                ^^^^^ meta.where-clause.julia keyword.control.julia
#                      ^ meta.where-clause.julia meta.generic-name.julia
#                        ^ keyword.operator.assignment.julia

  foo(1::Int, abc::Mytype)
# ^^^ meta.function-call.julia variable.function.julia meta.generic-name.julia
#     ^^^^^^^^^^^^^^^^^^^^ meta.function-call.julia meta.function-call.arguments.julia
  foo(x = 1::Int, y = abc::Mytype)
# ^^^ meta.function-call.julia variable.function.julia meta.generic-name.julia
#         ^ meta.function-call.julia meta.function-call.arguments.julia constant.numeric.julia
#                     ^^^ meta.function-call.julia meta.function-call.arguments.julia meta.generic-name.julia
  (foo).(1::Int, abc::Mytype)
# ^ meta.group.julia punctuation.section.group.begin.julia
#  ^^^ meta.group.julia meta.function-call.julia variable.function.julia meta.generic-name.julia
#     ^ meta.group.julia punctuation.section.group.end.julia
#      ^ keyword.operator.broadcast.julia
#       ^^^^^^^^^^^^^^^^^^^^^ meta.function-call.julia
  foo(x::Int, y::Mytype) = x
# ^^^ meta.function.inline.julia entity.name.function.julia meta.generic-name.julia
#     ^^^^^^^^^^^^^^^^^^ meta.function.inline.julia meta.function.parameters.julia
#                        ^ keyword.operator.assignment.julia
  foo(::Int, ::Mytype) = x
# ^^^ meta.function.inline.julia entity.name.function.julia meta.generic-name.julia
#     ^^ meta.function.inline.julia meta.function.parameters.julia keyword.operator.colons.julia
  foo(::Int = 1, ::Mytype = 1.0) = x
# ^^^ meta.function.inline.julia entity.name.function.julia meta.generic-name.julia
#     ^^ meta.function.inline.julia meta.function.parameters.julia keyword.operator.colons.julia
#           ^ meta.function.inline.julia meta.function.parameters.julia keyword.operator.assignment.julia
#             ^ meta.function.inline.julia meta.function.parameters.julia constant.numeric.julia
#              ^ meta.function.inline.julia meta.function.parameters.julia punctuation.separator.arguments.julia
#                  ^^^^^^ meta.function.inline.julia meta.function.parameters.julia meta.generic-name.julia
#                         ^ meta.function.inline.julia meta.function.parameters.julia keyword.operator.assignment.julia

(
# <- meta.group
  (foo::Int)(x::Int = 1, y::Mytype = 1.0) = x
#  ^^^ meta.function.inline.julia meta.generic-name.julia
#       ^^^ meta.function.inline.julia support.type.julia
#            ^^^^^^^^^^^^^^^^^^^^^^^^^^^^ meta.function.inline.julia meta.function.parameters.julia
#                     ^ meta.function.inline.julia meta.function.parameters.julia constant.numeric.julia
#
# <- meta.group
)
##
## Macros ====
##

  @eval hello()
# ^ punctuation.definition.macro.julia
#  ^^^^ variable.macro.julia support.function.macro.julia

  macro foo(x, y)
#       ^^^ meta.macro.julia entity.name.macro.julia meta.generic-name.julia
#           ^^^^ meta.macro.julia meta.macro.parameters.julia
  end

  macro (foo)(x, y)
#        ^^^ meta.macro.julia entity.name.macro.julia meta.generic-name.julia
#             ^^^^ meta.macro.julia meta.macro.parameters.julia
  end

  @foo x y
# ^ punctuation.definition.macro.julia
#  ^^^ variable.macro.julia meta.generic-name.julia
  @foo(x, y)
# ^ punctuation.definition.macro.julia
#  ^^^ variable.macro.julia meta.generic-name.julia

##
## Quotes ====
##

quote
  foo
# ^^^ meta.quote.julia meta.generic-name.julia
  Array
# ^^^^^ meta.quote.julia support.type.julia
  Base
# ^^^^ meta.quote.julia support.module.julia
  Foo{}
# ^^^^^ meta.quote.julia meta.parametric-type.julia
  function $(abc)(x)
    # <- - meta.function
    #        ^^^ meta.interpolation.julia meta.generic-name.julia
    #             ^ meta.quote.julia meta.generic-name.julia
  end
#
# <- meta.quote.julia
end

##
## Symbols ====
##

  foo(a = 1, :abc => foo())
# ^^^ meta.function-call.julia variable.function.julia meta.generic-name.julia
#            ^^^^ meta.function-call.julia meta.function-call.arguments.julia constant.other.symbol.julia
#                 ^^ meta.function-call.julia meta.function-call.arguments.julia keyword.operator.pair.julia

  :a.b
# ^ keyword.operator
#  ^ constant.other.symbol
#    ^ meta.generic-name
# (issue 3)
  ,:βa
#  ^ keyword.operator
#   ^^ constant.other.symbol
  [:+]
#  ^ keyword.operator
#   ^ constant.other.symbol
  (:∘)
#  ^ keyword.operator
#   ^ constant.other.symbol
  :a!
# ^ keyword.operator
#  ^^ constant.other.symbol
  :(:)
# ^ keyword.operator
#  ^ source.julia
#    ^ source.julia
  :(a)
# ^ keyword.operator
#  ^ source.julia
#   ^ meta.generic-name
#      ^ source.julia
  :++a
# ^ keyword.operator
#  ^^ constant.other.symbol
#    ^ meta.generic-name
  :+a
# ^ keyword.operator
#  ^ constant.other.symbol
#   ^ meta.generic-name
  :∘+a # Yes, this is correct, equivalent to +(:∘, a)
# ^ keyword.operator
#  ^ constant.other.symbol
#   ^ keyword.operator
#    ^ meta.generic-name
  :.///a
# ^ keyword.operator
#  ^^^ constant.other.symbol
#     ^ keyword.operator
#      ^ meta.generic-name
# (issue 43)
  :function
# ^ keyword.operator
#  ^^^^^^^^ constant.other.symbol

##
## Chars ====
##

  s = '\''
#      ^^ meta.string.julia string.quoted.single.julia constant.character.escape.julia
  s = 'ab'
#      ^^ meta.string.julia string.quoted.single.julia invalid.string.julia

##
## Strings ====
##

  s = "\U03b1\U03b2"
#      ^^^^^^^^^^^^ meta.string.julia string.quoted.double.julia constant.character.escape.unicode.16-bit-hex.julia
  s =  "foo + $apple"
#             ^ meta.string.julia meta.interpolation.julia keyword.operator.interpolation.julia
#              ^^^^^ meta.string.julia meta.interpolation.julia meta.generic-name.julia
  s = """
  apple
  orange
  banana
# ^^^^^^^ meta.string.julia string.quoted.double.block.julia
  """

  s = "$('\U03b1')\U03b2"
#       ^^^^^^^^^^ meta.string.julia meta.interpolation.julia
  s = "$(apple + 1)\U03b2"
#        ^^^^^^^^^ meta.string.julia meta.interpolation.julia
#                  ^^^^^^ meta.string.julia string.quoted.double.julia constant.character.escape.unicode.16-bit-hex.julia
  r = r"[a-z0-9]{2, 3}"
#       ^^^^^^^^^^^^^^ meta.string.julia string.quoted.double.julia source.regexp
  r"""
  [a-z0-9]{2, 3}"""
# ^^^^^^^^^^^^^^ meta.string.julia string.quoted.double.block.julia source.regexp

  r"\""
#   ^^ meta.string.julia string.quoted.double.julia source.regexp.python constant.character.escape.backslash.regexp
  r"""a"b"""
#     ^^^ meta.string.julia string.quoted.double.block.julia source.regexp.python

  b"DATA\xff\u2200"
#       ^^^^ meta.string.julia string.quoted.double.julia constant.character.escape.hex.julia
#           ^^^^^^ meta.string.julia string.quoted.double.julia constant.character.escape.unicode.32-bit-hex.julia
  b"""DATA\xff\u2200"""
#         ^^^^ meta.string.julia string.quoted.double.block.julia constant.character.escape.hex.julia
#             ^^^^^^ meta.string.julia string.quoted.double.block.julia constant.character.escape.unicode.32-bit-hex.julia
  v"3.0-"
# ^ storage.type.string.julia

  raw"abc\"few"
#        ^^ meta.string.julia string.quoted.double.julia constant.character.escape.julia

##
## List comprehension ====
##

  (s for s in S)
#  ^^^^^^^^^^^^^ meta.group.julia
#    ^^^ meta.group.julia keyword.control.julia
#                ^ - meta.group.julia
  (s for s in S if s > 0)
#  ^^^^^^^^^^^^^^^^^^^^^ meta.group.julia
#    ^^^ meta.group.julia keyword.control.julia
#               ^^ meta.group.julia keyword.control.julia
#                        ^ - meta.group.julia

  [s for s in S]
#  ^^^^^^^^^^^^^ meta.sequence.julia
#    ^^^ meta.sequence.julia keyword.control.julia
#                ^ - meta.sequence.julia
  [s for s in S if s > 0]
#  ^^^^^^^^^^^^^^^^^^^^^ meta.sequence.julia
#    ^^^ meta.sequence.julia keyword.control.julia
#               ^^ meta.sequence.julia keyword.control.julia
#

max_stars_repo_name: krrutkow/Julia-sublime

max_issues_repo_name: krrutkow/Julia-sublime

max_forks_repo_name: krrutkow/Julia-sublime


text:
@testset "TableDataset" begin
    @testset "TableDataset from rowaccess table" begin
        Tables.columnaccess(::Type{<:Tables.MatrixTable}) = false
        Tables.rowaccess(::Type{<:Tables.MatrixTable}) = true

        testtable = Tables.table([1 4.0 "7"; 2 5.0 "8"; 3 6.0 "9"])
        td = TableDataset(testtable)

        @test collect(@inferred(getobs(td, 1))) == [1, 4.0, "7"]
        @test numobs(td) == 3
    end

    @testset "TableDataset from columnaccess table" begin
        Tables.columnaccess(::Type{<:Tables.MatrixTable}) = true
        Tables.rowaccess(::Type{<:Tables.MatrixTable}) = false

        testtable = Tables.table([1 4.0 "7"; 2 5.0 "8"; 3 6.0 "9"])
        td = TableDataset(testtable)

        @test collect(@inferred(NamedTuple, getobs(td, 2))) == [2, 5.0, "8"]
        @test numobs(td) == 3

        @test getobs(td, 1) isa NamedTuple
    end

    @testset "TableDataset from DataFrames" begin
        testtable = DataFrame(
            col1 = [1, 2, 3, 4, 5],
            col2 = ["a", "b", "c", "d", "e"],
            col3 = [10, 20, 30, 40, 50],
            col4 = ["A", "B", "C", "D", "E"],
            col5 = [100.0, 200.0, 300.0, 400.0, 500.0],
            split = ["train", "train", "train", "valid", "valid"],
        )
        td = TableDataset(testtable)
        @test td isa TableDataset{<:DataFrame}

        @test collect(@inferred(getobs(td, 1))) == [1, "a", 10, "A", 100.0, "train"]
        @test numobs(td) == 5
    end

    @testset "TableDataset from CSV" begin
        open("test.csv", "w") do io
            write(io, "col1,col2,col3,col4,col5, split\n1,a,10,A,100.,train")
        end
        testtable = CSV.File("test.csv")
        td = TableDataset(testtable)
        @test td isa TableDataset{<:CSV.File}
        @test collect(@inferred(getobs(td, 1))) == [1, "a", 10, "A", 100.0, "train"]
        @test numobs(td) == 1
        rm("test.csv")
    end

    @testset "TableDataset is a table" begin
        testtable = DataFrame(
            col1 = [1, 2, 3, 4, 5],
            col2 = ["a", "b", "c", "d", "e"],
            col3 = [10, 20, 30, 40, 50],
            col4 = ["A", "B", "C", "D", "E"],
            col5 = [100.0, 200.0, 300.0, 400.0, 500.0],
            split = ["train", "train", "train", "valid", "valid"],
        )
        td = TableDataset(testtable)
        @testset for fn in (Tables.istable,
                            Tables.rowaccess, Tables.rows,
                            Tables.columnaccess, Tables.columns,
                            Tables.schema, Tables.materializer)
            @test fn(td) == fn(testtable)
        end
    end
end


text:
# Literals treated as constants
function Base.convert(::Type{Expression}, n::Number)
    if !(typeof(n) <: Expression)
        return Constant(n)
    else
        return n
    end
end

# Register functions and handle literals
macro register(sig)
    splitsig = splitdef(:($sig = nothing))
    name = splitsig[:name]
    args = splitsig[:args]
    typargs = typed_args(args)
    defs = :()
    for typarg in typargs
        splitsig[:args] = typarg
        splitsig[:body] = :(Operation($name, Expression[$(args...)]))
        defs = :($defs; $(combinedef(splitsig)))
    end
    esc(defs)
end
# Create all valid combinations of Expression,Number for function signature
function typed_args(args)
    cases = product(Tuple((:Expression, :Number) for i = 1:length(args))...)
    typargs = Vector{Any}[]
    for case in cases
        typarg = [Expr(:(::), args[i], case[i]) for i = 1:length(args)]
        push!(typargs, typarg)
    end
    typargs[1:end-1] # Last case is all Numbers but that should be specialized
end

# Binary & unary operators and functions
import DiffRules, SpecialFunctions, NaNMath
for (M, f, arity) in DiffRules.diffrules()
    fun = :($M.$f)
    sig = arity == 1 ? :($fun(x)) :
          arity == 2 && :($fun(x,y))
    @eval @register $sig
end

for fun = (:<, :>, :(==), :~, :!, :&, :|, :div, :max, :min)
    basefun = Expr(:., Base, QuoteNode(fun))
    sig = :($basefun(x,y))
    @eval @register $sig
end

# ifelse
@register Base.ifelse(cond,t,f)


text:
using IonSim
using QuantumOptics
using DifferentialEquations
using Distributions
using DataStructures

export simulate_with_ion_sim

# Note: IonSim seems to have problems with timescales other than 1e-6
timescale = 1e-6

function simulate_with_ion_sim(parameters, pulses, num_ions, b_field)
    #############################################
    # This function must return a dictionary of result values. The keys
    #   must represent each possible state, and the values must
    #   represent the probability of each state.
    # By convention, "S" and "D" are used for the bright and dark states,
    #   respectively. So the keys of the dictionary should look like:
    #   --> for num_ions == 1: "S", "D"
    #           e.g., Dict("S" => 0.1, "D" => 0.9)
    #   --> for num_ions == 2: "SS", "SD", "DS", "DD"
    #           e.g., Dict("SS" => 0.1, "SD" => 0.2, "DS" => 0.3, "DD" => 0.4)
    #   --> for num_ions == 3: "SSS", "SSD", "SDS", "SDD", etc.

    #############################################
    # Create and load the ions
    S = ["S-1/2", "S+1/2"]
    D = ["D-5/2", "D-3/2", "D-1/2" ,"D+1/2", "D+3/2", "D+5/2"]
    ions = Array{Ca40}(undef, num_ions)
    for i = 1:num_ions
        ions[i] = Ca40([S; D])
    end

    axial_frequency = parameters["TrapFrequencies.axial_frequency"]
    radial_frequency_1 = parameters["TrapFrequencies.radial_frequency_1"]
    radial_frequency_2 = parameters["TrapFrequencies.radial_frequency_2"]
    chain = LinearChain(
        ions=ions,
        com_frequencies=(x=radial_frequency_1, y=radial_frequency_2, z=axial_frequency),
        vibrational_modes=(;z=[1]))

    #############################################
    # Set up the lasers and the trap
    laser_pulses = Dict()
    for pulse_index in 1:length(pulses)
        pulse = pulses[pulse_index]
        if occursin("729G", pulse["dds_name"])
            laser = Laser(
                k=(x̂ + ẑ)/√2,
                ϵ=(x̂ - ẑ)/√2,
                Δ=pulse["freq"],
            )
            laser_pulses[laser] = pulse
        end
        # TODO: Add 729L1 and 729L2
    end

    lasers = Array{Laser}(undef, length(laser_pulses))
    for (i, (laser, pulse)) in enumerate(laser_pulses)
        lasers[i] = laser
    end

    trap = Trap(configuration=chain, B=b_field*1e-4, Bhat=ẑ, δB=0, lasers=lasers)
    mode = trap.configuration.vibrational_modes.z[1]

    for laser in lasers
        global_beam!(trap, laser)
    end

    #############################################
    # Determine the simulation time when the relevant lasers are on
    simulation_start_time = Inf
    simulation_stop_time = 0
    for laser in lasers
        pulse = laser_pulses[laser]
        simulation_start_time = min(simulation_start_time, pulse["time_on"] / timescale)
        simulation_stop_time = max(simulation_stop_time, pulse["time_off"] / timescale)
    end
    if isinf(simulation_start_time)
        simulation_start_time = 0
        simulation_stop_time = 0
    end
    simulation_total_time = simulation_stop_time - simulation_start_time
    simulation_tspan = range(simulation_start_time - (1e-3*timescale), simulation_stop_time + (1e-3*timescale), length=2)

    println("Total simulation time is $(simulation_total_time*timescale*1e6) μs")
    println("(start time = $(simulation_start_time*timescale*1e6) μs, stop time = $(simulation_stop_time*timescale*1e6) μs)")

    #############################################
    # Helpful functions
    function step_interval(t, t_begin, t_end)
        t >= t_begin && t < t_end ? 1 : 0
    end

    function project(solution, states...)
        real.(expect(ionprojector(trap, states...), solution))
    end

    #############################################
    # Set up the time-dependent E-field
    intensity_factor = 10^(-0.5) # this corresponds to pi time 3 μs
    pi_min = 3e-6
    for laser in lasers
        pulse = laser_pulses[laser]
        t_pi = pi_min * sqrt(intensity_factor / (pulse["amp"] * 10^(-1.5 * pulse["att"] / 10)))
        E = Efield_from_pi_time(t_pi, trap.Bhat, laser, ions[1], ("S-1/2", "D-1/2"))
        laser.E = t -> E * step_interval(t, pulse["time_on"] / timescale, pulse["time_off"] / timescale)
    end

    #############################################
    # Set up the time-dependent phase
    for laser in lasers
        pulse = laser_pulses[laser]
        laser.ϕ = t -> 2π * pulse["phase"]
    end

    #############################################
    # Run the simulation
    initial_state = undef
    if length(ions) == 1
        initial_state = ionstate(trap, "S-1/2")
    elseif length(ions) == 2
        initial_state = ionstate(trap, "S-1/2", "S-1/2")
    elseif length(ions) == 3
        initial_state = ionstate(trap, "S-1/2", "S-1/2", "S-1/2")
    end

    simulation_tstops = SortedSet{Float64}()
    for laser in lasers
        pulse = laser_pulses[laser]
        push!(simulation_tstops, pulse["time_on"] / timescale)
        push!(simulation_tstops, pulse["time_off"] / timescale)
    end
    simulation_tstops = collect(simulation_tstops)
    println("Calculated pulse start/stop times as $simulation_tstops")

    h = hamiltonian(trap, lamb_dicke_order=1, timescale=timescale, rwa_cutoff=1e5)
    @time tout, solution = timeevolution.schroedinger_dynamic(
        simulation_tspan,
        tensor(initial_state, fockstate(mode, 0)),
        h,
        callback=PresetTimeCallback(simulation_tstops, integrator -> return));
    solution = solution[end]

    #############################################
    # Measure the expectation values
    result = undef
    if length(ions) == 1
        result = Dict(
            "S" => sum([project(solution, state) for state=S]),
            "D" => sum([project(solution, state) for state=D]),
            )
    elseif length(ions) == 2
        result = Dict(
            "SS" => sum([project(solution, state1, state2) for state1=S, state2=S]),
            "SD" => sum([project(solution, state1, state2) for state1=S, state2=D]),
            "DS" => sum([project(solution, state1, state2) for state1=D, state2=S]),
            "DD" => sum([project(solution, state1, state2) for state1=D, state2=D]),
            )
    elseif length(ions) == 3
        result = Dict(
            "SSS" => sum([project(solution, state1, state2, state3) for state1=S, state2=S, state3=S]),
            "SSD" => sum([project(solution, state1, state2, state3) for state1=S, state2=S, state3=D]),
            "SDS" => sum([project(solution, state1, state2, state3) for state1=S, state2=D, state3=S]),
            "SDD" => sum([project(solution, state1, state2, state3) for state1=S, state2=D, state3=D]),
            "DSS" => sum([project(solution, state1, state2, state3) for state1=D, state2=S, state3=S]),
            "DSD" => sum([project(solution, state1, state2, state3) for state1=D, state2=S, state3=D]),
            "DDS" => sum([project(solution, state1, state2, state3) for state1=D, state2=D, state3=S]),
            "DDD" => sum([project(solution, state1, state2, state3) for state1=D, state2=D, state3=D]),
            )
    end

    #############################################
    # Apply projection noise and renormalize
    total_probability = 0
    for (state, probability) in result
        result[state] = mean(rand(Binomial(1, max(0, min(probability, 1))), 100))
        total_probability += result[state]
    end
    for (state, probability) in result
        result[state] = probability / total_probability
    end

    #############################################
    # Return the result
    println("Simulation results: $result")
    return result

end

max_stars_repo_name: rmshaffer/pulse-sequence-simulator

max_issues_repo_name: rmshaffer/pulse-sequence-simulator

max_forks_repo_name: rmshaffer/pulse-sequence-simulator


text:
###
# MatrixFactorizations.QRPackedQ
###

struct BandedFactorizationStyle <: ApplyStyle end

factorizestyle(::AbstractBandedLayout) = BandedFactorizationStyle()

copy(A::Applied{BandedFactorizationStyle,typeof(qr)}) = banded_qr(A.args...)
copy(A::Applied{BandedFactorizationStyle,typeof(factorize)}) = banded_qr(A.args...)

banded_lmul!(A::QRPackedQ, B::AbstractVecOrMat) = banded_qr_lmul!(A, B)
banded_lmul!(adjA::Adjoint{<:Any,<:QRPackedQ}, B::AbstractVecOrMat) = banded_qr_lmul!(adjA, B)
banded_rmul!(A::AbstractMatrix, Q::QRPackedQ) = banded_qr_rmul!(A, Q)
banded_rmul!(A::AbstractMatrix, adjQ::Adjoint{<:Any,<:QRPackedQ}) = banded_qr_rmul!(A, adjQ)

lmul!(A::QRPackedQ{<:Any,<:AbstractBandedMatrix}, B::AbstractVecOrMat) = banded_lmul!(A,B)
lmul!(adjA::Adjoint{<:Any,<:QRPackedQ{<:Any,<:AbstractBandedMatrix}}, B::AbstractVecOrMat) = banded_lmul!(adjA,B)
lmul!(A::QRPackedQ{<:Any,BandedSubBandedMatrix{T,C,R,I1,I2}}, B::AbstractVecOrMat) where {T,C,R,I1<:AbstractUnitRange,I2<:AbstractUnitRange} =
    banded_lmul!(A,B)
lmul!(adjA::Adjoint{T,<:QRPackedQ{T,<:BandedSubBandedMatrix{T,C,R,I1,I2,t}}}, B::AbstractVecOrMat) where {T,C,R,I1<:AbstractUnitRange,I2<:AbstractUnitRange,t} =
    banded_lmul!(adjA,B)
# rmul!(A::AbstractMatrix, adjQ::Adjoint{<:Any,<:QRPackedQ{<:Any,<:AbstractBandedMatrix}}) = banded_rmul!(A, adjA)
# rmul!(A::StridedMatrix, adjQ::Adjoint{<:Any,<:QRPackedQ{<:Any,<:AbstractBandedMatrix}}) = banded_rmul!(A, adjA)
rmul!(A::StridedVecOrMat{T}, Q::QRPackedQ{T,B}) where {T<:BlasFloat,B<:AbstractBandedMatrix{T}} = banded_rmul!(A, Q)
rmul!(A::StridedVecOrMat{T}, adjQ::Adjoint{<:Any,QRPackedQ{T,B}}) where {T<:BlasComplex,B<:AbstractBandedMatrix{T}} = banded_rmul!(A, adjQ)
rmul!(A::StridedVecOrMat{T}, adjQ::Adjoint{<:Any,QRPackedQ{T,B}}) where {T<:BlasReal,B<:AbstractBandedMatrix{T}} = banded_rmul!(A, adjQ)

###
# QL
###

ql(A::BandedMatrix{T}) where T = ql!(BandedMatrix{float(T)}(A, (max(bandwidth(A,1),bandwidth(A,1)+bandwidth(A,2)+size(A,1)-size(A,2)),bandwidth(A,2))))
ql(A::Tridiagonal{T}) where T = ql!(BandedMatrix{float(T)}(A, (2,1)))

ql!(A::BandedMatrix) = banded_ql!(A)

function banded_ql!(L::BandedMatrix{T}) where T
    D = bandeddata(L)
    l,u = bandwidths(L)
    ν = l+u+1
    m,n=size(L)
    τ = zeros(T, min(m,n))

    for k = n:-1:max((n - m + 1 + (T<:Real)),1)
        μ = m+k-n
        x = view(D,u+1+μ-k:-1:max(1,u-k+2), k)
        τk = reflector!(x)
        τ[k-n+min(m,n)] = τk
        N = length(x)
        for j = max(1,μ-l):k-1
            reflectorApply!(x, τk, view(D, u+1+μ-j:-1:u+2+μ-j-N,j))
        end
    end
    QL(L, τ)
end

function lmul!(A::QLPackedQ{<:Any,<:BandedMatrix}, B::AbstractVecOrMat)
    require_one_based_indexing(B)
    mA, nA = size(A.factors)
    mB, nB = size(B,1), size(B,2)
    if mA != mB
        throw(DimensionMismatch("matrix A has dimensions ($mA,$nA) but B has dimensions ($mB, $nB)"))
    end
    Afactors = A.factors
    l,u = bandwidths(Afactors)
    D = Afactors.data
    begin
        for k = max(nA - mA + 1,1):nA
            μ = mA+k-nA
            for j = 1:nB
                vBj = B[μ,j]
                for i = max(1,k-u):μ-1
                    vBj += conj(D[i-k+u+1,k])*B[i,j]
                end
                vBj = A.τ[k-nA+min(mA,nA)]*vBj
                B[μ,j] -= vBj
                for i = max(1,k-u):μ-1
                    B[i,j] -= D[i-k+u+1,k]*vBj
                end
            end
        end
    end
    B
end

function lmul!(adjA::Adjoint{<:Any,<:QLPackedQ{<:Any,<:BandedMatrix}}, B::AbstractVecOrMat)
    require_one_based_indexing(B)
    A = adjA.parent
    mA, nA = size(A.factors)
    mB, nB = size(B,1), size(B,2)
    if mA != mB
        throw(DimensionMismatch("matrix A has dimensions ($mA,$nA) but B has dimensions ($mB, $nB)"))
    end
    Afactors = A.factors
    l,u = bandwidths(Afactors)
    D = Afactors.data
    @inbounds begin
        for k = nA:-1:max(nA - mA + 1,1)
            μ = mA+k-nA
            for j = 1:nB
                vBj = B[μ,j]
                for i = max(1,k-u):μ-1
                    vBj += conj(D[i-k+u+1,k])*B[i,j]
                end
                vBj = conj(A.τ[k-nA+min(mA,nA)])*vBj
                B[μ,j] -= vBj
                for i = max(1,k-u):μ-1
                    B[i,j] -= D[i-k+u+1,k]*vBj
                end
            end
        end
    end
    B
end

### QBc/QcBc
function rmul!(A::AbstractMatrix,Q::QLPackedQ{<:Any,<:BandedMatrix})
    mQ, nQ = size(Q.factors)
    mA, nA = size(A,1), size(A,2)
    if nA != mQ
        throw(DimensionMismatch("matrix A has dimensions ($mA,$nA) but matrix Q has dimensions ($mQ, $nQ)"))
    end
    Qfactors = Q.factors
    l,u = bandwidths(Qfactors)
    D = Qfactors.data
    @inbounds begin
        for k = nQ:-1:max(nQ - mQ + 1,1)
            μ = mQ+k-nQ
            for i = 1:mA
                vAi = A[i,μ]
                for j = max(1,k-u):μ-1
                    vAi += A[i,j]*D[j-k+u+1,k]
                end
                vAi = vAi*Q.τ[k-nQ+min(mQ,nQ)]
                A[i,μ] -= vAi
                for j = max(1,k-u):μ-1
                    A[i,j] -= vAi*conj(D[j-k+u+1,k])
                end
            end
        end
    end
    A
end

### AQc
function rmul!(A::AbstractMatrix, adjQ::Adjoint{<:Any,<:QLPackedQ{<:Any,<:BandedMatrix}})
    Q = adjQ.parent
    mQ, nQ = size(Q.factors)
    mA, nA = size(A,1), size(A,2)
    if nA != mQ
        throw(DimensionMismatch("matrix A has dimensions ($mA,$nA) but matrix Q has dimensions ($mQ, $nQ)"))
    end
    Qfactors = Q.factors
    l,u = bandwidths(Qfactors)
    D = Qfactors.data
    @inbounds begin
        for k = max(nQ - mQ + 1,1):nQ
            μ = mQ+k-nQ
            for i = 1:mA
                vAi = A[i,μ]
                for j = max(1,k-u):μ-1
                    vAi += A[i,j]*D[j-k+u+1,k]
                end
                vAi = vAi*conj(Q.τ[k-nQ+min(mQ,nQ)])
                A[i,μ] -= vAi
                for j = max(1,k-u):μ-1
                    A[i,j] -= vAi*conj(D[j-k+u+1,k])
                end
            end
        end
    end
    A
end

function _banded_widerect_ldiv!(A::QL, B)
    error("Not implemented")
end
function _banded_longrect_ldiv!(A::QL, B)
    error("Not implemented")
end
function _banded_square_ldiv!(A::QL, B)
    L = A.factors
    lmul!(adjoint(A.Q), B)
    B .= Ldiv(LowerTriangular(L), B)
    B
end

for Typ in (:StridedVector, :StridedMatrix, :AbstractVecOrMat)
    @eval function ldiv!(A::QL{T,<:BandedMatrix}, B::$Typ{T}) where T
        m, n = size(A)
        if m == n
            _banded_square_ldiv!(A, B)
        elseif n > m
            _banded_widerect_ldiv!(A, B)
        else
            _banded_longrect_ldiv!(A, B)
        end
    end
end


text:
using MySQL
# using DataStreams
using DataFrames

"""
    init_mysql_database(;host = "127.0.0.1", dbname="test",
    username="root", pswd="", mysql_code=nothing, overwrite=false)

Create a MySQL database using the code inside mysql_code

### Arguments

* `host`, `dbname`, `user`, `pswd`
* `mysql_code::String`: String with MySQL code that crates all default tables
* `overwrite::Bool` : Flag, if true and dbname exists, drops all database and re-creates it

### Output

* `con`: Database connection and table-column names map

"""
function init_mysql_database(host="127.0.0.1",
    user="root", pwd="", dbname="test", overwrite=false)

    opts = Dict(MySQL.API.MYSQL_SET_CHARSET_NAME=>"utf8")

    con = MySQL.connect(host, user, pwd, opts=opts)

    init_mysql_database(con, dbname, overwrite)
end

function init_mysql_database(con::MySQL.Connection, dbname="test", overwrite=false)

    q = MySQL.query(con, "SHOW DATABASES LIKE '$dbname'; ")
    if length(q[1])>0
        if overwrite
            println("Set to overwrite MySQL database $dbname")
            MySQL.execute!(con, "DROP DATABASE IF EXISTS $dbname;")
            MySQL.execute!(con, "CREATE DATABASE $dbname
                CHARACTER SET utf8 COLLATE utf8_unicode_ci;")
        end
    else
        MySQL.execute!(con, "CREATE DATABASE $dbname
            CHARACTER SET utf8 COLLATE utf8_unicode_ci;")
    end

    MySQL.execute!(con, "USE $dbname;")

    return con
end

"""
    select_columns_mysql(con, table)

For a MySQL database, return an array of all columns in the given table
"""
function select_columns(con::MySQL.MySQLHandle, table)
    cols_query = MySQL.query(con, "SHOW COLUMNS FROM $table;", DataFrame)
    cols_query[1]
end

"""
    select_all_tables_mysql(con)

Return an array of all tables in a given MySQL database
"""
function select_all_tables(con::MySQL.MySQLHandle)
    tables_query = MySQL.query(con, "SHOW TABLES;", DataFrame)
    tables_query[1]
end

function print_error(con::MySQL.MySQLHandle)
    Base.showerror(stdout, MySQL.MySQLInternalError(con))
    println("\n")
end

"""
    query_mysql(con, query_code)

Execute a mysql command
"""
function db_query(con::MySQLHandle, query_code)
    try
        sel = MySQL.query(con, query_code, DataFrame)
        return sel
    catch
        #error("There was an error with MySQL")
        throw(MySQL.MySQLInternalError(con))
    end
end

function insert_row!(con::MySQL.MySQLHandle, tablename, data_values::Dict{Symbol, T},
    colname_dict::Dict{String, Array{String, 1}}, verbose = false) where T

    table_cols = colname_dict[symbol(tablename)]
    table_cols_backticks = [string("`", x, "`") for x in table_cols]
    cols_string = join(table_cols_backticks, ", ")
    vals_string = assemble_vals(data_values, table_cols)
    lastid = -1

    try
        MySQL.execute!(con, "INSERT INTO `$tablename` ($cols_string) values $vals_string;")
        lastid = MySQL.insertid(con)
    catch e
        # Base.showerror(STDOUT, MySQLInternalError(con))
        @warn "Warning: Row with values $vals_string not inserted into the table: $tablename"
        # throw(MySQLInternalError(con))
    end

    return lastid
end

"""
    insert_row!(db, tablename, values)
Insert a row of values into the specified table for a given a MySQL database handle

### Arguments:

* `db::MySQLDB`: Database object (connection and map)
* `data_values::Dict{String, Any}`: Array of (string) values
* `verbose`: Print debugging info
"""
function insert_row!(con::MySQL.MySQLHandle, tablename, data_values::Dict{Symbol, T},
    verbose = false) where T

    cols_string, vals_string = assemble_cols_and_vals(data_values)
    lastid = -1
    try
        MySQL.execute!(con, "INSERT INTO `$tablename` ($cols_string) values $vals_string;")
        lastid = MySQL.insertid(con)
    catch
        if verbose
            println("Warning! Row with values $vals_string not inserted into the table: $tablename")
            Base.showerror(stdout, MySQL.MySQLInternalError(con))
            println("\n")
        end
        # throw(MySQLInternalError(con))
    end
    return lastid
end

function create_server(con::MySQL.MySQLHandle, dbname; linkname = "fedlink", user="root", psswd="", host="127.0.0.1", port=3306)

    query_str = "CREATE SERVER $linkname
                 FOREIGN DATA WRAPPER mysql
                 OPTIONS (USER '$user', PASSWORD '$psswd', HOST '$host', PORT $port, DATABASE '$dbname');"

    MySQL.query(con, query_str, DataFrame)
end

"""
    disable_foreign_checks(con::MySQL.MySQLHandle)
Disables foreign checks for MySQL database
"""
function disable_foreign_checks(conn::MySQL.MySQLHandle)
    MySQL.execute!(conn, "SET FOREIGN_KEY_CHECKS = 0")
    return nothing
end

"""
    enable_foreign_checks(con::MySQL.MySQLHandle)
Enables foreign checks for MySQL database
"""
function enable_foreign_checks(conn::MySQL.MySQLHandle)
    MySQL.execute!(conn, "SET FOREIGN_KEY_CHECKS = 1")
    return nothing
end

"""
    set_innodb_checks(conn, autocommit = 1, foreign_keys = 1, unique = 1)
"""
function set_innodb_checks!(conn::MySQL.Connection, autocommit::Int = 1, foreign_keys::Int = 1, unique::Int = 1)
    MySQL.execute!(conn, "SET FOREIGN_KEY_CHECKS = $foreign_keys")
    MySQL.execute!(conn, "SET AUTOCOMMIT = $autocommit")
    MySQL.execute!(conn, "SET UNIQUE_CHECKS = $unique")
    return nothing
end


text:
struct Disk <: Shape
    id::Int64
    rmin::Float64
    rmax::Float64
    inverted::Bool
    obj_to_world::Transformation
    world_to_obj::Transformation
end

Disk(T::Transformation) = Disk(1, 0.0, 1.0, false, T, inv(T))
Disk() = Disk(1, 0.0, 1.0, false, Transformation(), Transformation())

can_intersect(D::Disk) = true
area(D::Disk) = π * (D.rmax^2 - D.rmin^2)

function obj_bounds(D::Disk)
    pMin = Point3(-D.rmax, -D.rmax, 0.0)
    pMax = Point3( D.rmax,  D.rmax, 0.0)
    return BoundingBox(pMin, pMax)
end

function world_bounds(D::Disk)
    D.obj_to_world(obj_bounds(D))
end

function (T::Transformation)(D::Disk)
    T2 = T * D.obj_to_world
    Disk(D.id, D.rmin, D.rmax, D.inverted, T2, inv(T2))
end

function shape_intersect(R::Ray, D::Disk)
    ray = D.world_to_obj(R)
    if ray.direction.z ≈ 0.0
        return nothing, NaN, NaN
    end
    t = -ray.origin.z / ray.direction.z
    if t < ray.tmin || t > ray.tmax
        return nothing, NaN, NaN
    end
    P = ray(t)
    r2 = P.x^2 + P.y^2
    if r2 > D.rmax^2 || r2 < D.rmin^2
        return nothing, NaN, NaN
    end

    dpdu = Vector3(-P.y, P.x, 0.0)

    n =Normal3(0, 0, sign(ray.origin.z))
    s =normalize(dpdu)

    DG = DifferentialGeometry(D.obj_to_world(P), D.obj_to_world(n), D.obj_to_world(s))

    n = flip_normal(n, ray.direction)

    return DG, t, 5e-4 * t
end

function intersectP(R::Ray, D::Disk)
    ray = D.world_to_obj(R)
    if ray.direction.z ≈ 0.0
        return false
    end
    t = -ray.origin.z / ray.direction.z
    if t < ray.tmin || t > ray.tmax
        return false
    end
    P = ray(t)
    r2 = P.x^2 + P.y^2
    if r2 > D.rmax^2 || r2 < D.rmin^2
        return false
    end
    return true
end

max_stars_repo_name: dronir/GRAYTR

max_issues_repo_name: dronir/GRAYTR

max_forks_repo_name: dronir/GRAYTR


text:
import EvalMetrics.Encodings: positive_label, negative_label
import EvalMetrics: apply, thresholds

using Base.Iterators: product

targets = [
    collect(1:10 .>= 6),
    collect(1:10 .>= 8),
    collect(1:10 .>= 3)
]

scores = [
    collect(range(0, 1.0; length=10)),
    collect(range(1.0, 0; length=10)),
    ones(10),
    zeros(10),
    [0.98, 0.26, 0.39, 0.13, 0.3, 0.84, 0.78, 0.48, 0.13, 0.31],
    [0.74, 0.48, 0.23, 0.91, 0.33, 0.92, 0.83, 0.61, 0.68, 0.09]
]

auroc_oracle = [1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5,
                0.62, 0.35714285714285715, 0.375, 0.6, 0.2857142857142857, 0.5]

auprc_oracle = [1.0, 1.0, 1.0, 0.30436507936507934, 0.16574074074074074,
                0.5927579365079365, 0.75, 0.65, 0.9, 0.75, 0.65, 0.9,
                0.5349999999999999, 0.22222222222222224, 0.6910714285714286,
                0.6477777777777778, 0.20925925925925926, 0.8595734126984127]

set_encoding(OneZero())

encs = [
    OneZero(),
    OneMinusOne(),
    OneTwo(),
    OneVsOne(3,4),
    OneVsOne(:three,:four),
    OneVsOne("three","four"),
    OneVsRest(1, [2,3,4]),
    OneVsRest(:one, [:two, :three]),
    OneVsRest("one", ["two", "three"]),
    RestVsOne([1,2,3], 4),
    RestVsOne([:one, :two], :three),
    RestVsOne(["one", "two"], "three")
]

@testset "curves for $enc encoding" for enc in encs
    global targets

    targets = [recode.(current_encoding(), enc, y) for y in targets]
    set_encoding(enc)

    @testset "check target" begin
        all_neg = fill(negative_label(enc), 10)
        all_pos = fill(positive_label(enc), 10)

        @test_throws ArgumentError au_roccurve(all_neg, rand(10))
        @test_throws ArgumentError au_roccurve(all_pos, rand(10))
        @test_throws ArgumentError au_prcurve(all_neg, rand(10))
        @test_throws ArgumentError au_prcurve(all_pos, rand(10))
        @test_throws ArgumentError au_roccurve(enc, all_neg, rand(10))
        @test_throws ArgumentError au_roccurve(enc, all_pos, rand(10))
        @test_throws ArgumentError au_prcurve(enc, all_neg, rand(10))
        @test_throws ArgumentError au_prcurve(enc, all_pos, rand(10))
    end

    iters = zip(auroc_oracle, auprc_oracle, Iterators.product(targets, scores))

    for (auroc_o, auprc_o, (y, s)) in iters
        thres = thresholds(s)

        fpr = false_positive_rate(y, s, thres)
        tpr = true_positive_rate(y, s, thres)
        rec = recall(y, s, thres)
        prec = precision(y, s, thres)
        cms = ConfusionMatrix(y, s, thres)

        @testset "apply" begin
            @test apply(ROCCurve, cms) == (fpr, tpr)
            @test apply(PRCurve, cms) == (rec, prec)
        end

        @testset "auc for ROCCurve" begin
            @test auc(ROCCurve, cms) ≈ auroc_o
            @test auc(ROCCurve, y, s) ≈ auroc_o
            @test auc(ROCCurve, enc, y, s)≈ auroc_o
            @test auc(ROCCurve, y, s, thres) ≈ auroc_o
            @test auc(ROCCurve, enc, y, s, thres) ≈ auroc_o

            @test au_roccurve(cms) ≈ auroc_o
            @test au_roccurve(y, s) ≈ auroc_o
            @test au_roccurve(enc, y, s) ≈ auroc_o
            @test au_roccurve(y, s, thres) ≈ auroc_o
            @test au_roccurve(enc, y, s, thres) ≈ auroc_o
        end

        @testset "auc for PRCurve" begin
            @test auc(PRCurve, cms) ≈ auprc_o
            @test auc(PRCurve, y, s) ≈ auprc_o
            @test auc(PRCurve, enc, y, s) ≈ auprc_o
            @test auc(PRCurve, y, s, thres) ≈ auprc_o
            @test auc(PRCurve, enc, y, s, thres) ≈ auprc_o

            @test au_prcurve(cms) ≈ auprc_o
            @test au_prcurve(y, s) ≈ auprc_o
            @test au_prcurve(enc, y, s) ≈ auprc_o
            @test au_prcurve(y, s, thres) ≈ auprc_o
            @test au_prcurve(enc, y, s, thres) ≈ auprc_o
        end
    end
end

using Plots

set_encoding(OneZero())

@testset "curves plotting" begin
    for trial in 1:10
        for n in [2, 5, 10, 10^3, 10^5]
            targets = rand(Bool, n)
            targets[1] = true
            targets[2] = false
            scores = rand(n)
            thres = sort(rand(max(2, n ÷ 2)))

            enc_args = [(OneZero(),), ()]
            target_scores_args = [(targets, scores), (fill(targets, 3), fill(scores, 3))]
            threshold_args = [(thres,), ()]

            for (e_args, ts_args, t_args) in product(enc_args, target_scores_args, threshold_args)
                # the following shouldn't fail
                prplot(e_args..., ts_args..., t_args...);
                prplot!(e_args..., ts_args..., t_args...);
                rocplot(e_args..., ts_args..., t_args...);
                rocplot!(e_args..., ts_args..., t_args...);
            end
        end
    end
end


text:
using Plots
using DynamicBounds
using DiffEqBase: ODEFunction
using DataFrames

function f(du,u,p,t)
    du[1] = 0.0
    for i in 2:20
        du[i] = p[1]*(u[i+1] - 2.0*u[i] + u[i-1])
    end
    du[21] = 0.0
    return
end
function fj(Jout,u,p,t)

end
function x0(p)
    x = collect(range(-1,1,step=0.1))
    return cos.(pi*x./2.0)
end
pL = [0.9]
pU = [1.1]
nt = 21
nx = 21

kwargs = Dict{Symbol,Any}()
prob = ODERelaxProb(f, (0,1), collect(1:(nt-1)), x0, 0.5*(pL+pU), pL, pU, nt,
                    xLconst = zeros(nx), xUconst = zeros(nx), kwargs...)
#integrator = make(prob, Wilhelm2019(AM2()))
npnts = 20
value = zeros(nx,nt)
lower_bound = zeros(nx,nt)
upper_bound = zeros(nx,nt)
convex = zeros(nx,nt)
concave = zeros(nx,nt)
for p in range(pL, pU, length = npnts)
    set!(integrator, ParameterValue(1), p)
    integrate!(integrator)
    getall!(values, integrator, Value())
    getall!(lower_bound, integrator, Bound{LOWER}())
    getall!(upper_bound, integrator, Bound{UPPER}())
    getall!(convex, integrator, Relaxation{LOWER}())
    getall!(concave, integrator, Relaxation{UPPER}())
end

#=
4x4 plot
1 surface
2 interval/convex
3 convex with
End Point... Relative Trajectory
=#


text:
using TC
using Base.Test

immutable Iterable{T} <: TypeClass
  start::Arrow{Tuple{T},State{T}}
  next::Arrow{Tuple{T,State{T}},Tuple{Value{T},State{T}}}
  done::Arrow{Tuple{T,State{T}},Bool}
end
@generated iterable{T<:Iterable}( ::Type{T} ) = TC._is_of_type( T )?:Iterable: :T

immutable Enumerable{T} <: TypeClass
  _super::Tuple{Iterable{T}}
  length::Arrow{Tuple{T},Int64}
end

function foo{T}( arg::T, _is::Type{Iterable} = iterable( Iterable{T} ) )
  first( arg )
end

@test is_of_type( Iterable{Vector{Int64}}) == true
@test is_of_type( Iterable{Symbol}) == false
@test is_of_type( Enumerable{Vector{Int64}}) == true

function fails( x )
  try
    foo( x )
    return false
  catch e
    return true
  end
end

@test foo( [ 1,2,3]) == 1
@test fails( :symbol ) == true


text:
using PyCall
using Base.Test

# Try your first TensorFlow program
# https://github.com/tensorflow/tensorflow

@pyimport tensorflow as tf

hello = tf.constant("Hello, TensorFlow!")
@test isa(hello, PyObject)

sess = tf.Session()
result = sess[:run](hello)
@test isa(result, String)
@test "Hello, TensorFlow!" == result

a = tf.constant(10)
b = tf.constant(32)
result = sess[:run](a[:__add__](b))
@test isa(result, Array{Int32,0})
@test 42 == result[1]


text:
# Windows
# =======
#
# Sliding window iterator.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

module Windows

export
    eachwindow,
    EachWindowIterator,
    EachWinItr,
    SeqWinItr,
    missed

import Bio

typealias ArrayOrStringOrSeq Union{AbstractArray, AbstractString, Bio.Seq.BioSequence}

"""
An iterator which moves across a container, as a sliding window.

Constructor requires the data to move the window across, the width of the
window, and the step of the window.
"""
immutable EachWindowIterator{T <: ArrayOrStringOrSeq}
    "A reference to the collection to iterate over."
    data::T
    "The width of the sliding window."
    width::Int
    "How many positions the sliding window moves along."
    step::Int
    function EachWindowIterator{T}(data::T, width::Int, step::Int = 1)
        @assert width >= 1 "Window width must be ≥ 1."
        @assert step >= 1 "step must be ≥ 1."
        @assert width <= length(data) "The window size cannot be greater than number of data elements."
        return new(data, width, step)
    end
end

typealias EachWinItr EachWindowIterator
typealias SeqWinItr EachWindowIterator{Bio.Seq.BioSequence}

"""
Calculate the number of windows that will result from iterating across the container.

Accepts one variable of type EachWindowIterator.
"""
function Base.size(winitr::EachWindowIterator)
    return length(StepRange(winitr.width, winitr.step, length(winitr.data)))
end

Base.length(winitr::EachWindowIterator) = size(winitr)

"""
Calculate the number of elements in the container that are missed during iteration.
Typically because of the width and step size of the window: some elements near the end of the
container are missed if window sizes and step sizes are too large.
"""
function missed(winitr::EachWindowIterator)
    l = length(winitr.data)
    return l - StepRange(winitr.width, winitr.step, l).stop
end

"""
A sliding window iterator is considered equal to another
sliding window iterator if the data is considered equal, and the
width and step of the iterator is also equivalent.
"""
function Base.:(==)(x::EachWindowIterator, y::EachWindowIterator)
    return (x.data == y.data) && (x.width == y.width) && (x.step == y.step)
end

"""
Convienient function for constructing an EachWindowIterator.

Accepts the same arguments as the EachWindowIterator constructor.
"""
function eachwindow{T <: ArrayOrStringOrSeq}(data::T, width::Int, step::Int = 1)
    EachWindowIterator{T}(data, width, step)
end

endpoint(state::Int, width::Int) = state + width - 1
winrange(state::Int, width::Int) = state:endpoint(state, width)
nextstate(state::Int, step::Int) = state + step

@inline function Base.start(it::EachWindowIterator)
    return 1
end

@inline function Base.next(it::EachWindowIterator, state::Integer)
    rng = winrange(state, it.width)
    return (rng, view(it.data, rng)), nextstate(state, it.step)
end

# Extra next method to account for fact than strings don't have sub method
# like arrays and NucleicAcid sequences do.
# The operation is an indexing operation, rather than a "proper"
# substring or subarray operation.
@inline function Base.next{T <: AbstractString}(it::EachWindowIterator{T}, state::Integer)
    rng = winrange(state, it.width)
    return (rng, it.data[rng]), nextstate(state, it.step)
end

@inline function Base.done(it::EachWindowIterator, state::Integer)
    return endpoint(state, it.width) > length(it.data)
end

function Base.show(io::IO, it::EachWindowIterator)
    print(io, "Sliding window iterator.\nContainer size: ", length(it.data), " elements.",
    "\nWindow width: ", it.width, "\nStep distance: ", it.step,
    "\nTotal windows: ", size(it), "\nMissed elements: ",
    missed(it))
end

end