text:
import ..InterpolationAndGrids:
    linear_interpolation,
    InterpolationGrid

import ..UncertainDatasets:
    UncertainIndexValueDataset,
    AbstractUncertainValueDataset

import ..SamplingConstraints:
    TruncateQuantiles

"""
    resample(udata::UncertainIndexValueDataset,
        grid::InterpolationGrid;
        trunc::TruncateQuantiles = TruncateQuantiles(0.001, 0.999))

Draw a realization of `udata`, then interpolate the data values to `grid`.

To avoid very large spans of interpolation, the uncertain indices are truncated to some
large quantile range. Values are not truncated.
"""
function resample(udata::UncertainIndexValueDataset,
        grid::InterpolationGrid;
        trunc::TruncateQuantiles = TruncateQuantiles(0.001, 0.999))

    # Constrain data so that all furnishing distributions for the indices have finite
    # supports.
    constrained_udata = UncertainIndexValueDataset(
        constrain(udata.indices, trunc),
        udata.values
        )

    # Resample the dataset
    inds, vals = resample(constrained_udata)

    # Interpolate to the provided grid
    intp = linear_interpolation(inds, vals, extrapolation_bc = grid.extrapolation_bc)

    # Return grid indices and the interpolated points.
    range(grid), intp(grid.min:grid.step:grid.max)
end

"""
    resample(udata::UncertainIndexValueDataset,
        sequential_constraint::SequentialSamplingConstraint,
        grid::InterpolationGrid;
        trunc::TruncateQuantiles = TruncateQuantiles(0.001, 0.999))

Draw a realization of `udata`, enforcing a `sequential_constraint` on the indices. Then,
interpolate the values of the realization to the provided grid of indices (`grid`).

To avoid very large spans of interpolation, the uncertain indices are truncated to some
large quantile range. Values are not truncated.
"""
function resample(udata::UncertainIndexValueDataset,
        sequential_constraint::SequentialSamplingConstraint,
        grid::InterpolationGrid;
        trunc::TruncateQuantiles = TruncateQuantiles(0.001, 0.999))

    # Constrain data so that all furnishing distributions for the indices have finite
    # supports.
    constrained_udata = UncertainIndexValueDataset(
        constrain(udata.indices, trunc),
        udata.values
        )

    # Resample the dataset with the sequential constraint.
    inds, vals = resample(constrained_udata, sequential_constraint)

    # Interpolate to the desired grid.
    intp = linear_interpolation(inds, vals, extrapolation_bc = grid.extrapolation_bc)

    # Return grid indices and the interpolated points.
    range(grid), intp(grid.min:grid.step:grid.max)

end

function resample(udata::UncertainIndexValueDataset,
        sequential_constraint::SequentialSamplingConstraint,
        grid::InterpolationGrid,
        n::Int;
        trunc::TruncateQuantiles = TruncateQuantiles(0.001, 0.999))

    # Constrain data so that all furnishing distributions for the indices have finite
    # supports.
    constrained_udata = UncertainIndexValueDataset(
        constrain(udata.indices, trunc),
        udata.values
        )

    # Interpolate
    grid = range(grid)
    resampled_vals = zeros(Float64, length(grid), n)

    for i = 1:n
        # Resample using the sequential constraint, then interpolate
        # the realization to the provided grid.
        inds, vals = resample(constrained_udata, sequential_constraint)

        intp = linear_interpolation(inds, vals,
            extrapolation_bc = grid.extrapolation_bc)

        # Each interpolated realization is a column in `resampled_vals`
        resampled_vals[:, i] = intp(grid.min:grid.step:grid.max)
    end

    # Return grid indices and the interpolated points. The interpolated
    # points now live in a matrix where each column is a realization.
    grid, resampled_vals
end

export resample


text:
using GraphBasedSystems
using LinearAlgebra
GBS = GraphBasedSystems

A = [
    0 1 0 1 1 0 1 0 1 0
    1 0 1 0 0 0 0 0 0 0
    0 1 0 0 0 0 0 0 0 0
    1 0 0 0 1 0 0 0 0 0
    1 0 0 1 0 1 0 0 0 0
    0 0 0 0 1 0 0 0 0 0
    1 0 0 0 0 0 0 1 0 0
    0 0 0 0 0 0 1 0 1 1
    1 0 0 0 0 0 0 1 0 0
    0 0 0 0 0 0 0 1 0 0
]

B = [
    0 1 0 0 0 1 0 0 0
    1 0 0 0 0 1 0 0 0
    0 0 0 0 0 1 1 0 1
    0 0 0 0 0 0 1 1 0
    0 0 0 0 0 0 0 1 1
    1 1 1 0 0 0 1 0 0
    0 0 1 1 0 1 0 1 1
    0 0 0 1 1 0 1 0 1
    0 0 1 0 1 0 1 1 0
]

C = [
    0 1 1 0 0 0
    1 0 0 1 0 0
    1 0 0 1 1 0
    0 1 1 0 0 1
    0 0 1 0 0 1
    0 0 0 1 1 0
]

ZAA = zeros(Int64,10,10)
ZAB = zeros(Int64,10,9)
ZBA = ZAB'
ZAC = zeros(Int64,10,6)
ZCA = ZAC'
ZBC = zeros(Int64,9,6)
ZCB = ZBC'

# Graph 1 is disconnected, 2-3-4-5 are connected, 6 is disconnected

A = [
    A   ZAA ZAB ZAC ZAA ZAA
    ZAA A   ZAB ZAC ZAA ZAA
    ZBA ZBA B   ZBC ZBA ZBA
    ZCA ZCA ZCB C   ZCA ZCA
    ZAA ZAA ZAB ZAC A   ZAA
    ZAA ZAA ZAB ZAC ZAA A
]

A[13,21] = A[21,13] = 1
A[16,30] = A[30,16] = 1
A[20,36] = A[36,20] = 1

system = System{Float64}(A, rand(0:3,size(A)[1]))

for i=1:10
    for entry in system.matrix_entries.nzval
        GBS.randomize!(entry)
    end
    for entry in system.vector_entries
        GBS.randomize!(entry)
    end

    F = full_matrix(system)
    f = full_vector(system)
    ldu_solve!(system)
    @test maximum(abs.(full_vector(system)-F\f)) < 1e-3
end

ext: jl


text:
using CUDArt
import CUDAdrv
using Base.Test

#########################
# Device init and close #
#########################
@test length(CUDArt.ptxdict) == 0
CUDArt.init(0)
@test in(0, keys(CUDArt.ptxdict))
CUDArt.close(0)
@test length(CUDArt.ptxdict) == 0

#########################
# Test CUDA array types #
#########################
# Test function for issue #41
function func41(a, dims)
    b = reinterpret(eltype(a), a, dims)
    Cint[size(b)...]
end

result = devices(dev->capability(dev)[1] >= 2, nmax=1) do devlist
    # Copying memory to and from device
    @test length(devlist) == 1
    device(devlist[1])
    # Test CudaArrays and CudaPitchedArrays
    for AT in (CudaArray, CudaPitchedArray)
        g = AT(Float64, (5,3))
        if AT == CudaPitchedArray
            p = pointer(g)
            @test p.xsize == 5
            @test p.ysize == 3
        end
        h = to_host(g)
        @test typeof(h) == Array{Float64,2}
        @test size(h) == (5,3)
        g2 = AT(1:5)
        g3 = similar(g2)
        @test eltype(g3) == eltype(g2)
        @test size(g3) == size(g2)
        g3 = similar(g2, Float32)
        @test eltype(g3) == Float32
        g3 = similar(g2, (7,))
        @test eltype(g3) == Int
        @test size(g3) == (7,)
        @test typeof(g3) == AT{Int, 1}
        h2 = to_host(g2)
        @test typeof(h2) == Array{Int,1}
        @test h2 == [1:5;]
        if AT == CudaPitchedArray
            p = pointer(g2)
            @test p.xsize == 5
            @test p.ysize == 1
        end
        # fill!
        fill!(g2, 0x00)
        h2 = to_host(g2)
        @test h2 == zeros(Int, 5)
        fill!(g2, 0)
        h2 = to_host(g2)
        @test h2 == zeros(Int, 5)
        fill!(g, 0.0)
        h = to_host(g)
        @test all(h .== 0)
        fill!(g, 0)
        h = to_host(g)
        @test all(h .== 0)
        fill!(g, 8.4)
        h = to_host(g)
        @test all(h .== 8.4)
        fill!(g, NaN)
        h = to_host(g)
        @test isequal(h, fill(NaN, size(g)))
        gbig = AT(Int32, (2000,1000))
        fill!(gbig, -17)
        hbig = to_host(gbig)
        @test all(hbig .== -17) && size(hbig) == (2000,1000)
        # Element-type conversions
        h32 = rand(Float32, (5,3))
        g64 = AT(Float64, (5,3))
        copy!(g64, h32)
        # AbstractArray fallbacks
        s32 = view(h32, 5:-1:1, :)
        copy!(g64, s32)
    end
    # Getting portions of an array
    h_src = reshape(1.0:15.0, 3, 5)
    d_src = CudaPitchedArray(h_src)
    d_dest = CudaPitchedArray(Float64, (1, 2))
    get!(d_dest, d_src, (2,2:3), NaN)
    h_dest = to_host(d_dest)
    @test h_dest == [5 8]
    get!(d_dest, d_src, (2,0:1), NaN)
    h_dest = to_host(d_dest)
    @test isequal(h_dest, [NaN 2])
    get!(d_dest, d_src, (3,5:6), 7.3)
    h_dest = to_host(d_dest)
    @test isequal(h_dest, [15 7.3])
    # Copies between CudaPitchedArrays and CudaArrays
    A = rand(map(UInt16, 5:11), 7, 2)
    GA = CudaArray(A)
    GP = CudaPitchedArray(eltype(A), size(A))
    copy!(GP, GA)
    H = to_host(GP)
    @test H == A
    A[3,2] = 0
    copy!(GP, A)
    copy!(GA, GP)
    H = to_host(GA)
    @test H == A
    # HostArray (pinned memory)
    A = HostArray(Int32, (6,4))
    A[2,2] = 7
    @test A[2,2] == 7
    G = CudaArray(A)
    device_synchronize()
    A[2,2] = 3
    @test A[2,2] == 3
    copy!(A, G)
#     @test A[2,2] == 3  # asynchronous
    device_synchronize()
    @test A[2,2] == 7
    S = unsafe_wrap(Array, pointer(A)+48, (6,2), false)  # a "SubArray"
    B = rand(map(Int32, 1:15), 6, 2)
    GB = CudaArray(B)
    copy!(S, GB)
    device_synchronize()
    @test A[:,3:4] == B
    GB = CudaPitchedArray(2B)
    copy!(S, GB)
    device_synchronize()
    @test A[:,3:4] == 2B
    # test vec
    let
        m = 20
        n = 5
        Am = rand(m, n)
        Av = rand(n)
        d_Am = CudaArray(Am)
        d_Av = CudaArray(Av)
        d_Amvec = vec(d_Am)
        d_Avvec = vec(d_Av)
        @test to_host(d_Amvec) == vec(Am)
        @test to_host(d_Avvec) == Av
    end
    # Issue #41
    a = CudaArray(zeros(2))
    @test func41(a, (2,1,1,1)) == Cint[2,1,1,1]
    # Comment in PR #50
    cudasleep(1)
end

gc()  # check for finalizer errors

#########################################
# Streams, multiple devices, and wait() #
#########################################
devices(dev->CUDArt.capability(dev)[1] >= 2, nmax=1) do devlist
    sleeptime = 0.5
    dev = first(devlist)
    device(dev)
    stream = Stream()
    # Force one run to precompile
    cudasleep(sleeptime; dev=dev, stream=stream)
    wait(stream)
    # Since the frequency is dynamic, cudasleep can be unreliable (at
    # least on some GPUs), so do a median-of-3
    sleeptimes = Float64[0, 0, 0]
    for i = 1:3
        tstart = time()
        cudasleep(sleeptime; dev=dev, stream=stream)
        wait(stream)
        tstop = time()
        sleeptimes[i] = tstop - tstart
    end
    @test 0.75*sleeptime <= median(sleeptimes) <= 1.25*sleeptime
end

if devcount() > 1
    devices(dev->true, nmax=2) do devlist
        sleeptime = 0.5
        results = Array{Any}(ceil(Int, 2.5*length(devlist)))
        streams = [(device(dev); Stream()) for dev in devlist]
        wait(streams[1])
        i = 1
        nextidx() = (idx=i; i+=1; idx)
        @sync begin
            for idev = 1:length(devlist)
                sleep(0.01)
                @async begin
                    while true
                        idx = nextidx()
                        if idx > length(results)
                            break
                        end
                        tstart = time()
                        dev = devlist[idev]
                        stream = streams[idev]
                        cudasleep(sleeptime; dev=dev, stream=stream)
                        wait(stream)
                        tstop = time()
                        results[idx] = (tstart, tstop, dev)
                    end
                end
            end
        end
        for i = 1:length(results)
            r = results[i]
            @test sleeptime <= (r[2]-r[1]) <= 1.1*sleeptime
            @test r[3] == devlist[mod1(i, 2)]
        end
        @test results[end][2]-results[1][1] < 1.1*sleeptime*ceil(Int, length(results)/length(devlist))
        nothing
    end
end


text:
using Deconvolution
using Base.Test

srand(42) # Fixed random seed

##### Wiener deconvolution

# Without blurring
x = linspace(0, 10, 15)
s = sinpi(x) - 1.5cos(x)
n = rand(length(s))
@test wiener(s + n, s, 0.5) ≈
    [-1.5300909616077885,-0.4607553796475224,-1.73313390900302,0.8555916473930066,
     2.2692599396286277,0.7989152977213652,1.1431959877178226,-0.8118586562136341,
     -1.6519310564150556,-0.4882369786466372,-1.0387206070662094,
     -0.8324097927723406,1.8208048770862535,0.7427114498907227,1.0587917379150844]

# With blurring
blurring_ft  = exp(-0.001 * (range(-div(length(x), 2), length(x)) .^ 2) .^ (5 // 6))
blurred_s_ft = fftshift(blurring_ft) .* fft(s) + fft(n)
blurred_s    = real(ifft(blurred_s_ft))
blurring     = ifft(fftshift(blurring_ft))
@test wiener(blurred_s, s, 0.5, blurring) ≈
    [-1.5288865581538547,-0.4592713125163481,-1.7360819833552192,0.8525544925861452,
     2.272966823544607,0.8021732396525801,1.14091217864045,-0.8166101187683248,
     -1.6455701546601862,-0.4926325667795202,-1.0323127293402468,
     -0.8380228416032021,1.8214493186523952,0.7447075584686396,1.0566048110361044]


text:
"""
Discharge estimates are biased, and this bias corrupts
material budgets based on them. We adjust for this
bias by assuming that the water balance is exactly satisfied.
"""
function adjustbias(Q::Discharge,T::Quantity,M::Mask)
    F = T*Q # Observed instantaneous flux
    t,Fo = tidalbalance(F,M) # Observed material balance
    _,To = tidalbalance(T,M) # Tidal integral of T balance
    _,μ = tidalaverage(Q,M) # Observed average water balance
    t,Fo.-μ.*To
end

function adjustbias_average(Q::Discharge,T::Quantity,M::Mask)
    F = T*Q
    t,Fo = tidalaverage(F,M)
    _,To = tidalaverage(T,M)
    _,μ  = tidalaverage(Q,M)
    t,Fo.-μ.*To
end

function adjustbias_discharge(Q::Discharge,M::Mask)
    t,Qs = segment(Q,M)
    _,μ = tidalaverage(Q,M)
    t,Qs.-μ
end

# We want a way to concatenate discharges that have been adjusted,
# while filling (with zeros?) in outside the mask
function adjustbias_dischargefill(Q::Discharge,M::Mask)
    irise,ifall = edgeindices(M)
    t,Qs = segment(Q,M)
    _,μ = tidalaverage(Q,M)
    Qb = Qs.-μ
    B = zeros(length(Q))
    for i in 1:length(Qb)
        B[irise[i]:ifall[i]] = Qb[i]
    end
    Discharge(times(Q),B)
end


text:
module helpers

using Dates

"""Helper methods for Hijri conversion."""

_DAYS_BEFORE_MONTH = [-1]
dbm = 0
for dim in Dates.DAYSINMONTH
    global dbm
    append!(_DAYS_BEFORE_MONTH, dbm)
    dbm += dim
end

function _days_before_year(year)
    "year -> number of days before January 1st of year."
    y = year - 1
    return y*365 + y÷4 - y÷100 + y÷400
end

function _days_before_month(year, month)
    "year, month -> number of days in year preceding first day of month."
    @assert 1 <= month && month <= 12 "month must be in 1..12"
    return _DAYS_BEFORE_MONTH[month + 1] + (month > 2 && _is_leap(year))
end
#TODO: The functions above should probably be implemented in stdlib

_DI400Y = _days_before_year(401)    # number of days in 400 years
_DI100Y = _days_before_year(101)    #    "    "   "   " 100   "
_DI4Y   = _days_before_year(5)      #    "    "   "   "   4   "

# A 4-year cycle has an extra leap day over what we'd get from pasting
# together 4 single years.
@assert _DI4Y == 4 * 365 + 1 _DI4Y

# Similarly, a 400-year cycle has an extra leap day over what we'd get from
# pasting together 4 100-year cycles.
@assert _DI400Y == 4 * _DI100Y + 1 _DI400Y

# OTOH, a 100-year cycle has one fewer leap day than we'd get from
# pasting together 25 4-year cycles.
@assert _DI100Y == 25 * _DI4Y - 1 _DI100Y

function bisect_left(a, x, lo = 1, hi = nothing)
    if lo < 1
        throw(BoundsError(a, lo))
    end
    if hi === nothing
        hi = length(a) + 1  # It's not `length(a)`!
    end
    while lo < hi
        mid = (lo + hi) ÷ 2
        a[mid] < x ? lo = mid + 1 : hi = mid
    end
    return lo
end

@doc  """Convert Julian day number (JDN) to date ordinal number.

      :param jdn: Julian day number (JDN).
      :type jdn: int
      """
function jdn_to_ordinal(jdn::Int)

    return jdn - 1721425
end

@doc  """Convert date ordinal number to Julian day number (JDN).

      :param n: Date ordinal number.
      :type n: int
      """
function ordinal_to_jdn(n::Int)

    return n + 1721425
end

@doc  """Return Reduced Julian Day (RJD) number from Julian day number (JDN).

      :param jdn: Julian day number (JDN).
      :type jdn: int
      """
function jdn_to_rjd(jdn::Int)

    return jdn - 2400000
end

@doc  """Return Julian day number (JDN) from Reduced Julian Day (RJD) number.

      :param rjd: Reduced Julian Day (RJD) number.
      :type rjd: int
      """
function rjd_to_jdn(rjd::Int)

    return rjd + 2400000
end

@doc   "year, month, day -> ordinal, considering 01-Jan-0001 as day 1."
function _ymd2ord(year, month, day)
    @assert 1 <= month && month <= 12 "month must be in 1..12"
    dim = Dates.daysinmonth(Date(year, month))
    @assert 1 <= day && day <= dim "day must be in 1..%d' % dim"
    return floor(Int, _days_before_year(year) +
            _days_before_month(year, month) +
            day) - 1 # TODO: -1 is not necessary; without it, we should get .5s
end

function _is_leap(year)
    "year -> 1 if leap year, else 0."
    return year % 4 == 0 && (year % 100 != 0 || year % 400 == 0)
end

@doc   "ordinal -> (year, month, day), considering 01-Jan-0001 as day 1."
function _ord2ymd(n)

    # n is a 1-based index, starting at 1-Jan-1.  The pattern of leap years
    # repeats exactly every 400 years.  The basic strategy is to find the
    # closest 400-year boundary at or before n, then work with the offset
    # from that boundary to n.  Life is much clearer if we subtract 1 from
    # n first -- then the values of n at 400-year boundaries are exactly
    # those divisible by _DI400Y:
    #
    #     D  M   Y            n              n-1
    #     -- --- ----        ----------     ----------------
    #     31 Dec -400        -_DI400Y       -_DI400Y -1
    #      1 Jan -399         -_DI400Y +1   -_DI400Y      400-year boundary
    #     ...
    #     30 Dec  000        -1             -2
    #     31 Dec  000         0             -1
    #      1 Jan  001         1              0            400-year boundary
    #      2 Jan  001         2              1
    #      3 Jan  001         3              2
    #     ...
    #     31 Dec  400         _DI400Y        _DI400Y -1
    #      1 Jan  401         _DI400Y +1     _DI400Y      400-year boundary
    n -= 1
    n400, n = divrem(n, _DI400Y)
    year = n400 * 400 + 1   # ..., -399, 1, 401, ...

    # Now n is the (non-negative) offset, in days, from January 1 of year, to
    # the desired date.  Now compute how many 100-year cycles precede n.
    # Note that it's possible for n100 to equal 4!  In that case 4 full
    # 100-year cycles precede the desired day, which implies the desired
    # day is December 31 at the end of a 400-year cycle.
    n100, n = divrem(n, _DI100Y)

    # Now compute how many 4-year cycles precede it.
    n4, n = divrem(n, _DI4Y)

    # And now how many single years.  Again n1 can be 4, and again meaning
    # that the desired day is December 31 at the end of the 4-year cycle.
    n1, n = divrem(n, 365)

    year += n100 * 100 + n4 * 4 + n1
    if n1 == 4 || n100 == 4
        @assert n == 0
        return year-1, 12, 31
    end

    # Now the year is correct, and n is the offset from January 1.  We find
    # the month via an estimate that's either exact or one too large.
    leapyear = n1 == 3 && (n4 != 24 || n100 == 3)
    @assert leapyear == _is_leap(year)
    month = (n + 50) >> 5
    preceding = _DAYS_BEFORE_MONTH[month+1] + (month > 2 && leapyear)
    if preceding > n  # estimate is too large
        month -= 1
        preceding -= Dates.DAYSINMONTH[month] + (month == 2 && leapyear)
    end

    n -= preceding
    @assert 0 <= n && n < Dates.daysinmonth(Date(year, month)) string("n: ", n, ", days in month: ", Dates.daysinmonth((year, month)))

    # Now the year and month are correct, and n is the offset from the
    # start of that month:  we're done!
    return year, month, n+1
end

end

max_stars_repo_name: ceferisbarov/hijri-converter

max_issues_repo_name: ceferisbarov/hijri-converter


text:
using WavePropBase
using WavePropBase.Geometry
using WavePropBase.IO
using WavePropBase.Mesh
using WriteVTK

@testset "VTK export" begin
    # This test should simply not throw an error
    Geometry.clear_entities!()
    Ω, M = WavePropBase.IO.gmsh_sphere()
    vtk_mesh_file(M, joinpath(@__DIR__,"ball")) |> vtk_save
    rm(joinpath(@__DIR__,"ball.vtu"))
    vtk_mesh_file(M, Ω, joinpath(@__DIR__,"ball")) |> vtk_save
    rm(joinpath(@__DIR__,"ball.vtu"))
    vtk_mesh_file(M, external_boundary(Ω), joinpath(@__DIR__,"sphere")) |> vtk_save
    rm(joinpath(@__DIR__,"sphere.vtu"))
end

max_stars_repo_name: WaveProp/WavePropBase

max_issues_repo_name: WaveProp/WavePropBase

max_forks_repo_name: WaveProp/WavePropBase


text:
"""A supertype for all tensor operations."""
abstract type Operation{rank} <: AbstractTensor{rank} end

"""Symbolises summing over given indices."""
struct IndexSumOperation{rank} <: Operation{rank}
    shape::Tuple{Vararg{AbstractVectorSpace}}
    children::Tuple{AbstractTensor, Indices}
    freeindices::Tuple{Vararg{FreeIndex}}
    function IndexSumOperation(A::AbstractTensor, indices::Indices, freeindices::Vararg{FreeIndex})
        new{length(A.shape)}(A.shape, (A, indices), freeindices)
    end
end

"""Checks if we have have an upper and lower index. If there is, we create
an IndexSumOperation node."""
function contractioncheck(A::AbstractTensor)
    contractionindices = []
    for (i, index) in enumerate(A.freeindices)
        if index' in A.freeindices[(i+1):end]
            push!(contractionindices, index)
        end
    end
    if isempty(contractionindices)
        return A
    end
    freeindices = tuple(setdiff(A.freeindices, [i for i in contractionindices], [i' for i in contractionindices])...)
    return IndexSumOperation(A, Indices(contractionindices...), freeindices...)
end

"""A node symbolising addition of tensors."""
struct AddOperation{rank} <: Operation{rank}
    shape::Tuple{Vararg{AbstractVectorSpace}}
    children::Tuple{Vararg{AbstractTensor}}
    freeindices::Tuple{Vararg{FreeIndex}}
    function AddOperation(children::Tuple{Vararg{AbstractTensor{rank}}}, freeindices::Tuple{Vararg{FreeIndex}}) where rank
        if length(children) > 1
            for child in children[2:length(children)]
                if child.shape != children[1].shape
                    throw(DimensionMismatch(string("Shapes ", children[1].shape, " and ", child.shape, " don't match")))
                end
            end
        end
        newchildren = tuple(filter!(x -> !(x isa ZeroTensor), [children...])...)
        if length(newchildren) == 1
            return newchildren[1]
        end
        if length(newchildren) == 0
            return children[1]
        end
        return contractioncheck(new{rank}(children[1].shape, newchildren, freeindices))
    end
end

function Base.:+(nodes::Vararg{Node})
    return AddOperation(nodes, tuple(union([node.freeindices for node=nodes]...)...))
end

"""A node symbolising tensor outer product."""
struct OuterProductOperation{rank} <: Operation{rank}
    shape::Tuple{Vararg{AbstractVectorSpace}}
    children::Tuple{AbstractTensor, AbstractTensor}
    freeindices::Tuple{Vararg{FreeIndex}}
    function OuterProductOperation(shape::Tuple{Vararg{AbstractVectorSpace}}, children::Tuple{AbstractTensor, AbstractTensor}, freeindices::Tuple{Vararg{FreeIndex}})
        return contractioncheck(new{length(shape)}(shape, children, freeindices))
    end
end

function ⊗(x::AbstractTensor, y::AbstractTensor)
    return OuterProductOperation((x.shape..., y.shape...), (x, y),  (x.freeindices..., y.freeindices...))
end

"""Shorthand for multiplying a tensor by a scalar.
If multiplying a number by a variable, the number must be first."""
function Base.:*(x::Scalar, A::AbstractTensor)
    if x == 1 || x == ConstantTensor(1)
        return A
    elseif A == ConstantTensor(1)
        return Tensor(x)
    elseif x == 0 || x isa ZeroTensor || A isa ZeroTensor
        return ZeroTensor(A.shape...)
    else
        return Tensor(x) ⊗ A
    end
end

function Base.:-(A::AbstractTensor)
    return -1*A
end

function Base.:-(A::AbstractTensor{rank}, B::AbstractTensor{rank}) where rank
    return A + (-1*B)
end

struct PowerOperation{rank} <: Operation{rank}
    shape::Tuple{}
    children::Tuple{AbstractTensor{0}, AbstractTensor{0}}
    freeindices::Tuple{Vararg{FreeIndex}}
    function PowerOperation(x::AbstractTensor{0}, y::AbstractTensor{0})
        new{0}((), (x, y), (x.freeindices..., y.freeindices...))
    end
end

function Base.:^(x::AbstractTensor{0}, y::AbstractTensor{0})
    return PowerOperation(x, y)
end

function Base.:^(x::AbstractTensor{0}, y::Number)
    return x^Tensor(y)
end

function Base.sqrt(x::AbstractTensor{0})
    return PowerOperation(x, Tensor(1//2))
end

struct DivisionOperation{rank} <: Operation{rank}
    shape::Tuple{Vararg{AbstractVectorSpace}}
    children::Tuple{AbstractTensor, AbstractTensor{0}}
    freeindices::Tuple{Vararg{FreeIndex}}
    function DivisionOperation(A::AbstractTensor, y::AbstractTensor{0})
        if y isa ZeroTensor
            throw(DivideError())
        elseif y == ConstantTensor(1)
            return A
        end
        return contractioncheck(new{length(A.shape)}((A.shape), (A, y), (A.freeindices..., y.freeindices...)))
    end
end

function Base.:/(A::AbstractTensor, y::AbstractTensor{0})
    return DivisionOperation(A, y)
end

function Base.:/(A::AbstractTensor, y::Number)
    return DivisionOperation(A, Tensor(y))
end

function Base.:/(y::Number, A::AbstractTensor{0})
    return DivisionOperation(Tensor(y), A)
end

"""A node symbolising indexing a tensor by its every dimension."""
struct IndexingOperation{rank} <: Operation{rank}
    shape::Tuple{}
    children::Tuple{AbstractTensor, Indices}
    freeindices::Tuple{Vararg{FreeIndex}}
    function IndexingOperation(x::AbstractTensor, indices::Indices)
        return contractioncheck(new{0}((),(x, indices), tuple(x.freeindices..., [i for i=indices.indices if i isa FreeIndex]...)))
    end
end

"""A node symbolising component tensor."""
struct ComponentTensorOperation{rank} <: Operation{rank}
    shape::Tuple{Vararg{AbstractVectorSpace}}
    children::Tuple{AbstractTensor, Indices}
    freeindices::Tuple{Vararg{FreeIndex}}
    function ComponentTensorOperation(shape::Tuple{Vararg{AbstractVectorSpace}}, children::Tuple{AbstractTensor, Indices}, freeindices::Tuple{Vararg{FreeIndex}})
        new{length(shape)}(shape, children, freeindices)
    end
end

"""Creates a component tensor of A over indices. Indices must be a
subset of the free indices of A."""
function componenttensor(A::AbstractTensor, indices::Vararg{Index})
    if length(indices) == 0
        return A
    end
    if !(indices ⊆ A.freeindices)
        throw(DomainError(string("The free indices ", indices, " are not a subset of ", A.freeindices)))
    end
    shape = tuple(A.shape..., [i.V for i in indices]...)
    freeindices = tuple(setdiff(A.freeindices, indices, [i' for i in indices])...)
    return ComponentTensorOperation(shape, (A, Indices(indices...)), freeindices)
end

idcounter = 0

"""Creates an IndexingOperation of A indexed by ys. If every dimension of A is
not indexed, creates a ComponentTensorOperation over the unindexed dimensions."""
function Base.getindex(A::AbstractTensor, ys::Vararg{Index})
    if length(ys) == 0
        return A
    end
    if length(A.shape) < length(ys)
        throw(BoundsError(A.shape, ys))
    end
    for i in 1:length(ys)
        if A.shape[i] != ys[i].V
            throw(DimensionMismatch(string(ys[i], " is not in the vector space ", A.shape[i])))
        end
    end
    addedindices = []
    global idcounter
    for i in 1:(length(A.shape) - length(ys))
        push!(addedindices, FreeIndex(A.shape[length(ys) + i], "", idcounter))
        idcounter += 1
    end
    return componenttensor(IndexingOperation(A, Indices(ys..., addedindices...)), addedindices...)
end

"""A convenience function that allows indexing a tensor by an integer. The
function creates a corresponding FixedIndex object in the appropiate vector
space.

Technically the function also does the same for a string, but this is not
recommended, because the syntax makes no distinction between upper and
lower indices.
"""
function Base.getindex(A::AbstractTensor, ys::Vararg{Union{String, Int, Index}})
    indexarray = []
    for i in 1:length(ys)
        push!(indexarray, toindex(A.shape[i], ys[i]))
    end
    return Base.getindex(A, indexarray...)
end

function Base.transpose(A::AbstractTensor)
    indices = []
    for V in A.shape
        global idcounter
        push!(indices, FreeIndex(V, "tranpose", idcounter))
        idcounter += 1
    end
    return componenttensor(A[indices...], reverse(indices)...)
end

function trace(A::AbstractTensor{2})
    if A.shape[1] != dual(A.shape[2])
        throw(DomainError(A, "Cannot get a trace from a tensor whose shape is not of the form (V, V')"))
    end
    global counter
    counter += 1
    return A[FreeIndex(A.shape[1], "trace", counter), FreeIndex(A.shape[2], "trace", counter)]
end

function Base.show(io::IO, op::Operation, depth::Int)
    println(io, ["\t" for i in 1:depth]..., typeof(op))
    for child in op.children
        if child isa Operation
            Base.show(io, child, depth + 1)
        else
            println(io, ["\t" for i in 1:(depth+1)]..., child)
        end
    end
end

Base.show(io::IO, op::Operation) = Base.show(io, op, 0)

max_stars_repo_name: Sohva/TensorDSL

max_issues_repo_name: Sohva/TensorDSL

max_forks_repo_name: Sohva/TensorDSL


text:
using Test

@testset "Functions" begin
    include("functions.jl")
end
@testset "Mutable Arithmetics" begin
    include("mutable_arithmetics.jl")
end
@testset "Sets" begin
    include("sets.jl")
end
@testset "Constraints" begin
    include("constraints.jl")
end
@testset "Variables" begin
    include("variables.jl")
end
@testset "Model" begin
    include("model.jl")
end
@testset "Universal Fallback" begin
    include("universalfallback.jl")
end
@testset "Parser" begin
    include("parser.jl")
end
@testset "Mock Optimizer" begin
    include("mockoptimizer.jl")
end
@testset "Caching Optimizer" begin
    include("cachingoptimizer.jl")
end
@testset "DenseDict" begin
    include("dense_dict.jl")
end
@testset "Copy" begin
    include("copy.jl")
end

@testset "CleverDicts" begin
    include("CleverDicts.jl")
end
@testset "DoubleDicts" begin
    include("DoubleDicts.jl")
end
@testset "Lazy iterators" begin
    include("lazy_iterators.jl")
end

@testset "Print with acronym" begin
    include("print_with_acronym.jl")
end


text:
module Tensorial

using LinearAlgebra, Statistics
# re-exports from LinearAlgebra and Statistics
export ⋅, ×, dot, tr, det, norm, normalize, mean, I, cross, eigen, eigvals, eigvecs

using StaticArrays
using Base: @pure, @_inline_meta, @_propagate_inbounds_meta
using ForwardDiff: Dual, value, partials
import StaticArrays: qr, lu, svd, diag, diagm # defined in LinearAlgebra, but call methods in StaticArrays
# re-exports from StaticArrays
export qr, lu, svd, diag, diagm

import Base: transpose, inv
import LinearAlgebra: dot, norm, normalize, tr, adjoint, det, cross, eigen, eigvals, eigvecs
import Statistics: mean

export
# Symmetry/Space
    Symmetry,
    @Symmetry,
    Space,
# AbstractTensor
    AbstractTensor,
    AbstractSecondOrderTensor,
    AbstractFourthOrderTensor,
    AbstractSymmetricSecondOrderTensor,
    AbstractSymmetricFourthOrderTensor,
    AbstractVec,
    AbstractMat,
# Tensor
    Tensor,
    SecondOrderTensor,
    FourthOrderTensor,
    SymmetricSecondOrderTensor,
    SymmetricFourthOrderTensor,
    Vec,
    Mat,
# macros
    @Vec,
    @Mat,
    @Tensor,
    @einsum,
# operations
    contraction,
    otimes,
    dotdot,
    symmetric,
    skew,
    ⊗,
    ⊡,
    rotmat,
    rotmatx,
    rotmaty,
    rotmatz,
    rotate,
# continuum mechanics
    vol,
    dev,
    stress_invariants,
    deviatoric_stress_invariants,
# voigt
    tovoigt,
    fromvoigt,
    tomandel,
    frommandel,
# ad
    gradient,
    hessian,
# quaternion
    Quaternion,
    quaternion,
    angleaxis

include("utils.jl")
include("Symmetry.jl")
include("Space.jl")
include("indexing.jl")
include("AbstractTensor.jl")
include("einsum.jl")
include("Tensor.jl")
include("permute.jl")
include("ops.jl")
include("continuum_mechanics.jl")
include("inv.jl")
include("voigt.jl")
include("ad.jl")
# include("simd.jl")
include("broadcast.jl")

include("quaternion.jl")

const ⊗ = otimes
const ⊡ = double_contraction

end # module


text:
function from_file(filename, sep = '\n')
    s = open(filename) do file
        read(file, String)
    end
    str_sudokus = split(strip(s), sep)
    grids = AbstractArray[]
    for str_sudoku in str_sudokus
        str_sudoku = replace(str_sudoku, "." => "0")
        one_line_grid = parse.(Int, split(str_sudoku, ""))
        grid = reshape(one_line_grid, 9, 9)
        push!(grids, grid)
    end
    return grids
end

function solve_sudoku(grid)
    m = CS.Optimizer(logging = [])

    x = [[MOI.add_constrained_variable(m, MOI.Integer()) for i in 1:9] for j in 1:9]
    for r in 1:9, c in 1:9
        MOI.add_constraint(m, x[r][c][1], MOI.GreaterThan(1.0))
        MOI.add_constraint(m, x[r][c][1], MOI.LessThan(9.0))
    end

    # set variables
    for r in 1:9, c in 1:9
        if grid[r, c] != 0
            sat = [MOI.ScalarAffineTerm(1.0, x[r][c][1])]
            MOI.add_constraint(
                m,
                MOI.ScalarAffineFunction{Float64}(sat, 0.0),
                MOI.EqualTo(convert(Float64, grid[r, c])),
            )
        end
    end
    # sudoku constraints
    for r in 1:9
        MOI.add_constraint(
            m,
            MOI.VectorOfVariables([x[r][c][1] for c in 1:9]),
            CS.AllDifferentSetInternal(9),
        )
    end
    for c in 1:9
        MOI.add_constraint(
            m,
            MOI.VectorOfVariables([x[r][c][1] for r in 1:9]),
            CS.AllDifferentSetInternal(9),
        )
    end
    variables = [MOI.VariableIndex(0) for _ in 1:9]
    for br in 0:2
        for bc in 0:2
            variables_i = 1
            for i in (br * 3 + 1):((br + 1) * 3), j in (bc * 3 + 1):((bc + 1) * 3)
                variables[variables_i] = x[i][j][1]
                variables_i += 1
            end
            MOI.add_constraint(m, variables, CS.AllDifferentSetInternal(9))
        end
    end

    MOI.optimize!(m)
end


text:
#####
##### Inception
#####

struct InceptionBlock
    path_1
    path_2
    path_3
    path_4
end

Flux.@treelike InceptionBlock

function InceptionBlock(
        in_chs,
        chs_1x1,
        chs_3x3_reduce,
        chs_3x3,
        chs_5x5_reduce,
        chs_5x5,
        pool_proj
    )
    path_1 = CrossCor((1, 1), in_chs=>chs_1x1, relu)

    path_2 = (CrossCor((1, 1), in_chs=>chs_3x3_reduce, relu),
            CrossCor((3, 3), chs_3x3_reduce=>chs_3x3, relu, pad = (1, 1, 1, 1)))

    path_3 = (CrossCor((1, 1), in_chs=>chs_5x5_reduce, relu),
            CrossCor((5, 5), chs_5x5_reduce=>chs_5x5, relu, pad = (2, 2, 2, 2)))

    path_4 = (MaxPool((3,3), stride = (1, 1), pad = (1, 1, 1, 1)),
            CrossCor((1, 1), in_chs=>pool_proj, relu))

    InceptionBlock(path_1, path_2, path_3, path_4)
end

function (m::InceptionBlock)(x)
    cat(
        m.path_1(x),
        m.path_2[2](m.path_2[1](x)),
        m.path_3[2](m.path_3[1](x)),
        m.path_4[2](m.path_4[1](x)),
        dims = 3
    )
end

@testset "Testing Inception" begin

    _googlenet() = Chain(CrossCor((7, 7), 3=>64, stride = (2, 2), relu, pad = (3, 3, 3, 3)),
          MaxPool((3, 3), stride = (2, 2), pad = (1, 1, 1, 1)),
          CrossCor((1, 1), 64=>64, relu),
          CrossCor((3, 3), 64=>192, relu, pad = (1, 1, 1, 1)),
          MaxPool((3, 3), stride = (2, 2), pad = (1, 1)),
          InceptionBlock(192, 64, 96, 128, 16, 32, 32),
          InceptionBlock(256, 128, 128, 192, 32, 96, 64),
          MaxPool((3, 3), stride = (2, 2), pad = (1, 1, 1, 1)),
          InceptionBlock(480, 192, 96, 208, 16, 48, 64),
          InceptionBlock(512, 160, 112, 224, 24, 64, 64),
          InceptionBlock(512, 128, 128, 256, 24, 64, 64),
          InceptionBlock(512, 112, 144, 288, 32, 64, 64),
          InceptionBlock(528, 256, 160, 320, 32, 128, 128),
          MaxPool((3, 3), stride = (2, 2), pad = (1, 1, 1, 1)),
          InceptionBlock(832, 256, 160, 320, 32, 128, 128),
          InceptionBlock(832, 384, 192, 384, 48, 128, 128),
          MaxPool((7, 7), stride = (1, 1), pad = (0, 0, 0, 0)),
          x -> reshape(x, :, size(x, 4)),
          #Dropout(0.4),
          Dense(1024, 1000, relu),
          x -> log.(max.(x, Float32(1e-9))),
          softmax,
    )

    model = _googlenet()

    x = rand(Float32, 224, 224, 3, 16)

    backend = nGraph.Backend()
    z = model(x)
    f = nGraph.compile(backend, model, x)

    @test isapprox(z, parent(f()))

    @info "Compiling Training Pass for Inception"
    y = similar(model(x))
    y .= zero(eltype(y))

    model = _googlenet()
    loss(x, y) = sum(model(x) .- y)
    l = loss(x, y)
    @show l

    # @info "Testing nGraph Gradients"
    # h = nGraph.compile(backend, loss, x, y; optimizer = nGraph.Gradient)
    # @show parent(h())

    # # Run the Flux back-propagation
    # @info "Taking Flux Gradients"
    # ps = Flux.params(model)
    # @show length(ps)
    # Flux.back!(l)

    # # Now, we compare.
    # # The `Gradient` optimizer has a dict that maps the original tracked array to the
    # # tensors that were built from them.
    # #
    # # We leverage this to check that
    # #
    # # - The data was captured correctly when parameters were constructed
    # # - The gradients computed by the two frameworks are rougly equal
    # gradient_map = h.optimizer._id

    # for p in ps
    #     @test haskey(gradient_map, p)
    # end

    #=
    for p in ps
        # If this is a convolution weight, we need to flip it to compare data and gradients
        if ndims(p) == 4
            # Get the data tensor from the gradient map. Make sure we copied it correctly.
            x = copy(p.data)
            nGraph.flip!(x)
            @test isapprox(parent(gradient_map[p][1]), x)

            # Check that the gradients are the same as well
            x = copy(p.grad)
            nGraph.flip!(x)
            @test isapprox(parent(gradient_map[p][2]), x)
        else
            @test isapprox(parent(gradient_map[p][1]), p.data)
            @test isapprox(parent(gradient_map[p][2]), p.grad)

        end
    end
    =#
end


text:
using Base: Ryu, IEEEFloat

function formatinfo(s::Spec, x::IEEEFloat)
    # The cost of computing the exact size is high, so we estimate an upper bound.
    if s.type == 'f' || s.type == 'F' || s.type == '%'
        width = Ryu.neededdigits(typeof(x))
        if s.type == '%'
            width += 1
        end
    else
        # Otherwise, the scientific notation of Float64 will be an upper bound.
        # sign (1) + fraction (17) + point (1) + e/E (1) + sign (1) + exponent (3)
        width = 24
    end
    if isspecified(s.precision)
        width += s.precision
    end
    size = width
    if isspecified(s.width)
        size += paddingsize(s, width)
    end
    return size, nothing
end

@inline function formatfield(data::Vector{UInt8}, p::Int, s::Spec, x::IEEEFloat, ::Nothing)
    start = p
    p, signed = sign(data, p, x, s.sign)
    x = abs(x)

    uppercase = s.type == 'F' || s.type == 'E' || s.type == 'A' || s.type == 'G'
    plus = false
    space = false
    hash = s.altform
    expchar = uppercase ? UInt8('E') : UInt8('e')
    if isinf(x)
        if uppercase
            p = @copy data p "INF"
        else
            p = @copy data p "inf"
        end
    elseif isnan(x)
        if uppercase
            p = @copy data p "NAN"
        else
            p = @copy data p "nan"
        end
    elseif s.type == 'F' || s.type == 'f'
        precision = default(s.precision, 6)
        p = writefixed(data, p, x, precision, plus, space, hash)
    elseif s.type == 'E' || s.type == 'e'
        precision = default(s.precision, 6)
        p = writeexp(data, p, x, precision, plus, space, hash, expchar)
    elseif s.type == 'G' || s.type == 'g'
        prec = max(default(s.precision, 6), 1)
        exp = parseexp(data, writeexp(data, p, x, prec))
        if -4 ≤ exp < prec
            p = writefixed(data, p, x, prec - (exp + 1), plus, space, hash; trimtrailingzeros = !hash)
        else
            p = writeexp(data, p, x, prec - 1, plus, space, hash, expchar; trimtrailingzeros = !hash)
        end
    elseif s.type == 'A' || s.type == 'a'
        if uppercase
            p = @copy data p "0X"
        else
            p = @copy data p "0x"
        end
        precision = default(s.precision, -1)
        p = hexadecimal(data, p, x, precision, uppercase)
    elseif s.type == '%'
        precision = default(s.precision, 6)
        p = writefixed(data, p, 100x, precision, plus, space, hash)
        data[p] = UInt8('%')
        p += 1
    else
        @assert !isspecified(s.type)
        if isspecified(s.precision)
            prec = max(s.precision, 1)
            exp = parseexp(data, writeexp(data, p, x, prec))
            if -4 ≤ exp < prec
                hash = true
                p = writefixed(data, p, x, prec - (exp + 1), plus, space, hash; trimtrailingzeros = true)
            else
                p = writeexp(data, p, x, prec - 1, plus, space, hash, expchar; trimtrailingzeros = true)
            end
        else
            hash = true
            p = writeshortest(data, p, x, plus, space, hash, -1, expchar)
        end
    end

    return processfloat(data, start, p, s, signed)
end

# This is to avoid excessive inlining.
@noinline writefixed(data, p, x, precision, plus, space, hash; trimtrailingzeros = false) =
    Ryu.writefixed(data, p, x, precision, plus, space, hash, UInt8('.'), trimtrailingzeros)
@noinline writeexp(data, p, x, precision, plus = false, space = false, hash = false, expchar = UInt8('e'); trimtrailingzeros = false) =
    Ryu.writeexp(data, p, x, precision, plus, space, hash, expchar, UInt8('.'), trimtrailingzeros)
@noinline writeshortest(data, p, x, plus, space, hash, precision, expchar) =
    Ryu.writeshortest(data, p, x, plus, space, hash, precision, expchar)

function parseexp(data, p)
    b1, b2, b3, b4 = data[p-4], data[p-3], data[p-2], data[p-1]
    exp = 10 * (b3 - Z) + (b4 - Z)
    if b1 == UInt8('e')
        # two-digit exponent
        sign = b2
    else
        # three-digit exponent
        sign = b1
        exp += 100 * (b2 - Z)
    end
    return flipsign(exp, sign == UInt8('-') ? -1 : 1)
end

# NOTE: the sign of `x` is ignored
function hexadecimal(data::Vector{UInt8}, p::Int, x::IEEEFloat, precision::Int, uppercase::Bool)
    if iszero(x)
        data[p] = Z
        p += 1
        if precision > 0
            data[p] = UInt8('.')
            p += 1
            p = pad(data, p, '0', precision)
        end
        exp = 0
    else
        # make fr ∈ [1, 2) so that the digit before point is always '1'
        fr, exp = frexp(abs(x))
        fr *= 2
        exp -= 1
        data[p] = UInt8('1')
        p += 1
        if isone(fr)
            if precision > 0
                data[p] = UInt8('.')
                p += 1
                p = pad(data, p, '0', precision)
            end
        else
            # fr ∈ (1, 2)
            if precision != 0
                data[p] = UInt8('.')
                p += 1
            end
            p, carry = hexadecimal_fraction(data, p, Float64(fr), precision, uppercase)
            exp += carry
        end
    end
    data[p] = uppercase ? UInt8('P') : UInt8('p')
    p += 1
    data[p] = exp < 0 ? UInt8('-') : UInt8('+')
    p += 1
    uexp = unsigned(abs(exp))
    return decimal(data, p, uexp, ndigits_decimal(uexp))
end

# fr must be in (1, 2)
function hexadecimal_fraction(data::Vector{UInt8}, p::Int, fr::Float64, precision::Int, uppercase::Bool)
    @assert 1 < fr < 2
    m = 13  # the number of nibbles in significand (excluding implicit bit)
    u = reinterpret(UInt64, fr) & Base.significand_mask(Float64)
    if precision < 0
        # trim trailing zero nibbles
        tz = trailing_zeros(u) ÷ 4
        return hexadecimal(data, p, u >> 4tz, m - tz, uppercase), 0
    elseif precision ≥ m
        # emit all digits and padding zeros
        p = hexadecimal(data, p, u, m, uppercase)
        return pad(data, p, '0', precision - m), 0
    else
        # half-to-even rounding
        k = m - precision
        u1 = u >> 4k
        u2 = u & ((one(u) << 4k) - 0x1)
        half = oftype(u, 0x8) << 4(k - 1)
        carry = 0
        if u2 > half || u2 == half && (u1 & 0x1) > 0
            u1 += 1
            carry = Int(u1 ≥ (one(u) << 4precision))
        end
        return hexadecimal(data, p, u1, precision, uppercase), carry
    end
end

# group digits in integral part by inserting separators (and leading zeros)
function groupfloat(data::Vector{UInt8}, start::Int, p::Int, zero::Bool, minwidth::Int, sep::UInt8)
    # determine integral part
    @assert Z ≤ data[start] ≤ Z + 0x9
    i = start + 1
    while i < p && Z ≤ data[i] ≤ Z + 0x9
        i += 1
    end

    # insert leading zeros
    k = 3           # number of digits between separators
    m = i - start   # number of digits in integral part
    width = p - start  # current content width
    if zero
        z = number_of_leading_zeros(m, k, m + minwidth - width)
        if z > 0
            insert_zeros(data, start, p, z)
            p += z; m += z; i += z
        end
    end

    # insert separators
    n = div(m - 1, k)
    if n > 0
        s = i - k
        copyto!(data, s + n, data, s, p - s)
        data[s+n-1] = sep
        p += n
        s -= k
        n -= 1
        while n > 0
            copyto!(data, s + n, data, s, k)
            data[s+n-1] = sep
            s -= k
            n -= 1
        end
    end

    # return the next position after the last byte
    return p
end

isexpchar(b) = b == UInt8('e') || b == UInt8('E')

# post-processing of formatted floating-point numbers
function processfloat(data, start, p, s, signed)
    # decimal point and fractional zero
    if !isspecified(s.type) && data[p-1] == UInt8('.')
        data[p] = Z
        p += 1
    end

    # grouping and sign-aware zero padding
    if isspecified(s.grouping)
        minwidth = isspecified(s.width) ? s.width - signed : 0
        sep = s.grouping == GROUPING_COMMA ? UInt8(',') : UInt8('_')
        p = groupfloat(data, start + signed, p, s.zero, minwidth, sep)
    elseif s.zero && isspecified(s.width)
        p = insert_zeros(data, start + signed, p, s.width - (p - start))
    end

    # alignment
    if isspecified(s.width)
        p = aligncontent(data, p, start, p - start, s.fill, default(s.align, ALIGN_RIGHT), s.width)
    end
    return p
end


text:
# fitting noisy data to an exponential model
model(xpts, p) = p[1]*exp(-xpts.*p[2])

# some example data
srand(12345)
xpts = linspace(0,10,20)
data = model(xpts, [1.0 2.0]) + 0.01*randn(length(xpts))

beta, r, J = curve_fit(model, xpts, data, [0.5, 0.5])
@assert norm(beta - [1.0, 2.0]) < 0.05

# can also get error estimates on the fit parameters
errors = estimate_errors(beta, r, J)
@assert norm(errors - [0.016, 0.075]) < 0.01


text:
struct Polynomial{DT<:AbstractFloat, BT <: Basis}
    b::BT
    c::Vector{DT}

    function Polynomial{DT,BT}(b::BT, c) where {DT, BT <: Basis{DT}}
        @assert nnodes(b) == length(c)
        new(b, c)
    end
end

function Polynomial(b::Basis{DT}, c::Vector{DT}) where {DT}
    Polynomial{DT, typeof(b)}(b, c)
end

function LagrangePolynomial(x::Vector{DT}, c::Vector{DT}) where {DT}
    b = LagrangeBasis(x)
    Polynomial{DT, typeof(b)}(b, c)
end

function BernsteinPolynomial(x::Vector{DT}, c::Vector{DT}) where {DT}
    b = BernsteinBasis(x)
    Polynomial{DT, typeof(b)}(b, c)
end

function Base.similar(lag::Polynomial{T}, c::Vector{T}) where {T}
    @assert length(c) == length(lag.c)
    Polynomial(lag.b, c)
end

function evaluate!(pol::Polynomial{T}, x::Vector{T}, y::Vector{T}) where {T}
    @assert length(x) == length(y)

    local tx::T
    local ty::T

    for k in eachindex(x,y)
        tx = x[k]
        ty = 0
        for j in eachindex(pol.c)
            ty += pol.c[j] * eval_basis(pol.b, j, tx)
        end
        y[k] = ty
    end

    return nothing
end

function evaluate!(b::LagrangeBasis{T}, c::Vector{T}, x::Vector{T}, y::Vector{T}) where {T}
    evaluate!(Polynomial(b, c), x, y)
end


text:
export Query,
    NumericFilter,
    GeoFilter,
    SortByField,
    limit_ids!,
    return_fields!,
    return_field!,
    summarize,
    highlight!,
    get_args

abstract type AbstractFilter end

"""
    Filter
    Filter{AbstractFilter}

Filter object that can be used in RediSearch

# Fields
- `args::AbstractArray`: Redis args that will be used when searching
"""
mutable struct Filter <: AbstractFilter
    args::AbstractArray
    Filter(keyword, field, args...) = new([keyword, field, args...])
end

function NumericFilter(
    field::String; minval="-inf", maxval="+inf", min_exclusive=true, max_exclusive=true
)
    args = [min_exclusive ? "($(minval)" : minval, max_exclusive ? "$(minval))" : maxval]

    return Filter("FILTER", field, args...)
end

"""
    GeoFilterUnit

BareModule for Geo units that can be used in a filter object.

# Members
- `METERS`: Meters unit of measurement.
- `KILOMETERS`: Kilometers unit of measurement.
- `FEET`: Feet unit of measurement.
- `MILES`: Miles unit of measurement.

# Returns
- `String`
"""
baremodule GeoFilterUnit
    const METERS = "m"
    const KILOMETERS = "km"
    const FEET = "ft"
    const MILES = "mi"
end

function GeoFilter(field, lon, lat, radius; unit=GeoFilterUnit.KILOMETERS)
    return Filter("GEOFILTER", field, lon, lat, radius, unit)
end

"""
    SortByField(field::AbstractString; asc::Bool=true)

Create a `Fitler` object for sorting specific fields.

# Arguments
- `field::AbstractString`: Field name that is being sorted

# KeyWords
- `asc::Bool`: Sort by ascending, default as `true`

# Returns
- `Fitler::AbstractFilter`
"""
function SortByField(field::AbstractString; asc::Bool=true)
    return Filter(field, asc ? "ASC" : "DESC")
end

abstract type AbstractQuery end

"""
    Query(
        query_string::AbstractString;
        offset = 0,
        num = 10,
        no_content = false,
        no_stopwords = false,
        fields = [],
        verbatim = false,
        with_payloads = false,
        with_scores = false,
        scorer = nothing,
        filters = [],
        ids = [],
        slop = -1,
        in_order = false,
        sortby = nothing,
        return_fields = [],
        summarize_fields = [],
        highlight_fields = [],
        language = nothing
    )::Query

Creates a query instance that can be used to search a RediSearch instance

# Arguments
- `query_string::String`

# KeyWords
- `offset::Integer`
- `num::Integer`
- `no_content::Bool`
- `no_stopwords::Bool`
- `fields::AbstractArray`
- `verbatim::Bool`verbatim setting for
- `with_payloads::Bool`
- `with_scores::Bool`
- `scorer::Union{String, Nothing}`
- `filters::AbstractArray{AbstractFilter}`
- `ids::AbstractArray`
- `slop::Integer`
- `in_order::Bool`
- `sortby::Union{AbstractFilter, Nothing}`
- `return_fields::AbstractArray`
- `summarize_fields::AbstractArray`
- `highlight_fields::AbstractArray`
- `language::Union{String, Nothing}`

# Returns
- `Query::AbstractQuery`
"""
mutable struct Query <: AbstractQuery
    query_string::String
    offset::Integer
    num::Integer
    no_content::Bool
    no_stopwords::Bool
    fields::AbstractArray
    verbatim::Bool
    with_payloads::Bool
    with_scores::Bool
    scorer::Union{String,Nothing}
    filters::AbstractArray{AbstractFilter}
    ids::AbstractArray
    slop::Integer
    in_order::Bool
    sortby::Union{AbstractFilter,Nothing}
    return_fields::AbstractArray
    summarize_fields::AbstractArray
    highlight_fields::AbstractArray
    language::Union{String,Nothing}
end

function Query(
    query_string::AbstractString;
    offset=0,
    num=10,
    no_content=false,
    no_stopwords=false,
    fields=[],
    verbatim=false,
    with_payloads=false,
    with_scores=false,
    scorer=nothing,
    filters=[],
    ids=[],
    slop=-1,
    in_order=false,
    sortby=nothing,
    return_fields=[],
    summarize_fields=[],
    highlight_fields=[],
    language=nothing,
)::Query
    return Query(
        query_string,
        offset,
        num,
        no_content,
        no_stopwords,
        fields,
        verbatim,
        with_payloads,
        with_scores,
        scorer,
        filters,
        ids,
        slop,
        in_order,
        sortby,
        return_fields,
        summarize_fields,
        highlight_fields,
        language,
    )
end

"""
    limit_ids!(q::AbstractQuery, ids::AbstractArray)

Alter a query to only return specific IDs
"""
function limit_ids!(q::AbstractQuery, ids::AbstractArray)
    q.ids = ids
    return q
end

"""
    return_fields!(q::AbstractQuery, fields::AbstractArray)

Alter the query object to only return specific fields

# Arguments
- `q::AbstractQuery`: Query object
- `fields::AbstractArray`: List of fields to be returned
"""
function return_fields!(q::AbstractQuery, fields::AbstractArray)
    push!(q.return_fields, fields...)
    return q
end
return_fields!(q::AbstractQuery, field) = return_fields!(q, [field])

"""
    return_field!(q::AbstractQuery, field::AbstractString, as_field::AbstractString)

Map a return field to a new name

# Arguments
- `q::AbstractQuery`: Query object
- `field::AbstractString`: Field to be mapped
- `as_field::AbstractString`: Name for mapping to be made.
"""
function return_field!(q::AbstractQuery, field::AbstractString, as_field::AbstractString)
    push!(q.return_fields, [field, "AS", as_field]...)
    return q
end

"""
    summarize!(
        q::AbstractQuery;
        fields::Union{Nothing,AbstractArray}=nothing,
        context_len::Union{Nothing,Number}=nothing,
        num_frags::Union{Nothing,Number}=nothing,
        sep::Union{Nothing,AbstractString}=nothing,
    )

Create a summary summary object to be used in the query.

# Arguments
- `q::AbstractQuery`: Query object

# KeyWords

- `fields::Union{Nothing,AbstractArray}`:Array of fields to summarize
- `context_len::Union{Nothing,Number}`:Context length
- `num_frags::Union{Nothing,Number}`:Number of Fragments to include
- `sep::Union{Nothing,AbstractString}`:Separater string.
"""
function summarize!(
    q::AbstractQuery;
    fields::Union{Nothing,AbstractArray}=nothing,
    context_len::Union{Nothing,Number}=nothing,
    num_frags::Union{Nothing,Number}=nothing,
    sep::Union{Nothing,AbstractString}=nothing,
)
    args = ["SUMMARIZE"]

    if !isempty(fields)
        push!(args, "LEN", string(length(fields)), fields...)
    end

    if !isnothing(context_len)
        push!(args, "LEN", string(context_len))
    end

    if !isnothing(num_frags)
        push!(args, "FRAGS", string(num_frags))
    end

    if !isnothing(sep)
        push!(args, "SPEARATOR", sep)
    end

    q.summarize_fields = args

    return q
end

"""
    highlight!(
        q::AbstractQuery; fields::AbstractArray=nothing, tags::AbstractArray=nothing
    )

    Highlighting will highlight the found term (and its variants) with a user-defined tag.
    This may be used to display the matched text in a different typeface using a markup
    language, or to otherwise make the text appear differently.

# Arguments
- `q::AbstractQuery`: Query object.

# KeyWords
- `fields::AbstractArray`: List of fields
- `tags::AbstractArray`: User Tag objects.
"""
function highlight!(
    q; fields::AbstractArray=Vector{AbstractFilter}(), tags::AbstractArray=[]
)
    args = ["HIGHLIGHT"]

    if !isempty(fields)
        push!(args, "LEN", string(length(fields)), fields...)
    end

    if !isnothing(tags) && !isempty(tags)
        push!(args, "TAGS", tags...)
    end

    q.highlight_fields = args
    return q
end

"""
    get_args(q::AbstractQuery) -> Vector{Union{String,Number}}

For a query object, compile all the redis arguemnts that will be used when searching in
RediSearch.

# Arguments
- `q::AbstractQuery`: Query object.
"""
function get_args(q::AbstractQuery)
    args = Vector{Union{String,Number}}()
    push!(args, q.query_string)

    if q.no_content
        push!(args, "NOCONTENT")
    end

    if !isempty(q.fields)
        push!(args, "INFIELDS", string(length(q.fields)), q.fields...)
    end

    if q.verbatim
        push!(args, "VERBATIM")
    end

    if q.no_stopwords
        push!(args, "NOSTOPWORDS")
    end

    if !isempty(q.filters)
        for filter in q.filters
            if filter <: AbstractFilter
                push!(args, filter.args...)
            end
        end
    end

    if q.with_payloads
        push!(args, "WITHPAYLOADS")
    end

    if !isnothing(q.scorer)
        push!(args, "SCORER", self.scorer)
    end

    if q.with_scores
        push!(args, "WITHSCORES")
    end

    if !isempty(q.ids)
        push!(args, "INKEYS", string(length(q.ids), q.ids...))
    end

    if q.slop >= 0
        push!(args, "SLOP", q.slop)
    end

    if q.in_order
        push!(args, "INORDER")
    end

    if !isempty(q.return_fields)
        push!(args, "RETURN", string(length(q.return_fields), q.return_fields...))
    end

    if !isnothing(q.sortby)
        if !isempty(q.sortby.args)
            push!(args, "SORTBY", q.sortby.args...)
        end
    end

    if !isnothing(q.language)
        push!(args, "LANGUAGE", q.language)
    end

    push!(args, q.summarize_fields..., q.highlight_fields...)
    push!(args, "LIMIT", q.offset, q.num)

    return args
end


text:
for (T, funcpairs) in (
    (
        MagLike,
        (
            (:(Base.isfinite), :is_finite),
            (:(Base.isinf), :is_inf),
            (:isspecial, :is_special),
            (:(Base.iszero), :is_zero),
        ),
    ),
    (
        ArfLike,
        (
            (:(Base.isfinite), :is_finite),
            (:(Base.isinf), :is_inf),
            (:(Base.isinteger), :is_int),
            (:(Base.isnan), :is_nan),
            (:isneginf, :is_neg_inf),
            (:isnormal, :is_normal),
            (:(Base.isone), :is_one),
            (:isposinf, :is_pos_inf),
            (:isspecial, :is_special),
            (:(Base.iszero), :is_zero),
        ),
    ),
    (
        ArbLike,
        (
            (:isexact, :is_exact),
            (:(Base.isfinite), :is_finite),
            (:(Base.isinteger), :is_int),
            (:isnegative, :is_negative),
            (:isnonnegative, :is_nonnegative),
            (:isnonpositive, :is_nonpositive),
            (:isnonzero, :is_nonzero),
            (:(Base.isone), :is_one),
            (:ispositive, :is_positive),
            (:(Base.iszero), :is_zero),
        ),
    ),
    (
        AcbLike,
        (
            (:isexact, :is_exact),
            (:(Base.isfinite), :is_finite),
            (:(Base.isinteger), :is_int),
            (:(Base.isone), :is_one),
            (:(Base.isreal), :is_real),
            (:(Base.iszero), :is_zero),
            (:containszero, :contains_zero),
            (:containsint, :contains_int),
        ),
    ),
    (ArbVectorLike, ((:(Base.isfinite), :is_finite), (:(Base.iszero), :is_zero))),
    (
        AcbVectorLike,
        (
            (:(Base.isfinite), :is_finite),
            (:(Base.isreal), :is_real),
            (:(Base.iszero), :is_zero),
        ),
    ),
    (ArbPoly, ((:(Base.isone), :is_one), (:(isx), :is_x), (:(Base.iszero), :is_zero))),
    (ArbSeries, ((:(Base.isone), :is_one), (:(isx), :is_x), (:(Base.iszero), :is_zero))),
    (
        AcbPoly,
        (
            (:(Base.isone), :is_one),
            (:(isx), :is_x),
            (:(Base.iszero), :is_zero),
            (:(Base.isreal), :is_real),
        ),
    ),
    (
        AcbSeries,
        (
            (:(Base.isone), :is_one),
            (:(isx), :is_x),
            (:(Base.iszero), :is_zero),
            (:(Base.isreal), :is_real),
        ),
    ),
    (
        ArbMatrixLike,
        (
            (:isexact, :is_exact),
            (:(Base.isfinite), :is_finite),
            (:(Base.isone), :is_one),
            (:(Base.iszero), :is_zero),
        ),
    ),
    (
        AcbMatrixLike,
        (
            (:isexact, :is_exact),
            (:(Base.isfinite), :is_finite),
            (:(Base.isone), :is_one),
            (:(Base.isreal), :is_real),
            (:(Base.iszero), :is_zero),
        ),
    ),
)
    for (jlf, arbf) in funcpairs
        @eval $jlf(x::$T) = !iszero($arbf(x))
    end
end

for ArbT in (ArfLike, ArbLike, AcbLike)
    @eval begin
        Base.isequal(y::$ArbT, x::$ArbT) = !iszero(equal(x, y))
        # Comparison of non-floating point values should use ==
        Base.:(==)(y::Integer, x::$ArbT) = !iszero(equal(x, y))
        Base.:(==)(x::$ArbT, y::Integer) = !iszero(equal(x, y))
    end
end
Base.:(==)(x::MagLike, y::MagLike) = !iszero(equal(x, y))
Base.isless(x::MagLike, y::MagLike) = cmp(x, y) < 0
Base.:(<)(x::MagLike, y::MagLike) = cmp(x, y) < 0
Base.:(<=)(x::MagLike, y::MagLike) = cmp(x, y) <= 0

for jltype in (Arf, Integer, Unsigned, Base.GMP.CdoubleMax)
    @eval begin
        Base.isless(x::ArfLike, y::$jltype) = (isnan(y) && !isnan(x)) || cmp(x, y) < 0
        Base.:(<)(x::ArfLike, y::$jltype) = !isnan(x) && !isnan(y) && cmp(x, y) < 0
        Base.:(<=)(x::ArfLike, y::$jltype) = (x < y) || isequal(x, y)
    end
end

for (ArbT, args) in (
    (ArfLike, ((:(==), :equal),)),
    (
        ArbLike,
        ((:(==), :eq), (:(!=), :ne), (:(<), :lt), (:(<=), :le), (:(>), :gt), (:(>=), :ge)),
    ),
    (AcbLike, ((:(==), :eq), (:(!=), :ne))),
    (ArbMatrixLike, ((:(==), :eq), (:(!=), :ne))),
    (AcbMatrixLike, ((:(==), :eq), (:(!=), :ne))),
)
    for (jlf, arbf) in args
        @eval begin
            Base.$jlf(x::$ArbT, y::$ArbT) = !iszero($arbf(x, y))
        end
    end
end

Base.isequal(x::T, y::T) where {T<:Union{ArbPoly,AcbPoly}} = !iszero(equal(x, y))
Base.isequal(x::T, y::T) where {T<:Union{ArbSeries,AcbSeries}} =
    degree(x) == degree(y) && !iszero(equal(x, y))

function Base.:(==)(x::T, y::T) where {T<:Union{ArbPoly,ArbSeries,AcbPoly,AcbSeries}}
    degree(x) == degree(y) || return false
    for i = 0:degree(x)
        x[i] == y[i] || return false
    end
    return true
end

function Base.:(!=)(x::T, y::T) where {T<:Union{ArbPoly,AcbPoly}}
    return iszero(overlaps(x, y))
end

function Base.:(!=)(x::T, y::T) where {T<:Union{ArbSeries,AcbSeries}}
    return (degree(x) != degree(y)) || iszero(overlaps(x, y))
end


text:
@testset "EDF.Signal label handling" begin
    signal_names = ["eeg", "eog", "test"]
    canonical_names = OndaEDF.STANDARD_LABELS[["eeg"]]
    @test OndaEDF.match_edf_label("EEG C3-(M1 +A2)/2 - rEf3", signal_names, "c3", canonical_names) == "c3-m1_plus_a2_over_2"
    @test OndaEDF.match_edf_label("EEG C3-(M1 +A2)/2 - rEf3", ["ecg"], "c3", canonical_names) == nothing
    @test OndaEDF.match_edf_label("EEG C3-(M1 +A2)/2 - rEf3", signal_names, "c4", canonical_names) == nothing
    @test OndaEDF.match_edf_label(" TEsT   -Fpz  -REF-cpz", signal_names, "fpz", canonical_names) == "-fpz-ref-cpz"
    @test OndaEDF.match_edf_label(" TEsT   -Fpz  -REF-cpz", signal_names, "fp", canonical_names) == nothing
    @test OndaEDF.match_edf_label("  -Fpz  -REF-cpz", signal_names, "fpz", canonical_names) == "-fpz-ref-cpz"
    @test OndaEDF.match_edf_label("EOG L", signal_names, "left", OndaEDF.STANDARD_LABELS[["eog", "eeg"]]) == "left"
    @test OndaEDF.match_edf_label("EOG R", signal_names, "right", OndaEDF.STANDARD_LABELS[["eog", "eeg"]]) == "right"
    for (signal_names, channel_names) in OndaEDF.STANDARD_LABELS
        for channel_name in channel_names
            name = channel_name isa Pair ? first(channel_name) : channel_name
            x = OndaEDF.match_edf_label("EKGR", signal_names, name, channel_names)
            @test name == "avr" ? x == "avr" : x == nothing
        end
    end
    @test OndaEDF.export_edf_label("eeg", "t4") == "EEG T4-Ref"
    @test OndaEDF.export_edf_label("eeg", "t4-a1") == "EEG T4-A1"
    @test OndaEDF.export_edf_label("emg", "lat") == "EMG LAT"
end


text:
abstract type AbstractGaussNewton{T} <: Direction{T} end

function objective!(GN::AbstractGaussNewton, y::AbstractVector, x::AbstractVector)
	norm(GN.f(y, x))
end
# WARNING overwrites value vector
function objective!(GN::AbstractGaussNewton, x::AbstractVector)
	y = DiffResults.value(GN.r)
    norm(GN.f(y, x))
end
function objective(GN::AbstractGaussNewton, x::AbstractVector)
	y = similar(DiffResults.value(GN.r))
	norm(GN.f(y, x))
end
function get_value_jacobian(GN::AbstractGaussNewton)
    return DiffResults.value(GN.r), DiffResults.jacobian(GN.r)
end

function update_jacobian!(GN::AbstractGaussNewton, x::AbstractVector)
    y, J = get_value_jacobian(GN)
    jacobian!(GN.r, GN.f, y, x, GN.cfg) # constant allocation
    return y, J
end

# update_jacobian! needs to be called prior to update_direction for correctness
function update_direction!(GN::AbstractGaussNewton, λ::Real, min_diagonal::Real = 1e-12)
    !(isinf(λ) || isnan(λ)) || throw(DomainError("$λ: has to be finite"))
    y, J = get_value_jacobian(GN)
    mul!(GN.d, J', y)
    mul!(GN.JJ, J', J)
    lm_scaling!(GN.JJ, λ, min_diagonal)
    C = cholesky!(GN.JJ) # constant allocation
    ldiv!(C, GN.d)
    GN.d .*= -1
    return y, GN.d, J, C
end

function lm_scaling!(JJ::AbstractMatrix, λ::Real, min_diagonal::Real = 1e-12)
    for i in diagind(JJ)
        JJ[i] = max(JJ[i] * (1+λ), min_diagonal)
    end
    JJ
end

################################################################################
# Gauss-Newton algorithm, minimizes sum of squares of vector-valued function f w.r.t. x
struct GaussNewton{T, F, D, R, C, JT} <: AbstractGaussNewton{T}
    f::F
    d::D # pre-allocation for parameter direction vector
    # y::Y # pre-allocation for output / residual vector
    r::R # diff results
    cfg::C # Jacobian config
    JJ::JT # pre-allocation for Hessian approximation
    function GaussNewton(f::F, x::U, y::V) where {T, U<:AbstractVector{T},
                                                V<:AbstractVector, F<:Function}
        d = similar(x)
        r = DiffResults.JacobianResult(y, x)
        cfg = JacobianConfig(f, y, x)
        J = DiffResults.jacobian(r)
        JJ = J'J
        new{T, F, U, typeof(r), typeof(cfg), typeof(JJ)}(f, d, r, cfg, JJ)
    end
end

# using qr factorization instead of cholesky of normal equations
function valdir(GN::GaussNewton, x::AbstractVector, qr::Val{true})
    y, J = update_jacobian!(GN, x)
    ldiv!(GN.d, qr!(J), y)
    GN.d .*= -1
    return norm(y), GN.d
end

function valdir(GN::GaussNewton, x::AbstractVector, min_diagonal::Real = 1e-12)
    update_jacobian!(GN, x)
    y, d, _, _ = update_direction!(GN, 0, min_diagonal)
	return norm(y), GN.d
end

################################################################################
# TODO: geodesic acceleration (see https://en.wikipedia.org/wiki/Levenberg–Marquardt_algorithm)
struct LevenbergMarquart{T, F, D, R, C, JT} <: AbstractGaussNewton{T}
    f::F
    d::D # pre-allocation for parameter direction vector
    r::R # DiffResult
    cfg::C # Jacobian config
    JJ::JT # pre-allocation for Hessian approximation
    function LevenbergMarquart(f::F, x::U, y::V) where {T, U<:AbstractVector{T},
                                                V<:AbstractVector, F<:Function}
		d = similar(x)
        r = DiffResults.JacobianResult(y, x)
        cfg = JacobianConfig(f, y, x)
        J = DiffResults.jacobian(r)
        JJ = J'J
        new{T, F, U, typeof(r), typeof(cfg), typeof(JJ)}(f, d, r, cfg, JJ)
    end
end

function valdir(LM::LevenbergMarquart, x::AbstractVector, λ::Real, min_diagonal::Real = 1e-12)
    update_jacobian!(LM, x)
    update_direction!(LM, λ, min_diagonal)
end

################################################################################
# settings for LM optimization algorithm
struct LevenbergMarquartSettings{T<:Real}
	max_iter::Int # maximum number of iterations
	min_iter::Int # see min_decrease below
	max_backtrack::Int # maximum number of consecutive backtracks (λ increases)

	min_resnorm::T # minimum norm of residual vector for termination
	min_res::T # elementwise minimum of residual vector for termination

	min_decrease::T # minimum decrease in residual norm for termination
	 # ... after at least miniter iterations have been executed

	max_step::T # maximum norm of step proposal
	min_diagonal::T # minimum diagonal entries of Hessian approximation

	min_λ::T
	max_λ::T

	increase_factor::T # factor by which λ is increased
	decrease_factor::T
end

function LevenbergMarquartSettings(;max_iter::Int = 128, max_backtrack::Int = 32,
	min_resnorm::Real = 1e-6, min_res::Real = 1e-6,
	min_decrease::Real = 1e-12, min_iter::Int = 1,
	max_step::Real = Inf, min_diagonal::Real = 1e-12,
	min_lambda::Real = 1e-12, max_lambda::Real = 1e12,
	increase_factor::Real = 10, decrease_factor::Real = 7)

	LevenbergMarquartSettings(max_iter, min_iter, max_backtrack,
								promote(min_resnorm, min_res,
								min_decrease, max_step, min_diagonal,
								min_lambda, max_lambda,
								increase_factor, decrease_factor)...)
end

################################################################################

function optimize!(LM::LevenbergMarquart{<:Real},
			x::AbstractVector, y::AbstractVector, geodesic::Val{false},
			stn::LevenbergMarquartSettings = LevenbergMarquartSettings(),
			λ::Real = 1e-6, verbose::Union{Val{true}, Val{false}} = Val(false)
			)
    oldx = copy(x)
	val, newval = objective!(LM, y, x), Inf
    for i in 1:stn.max_iter
		update_jacobian!(LM, oldx)
		newval, λ = lm_backtrack!(LM, x, oldx, y, λ, stn)
		verbose isa Val{false} || lm_verbose(i, newval, val, λ, y)
		if converged(stn, newval, val, y, i)
			return x, i
		end
		val = newval
        copy!(oldx, x)
	end
	return x, stn.max_iter
end

function optimize!(LM::LevenbergMarquart{<:Real}, x::AbstractVector, y::AbstractVector, geodesic::Val{true},
	stn::LevenbergMarquartSettings = LevenbergMarquartSettings(),
	λ::Real = 1e-6, h::Float64=0.01, α::Float64=0.75, verbose::Union{Val{true}, Val{false}} = Val(false))

    oldx = copy(x)
	val, newval = objective!(LM, y, x), Inf
	temp = zero(y)
	temp2 = zero(y)
	r = zero(y)
	dx2 = zero(x)

    for i in 1:stn.max_iter
		update_jacobian!(LM, oldx)
		newval, λ = lm_geo_backtrack!(LM, x, oldx, y, temp, temp2, r, dx2, λ, h, α, stn)
		verbose isa Val{false} || lm_verbose(i, newval, val, λ, y)
		if converged(stn, newval, val, y, i)
			return x, i
		end
		val = newval
        copy!(oldx, x)
	end
	return x, stn.max_iter
end

function lm_backtrack!(LM::LevenbergMarquart,
						x::AbstractVector, oldx::AbstractVector, y::AbstractVector,
						λ::Real, stn::LevenbergMarquartSettings)
	val = 0.
	for j in 1:stn.max_backtrack
		w, dx, _, _ = update_direction!(LM, λ, stn.min_diagonal)
		val = norm(w)
		@. x = oldx + dx
		newval = objective!(LM, y, x) # norm(LM.f(y, x))
		if !(newval ≥ val || isnan(newval) || any(x->abs(x) > stn.max_step, dx))
			λ = max(λ / stn.decrease_factor, stn.min_λ)
			return newval, λ
		end
		λ = min(λ * stn.increase_factor, stn.max_λ)
	end
	copy!(x, oldx) # if backtracking wasn't successful
	return val, λ
end

function lm_geo_backtrack!(LM::LevenbergMarquart,
	x::AbstractVector, oldx::AbstractVector, y::AbstractVector,
	temp::AbstractVector, temp2::AbstractVector, r::AbstractVector, dx2::AbstractVector,
	λ::Real, h::Float64=0.01, α::Float64=0.75, stn::LevenbergMarquartSettings = LevenbergMarquartSettings())

	val = 0.
	for j in 1:stn.max_backtrack
		w, dx, J, C = update_direction!(LM, λ, stn.min_diagonal)
		val = norm(w)
		r .= 2/h * ( (LM.f(temp2, oldx + h.*dx)- w)/h - mul!(temp, J, dx) )
		mul!(dx2, J', r)
		ldiv!(C, dx2)
		if 2(norm(dx2))/norm(dx) < α
			@. x = oldx + dx - 0.5dx2
			newval = objective!(LM, y, x) # norm(LM.f(y, x))
			if !(newval ≥ val || isnan(newval) || any(x->abs(x) > stn.max_step, dx - 0.5dx2))
				λ = max(λ / stn.decrease_factor, stn.min_λ)
				return newval, λ
			end
		end
		λ = min(λ * stn.increase_factor, stn.max_λ)
	end
	copy!(x, oldx) #
end

function converged(stn::LevenbergMarquartSettings, newval::Real, val::Real,
						y::AbstractVector, i::Int)
	newval ≤ stn.min_resnorm && return true
	(abs(newval-val) < stn.min_decrease && i ≥ stn.min_iter) && return true
	all(x -> abs(x) < stn.min_res, y) && return true
	return false
end

function lm_verbose(i::Int, newval::Real, val::Real, λ::Real, y::AbstractVector)
	println("iteration $i:
		residual norm = $(norm(y)),
		maximum residual = $(maximum(y)),
		λ = $λ,
		improvement = $(val-newval)")
end


text:
using Test
using ShuffleProofs: Leaf, Tree, Node, encode, decode, Hash, PRG

h = Hash("sha256")

version = "3.0.4"
sid = "SessionID"
auxsid = "default"
nr = UInt32(100)
nv = UInt32(256)
ne = UInt32(256)

s_H = "SHA-256"
s_PRG = "SHA-256"
s_Gq = "ModPGroup(safe-prime modulus=2*order+1. order bit-length = 511)::00000000020100000020636f6d2e766572696669636174756d2e61726974686d2e4d6f645047726f757000000000040100000041009a91c3b704e382e0c772fa7cf0e5d6363edc53d156e841555702c5b6f906574204bf49a551b695bed292e0218337c0861ee649d2fe4039174514fe2c23c10f6701000000404d48e1db8271c17063b97d3e7872eb1b1f6e29e8ab7420aaab8162db7c832ba1025fa4d2a8db4adf69497010c19be0430f7324e97f201c8ba28a7f1611e087b3010000004100300763b0150525252e4989f51e33c4e6462091152ef2291e45699374a3aa8acea714ff30260338bddbb48fc7446b273aaada90e3ee8326f388b582ea8a073502010000000400000001"

data = (version, sid * "." * auxsid, nr, nv, ne, s_PRG, s_Gq, s_H)

tree = Tree(data)
binary = encode(tree)

ρ = "15e6c97600bbe30125cbc08598dcde01a769c15c8afe08fe5b7f5542533159e9"
@test bytes2hex(h(binary)) == ρ


text:
module ProblemTest
using LinearAlgebra
P1 = (sol = [[1,0], [0, 1]], vsol = -1,  A = [1 1], b = [1], c = [-1, -1], isoptimal = true, isbounded = true, isfeasible = true)
P2 = (sol = [[1,0], [0, 1]], vsol = -1, A = [1 1; 1 1], b = [1, 1], c = [-1, -1], isoptimal = true, isbounded = true, isfeasible = true)
function kmstar(n::Int, ::Type{T} = Int) where T
    xstar = zeros(T, 2*n)
    xstar[n] = 100^(n - 1)
    rest = [100^(k - 1) for k in 1:n]
    rest[end] = 0
    xstar[n+1:end] = rest
    return xstar
end
function km(n::Int, ::Type{T} = Int) where T
    c = [10^(n - i) for i in 1:n]
    b = [100^(i - 1) for i in 1:n]
    A = zeros(T, n, n)
    for i in 1:n
        A[i, 1:i] = 2*[10^(i - j) for j in 1:i]
    end
    A[:, :] -= I
    A = [A I]
    c = [c; zeros(T, n)]
    return A, b, -c, -100^(n - 1), kmstar(n)
end
PKM = [
    ((A, b, c, vstar, xstar) = km(i);
        (sol = [xstar],
        vsol = vstar,
            A = A,
            b = b,
            c = c,
            isoptimal = true,
            isbounded = true,
            isfeasible = true
        )
    )
    for i in 1:10
]
#Problem List
PL = [P1, P2, PKM...]
end

max_stars_repo_name: fbastin/Jasmin

max_issues_repo_name: fbastin/Jasmin

max_forks_repo_name: fbastin/Jasmin