text:
# Copyright (c) 2016, Jernej Kovacic
# Licensed under the MIT "Expat" License.
#
# See LICENSE.md for more details.

#
# Note: code in this file is based on the C++ library "Math", available at:
# https://github.com/jkovacic/math and
# https://github.com/jkovacic/math/blob/master/lib/specfun/SpecFunGeneric.cpp
#

include("ctdfrac.jl")
include("settings.jl")
#include("exception.jl")

#
# The following two publications are often referred in this file:
#
# - [Numerical Recipes] or shorter [NR]:
#      William H. Press, Saul A. Teukolsky, William T. Vetterling, Brian P. Flannery
#      Numerical Recipes, The Art of Scientific Computing, 3rd Edition,
#      Cambridge University Press, 2007
#      http://www.nr.com
#
# - [Abramowitz & Stegun] or shorter [AS]:
#      Milton Abramowitz, Irene A. Stegun
#      Handbook of Mathematical Functions with Formulas, Graphs and Mathematical Tables
#      National Bureau of Standards,
#      Applied Mathematics Series 55
#      Issued in June 1964
#
#      The book is copyright free and can be downloaded from:
#      http://www.cs.bham.ac.uk/~aps/research/projects/as/book.php
#

#     betac(x, y)
#
# Unlike Julia 0.5, Julia 0.4 does not support the beta function
# for complex arguments.
# To support regularized beta function on complex numbers as well,
# this convenience function has been introduced.
# TODO: the function may be removed when Julia 0.4 is no longer supported.
function betac{T<:AbstractFloat}(x::Complex{T}, y::Complex{T})

    #
    # It can be shown that B(x,y) can be expressed with gamma
    # functions:
    #
    #             G(x) * G(y)
    #   B(x,y) = -------------
    #               G(x+y)
    #

    const EPS2 = eps(T) * eps(T)
    const gxy = gamma(x + y);

    # handle a very unlikely occassion that G(x,y) gets very close to 0:
    if ( abs2(gxy) < EPS2 )
        throw(SpecfunUndefinedError)
    end

    return gamma(x) * gamma(y) / gxy
end

#     inc_gamma_ctdfrac_closure(a, x)
#
# Return a tupple of parameterized functions `fa(i)` and `fb(i)` that
# evaluate coefficients of the following continued fraction, necessary
# to evaluate the incomplete gamma function:
#
# When x > (1+a), the following continued fraction will be evaluated:
#
#                            1*(1-a)
#   cf = (x-a+1) - ---------------------------
#                                  2*(2-a)
#                    (x-a+3) - ---------------
#                               (x-a+5) - ...
#
# `a(i)` and `b(i)` depend on parameters `a` and `x` and
# are defined as follows:
#
# * a(x,i) = - i * (i-a)
# * b(x,i) = x - a + 1 + 2*i
#
# `a` and `x` are the returned functions' parameters.
#

function inc_gamma_ctdfrac_closure(a, x)
    return i::Integer -> (-i) * (i-a),
           i::Integer -> x - a + 1 + 2*i
end

#     inc_beta_ctdfrac_closure(a, b, x)
#
# Return a tupple of parameterized functions `fa(i)` and `fb(i)` that
# evaluate coefficients of the following continued fraction, necessary
# to evaluate the incomplete beta function:
#
#                     a1
#   cf = 1 + ---------------------
#                       a2
#             1 + ---------------
#                         a3
#                  1 + ---------
#                       1 + ...
# `a(i)` and `b(i)` depend on parameters `a`, `b` and `x` and
# are defined as follows:
#
# * b(x,i) = 1
# *
# *           /    (a+m) * (a+b+m) * x
# *           | - -------------------------   <== i = 2*m+1
# *           /    (a+2*m) * (a + 2*m + 1)
# * a(x,i) = {
# *           \       m * (b-m) * x
# *           |   ---------------------       <== i = 2*m
# *           \    (a+2*m-1) * (a+2*m)
#
# `a`, `b` and `x` are the returned functions' parameters.
#

function inc_beta_ctdfrac_closure(a, b, x)
    T = typeof(x)

    fa =
        function (i::Integer)
            m = div(i, 2)

            return ( 1 == rem(i, 2) ?
                     -((a+m) * (a+b+m) * x) / ((a+i-1) * (a+i)) :
                     (m * (b-m) * x) / ((a+i-1) * (a+i))  )
        end

    fb = i::Integer -> T(1)

    return fa, fb
end

#     incGamma(a, x, upper, reg)
#
# Evaluate an incomplete gamma function with real arguments.
# The exact kind of the returned value depends on parameters
# `upper` and `reg`.
#
# # Arguments:
# * `a`: parameter of the incomplete gamma function
# * `x`: second input argument, the integration limit
# * `upper::Bool`: should the upper (if 'true') or the lower (if 'false') inc. gamma function be returned
# * `reg::Bool`: if 'true', the regularized gamma function is returned, i.e. divided by `gamma(a)`
#
# Note that both `a` and `x` must be strictly greater than 0.
#
# The function throws `SpecfunUndefinedError` when the inc. gamma is
# not defined for any input argument, or `NoConvergenceError`
# when the internal algorithm does not converge.
#
function incGamma{T<:AbstractFloat}(a::T, x::T, upper::Bool, reg::Bool)

    const EPS = eps(T)
    # sanity check:
    if ( a < EPS || x < EPS )
        throw(SpecfunUndefinedError)
    end

    #
    # The algorithm for numerical approximation of the incomplete gamma function
    # as proposed by [Numerical Recipes], section 6.2:
    #
    # When x > (a+1), the upper gamma function can be evaluated as
    #
    #                 -x    a
    #                e   * x
    #   G(a,x) ~= --------------
    #                cf(a,x)
    #
    # where 'cf(a,x) is the continued fraction defined above, its coefficients
    # 'a(i)' and 'b(i)' are implemented in 'inc_gamma_ctdfrac_closure'.
    #
    # When x < (a+1), it is more convenient to apply the following Taylor series
    # that evaluates the lower incomplete gamma function:
    #
    #                          inf
    #                         -----
    #              -x    a    \        G(a)       i
    #   g(a,x) ~= e   * x  *   >    ---------- * x
    #                         /      G(a+1+i)
    #                         -----
    #                          i=0
    #
    # Applying the following property of the gamma function:
    #
    #   G(a+1) = a * G(a)
    #
    # The Taylor series above can be further simplified to:
    #
    #                          inf
    #                         -----              i
    #              -x    a    \                 x
    #   g(a,x) ~= e   * x  *   >    -------------------------
    #                         /      a * (a+1) * ... * (a+i)
    #                         -----
    #                          i=0
    #
    # Once either a lower or an upper incomplete gamma function is evaluated,
    # the other value may be quickly obtained by applying the following
    # property of the incomplete gamma function:
    #
    #   G(a,x) + g(a,x) = G(a)
    #
    # A value of a regularized incomplete gamma function is obtained
    # by dividing g(a,x) or G(a,x) by G(a).
    #

    # This factor is common to both algorithms described above:
    ginc = exp(-x) * (x^a)

    if ( x > (a + 1) )
        #
        # x > (a + 1)
        #
        # In this case evaluate the upper gamma function as described above.
        #

        fa, fb = inc_gamma_ctdfrac_closure(a, x)
        const G = ( true==upper && false==reg ? T(0) : gamma(a) )

        ginc /= ctdfrac_eval(fa, fb)

        #
        # Apply properties of the incomplete gamma function
        # if anything else except a generalized upper incomplete
        # gamma function is desired.
        #
        if ( false == upper )
            ginc = G - ginc
        end

        if ( true == reg )
            # Note: if a>0, gamma(a) is always greater than 0
            ginc /= G
        end
    else
        #
        # x < (a + 1)
        #
        # In this case evaluate the lower gamma function as described above.
        #
        const G = ( false==upper && false==reg ? T(0) : gamma(a) )

        # Initial term of the Taylor series at i=0:
        ginc /= a
        term = ginc

        # Proceed the Taylor series for i = 1, 2, 3... until it converges:
        const TOL = SPECFUN_ITER_TOL(T)
        at = a
        i = 1
        while ( abs(term) > TOL && i<SPECFUN_MAXITER )
            at += T(1)
            term *= (x / at)
            ginc += term
            i += 1
        end

        # has the series converged?
        if ( i >= SPECFUN_MAXITER )
            throw(NoConvergenceError)
        end

        #
        # Apply properties of the incomplete gamma function
        # if anything else except a generalized lower incomplete
        # gamma function is requested.
        #
        if ( true == upper )
            ginc = G - ginc
        end

        if ( true == reg )
            # Note: if a>0, gamma(a) is always greater than 0
            ginc /= G
        end
    end

    return ginc
end

#     incGamma(a, x, upper, reg)
#
# Evaluate an incomplete gamma function with complex arguments.
# The exact kind of the returned value depends on parameters
# `upper` and `reg`.
#
# # Arguments:
# * `a::Complex`: parameter of the incomplete gamma function
# * `x::Complex`: second input argument, the integration limit
# * `upper::Bool`: should the upper (if 'true') or the lower (if 'false') inc. gamma function be returned
# * `reg::Bool`: if 'true', the regularized gamma function is returned, i.e. divided by `gamma(a)`
#
# Note that unlike at real numbers, incomplete gamma function is defined
# virtually everywhere on the complex plane except at `a` = negative integer.
#
# The function throws `SpecfunUndefinedError` when the inc. gamma is
# not defined for any input argument, or `NoConvergenceError`
# when the internal algorithm does not converge.
#
function incGamma{T<:AbstractFloat}(a::Complex{T}, x::Complex{T}, upper::Bool, reg::Bool)

    #
    # The lower incomplete gamma function can be expanded into:
    #
    #                          inf
    #                         -----              i
    #              -x    a    \                 x
    #   g(a,x) ~= e   * x  *   >    -------------------------
    #                         /      a * (a+1) * ... * (a+i)
    #                         -----
    #                          i=0
    #
    # Once the lower incomplete gamma function is evaluated,
    # the upper one may be quickly obtained by applying the following
    # property of the incomplete gamma function:
    #
    #   G(a,x) + g(a,x) = G(a)
    #
    # A value of a regularized incomplete gamma function is obtained
    # by dividing g(a,x) or G(a,x) by G(a).
    #

    const EPS2 = eps(T) * eps(T)
    # If 'a' equals zero, division by zero would occur
    if ( abs2(a) < EPS2 )
        throw(SpecfunUndefinedError)
    end

    # G = gamma(a) or (0,0) when the value is not used
    const G = ( true==reg || true==upper ? gamma(a) : Complex{T}(0) )

    at = a

    # The first term of the series
    ginc = (x^a) * exp(-x) / a
    term = ginc

    # proceed the series until it converges
    const TOL = SPECFUN_ITER_TOL(T)
    const TOL2 = TOL * TOL
    i = 0
    while ( abs2(term) > TOL2 && i < SPECFUN_MAXITER )
        at += T(1)

        # if 'a' is a negative integer, sooner or later this exception will be thrown
        if ( abs2(at) < EPS2 )
            throw(SpecfunUndefinedError)
        end

        term *= x / at
        ginc += term

        i += 1
    end

    # has the series converged?
    if ( i >= SPECFUN_MAXITER )
        throw(NoConvergenceError)
    end

    #
    # Apply properties of the incomplete gamma function
    # if anything else except a non-regularized lower incomplete
    # gamma function is requested.
    #
    if ( true == upper )
        ginc = G - ginc
    end

    if ( true == reg )
        # very unlikely but check it anyway
        if ( abs2(G) < EPS2 )
            throw(SpecfunUndefinedError)
        end

        ginc /= G
    end

    return ginc
end

#     incBeta(a, b, x, lower, reg)
#
# Evaluate the incomplete beta function for real arguments. The exact
# kind of the returned value depends on parameters `upper` and `reg`.
#
# # Arguments
# * `a`: parameter `a` of the beta function
# * `b`: parameter `b` of the beta function
# * `x`: the integration limit
# * `lower::Bool`: should the lower (if 'true') or the upper (if 'false) inc. beta function be returned
# * `reg::Bool`: if 'true', the regularized beta function is returned, i.e. divided by `beta(a, b)`
#
# Note that `a`, `b` and `x` must be strictly greater than 0,
# `x` must be greater than 0 and less than 1
#
# The function throws `SpecfunUndefinedError` when the inc. gamma is
# not defined for any input argument, or `NoConvergenceError`
# when the internal algorithm does not converge.
#

function incBeta{T<:AbstractFloat}(a::T, b::T, x::T, lower::Bool, reg::Bool)

    const EPS = eps(T)

    # sanity check
    if ( a < EPS || b < EPS || x < T(0) || x > T(1) )
        throw(SpecfunUndefinedError)
    end

    #
    # The algorithm is described in detail in [Numerical Recipes], section 6.4
    #
    # If x < (a+1)/(a+b+2), the incomplete beta function will be evaluated as:
    #
    #                  a        b
    #                 x  * (1-x)
    #   Bx(a,b) ~= ------------------
    #                a * cf(x,a,b)
    #
    # where cf(x,a,b) is the continued fraction defined above, its
    # coefficients are implemented in 'inc_beta_ctdfrac_closure'.
    #
    # If x > (a+1)/(a+b+2), the algorithm will converge faster if the following
    # property of the incomplete beta function is applied:
    #
    #   B (a, b) + B   (b, a) = B(a, b)
    #    x          1-x
    #
    # and B_1-x(b,a) can be evaluated as described above.
    #
    # Once either a lower or an upper incomplete beta function is evaluated,
    # the other value may be quickly obtained by applying the following
    # property of the incomplete beta function:
    #
    #   Bx(a,b) + bx(a,b) = B(a,b)
    #
    # A value of a regularized incomplete beta function is obtained
    # by dividing Bx(a,b) or bx(a,b) by B(a,b).
    #

    # checks if x is "large", i.e. greater than the threshold defined above:
    const xlarge = (x > ( (a+T(1)) / (a+b+T(2)) ) )

    # If x is large, the parameters must be swapped
    const xn = ( false==xlarge ? x :  T(1)-x )
    const an = ( false==xlarge ? a : b )
    const bn = ( false==xlarge ? b : a )

    # B = B(a,b) or 0 when the value is not used
    const B = (
        true==reg || (false==xlarge && false==lower) || (true==xlarge  && true==lower) ?
        beta(a, b) : T(0) )

    binc = T(0)

    if ( abs(xn) < EPS )
        # Both boundary conditions are handled separately:
        #   B0(a,b) = 0  and  B1(a,b) = B(a,b)

        binc = ( true==xlarge ? B : T(0) )
    else
        # 'x' is somewhere between 0 and 1, apply the algorithm described above

        fa, fb = inc_beta_ctdfrac_closure(an, bn, xn)
        binc = (xn^an) * ((T(1)-xn)^bn) / an
        binc /= ctdfrac_eval(fa, fb)
    end

    #
    # When x is "large", the algorithm actually returns the
    # upper incomplete beta function!
    #
    # Depending on the requested result ('lower') adjust
    # 'binc' if necessary
    #
    if ( (false==xlarge && false==lower) ||
         (true==xlarge  && true==lower) )
        binc = B - binc
    end

    # Finally regularize the result if requested (via 'reg')
    if ( true == reg )
        # Just in case handle the very unlikely case
        if ( abs(B) < EPS )
            throw(SpecfunUndefinedError)
        end

        binc /= B
    end

    return binc
end

#     incBeta(a, b, x, lower, reg)
#
# Evaluate the incomplete beta function for real arguments. The exact
# kind of the returned value depends on parameters `upper` and `reg`.
#
# # Arguments
# * `a`: parameter `a` of the beta function
# * `b`: parameter `b` of the beta function
# * `x`: the integration limit
# * `lower::Bool`: should the lower (if 'true') or the upper (if 'false) inc. beta function be returned
# * `reg::Bool`: if 'true', the regularized beta function is returned, i.e. divided by `beta(a, b)`
#
# Note that `a`, `b` and `x` must be strictly greater than 0,
# `x` must be greater than 0 and less than 1
#
# The function throws `SpecfunUndefinedError` when the inc. gamma is
# not defined for any input argument, or `NoConvergenceError`
# when the internal algorithm does not converge.
#

#     incBeta(a, b, x, lower, reg)
#
# Evaluate the incomplete beta function for real arguments. The exact
# kind of the returned value depends on parameters `upper` and `reg`.
#
# Arguments
# * `a::Complex`: parameter `a` of the beta function
# * `b::Complex`: parameter `b` of the beta function
# * `x::Complex`: the integration limit
# * `lower::Bool`: should the lower (if 'true') or the upper (if 'false) inc. beta function be returned
# * `reg::Bool`: if 'true', the regularized beta function is returned, i.e. divided by `beta(a, b)`
#
# Note that unlike at real numbers, incomplete beta function is defined
# virtually everywhere except where 'a' is a negative integer, currently this
# function additionally requires that /x/ < 1 otherwise it may not converge.
#
# The function throws `SpecfunUndefinedError` when the inc. gamma is
# not defined for any input argument, or `NoConvergenceError`
# when the internal algorithm does not converge.
#
function incBeta{T<:AbstractFloat}(a::Complex{T}, b::Complex{T}, x::Complex{T}, lower::Bool, reg::Bool)

    #
    # Series expansion of the incomplete beta function as explained at:
    # http://functions.wolfram.com/GammaBetaErf/Beta3/06/01/03/01/01/
    #
    #                a    /                                   2        \
    #               x     |      a*(1-b)*x     a*(1-b)*(2-b)*x         |
    #   Bx(a,b) ~= ---- * | 1 + ----------- + ------------------ + ... | =
    #               a     |        (a+1)           2*(a+2)             |
    #                     \                                            /
    #
    #
    #                     /       inf                                         \
    #                a    |      -----                                     i  |
    #               x     |      \      a * (1-b) * (2-b) * ... * (i-b) * x   |
    #            = ---- * | 1 +   >    -------------------------------------- |
    #               a     |      /                  i! * (a+i)                |
    #                     |      -----                                        |
    #                     \       i=1                                         /
    #
    #

    const EPS2 = eps(T) * eps(T)

    # Division by zero will occur if a==0
    if ( abs2(a) < EPS2 )
        throw(SpecfunUndefinedError)
    end

    # Currently the function requires that /x/ < 1
    # when the series is guaranteed to converge
    if ( abs2(x) >= T(1) )
        throw(SpecfunUndefinedError)
    end

    # B = B(a,b) or (0,0) when the value is not used
    const B = ( true==reg || false==lower ? betac(a, b) : Complex{T}(0) )

    # The first term of the series
    binc = (x^a) / a
    term = binc * a

    at = a;
    bt = -b;

    # proceed the series until it converges
    const TOL = SPECFUN_ITER_TOL(T)
    const TOL2 = TOL * TOL
    i = 1
    while ( abs2(term) > TOL2 && i <= SPECFUN_MAXITER )
        at += T(1)
        bt += T(1)

        # if 'a' is a negative integer, sooner or later this exception will be thrown
        if ( abs2(at) < EPS2 )
            throw(SpecfunUndefinedError)
        end

        term *= x * bt / T(i)

        binc += term / at

        i += 1
    end

    # has the series converged?
    if ( i >= SPECFUN_MAXITER )
        throw(NoConvergenceError)
    end

    #
    # Apply properties of the incomplete beta function
    # if anything else except a non-regularized lower incomplete
    # beta function is desired.
    #
    if ( false == lower )
        binc = B - binc
    end

    if ( true == reg )
        # Very unlikely but check it anyway
        if ( abs2(B) < EPS2 )
            throw(SpecfunUndefinedError)
        end

        binc /= B
    end

    return binc
end

"""
    inc_gamma_upper(a, x)

Compute the upper incomplete gamma function, defined as:

           inf
            /
            |  a-1    -t
   G(a,x) = | t    * e   dt
            |
            /
            x

# Arguments:
* `a`: parameter of the incomplete gamma function
* `x`: second input argument, the lower integration limit

For real arguments, `a` and `x` must be greater than 0.
For complex arguments, the function is defined for virtually
any combination of `a` and `x` except when `a` is 0 or a
negative integer.

The function throws `SpecfunUndefinedError` when the function is
not defined for any input argument, or `NoConvergenceError`
when the internal algorithm does not converge.

# Examples:

```jldoctest
julia> inc_gamma_upper(2.0, 0.5)
0.9097959895689625

julia> inc_gamma_upper(2.0, 5.0)
0.040427681994512805

julia> inc_gamma_upper(2.0+im, 1.0-im)
0.294740622917453 + 1.1235452373416361im
```
"""
function inc_gamma_upper(a, x)
    return incGamma(a, x, true, false)
end

"""
    inc_gamma_lower(a, x)

Compute the lower incomplete gamma function, defined as:

            x
            /
            |  a-1    -t
   g(a,x) = | t    * e   dt
            |
            /
            0

# Arguments:
* `a`: parameter of the incomplete gamma function
* `x`: second input argument, the lower integration limit

For real arguments, `a` and `x` must be greater than 0.
For complex arguments, the function is defined for virtually
any combination of `a` and `x` except when `a` is 0 or a
negative integer.

The function throws `SpecfunUndefinedError` when the function is
not defined for any input argument, or `NoConvergenceError`
when the internal algorithm does not converge.

# Examples:

```jldoctest
julia> inc_gamma_lower(3.0, 2.5)
0.9123737682332481

julia> inc_gamma_lower(1.5, 4.0)
0.8454501129849447

julia> inc_gamma_lower(1.4-0.7im, 0.3+0.1im)
-0.016291214690464333 + 0.13052172604872464im
```
"""
function inc_gamma_lower(a, x)
    return incGamma(a, x, false, false)
end

"""
    inc_gamma_upper_reg(a, x)

Compute the regularized upper incomplete gamma function, defined as:

                  inf
                   /
              1    |  a-1    -t
   Q(a,x) = ------ | t    * e   dt
             G(a)  |
                   /
                   x

# Arguments:
* `a`: parameter of the incomplete gamma function
* `x`: second input argument, the lower integration limit

For real arguments, `a` and `x` must be greater than 0.
For complex arguments, the function is defined for virtually
any combination of `a` and `x` except when `a` is 0 or a
negative integer.

The function throws `SpecfunUndefinedError` when the function is
not defined for any input argument, or `NoConvergenceError`
when the internal algorithm does not converge.

# Examples:

```jldoctest
julia> inc_gamma_upper_reg(3.0, 0.5)
0.9856123220330416

julia> inc_gamma_upper_reg(-1.1+0.8im, 2.4-1.3im)
0.001971575402883159 + 0.016252537792016058im
```
"""
function inc_gamma_upper_reg(a, x)
    return incGamma(a, x, true, true)
end

"""
    inc_gamma_lower_reg(a, x)

Compute the regularized lower incomplete gamma function, defined as:

                   x
                   /
              1    |  a-1    -t
   P(a,x) = ------ | t    * e   dt
             G(a)  |
                   /
                   0

# Arguments:
* `a`: parameter of the incomplete gamma function
* `x`: second input argument, the lower integration limit

For real arguments, `a` and `x` must be greater than 0.
For complex arguments, the function is defined for virtually
any combination of `a` and `x` except when `a` is 0 or a
negative integer.

The function throws `SpecfunUndefinedError` when the function is
not defined for any input argument, or `NoConvergenceError`
when the internal algorithm does not converge.

# Examples:

```jldoctest
julia> inc_gamma_lower_reg(1.5, 3.2)
0.9063092095923061

julia> inc_gamma_lower_reg(-3.2-1.4im, -2.0-0.5im)
0.23971636700888468 + 0.4357835552823487im
```
"""
function inc_gamma_lower_reg(a, x)
    return incGamma(a, x, false, true)
end

"""
    inc_beta_lower(a, b, x)

Compute the lower incomplete beta function, defined as:

             x
             /
             |  a-1        b-1
   Bx(a,b) = | t    * (1-t)    dt
             |
             /
             0

# Arguments:
* `a`: parameter `a` of the beta function
* `b`: parameter `b` of the beta function
* `x`: integration limit

For real arguments, `a` and `b` must be greater than 0,
`x` must be greater than 0 and less than 1.
For complex arguments, `a` and `b` can be any complex values except that
`a` must not be 0 or a negative integer. Currently the `x` is restricted
to: /x/ < 1.

The function throws `SpecfunUndefinedError` when the function is
not defined for any input argument, or `NoConvergenceError`
when the internal algorithm does not converge.

# Examples:

```jldoctest
julia> inc_beta_lower(2.0, 5.0, 0.2)
0.01148800000000001

julia> inc_beta_lower(2.0, 5.0, 0.7)
0.032968833333333336

julia> inc_beta_lower(2.0-im, 3.0+im, 0.1+0.1im)
-0.014279979073458599 - 0.010562832806049189im
```
"""
function inc_beta_lower(a, b, x)
    return incBeta(a, b, x, true, false)
end

"""
    inc_beta_upper(a, b, x)

Compute the upper incomplete beta function, defined as:

             1
             /
             |  a-1        b-1
   bx(a,b) = | t    * (1-t)    dt = B(a,b) - Bx(a,b)
             |
             /
             x

# Arguments:
* `a`: parameter `a` of the beta function
* `b`: parameter `b` of the beta function
* `x`: integration limit

For real arguments, `a` and `b` must be greater than 0,
`x` must be greater than 0 and less than 1.
For complex arguments, `a` and `b` can be any complex values except that
`a` must not be 0 or a negative integer. Currently the `x` is restricted
to: /x/ < 1.

The function throws `SpecfunUndefinedError` when the function is
not defined for any input argument, or `NoConvergenceError`
when the internal algorithm does not converge.

# Examples:

```jldoctest
julia> inc_beta_upper(1.0, 2.0, 0.15)
0.36125

julia> inc_beta_upper(1.0, 2.0, 0.82)
0.01620000000000001

julia> inc_beta_upper(-1.2+im, -3.0-im, 0.1-0.2im)
35.013019655187826 + 12.573969337465748im
```
"""
function inc_beta_upper(a, b, x)
    return incBeta(a, b, x, false, false)
end

"""
    inc_beta_lower_reg(a, b, x)

Compute the regularized lower incomplete beta function, defined as:

                      x
                      /
                1     |  a-1        b-1
   Ix(a,b) = -------- | t    * (1-t)    dt
              B(a,b)  |
                      /
                      0

# Arguments:
* `a`: parameter `a` of the beta function
* `b`: parameter `b` of the beta function
* `x`: integration limit

For real arguments, `a` and `b` must be greater than 0,
`x` must be greater than 0 and less than 1.
For complex arguments, `a` and `b` can be any complex values except that
`a` must not be 0 or a negative integer. Currently the `x` is restricted
to: /x/ < 1.

The function throws `SpecfunUndefinedError` when the function is
not defined for any input argument, or `NoConvergenceError`
when the internal algorithm does not converge.

# Examples:

```jldoctest
julia> inc_beta_lower_reg(3.0, 2.0, 0.12)
0.00628992

julia> inc_beta_lower_reg(2.2-0.2im, 1.4+0.7im, -0.3+0.4im)
0.8221304290364263 - 0.40178180061491814im
```
"""
function inc_beta_lower_reg(a, b, x)
    return incBeta(a, b, x, true, true)
end

"""
    inc_beta_upper_reg(a, b, x)

Compute the regularized upper beta function, defined as:

                      1
                      /
                1     |  a-1        b-1
   ix(a,b) = -------- | t    * (1-t)    dt = 1 - Ix(a,b)
              B(a,b)  |
                      /
                      x

# Arguments:
* `a`: parameter `a` of the beta function
* `b`: parameter `b` of the beta function
* `x`: integration limit

For real arguments, `a` and `b` must be greater than 0,
`x` must be greater than 0 and less than 1.
For complex arguments, `a` and `b` can be any complex values except that
`a` must not be 0 or a negative integer. Currently the `x` is restricted
to: /x/ < 1.

The function throws `SpecfunUndefinedError` when the function is
not defined for any input argument, or `NoConvergenceError`
when the internal algorithm does not converge.

# Examples:

```jldoctest
julia> inc_beta_upper_reg(4.0, 1.0, 0.32)
0.98951424

julia> inc_beta_upper_reg(-3.8-0.7im, -3.4+2.7im, -0.6-0.1im)
0.0010989969710151117 + 0.0021282694661565148im
```
"""
function inc_beta_upper_reg(a, b, x)
    return incBeta(a, b, x, false, true);
end

#     as26_2_22(p)
#
# Implement the formula 26.2.22 in [Abramowitz & Stegun], i.e
# return such `x` that approximately satisfies Q(x) ~= p.
#
# # Arguments:
# * `p` - desired probability
#
# Unlike the original formula, this function accepts `p` being
# greater than 0.5.
#
# This "private" function expects that `p` is greater than 0
# and less than 1.
#
function as26_2_22(p)

    #
    # If 0 < p < 0.5:
    #
    #             +--------+
    #      -+    /      1               +------------+
    #   t =  \  /  ln -----   =    -+  / (-2) * ln(p)
    #         \/       p^2           \/
    #
    # Then initial approximation of 'x' can be calculated as:
    #
    #                 2.30753 + 0.27061 * t
    #   x ~= t - ---------------------------------
    #             1 + t * (0.99229 + 0.04481 * t)
    #
    # If p > 0.5, the expressions above are performed on its complementary
    # value (1-p) and the sign of the final 'x' is reversed.
    #

    T = typeof(p)

    const pp = ( p>= T(0.5) ? T(1) - p : p )

    const t = sqrt(T(-2) * log(pp) )

    # [Abramowitz & Stegun], section 26.2.22:
    const a0 = T(2.30753)
    const a1 = T(0.27061)
    const b0 = T(0.99229)
    const b1 = T(0.04481)

    x = t - (a0 + a1 * t) / (T(1) + t * (b0 + b1 * t))
    if ( p > T(0.5) )
        x = -x
    end

    return x
end

#     invIncGamma(a, g, upper, reg)
#
# Compute the inverse of the incomplete gamma function. The exact kind
# of the returned value depends on parameters `upper` and `reg`.
#
# # Arguments:
# * `a`: parameter of the incomplete gamma function
# * `g`: desired value of the incomplete gamma function
# * `upper::Bool`: if 'true', inverse of the upper incomplete gamma function will be evaluated
# * `reg::Bool`: if 'true', inverse of the regularized incomplete gamma function (divided by `gamma(a)`) will be evaluated
#
# Note that both `a` and `g` must be strictly greater than 0.
#
# The function is defined for real arguments only.
#
# The function throws `SpecfunUndefinedError` when the function is
# not defined for any input argument, or `NoConvergenceError`
# when the internal algorithm does not converge.
#
function invIncGamma{T<:AbstractFloat}(a::T, g::T, upper::Bool, reg::Bool)

    const EPS = eps(T)

    # sanity check
    if ( a < EPS || g < T(0) )
        throw(SpecfunUndefinedError)
    end

    #
    # The algorithm finds an inverse of the regularized lower
    # incomplete gamma function. If 'g' represents a result of
    # any other kind of the incomplete gamma function (specified)
    # by 'upper' and 'reg', it should be converted to a result of
    # the regularized lower incomplete gamma function by applying
    # well known properties of this function.
    #
    # The algorithm is proposed in [Numerical Recipes], section 6.2.1.
    #

    const G = gamma(a)

    p = g;
    if ( false == reg )
        p /= G
    end

    if ( true == upper )
        p = T(1) - p
    end

    # if 'p' equals 0 (or is very close to it),
    # the inverse inc. gamma function will also equal 0.
    if ( abs(p) < EPS )
        return T(0)
    end

    #
    # If 'p' is greater than 1, the incomplete gamma function
    # is not defined at all.
    # If 'p' equals 1 (or is very close to it), the result would
    # be infinity which is not supported by this algorithm.
    #
    if ( p >= (T(1) - EPS) )
        throw(SpecfunUndefinedError)
    end

    #
    # There is no direct series or continued fraction to evaluate
    # an inverse of the incomplete gamma function. On the other hand
    # there are known algorithms that estimate this value very well.
    # This algorithm implements approximation method proposed by the
    # [Numerical Recipes], section 6.2.1 and [Abramowitz & Stegun],
    # sections 26.2.22 and 26.4.17.
    #
    # Note: another method to estimate the initial value of the incomplete
    # gamma function is described in:
    #   Armido R. DiDonato, Alfred H. Morris
    #   Computation of the incomplete gamma function ratios and their inverse
    #   ACM Transactions on Mathematical Software
    #   Volume 12 Issue 4, Dec. 1986, pp. 377-393
    #   http://dl.acm.org/citation.cfm?id=23109&coll=portal&dl=ACM
    #
    # When an approximation is known, the Newton - Raphson method can be
    # applied to refine it further to the desired tolerance.
    #

    x = T(1)

    if ( a <= T(1) )
        #
        # a <= 1:
        #
        # First P(a,1) is approximated by:
        #
        #   Pa = P(a,1) ~= a* (0.253 + 0.12 * a)
        #
        # If 'p' is less than Pa's complement value, 'x' is
        # approximated as:
        #
        #         a      +------+
        #        --+    /   p
        #   x ~=    \  / ------
        #            \/   1-Pa
        #
        # If 'p' is greater than 1-Pa, 'x' is approximated as:
        #
        #                1 - p
        #   x ~= 1 - ln -------
        #                 Pa
        #
        # Note: as long as 0<a<1, Pa will never get even close to 1,
        #       hence the first expression is always defined.
        #

        const c1 = T(0.253)
        const c2 = T(0.12)
        const Pa = a * (c1 + c2 * a)
        const cpa = T(1) - Pa

        x = ( p < cpa ?
                 (p / cpa)^ (T(1) / a) :
                 T(1) - log((T(1) - p)/ Pa)  )
    else
        #
        # a > 1:
        #
        # The initial approximation of 'x' is obtained by the function
        # as26_2_22() that implements [Abramowitz & Stegun], section 26.2.22.
        #
        # Finally the approximated x0 is evaluated as:
        #
        #             /       1          x      \ 3
        #   x0 ~= a * | 1 - ----- + ----------- |
        #             \      9*a     3*sqrt(a)  /
        #

        x = -as26_2_22(p)

        # [Abramowitz & Stegun], section 26.4.17:
        x = T(1) - ( T(1) / (T(9) * a) ) + ( x / ( T(3) * sqrt(a) ) )

        # x = a * x^3:
        x  *= a * x * x
    end

    #
    # When the initial value of 'x' is obtained, a root finding
    # method is applied to determine the exact inverse. This algorithm
    # applies the slightly modified Newton - Raphson method:
    # if 'x' tries to go negative, its half value is assigned
    # to the new value.
    #

    #
    # The Newton - Raphson algorithm also requires the differentiation of
    # the regularized lower incomplete gamma function:
    #
    #                   a-1    -x
    #    d P(a,x)      x    * e
    #   ---------- = --------------
    #       dx         gamma(a)
    #
    # Verified by Maxima:
    # (%i1)  diff(gamma_incomplete_regularized(a, x), x);
    # (%o1)  (x^a-1*%e^-x)/gamma(a)
    #

    const TOL = SPECFUN_ITER_TOL(T)
    xn = T(0)
    f = inc_gamma_lower_reg(a, x) - p
    i = 0
    while ( abs(f) > TOL && i < SPECFUN_MAXITER )
        xn = x - f * G * exp(x) / (x^(a-T(1)))

        # x must not go negative!
        x = ( xn > EPS ? xn : x * T(0.5) )

        f = inc_gamma_lower_reg(a, x) - p

        i += 1
    end

    # Has the algorithm converged?
    if ( i >= SPECFUN_MAXITER )
        throw(NoConvergenceError)
    end

    return x
end

#     invIncBeta(a, b, y, lower, reg)
#
# Compute the inverse of the incomplete beta function. The exact kind
# of the returned value depends on parameters `lower` and `reg`.
#
# # Arguments:
# * `a`: parameter `a` of the incomplete beta function
# * `b`: parameter `b` of the incomplete beta function
# * `y`: desired value of the incomplete beta function
# * `lower::Bool`: if 'true', inverse of the lower incomplete beta function will be evaluated
# * `reg::Bool`: if 'true', inverse of the regularized incomplete beta function (divided by `beta(a,b)`) will be evaluated
#
# Note that `a`, `b` and `y` must be strictly greater than 0.
#
# The function is defined for real arguments only.
#
# The function throws `SpecfunUndefinedError` when the function is
# not defined for any input argument, or `NoConvergenceError`
# when the internal algorithm does not converge.
#
function invIncBeta{T<:AbstractFloat}(a::T, b::T, y::T, lower::Bool, reg::Bool)

    const EPS = eps(T)

    # sanity check
    if ( a < EPS || b < EPS || y < T(0) )
        throw(SpecfunUndefinedError)
    end

    #
    # The algorithm finds an inverse of the regularized lower
    # incomplete beta function. If 'y' represents a result of
    # any other kind of the incomplete beta function (specified
    # by 'lower' and 'reg'), it should be converted to a result of
    # the regularized lower incomplete beta function by applying
    # well known properties of this function.
    #
    # The algorithm is proposed in [Numerical Recipes], section 6.4.
    #

    const B = beta(a, b)

    p = y

    if ( false == reg )
        p /= B
    end

    if (false == lower)
        p = T(1) - p
    end

    # if 'p' equals 0 (or is very close to it),
    # the inverse inc. beta function will also equal 0.
    if ( abs(p) < EPS )
        return T(0)
    end

    #
    # If 'p' is greater than 1, the incomplete beta function
    # is not defined at all.
    # If 'p' equals 1 (or is very close to it), the result would
    # be infinity which is not supported by this algorithm.
    #
    if ( p >= (T(1) - EPS) )
        throw(SpecfunUndefinedError)
    end

    #
    # There is no known direct series or continued fraction to evaluate
    # an inverse of the incomplete beta function. On the other hand
    # there are known algorithms that estimate this value quite well.
    # This algorithm implements approximation method proposed by the
    # [Numerical Recipes], section 6.4 and [Abramowitz & Stegun],
    # section 26.5.22.
    #
    # When an approximation is known, the Newton - Raphson method can be
    # applied to refine it further to the desired tolerance.
    #

    x = T(0)

    if ( a>=T(1) && b>=T(1) )
        #
        # a >= 1  and  b >= 1:
        #
        # Initial 'x' is approximated as described in [Abramowitz & Stegun],
        # section 26.2.22 (implemented in a separate function).
        #
        # Then it is further refined as described in [Abramowitz & Stegun],
        # section 26.5.22:
        #
        #           1                   1
        #   ta = -------   and  tb = -------
        #         2*a-1               2*b-1
        #
        #              2
        #             x  - 3
        #   lambda = --------
        #               6
        #
        #           2
        #   h = ---------
        #        ta + tb
        #
        #        x * sqrt(h + lambda)                /           5      2   \
        #   w = ---------------------- - (tb - ta) * | lambda + --- - ----- |
        #                h                           \           6     3*h  /
        #
        #                   a
        #   x  ~=   ------------------
        #                       2*w
        #              a + b * e
        #

        # [Abramowitz & Stegun], section 26.2.22:
        x = as26_2_22(p)

        const lambda = (x*x - T(3)) / T(6)
        const ta = T(1) / ( T(2) * a - T(1) )
        const tb = T(1) / ( T(2) * b - T(1) )
        const h = T(2) / (ta + tb)

        const w = x * sqrt(h + lambda) / h - (tb - ta) *
                    (lambda + T(5) / T(6) - T(2) / (T(3) * h) )

        x = a / (a + b * exp(T(2) * w))
    else
        #
        # Either 'a' or 'b' or both are less than 1:
        #
        # In this case the initial 'x' is approximated as
        # described in [Numerical Recipes], section 6.4:
        #
        #
        #         1    /    a    \a              1    /    b    \b
        #   ta = --- * | ------- |    and  tb = --- * | ------- |
        #         a    \  a + b  /               b    \  a + b  /
        #
        #
        #   S = ta + tb
        #
        # If p < ta/S:
        #
        #         a    +----------+
        #   x ~=  -+  / p * S * a
        #           \/
        #
        # If p >= ta/S:
        #
        #        b    +----------+
        #   x ~= -+  / 1 - p*S*b
        #          \/
        #

        const ta = ( (a/(a+b))^a ) / a
        const tb = ( (b/(a+b))^b ) / b
        const S = ta + tb

        x = ( p < ta/S ?   (p*S*a)^(T(1)/a) :  (T(1)-p*S*b)^(T(1)/b) )
    end

    #
    # When the initial value of 'x' is obtained, a root finding
    # method is applied to determine the exact inverse. This algorithm
    # applies the slightly modified Newton - Raphson method:
    # if 'x' tries to go negative or beyond 1, it is bisected between
    # its current value and the upper/lower bondary.
    #

    #
    # The Newton - Raphson algorithm also requires the differentiation of
    # the regularized lower incomplete beta function:
    #
    #                   a-1        b-1
    #    d Ix(a,b)     x    * (1-x)
    #   ----------- = ------------------
    #      dx             beta(a,b)
    #
    # Verified by Maxima:
    # (%i1) diff(beta_incomplete_regularized(a, b, x), x);
    # (%o1) ((1-x)^b-1*x^a-1)/beta(a,b)
    #

    xn = T(0)
    f = inc_beta_lower_reg(a, b, x) - p

    const TOL = SPECFUN_ITER_TOL(T)
    i = 0
    while ( abs(f) > TOL && i < SPECFUN_MAXITER )
        xn = x - f * B / (x^(a-T(1)) * (T(1)-x)^(b-T(1)) )

        # x must not go negative or beyond 1!
        if ( xn < EPS )
            x *= T(0.5)
        elseif ( xn > (T(1)-EPS) )
            x = (T(1) + x) * T(0.5)
        else
            x = xn
        end

        f = inc_beta_lower_reg(a, b, x) - p

        i += 1
    end

    # Has the algorithm converged?
    if ( i >= SPECFUN_MAXITER )
        throw(NoConvergenceError)
    end

    return x
end

"""
    inc_gamma_lower_inv(a, g)

Compute the inverse of the lower incomplete gamma function,
i.e. return such `x` that satisfies:

    x
    /
    |  a-1    -t
    | t    * e   dt  =  g
    |
    /
    0

# Arguments:
* `a`: parameter of the incomplete gamma function
* `g`: desired value of the lower incomplete gamma function

`a` must be strictly greater than 0, `g` must be greater or equal
to 0 and less than `gamma(a)`.

The function is defined for real arguments only.

The function throws `SpecfunUndefinedError` when the function is
not defined for any input argument, or `NoConvergenceError`
when the internal algorithm does not converge.

# Examples:

```jldoctest
julia> inc_gamma_lower_inv(0.24, 0.94)
0.002024362637442518

julia> inc_gamma_lower_inv(2.8, 0.17)
0.9873918720287759
```
"""
function inc_gamma_lower_inv(a, g)
    return invIncGamma(a, g, false, false)
end

"""
    inc_gamma_upper_inv(a, g)

Compute the inverse of the upper incomplete gamma function,
i.e. return such `x` that satisfies:

   inf
    /
    |  a-1    -t
    | t    * e   dt  =  g
    |
    /
    x

# Arguments:
* `a`: parameter of the incomplete gamma function
* `g`: desired value of the upper incomplete gamma function

`a` must be strictly greater than 0, `g` must be greater or equal
to 0 and less than `gamma(a)`.

The function is defined for real arguments only.

The function throws `SpecfunUndefinedError` when the function is
not defined for any input argument, or `NoConvergenceError`
when the internal algorithm does not converge.

# Examples:

```jldoctest
julia> inc_gamma_upper_inv(0.65, 0.86)
0.21736651483074873

julia> inc_gamma_upper_inv(3.5, 0.43)
5.609013631489573
```
"""
function inc_gamma_upper_inv(a, g)
    return invIncGamma(a, g, true, false)
end

"""
    inc_gamma_lower_reg_inv(a, g)

Compute the inverse of the regularized lower incomplete gamma function,
i.e. return such `x` that satisfies:

          x
          /
     1    |  a-1    -t
   ------ | t    * e   dt  =  g
    G(a)  |
          /
          0

# Arguments:
* `a`: parameter of the incomplete gamma function
* `g`: desired value of the regularized lower incomplete gamma function

`a` must be strictly greater than 0, `g` must be greater or equal
to 0 and less than 1.

The function is defined for real arguments only.

The function throws `SpecfunUndefinedError` when the function is
not defined for any input argument, or `NoConvergenceError`
when the internal algorithm does not converge.

# Examples:

```jldoctest
julia> inc_gamma_lower_reg_inv(0.2, 0.3)
0.0015877907243440615

julia> inc_gamma_lower_reg_inv(3.0, 0.7)
3.6155676658661333
```
"""
function inc_gamma_lower_reg_inv(a, g)
    return invIncGamma(a, g, false, true)
end

"""
    inc_gamma_upper_reg_inv(a, g)

Compute the inverse of the upper incomplete gamma function,
i.e. return such `x` that satisfies:

         inf
          /
     1    |  a-1    -t
   ------ | t    * e   dt  =  g
    G(a)  |
          /
          x

# Arguments:
* `a`: parameter of the incomplete gamma function
* `g`: desired value of the regularized upper incomplete gamma function

`a` must be strictly greater than 0, `g` must be greater or equal
to 0 and less than 1.

The function is defined for real arguments only.

The function throws `SpecfunUndefinedError` when the function is
not defined for any input argument, or `NoConvergenceError`
when the internal algorithm does not converge.

# Examples:

```jldoctest
julia> inc_gamma_upper_reg_inv(0.3, 0.4)
0.1412525036310711

julia> inc_gamma_upper_reg_inv(5.2, 0.82)
3.1296773937114972
```
"""
function inc_gamma_upper_reg_inv(a, g)
    return invIncGamma(a, g, true, true)
end

"""
    inc_beta_lower_inv(a, b, y)

Compute the inverse of the lower incomplete beta function,
i.e. returns such `x` that satisfies:

             x
             /
             |  a-1        b-1
   Bx(a,b) = | t    * (1-t)    dt = y
             |
             /
             0

# Arguments:
* `a`: parameter `a` of the incomplete beta function
* `b`: parameter `b` of the incomplete beta function
* `y`: desired value of the lower incomplete beta function (`Bx(a,b)`)

`a` and `b` must be strictly greater than 0, `y` must be greater or equal
to 0 and less than `beta(a,b)`.

The function is defined for real arguments only.

The function throws `SpecfunUndefinedError` when the function is
not defined for any input argument, or `NoConvergenceError`
when the internal algorithm does not converge.

# Examples:

```jldoctest
julia> inc_beta_lower_inv(2.8, 0.3, 2.0)
0.9997335668051313

julia> inc_beta_lower_inv(1.1, 1.3, 0.4)
0.5190110112035177
```
"""
function inc_beta_lower_inv(a, b, y)
    return invIncBeta(a, b, y, true, false)
end

"""
    inc_beta_upper_inv(a, b, y)

Return the inverse of the upper incomplete beta function,
i.e. it returns such `x` that satisfies:

             1
             /
             |  a-1        b-1
   bx(a,b) = | t    * (1-t)    dt = y
             |
             /
             x

# Arguments:
* `a`: parameter `a` of the incomplete beta function
* `b`: parameter `b` of the incomplete beta function
* `y`: desired value of the upper incomplete beta function (`bx(a,b)`)

`a` and `b` must be strictly greater than 0, `y` must be greater or equal
to 0 and less than `beta(a,b)`.

The function is defined for real arguments only.

The function throws `SpecfunUndefinedError` when the function is
not defined for any input argument, or `NoConvergenceError`
when the internal algorithm does not converge.

# Examples:

```jldoctest
julia> inc_beta_upper_inv(0.4, 0.5, 1.8)
0.41086943388583214

julia> inc_beta_upper_inv(1.7, 1.1, 0.2)
0.7205525719752999
```
"""
function inc_beta_upper_inv(a, b, y)
    return invIncBeta(a, b, y, false, false)
end

"""
    inc_beta_lower_reg_inv(a, b, y)

Compute the inverse of the regularized lower incomplete beta function,
i.e. return such `x` that satisfies:

                       x
                       /
                1      |  a-1        b-1
   Ix(a,b) = --------  | t    * (1-t)    dt = y
              B(a,b)   |
                       /
                       0

# Arguments:
* `a`: parameter `a` of the incomplete beta function
* `b`: parameter `b` of the incomplete beta function
* `y`: desired value of the regularized lower incomplete beta function (`Ix(a,b)`)

`a` and `b` must be strictly greater than 0, `y` must be greater or equal
to 0 and less than 1.

The function is defined for real arguments only.

The function throws `SpecfunUndefinedError` when the function is
not defined for any input argument, or `NoConvergenceError`
when the internal algorithm does not converge.

# Examples:

```jldoctest
julia> inc_beta_lower_reg_inv(0.3, 0.2, 0.7)
0.9785534149620168

julia> inc_beta_lower_reg_inv(2.4, 3.5, 0.6)
0.44942828803639806
```
"""
function inc_beta_lower_reg_inv(a, b, y)
    return invIncBeta(a, b, y, true, true)
end

"""
    inc_beta_upper_reg_inv(a, b, y)

Compute the inverse of the regularized upper incomplete beta function,
i.e. return such `x` that satisfies:

                      1
                      /
                1     |  a-1        b-1
   ix(a,b) = -------  | t    * (1-t)    dt = y
              B(a,b)  |
                      /
                      x

# Arguments:
* `a`: parameter `a` of the incomplete beta function
* `b`: parameter `b` of the incomplete beta function
* `y`: desired value of the regularozed upper incomplete beta function (`ix(a,b)`)

`a` and `b` must be strictly greater than 0, `y` must be greater or equal
to 0 and less than 1.

The function is defined for real arguments only.

The function throws `SpecfunUndefinedError` when the function is
not defined for any input argument, or `NoConvergenceError`
when the internal algorithm does not converge.

# Examples:

```jldoctest
julia> inc_beta_upper_reg_inv(0.9, 1.5, 0.7)
0.18163313417789242

julia> inc_beta_upper_reg_inv(1.9, 2.7, 0.25)
0.5648299809468559
```
"""
function inc_beta_upper_reg_inv(a, b, y)
    return invIncBeta(a, b, y, false, true)
end


text:
function bern(t, order)
# construct bernstein polynomial of order ``order"
out = zeros(size(t,1),1)
for o = 0:order
   if o ==0
      out = b(t,order,o)
   else
      out = [out b(t,order,o)];
   end

end
return out
end

ext: jl


text:
__precompile__()

module FastGaussQuadrature

using Compat, SpecialFunctions

if VERSION ≥ v"0.7-"
    using LinearAlgebra
    cumprod(A::AbstractArray) = Base.cumprod(A, dims=1)
    cumprod(A::AbstractArray, d::Int) = Base.cumprod(A, dims=d)
    sum(A::AbstractArray, n::Int) = Base.sum(A, dims=n)
    sum(A) = Base.sum(A)
    flipdim(A, d) = reverse(A, dims=d)
else
    rmul!(A::AbstractArray, x::Number) = scale!(A, x)
    repeat(A, n, m) = repmat(A, n, m)
    repeat(A...; kwds...) = Base.repeat(A...; kwds...)
    eigen(A) = eig(A)
    const floatmax = realmax
end

export gausslegendre
export gausschebyshev
export gausslaguerre
export gausshermite
export gaussjacobi
export gausslobatto
export gaussradau
export besselroots

import SpecialFunctions: besselj, airyai, airyaiprime

include("gausslegendre.jl")
include("gausschebyshev.jl")
include("gausslaguerre.jl")
include("gausshermite.jl")
include("gaussjacobi.jl")
include("gausslobatto.jl")
include("gaussradau.jl")
include("besselroots.jl")

end


text:
"Causal Inductive Synthesis Corpus"
module CISC

f(x) = abs(x)

end


text:
import Observables
import Makie

export Twospinor
export update

"""
    Represents a two-component spinor.

fields: point, gauge1, gauge2, configuration, radius, segments1, segments2, color and observable.
"""
mutable struct Twospinor <: Sprite
    point::ComplexPlane
    gauge1::U1
    gauge2::U1
    configuration::Biquaternion
    radius::Float64
    segments1::Int
    segments2::Int
    color::Observables.Observable{Array{Makie.ColorTypes.RGBA{Float32},2}}
    observable::Tuple{Observables.Observable{Array{Float64,2}},
                      Observables.Observable{Array{Float64,2}},
                      Observables.Observable{Array{Float64,2}}}
end

"""
    Twospinor(scene, point; [rotation, [gauge1, [gauge2, [configuration, [radius, [segments1, [segments2, [color, [transparency]]]]]]]]])

Construct a 2-spinor with the given `scene`, `point` in the principal bundle,
`gauge1`, `gauge2`, `configuration`, `radius`, `segments1`, `segments2`, `color` and `transparency`.
"""
function Twospinor(scene::Makie.Scene,
                   point::ComplexPlane;
                   gauge1::U1 = U1(0),
                   gauge2::U1 = U1(2π),
                   configuration::Biquaternion = Biquaternion(ℝ³(0, 0, 0)),
                   radius::Float64 = 0.005,
                   segments1::Int = 36,
                   segments2::Int = 12,
                   color::Makie.RGBAf = Makie.RGBAf(1.0, 0.2705, 0.0, 0.5),
                   transparency::Bool = false)
    spinor = constructtwospinor(point, gauge1, gauge2, configuration, radius, segments1, segments2)
    colorarray = Observables.Observable(fill(color, size(spinor)...))
    observable = buildsurface(scene, spinor, colorarray, transparency = transparency)
    Twospinor(point, gauge1, gauge2, configuration, radius, segments1, segments2, colorarray, observable)
end

"""
    update(spinor, point)

Update a Twospinor by changing its observable with the given `spinor` and `point` in the principal bundle.
"""
function update(spinor::Twospinor, point::ComplexPlane)
    spinor.point = point
    value = constructtwospinor(spinor.point,
                               spinor.gauge1,
                               spinor.gauge2,
                               spinor.configuration,
                               spinor.radius,
                               spinor.segments1,
                               spinor.segments2)
    updatesurface(value, spinor.observable)
end

"""
    update(spinor, gauge1, gauge2)

Update a Twospinor by changing its observable with the given `spinor`, `gauge1` and `gauge2`.
"""
function update(spinor::Twospinor, gauge1::U1, gauge2::U1)
    spinor.gauge1 = gauge1
    spinor.gauge2 = gauge2
    value = constructtwospinor(spinor.point,
                               spinor.gauge1,
                               spinor.gauge2,
                               spinor.configuration,
                               spinor.radius,
                               spinor.segments1,
                               spinor.segments2)
    updatesurface(value, spinor.observable)
end

"""
    update(spinor, configuration)

Update a Twospinor by changing its observable with the given `spinor` and `configuration`.
"""
function update(spinor::Twospinor, configuration::Biquaternion)
    spinor.configuration = configuration
    value = constructtwospinor(spinor.point,
                               spinor.gauge1,
                               spinor.gauge2,
                               spinor.configuration,
                               spinor.radius,
                               spinor.segments1,
                               spinor.segments2)
    updatesurface(value, spinor.observable)
end

"""
    update(spinor, color)

Update a Twospinor by changing its observable with the given `spinor` and `color`.
"""
function update(spinor::Twospinor, color::Makie.RGBAf)
    spinor.color[] = fill(color, spinor.segments1, spinor.segments2)
end

"""
    update(spinor, point, gauge1, gauge2)

Update a Twospinor by constructing it from scratch with the given
`point` in the principal bundle, `gauge1` and `gauge2`.
"""
function update(spinor::Twospinor, point::ComplexPlane, gauge1::U1, gauge2::U1)
    spinor.point = point
    spinor.gauge1 = gauge1
    spinor.gauge2 = gauge2
    value = constructtwospinor(spinor.point,
                               spinor.gauge1,
                               spinor.gauge2,
                               spinor.configuration,
                               spinor.radius,
                               spinor.segments1,
                               spinor.segments2)
    updatesurface(value, spinor.observable)
end

max_stars_repo_name: iamazadi/geometrizer

max_issues_repo_name: iamazadi/Porta

max_forks_repo_name: iamazadi/Porta


text:
@inline cst_S(n, i, x) = binomial(n,i)*x^i*(1-x)^(n-i)
@inline cst_C(n1, n2, x) = x^n1*(1-x)^n2
@inline function cst_dS(n, i, x)
    if i == 0 # (1-x)^n
        -n*(1-x)^(n-1)
    elseif i == n # x^n
        n*x^(n-1)
    else
        binomial(n, i)*(i*x^(i-1)*(1-x)^(n-i)-(n-i)*x^i*(1-x)^(n-i-1))
    end
end

@inline function cst_dC(n1, n2, x)
    n1*x^(n1-1)*(1-x)^n2 - n2*x^n1*(1-x)^(n2-1)
end
function fit_cst(n::Int, dz::T,
    xu::AbstractVector{T},yu::AbstractVector{T},
    xl::AbstractVector{T},yl::AbstractVector{T},
    opt::OptimizerFactory;
    c1::T = 0.5,
    c2::T = 1.0 ) where {T}
    m = Model(opt)
    lb = oneunit(T)*-2
    up = oneunit(T)*2
    n_up = length(xu)
    n_lo = length(xl)
    @variable(m, lb <= A₀ <= up, start = 1.0)
    @variable(m, lb <= Aᵤ[1:n] <= up, start = 1.0)
    @variable(m, lb <= Aₗ[1:n] <= up, start = -1.0)
    @expression(m, yuc[j=1:n_up], cst_C(c1,c2,xu[j])*( A₀*cst_S(n,0,xu[j])+sum(Aᵤ[i]*cst_S(n,i,xu[j]) for i = 1:n))+xu[j]*dz)
    @expression(m, ylc[j=1:n_lo], cst_C(c1,c2,xl[j])*(-A₀*cst_S(n,0,xl[j])+sum(Aₗ[i]*cst_S(n,i,xl[j]) for i = 1:n))-xl[j]*dz)
    @expression(m, fup, sum((yuc[j]-yu[j])^2 for j=1:n_up))
    @expression(m, flo, sum((ylc[j]-yl[j])^2 for j=1:n_lo))
    @objective(m, Min, fup+flo)
    optimize!(m)
    status = termination_status(m)
    value(A₀), value.(Aᵤ), value.(Aₗ), status
end


text:
### A Pluto.jl notebook ###
# v0.17.2

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ c7952ee6-45b5-11eb-1158-5bb2ff274ce9
begin
	psd = "C:\\Users\\dusty\\Documents\\PlaymatProjects\\PlaymatSimulator"

	import Pkg
	Pkg.activate(psd)

	using JSON
	using HTTP
	using GameOne
	using PlutoUI
	using Serialization
	using Images: load, RGB
	using PlaymatSimulator
	using ImageTransformations: imresize

	gd = "C:/Users/dusty/Documents/PlaymatProjects/PlaymatGames/Mage"

	md"""
	*mtg_cards loaded from JSON*

	## EDH Deck: Vannifar's Circus
	"""
end

# ╔═╡ 0c145632-3927-11eb-19b9-877e05c1bcdc
begin
	# card db: https://scryfall.com/docs/api/bulk-data

	mtg_cards = JSON.parsefile("$gd/Base/json/oracle-cards-20211125220340.json")

	md"**Found $(length(mtg_cards)) cards in JSON db!**"
end

# ╔═╡ 621b08a4-384e-11eb-0109-61e9b9ecf125
if !(@isdefined deck)
deck = Dict{Symbol,Any}(
:name => "Vannifar's Circus",
:commander_names => ["Prime Speaker Vannifar"],
:card_names => [

"Acidic Slime",
"Augmenter Pugilist // Echoing Equation",
"Augur of Autumn",
"Alchemist's Refuge",
"Amphin Mutineer",
"Birds of Paradise",
"Barkchannel Pathway // Tidechannel Pathway",
"Brainstorm",
"Breeding Pool",
"Coiling Oracle",
"Command Tower",
"Counterspell",
"Crop Rotation",
"Cryptic Command",
"Cyclonic Rift",
"Deadeye Navigator",
"Devoted Druid",
"Dryad Arbor",
"Dryad of the Ilysian Grove",
"Elvish Mystic",
"Elvish Reclaimer",
"Eternal Witness",
"Experiment Kraj",
"Bala Ged Recovery // Bala Ged Sanctuary",
"Fatestitcher",
"Fblthp, the Lost",
"Forbidden Orchard",
"Forest",
"Glasspool Mimic // Glasspool Shore",
"Glen Elendra Archmage",
"Grand Architect",
"Green Sun's Zenith",
"Growth Spiral",
"Incubation Druid",
"Island",
"Kinnan, Bonder Prodigy",
"Kiora's Follower",
"Leech Bonder",
"Ley Weaver",
"Lightning Greaves",
"Llanowar Reborn",
"Lore Weaver",
"Magosi, the Waterveil",
"Mana Drain",
"Maze of Ith",
"Minamo, School at Water's Edge",
"Mirrorhall Mimic // Ghastly Mimicry",
"Misty Rainforest",
"Mulldrifter",
"Murkfiend Liege",
"Nykthos, Shrine to Nyx",
"Ornithopter of Paradise",
"Paradise Mantle",
"Parcelbeast",
"Pemmin's Aura",
"Pili-Pala",
"Ramunap Excavator",
"Rapid Hybridization",
"Reclamation Sage",
"Rejuvenating Springs",
"Rift Sower",
"Rimewood Falls",
"Rishkar, Peema Renegade",
"Tangled Florahedron // Tangled Vale",
"Thassa, Deep-Dwelling",
"Timeless Witness",
"Trickbind",
"Saryth, the Viper's Fang",
"Seedborn Muse",
"Simic Growth Chamber",
"Skullclamp",
"Snapcaster Mage",
"Sol Ring",
"Spellseeker",
"Spore Frog",
"Strip Mine",
"Swan Song",
"Teferi, Mage of Zhalfir",
"Temple of Mystery",
"Temporal Mastery",
"Thousand-Year Elixir",
"Tolaria West",
"Tribute Mage",
"Trinket Mage",
"Tropical Island",
"Vastwood Fortification // Vastwood Thicket",
"Venser, Shaper Savant",
"Scryb Ranger",
"Vizier of Tumbling Sands",
"Waterlogged Grove",
"Willbreaker",
"Wirewood Symbiote",
"Worldly Tutor",
"Yavimaya, Cradle of Growth",
"Yavimaya Coast",
"Young Wolf",
"Urza's Saga",
"Viridian Longbow",
]
)
else
	nothing
end

# ╔═╡ c61fa79e-4583-11eb-3b71-2d334aca843d
md"""
Look OK? Note: images that do not require in-game scaling will suffer less distortion.
"""

# ╔═╡ 7420cf10-45cc-11eb-2780-4f320bd8a2cf
begin
	cards = []
	card_names = mtg_cards |> keys |> collect

	for n in sort(vcat(deck[:card_names], deck[:commander_names]))

		for c in mtg_cards

			if n == c["name"]
				push!(cards, c)
				break
			end
		end

	end

	md"""
	**Found $(length(cards)) matching cards in mtg_cards!**
	"""
end

# ╔═╡ c31dd202-50c6-11eb-0631-13c70535635e
missing_cards = filter!(x->!(x in [ c["name"] for c in cards ]), vcat(deck[:card_names], deck[:commander_names]))

# ╔═╡ c5374766-4ef1-11eb-2555-c159dba953f0
md"""Card #: $(@bind i Slider(1:length(cards), default=50, show_value=true))
Shrink card by $(@bind card_ratio Slider(0.1:0.05:1.25, default=0.5, show_value=true))

or: $(@bind customw CheckBox(default=false)) specify the default card width: $(@bind card_width NumberField(1:10; default=270)) in pixels
"""

# ╔═╡ 0899b1d0-5972-11eb-0470-4d480cf95d53
md"""
##### Download card imgs? $(@bind download_imgs CheckBox())
"""

# ╔═╡ be2af776-5971-11eb-13ec-3d7982a01ea3
md"""
##### Save data to disk? $(@bind save_data CheckBox())
"""

# ╔═╡ 43cdcbdc-5f68-11eb-11bd-51fbbf610333
deck

# ╔═╡ ce216c54-468a-11eb-13b8-7f3dac7af44a
function get_face_img(img_uri::String)
	img_resp = HTTP.get(img_uri)
	card_img = img_resp.body |> IOBuffer |> load
end

# ╔═╡ 2ab53d00-50cd-11eb-1cd4-5bf94ce53692
function get_card_preview_img(c)
	if haskey(all_cards[i], "card_faces") && haskey(all_cards[i]["card_faces"][1], "image_uris")
		hcat([get_face_img(f["image_uris"]["border_crop"])
			for f in all_cards[i]["card_faces"] ]...)
	else
		get_face_img(all_cards[i]["image_uris"]["border_crop"])
	end
end

# ╔═╡ 97bc3768-50ce-11eb-3f74-95d4fefe3792
function get_card_faces(c, faces=[])
	if haskey(c, "image_uris")
		push!(faces, get_face_img(c["image_uris"]["border_crop"]))

	elseif haskey(c, "card_faces")
		for f in c["card_faces"]
			pushfirst!(faces, get_face_img(f["image_uris"]["border_crop"])) # MtG card faces come in 2nd face 1st
		end
	end

	return faces
end

# ╔═╡ 73d1cd18-4647-11eb-3994-7d4eb92eddca
if length(cards) > 0 && (@isdefined cards)

	if customw
		fs = get_card_faces(cards[i])
		sz = size(fs[end])
		imresize.(fs, ceil(Int, card_width * sz[1]/sz[2]), card_width)[end]
	else
		imresize.(get_card_faces(cards[i]), ratio=card_ratio)[end]
	end
else
	nothing
end

# ╔═╡ dfc9b56e-50ce-11eb-0e7f-83ec6e831901
if download_imgs
	CARD_FACES = []

	for c in cards
		faces = get_card_faces(c)

		if customw
			sz = size(faces[end])
			faces = imresize.(faces, ceil(Int, card_width * sz[1]/sz[2]), card_width)
		else
			faces = imresize.(faces, ratio=card_ratio)
		end

		push!(CARD_FACES, c["name"] => faces)
		sleep(0.1)
	end
end

# ╔═╡ 3a17fc4c-5f69-11eb-29a1-c9bb584fee32
if @isdefined CARD_FACES
	CARD_FACES[i][end][end]
end

# ╔═╡ d654ad1e-468a-11eb-2348-695621b7b9b0
function search_mtg_cards_by_keyword(q::String, mtg_cards::Array)
	[ n for n in [ c["name"] for c in mtg_cards ] if occursin(q, n) ]
end

# ╔═╡ 7cbdfda6-5eb5-11eb-2be6-0f93c3374a47
function hbox(x, y, gap=16; sy=size(y), sx=size(x))
	w,h = (max(sx[1], sy[1]),
		   gap + sx[2] + sy[2])

	slate = fill(RGB(1,1,1), w,h)
	slate[1:size(x,1), 1:size(x,2)] .= RGB.(x)
	slate[1:size(y,1), size(x,2) + gap .+ (1:size(y,2))] .= RGB.(y)
	slate
end

# ╔═╡ 5f7ebd78-3db7-11eb-0690-1b8ee4ebe7db
begin
	CARD_BACK_PATH = "$gd/Base/ui/cards/card_back.png"
	CARD_BACK_RAW = load(CARD_BACK_PATH)
	CARD_BACK_IMG = imresize(CARD_BACK_RAW, ratio=card_ratio)

	hbox(CARD_BACK_RAW, CARD_BACK_IMG)
end

# ╔═╡ 5e6f7046-4da5-11eb-0122-bd82397aab4f
if save_data
	deck[:CARD_BACK_IMG] = CARD_BACK_IMG
	deck[:CARD_FACES] = CARD_FACES
	deck[:CARD_WIDTH] = CARD_FACES[1][end][end] |> width |> Int32
	deck[:CARD_HEIGHT] = CARD_FACES[1][end][end] |> height |> Int32

	fn = "$gd/EDH/decks/$(deck[:name])/$(deck[:name]).jls"
	serialize(fn, deck)

	md"Deck data saved to $fn"
end

# ╔═╡ 85103516-5eb5-11eb-3abc-dfc9ae9caf96
vbox(x,y, gap=16) = hbox(x', y')'

# ╔═╡ b0ae3401-4c58-4832-9f50-5669949717d3
# page/cell width control
"""
<style>
	main {
	max-width: 1100px;
	align-left: flex-start;
	margin-left: 50px;
	}
""" |> HTML

# ╔═╡ Cell order:
# ╟─c7952ee6-45b5-11eb-1158-5bb2ff274ce9
# ╟─0c145632-3927-11eb-19b9-877e05c1bcdc
# ╟─621b08a4-384e-11eb-0109-61e9b9ecf125
# ╟─5f7ebd78-3db7-11eb-0690-1b8ee4ebe7db
# ╟─c61fa79e-4583-11eb-3b71-2d334aca843d
# ╟─7420cf10-45cc-11eb-2780-4f320bd8a2cf
# ╟─c31dd202-50c6-11eb-0631-13c70535635e
# ╟─c5374766-4ef1-11eb-2555-c159dba953f0
# ╟─73d1cd18-4647-11eb-3994-7d4eb92eddca
# ╟─0899b1d0-5972-11eb-0470-4d480cf95d53
# ╟─3a17fc4c-5f69-11eb-29a1-c9bb584fee32
# ╟─dfc9b56e-50ce-11eb-0e7f-83ec6e831901
# ╟─be2af776-5971-11eb-13ec-3d7982a01ea3
# ╟─5e6f7046-4da5-11eb-0122-bd82397aab4f
# ╟─43cdcbdc-5f68-11eb-11bd-51fbbf610333
# ╟─ce216c54-468a-11eb-13b8-7f3dac7af44a
# ╟─2ab53d00-50cd-11eb-1cd4-5bf94ce53692
# ╟─97bc3768-50ce-11eb-3f74-95d4fefe3792
# ╟─d654ad1e-468a-11eb-2348-695621b7b9b0
# ╟─7cbdfda6-5eb5-11eb-2be6-0f93c3374a47
# ╟─85103516-5eb5-11eb-3abc-dfc9ae9caf96
# ╟─b0ae3401-4c58-4832-9f50-5669949717d3


text:
@testset "41.first-missing-positive.jl" begin
    @test first_missing_positive([1, 2, 0]) == 3
    @test first_missing_positive([3, 4, -1, 1]) == 2
    @test first_missing_positive([7, 8, 9, 11, 12]) == 1
end


text:
import Distributions, Random, StatsFuns, DataFrames

"""
    simulate_data_logistic(coefs; n::Int=100, noise::Number=0.0, groupnames=:random)

Generate a DataFrame of variables related a binary dependent variable by specified regression `coefs`.

!!! warning

    This function, adapted from [this](https://stats.stackexchange.com/questions/46523/how-to-simulate-artificial-data-for-logistic-regression) thread, doesn't work as expected (See [#27](https://github.com/neuropsychology/Psycho.jl/issues/27)).

# Multiple Variables / Groups
- If `coefs` is a vector (*e.g., `[0.1, 0.2]`), the DataFrame will contain `length(coefs)` variables (`Var1, Var2, ...`). Altough uncorrelated between them, they are  correlated to the outcome (`y`) by the specified coefs.
- If `coefs` is a vector of vectors (*e.g., `[[0.1], [0.2]]`), it will create `length(coefs)` groups, *i.e.*, stacked DataFrames with a correlation between the variables and the outcome varying between groups. It is possible to specify the `groupnames`.

# Arguments
- `coefs`: Regression coefficients. Can be a number, a vector of numbers or a vector of vectors.
- `n::Int`: Number of observations.
- `noise::Number`: The SD of the random gaussian noise.
- `groupnames::Vector`: Vector of group names (default to `:random`).
- `kwargs...`: Arguments to pass to other functions.

!!! note

    **Ideas / help required:**
    - Different group sizes (See [#9](https://github.com/neuropsychology/Psycho.jl/issues/9))

# Examples
```julia
simulate_data_logistic(0.2)
```
"""
function simulate_data_logistic end

function simulate_data_logistic(coefs::Vector{<:Number}; n::Int=100, noise::Number=0.0)
    # https://stats.stackexchange.com/questions/46523/how-to-simulate-artificial-data-for-logistic-regression
    n_var = length(coefs)
    X = randn(n, n_var)
    Z = X * coefs .+ 1

    p = StatsFuns.logistic.(Z)
    y = Float64.(rand(length(p)) .< p)

    # Add noise
    if noise != 0
      X = X .+ Random.rand(Distributions.Normal(0, noise), n, n_var)
    end

    data = DataFrames.DataFrame(hcat(y, X),
      Symbol.(vcat(["y"], "Var" .* string.(1:n_var))))
    return data

end

function simulate_data_logistic(coefs::Number; n::Int=100, noise::Number=0.0)
  simulate_data_logistic([coefs], n=n, noise=noise)
end

function simulate_data_logistic(coefs::Vector{<:Vector}; n::Int=100, noise::Number=0.0, groupnames=:random, kwargs...)
  data = simulate_data_groups(simulate_data_logistic, coefs=coefs, n=n, noise=noise, groupnames=groupnames, kwargs...)
  return data
end


text: const CROSSVALIDATION_EXTRAKEY = :crossvalidation


text:
"""This module implements classical constraints and their gradients."""
# Author: Henri Gérard <hgerard.proy@gmail.com>
# License: MIT

function div_function(nameDiv::DivergenceConstraint)
    error("'error_loss' not defined for $(typedof(s))")
end

function div_function(nameDiv::KLConstraint)
    return (f = x->exp(x)-1)
end

function div_grad(nameDiv::DivergenceConstraint)
    error("'error_loss' not defined for $(typedof(s))")
end

function div_grad(nameDiv::KLConstraint)
    return (f = x->exp(x))
end

function fconstraint(x, Z, regressionModel, i, constraint::GeneralConstraint)
    error("'error_loss' not defined for $(typedof(s))")
end

function fconstraint(x, Z, regressionModel, i, constraint::PositiveConstraint)
    return -x[1]
end

function fconstraint(x, Z, regressionModel, i, constraint::DivergenceConstraint)
    lt = size(Z)[2]-1

    λ = x[1]
    μ = x[2]
    θ = x[3:3+lt-1]
    s = x[3+lt:end]

    # f = x->exp(x)-1
    # f = constraint.f
    f = div_function(constraint)

    # inside_number = (error_loss(θ, Z[i,:], regressionModel) - μ)/λ

    if (λ<=0)
        return error_loss(θ, Z[i,:], regressionModel)-μ
    else
        return λ*f((error_loss(θ, Z[i,:], regressionModel)-μ)/λ)-s[i]
    end
end

function fconstraint(x, Z, regressionModel, i, constraint::EntropicConstraint)
    lt = size(Z)[2]-1

    μ = x[1]
    θ = x[2:2+lt-1]
    s = x[2+lt:end]

    # f = x->exp(x)-1
    # f = constraint.f
    f = div_function(KLConstraint())

    return f((error_loss(θ, Z[i,:], regressionModel)-μ))-s[i]

end

function fconstraint(x, Z, regressionModel, i, constraint::DROConstraint)
    lt = size(Z)[2]-1

    τ = x[1]
    θ = x[2:2+lt-1]
    s = x[2+lt:end]

    N = length(s)

    l = div(i-1,N)+1
    j = i-(l-1)*N

    return error_loss(θ, Z[j,:], regressionModel)- τ*(Z[l,:]-Z[j,:])'*(Z[l,:]-Z[j,:])-s[l]
end

# function circle(x,i)
#      return (x[1]-Z[i,1])^2 + (x[2]-Z[i,2])^2 - Z[i,3]^2
# end
#
# function expo(x,i)
#     return exp(Z[i,1]*x[1]) - x[2]
# end
# if fname == "circle"
#          return [2*(x[1]-Z[i,1]); 2*(x[2]-Z[i,2])]
#     elseif fname == "expo"
#         return [Z[i,1]*exp(Z[i,1]*x[1]); -1]

#
# function subgradient_lambda(f, fp, λ, inside_number, loss_eval)
#     return f(inside_number) - loss_eval/λ*fp(inside_number)
# end
#
# function subgradient_mu(fp, inside_number)
#     return -fp(inside_number)
# end
#
# function subgradient_theta(fp, λ, inside_number, lossp)
#     grad = zeros(length(lossp))
#     for k in 1:length(lossp)
#         grad[k] = lossp[k]*fp(inside_number)
#     end
#     return grad
# end

function subgradient(x, Z, regressionModel, i, constraint::GeneralConstraint)
    error("'error_loss' not defined for $(typedof(s))")
end

function subgradient(x, Z, regressionModel, i, constraint::PositiveConstraint)
    return [-1; zeros(length(x)-1)]
end

function subgradient(x, Z, regressionModel, i, constraint::DivergenceConstraint)
    lt = size(Z)[2]-1

    λ = x[1]
    μ = x[2]
    θ = x[3:3+lt-1]
    s = x[3+lt:end]

    # f  = x->exp(x)-1
    # fp = x->exp(x)
    f  = div_function(constraint)
    fp = div_grad(constraint)

    inside_number = (error_loss(θ, Z[i,:], regressionModel) - μ)/λ
    lossp = zeros(lt)
    for k in 1:lt
        lossp[k] = error_prime(θ, Z[i,:], k, regressionModel)
    end

    if  λ<=0

        grad = zeros(length(x))
        grad[2] = -1
        grad[3:3+lt-1] = lossp
        grad[3+lt:end] = zeros(length(s))

        return grad
    else

        grad = zeros(length(x))
        # grad[1] = subgradient_lambda(f, fp, λ, inside_number, (error_loss(θ, Z[i,:], regressionModel) - μ))
        # grad[2] = subgradient_mu(fp, inside_number)
        # grad[3:3+lt-1] = subgradient_theta(fp, λ, inside_number, lossp)
        grad[1] = f(inside_number) - (error_loss(θ, Z[i,:], regressionModel) - μ)/λ*fp(inside_number)
        grad[2] = -fp(inside_number)
        for k in 1:length(lossp)
            grad[3+k-1] = lossp[k]*fp(inside_number)
        end
        grad[3+lt-1+i] = -1
        return grad
    end
end

function subgradient(x, Z, regressionModel, i, constraint::EntropicConstraint)
    lt = size(Z)[2]-1

    μ = x[1]
    θ = x[2:2+lt-1]
    s = x[2+lt:end]

    # f  = x->exp(x)-1
    # fp = x->exp(x)
    f  = div_function(KLConstraint())
    fp = div_grad(KLConstraint())

    inside_number = (error_loss(θ, Z[i,:], regressionModel) - μ)
    lossp = zeros(lt)
    for k in 1:lt
        lossp[k] = error_prime(θ, Z[i,:], k, regressionModel)
    end

    grad = zeros(length(x))
    grad[1] = -fp(inside_number)
    for k in 1:length(lossp)
        grad[2+k-1] = lossp[k]*fp(inside_number)
    end
    grad[2+lt-1+i] = -1
    return grad

end

function subgradient(x, Z, regressionModel, i, constraint::DROConstraint)
    lt = size(Z)[2]-1

    τ = x[1]
    θ = x[2:2+lt-1]
    s = x[2+lt:end]

    N = length(s)

    l = div(i-1,N)+1
    j = i-(l-1)*N

    lossp = zeros(lt)
    for k in 1:lt
        lossp[k] = error_prime(θ, Z[j,:], k, regressionModel)
    end

    grad = zeros(length(x))
    grad[1] = -(Z[l,:]-Z[j,:])'*(Z[l,:]-Z[j,:])
    grad[2:2+lt-1] = lossp
    grad[2+lt-1+l] = -1

    return grad
end

max_stars_repo_name: Henri-Gerard/robox

max_issues_repo_name: Henri-Gerard/robox

max_forks_repo_name: Henri-Gerard/robox


text:
# GPU implementation

# a port of the GPU kernels from Baidu's C++ warp-ctc package,
# which itself is Copyright 2015-2016 Baidu USA LLC
# and available under the Apache 2.0 license
#
# Apache 2.0 license: https://www.apache.org/licenses/LICENSE-2.0
# GitHub: https://github.com/baidu-research/warp-ctc/
# paper: https://arxiv.org/pdf/1512.02595.pdf

using Flux
using Statistics
using CUDA
using NNlib

const MAX_THREADS = 256

function log_plus_f(p1, p2)
  isinf(p1) && return p2
  isinf(p2) && return p1
  if p1 < p2
    p1, p2 = p2, p1
  end
  return p1 + log(1+exp(p2 - p1))
end

function count_repeats(A)
  repeats = 0
  for (i,elem) in enumerate(A)
    if i > 1 && A[i] == A[i-1]
      repeats += 1
    end
  end
  return repeats
end

function compute_alpha_kernel(probs, labelSize, uttLength, repeats, labelsWithoutBlanks, labelsWithBlanks, alpha, blankLabel)

  tid = threadIdx().x
  L = labelSize
  T = uttLength
  S = length(labelsWithBlanks)

  if L + repeats > T
    return nothing
  end
  labels = labelsWithBlanks

  # Corner-case checking
  start = (L + repeats <= T) ? 0 : 1
  last = S > 1 ? 2 : 1

  # Fill in first column (time step)
  i = tid
  while i <= last - start
    alpha[start+i, 1] = probs[labels[start+i], 1]
    i += blockDim().x
  end
  sync_threads()

  # Fill in coefficients for each time step
  for t=2:T
    # Corner-case checking
    if tid == 1 && !(1 < S - 2*(T-t) - 1)
      if start == 0
        alpha[1, t] = alpha[1, t-1] + probs[blankLabel, t]
      elseif start == 1
        alpha[1, t] = alpha[1, t-1]
      end
    end
    sync_threads()

    # Fill in coefficients for each label class in the target output sequence;
    # each thread will process the calculations for one class
    idx = tid+1
    while idx <= S
      prevSum = log_plus_f(alpha[idx, t-1], alpha[idx-1, t-1])
      if labels[idx] != blankLabel && idx != 2 && labels[idx] != labels[idx-2]
        prevSum = log_plus_f(prevSum, alpha[idx-2, t-1])
      end
      if idx < S - 2*(T-t) - 1
        alpha[idx, t] = -Inf32
      else
        alpha[idx, t] = prevSum + probs[labels[idx], t]
      end
      idx += blockDim().x
    end
    sync_threads()
  end
  return nothing
end

function compute_beta_and_grad_kernel(probs, labelSize, uttLength,
                  repeatsInLabel, labelsWithBlanks,
                  alphas, beta, output, accum,
                  grad, blankLabel, loss)

  tid = threadIdx().x
  L = labelSize
  T = uttLength
  S = 2*L + 1
  repeats = repeatsInLabel
  labels = labelsWithBlanks

  if (L+repeats) > T
    return nothing
  end

  # Corner-case checking
  start = S > 1 ? S-2 : 0
  last = L + repeats < T ? S : S-1
  sync_threads()
  i = tid

  # Calculate coefficients for last column (time step)
  # then determine alpha and beta product
  while i <= last - start
    beta[i+start, T] = 0
    output[i+start, T] = beta[i+start, T] + alphas[i+start, T]
    i += blockDim().x
  end
  sync_threads()

  # Fill in `accum` for last column (time step)
  if tid == 1
    for i=1:S
      labelIdx = labels[i]
      accum[labelIdx, T] = log_plus_f(accum[labelIdx, T], output[i, T])
    end
  end
  sync_threads()

  # Fill in `grad` for last column (time step)
  idx = tid
  while idx <= size(grad, 1)
    s = -Inf32
    for i=1:S
      s = log_plus_f(s, output[i, T])
    end

    # ∂L/∂a (where a is activation before logsoftmax)
    grad[idx, T] = exp(probs[idx, T]) - exp(accum[idx, T] - s)
    idx += blockDim().x
  end
  sync_threads()

  # Fill in the rest of the coefficients
  t = T-1
  while t >= 1
    if t < T
      idx = tid
      while idx <= S
        nextSum = probs[labels[idx], t+1] + beta[idx, t+1]
        if idx < S
          nextSum = log_plus_f(nextSum,
            probs[labels[idx+1], t+1] + beta[idx+1, t+1])
        end
        if labels[idx] != blankLabel && idx != S-1 && labels[idx] != labels[idx+2]
          nextSum = log_plus_f(nextSum,
            probs[labels[idx+2], t+1] + beta[idx + 2, t+1])
        end
        if idx > 2*t
          beta[idx, t] = -Inf32
        else
          beta[idx, t] = nextSum
        end
        idx += blockDim().x
      end
      sync_threads()
      idx = tid
      while idx <= S
        output[idx, t] = alphas[idx, t] + beta[idx, t]
        idx += blockDim().x
      end
      sync_threads()
    end
    sync_threads()

    # Calculate accumulated alpha-beta products for each label class for
    # each time step; used in calculating gradients
    if tid == 1
      for i=1:S
        labelIdx = labels[i]
        accum[labelIdx, t] = log_plus_f(accum[labelIdx, t], output[i, t])
      end
    end
    sync_threads()
    idx = tid

    # Calculate gradients
    while idx <= size(grad, 1)

      # ∂L/∂a (where a is activation before logsoftmax)
      grad[idx, t] = exp(probs[idx, t]) - exp(accum[idx, t] + loss)
      idx += blockDim().x
    end
    sync_threads()
    t -= 1
    sync_threads()
  end
  return nothing
end

function ctc_alpha(ŷ::CuArray, y)
  ŷ = logsoftmax(ŷ)
  blank = size(ŷ, 1)
  ycu = cu(y)
  z′ = CUDA.fill(blank, 2 * length(y) + 1)
  z′[eachindex(y) .* 2] .= ycu
  T = size(ŷ, 2)
  U′ = 2*length(y) + 1
  alphas = CUDA.fill(log(zero(eltype(ŷ))), U′,T)
  nRepeats = count_repeats(cpu(y))
  nThreads = min(U′, MAX_THREADS)
  @cuda blocks=1 threads=nThreads compute_alpha_kernel(ŷ, length(y), T, nRepeats, ycu, z′, alphas, blank)
  return (loss=-1 * logsumexp(alphas[end-1:end]), alpha=alphas, z′=z′, yhat=ŷ, nRepeats=nRepeats)
end

ctc_loss(ŷ::CuArray, y) = ctc_alpha(ŷ::CuArray, y).loss

function ∇ctc_loss(ŷ::CuArray, y, out)
  loss, alphas, z′, ŷ, nRepeats = out
  U′, T = size(alphas)
  blank = size(ŷ, 1)
  typed_zero = zero(eltype(ŷ))
  betas = CUDA.fill(log(typed_zero), U′, T)
  output = CUDA.fill(log(typed_zero), U′, T)
  nThreads = min(U′, MAX_THREADS)
  grads = CUDA.fill(log(typed_zero), size(ŷ))
  accum = CUDA.fill(log(typed_zero), size(ŷ))
  @cuda blocks=1 threads=nThreads compute_beta_and_grad_kernel(ŷ, length(y), T, nRepeats, CuArray(z′), alphas, betas, output, accum, grads, blank, loss)
  return grads
end


text:
# eml_em.jl
using MAT: matopen
using LinearMapsAA: LinearMapAA, LinearMapAO
using LinearAlgebra: norm
using SPECTrecon: RealU, SPECTplan, project!, backproject!
using MIRTjim: jim

"""
    xs = eml_em(x0, G, yi, ri, niter)
E-ML-EM algorithm for image reconstruction from Poisson emission data
with system model: y = Poisson(G * x + r)
* Input:
-`x0`: initial image guess
-`G`: system model
-`yi`: projections
-`ri`: scatters
* Output:
-`xs`: estimates each iteration
"""
function eml_em(x0::AbstractArray{<:RealU, 3},
                G::LinearMapAO,
                yi::AbstractArray{<:RealU, 3},
                ri::AbstractArray{<:RealU, 3},
                niter::Int,
                )
    nx, ny, nz = size(x0)
    xs = zeros(eltype(x0), nx, ny, nz, niter + 1)
    xs[:, :, :, 1] .= x0
    Asum = G' * ones(eltype(yi), size(yi))
    x = x0
    x_old = x0
    nrmsd(x, x_old) = norm(vec(x .- x_old)) / norm(vec(x_old)) * 100
    for iter = 1:niter
        x_old = x
        yp = (G * x) .+ ri
        yp[yp .== 0] .= Inf
        eterm = G' * (yi ./ yp)
        x = x .* eterm ./ Asum
        x = max.(x, 0)
        x = min.(x, Inf)
        xs[:, :, :, iter + 1] .= x
        if iter % 2 == 0
            print("iter: ", iter,
                  ", range: ", round.((minimum(x), maximum(x)); digits = 3),
                  ", relative diff: ", round(nrmsd(x, x_old); digits = 3), "%",
                  "\n"
                  )
        end
    end
    return xs
end

"""
    linearmap(mumap, psfs, dy; interpmeth, mode)
Return a linear map operator of SPECT plan
*Input:
-`mumap`: attenuation map
-`psfs`: point spread function
-`dy`: pixel width for `y` direction
-`interpmeth`: interpolation method, `:one` for 3-pass 1d, `:two` for 2d
-`mode`: `:fast` or `:mem`
"""
function linearmap(mumap::AbstractArray{<:RealU, 3},
                   psfs::AbstractArray{<:RealU, 4},
                   dy::RealU;
                   interpmeth::Symbol = :two,
                   mode::Symbol = :fast,
                   )
        T = eltype(mumap)
        nx, ny, nz = size(mumap)
        nview = size(psfs, 4)
        idim = (nx, ny, nz)
        odim = (nx, nz, nview)
        plan = SPECTplan(mumap, psfs, dy; interpmeth, mode)
        forw! = (y,x) -> project!(y, x, plan)
        back! = (x,y) -> backproject!(x, y, plan)
        return LinearMapAA(forw!, back!, (prod(odim),prod(idim)); T, idim, odim)
end

function test_eml_em()
    T = Float32
    path = "../data/"
    file = matopen(path*"mumap208.mat")
    mumap = read(file, "mumap208")
    close(file)

    file = matopen(path*"psf_208.mat")
    psfs = read(file, "psf_208")
    close(file)

    file = matopen(path*"xinit.mat")
    xinit = read(file, "xinit")
    close(file)

    file = matopen(path*"yi-80.mat")
    yi = read(file, "yi")
    close(file)

    file = matopen(path*"ri-80.mat")
    ri = read(file, "ri")
    close(file)

    dy = T(4.7952)
    G = linearmap(mumap, psfs, dy; interpmeth = :two, mode = :fast)
    xs = eml_em(xinit, G, yi, ri, 20)
    return xs
end

xs = test_eml_em()

path = "../data/"
file = matopen(path*"xtrue.mat")
xtrue = convert(Array{Float32, 3}, read(file, "xtrue"))
close(file)

file = matopen(path*"xrecon_hk.mat")
xrecon_hk = convert(Array{Float32, 4}, read(file, "xrecon_hk"))
close(file)

idx = 40
niter = 21
jim(jim(xtrue[:, :, idx], "xtrue"),
    jim(xrecon_hk[:, :, idx, niter], "recon_hk"),
    jim(xs[:, :, idx, niter], "recon_zl"),
    jim(xrecon_hk[:, :, idx, niter] - xs[:, :, idx, niter], "diff (hk-zl)"))


text:
@testset "Boykov Kolmogorov" begin
    # construct graph
    gg = DiGraph(3)
    add_edge!(gg, 1, 2)
    add_edge!(gg, 2, 3)

    # source and sink terminals
    source, target = 1, 3

    for g in testdigraphs(gg)
    # default capacity
      capacity_matrix = LightGraphs.DefaultCapacity(g)
      residual_graph = @inferred(LightGraphs.residual(g))
      T = eltype(g)
      flow_matrix = zeros(T, 3, 3)
      TREE = zeros(T, 3)
      TREE[source] = T(1)
      TREE[target] = T(2)
      PARENT = zeros(T, 3)
      A = [T(source), T(target)]
# see https://github.com/JuliaLang/julia/issues/21077
# @show("testing $g with eltype $T, residual_graph type is $(eltype(residual_graph)), flow_matrix type is $(eltype(flow_matrix)), capacity_matrix type is $(eltype(capacity_matrix))")
      path = LightGraphs.find_path!(
        residual_graph, source, target, flow_matrix,
        capacity_matrix, PARENT, TREE, A)

      @test path == [1, 2, 3]
  end
end


text:
include("SnapData/partsnapLib.jl")
#TestSnap()
 ###############para##############
folderB="E:/a/b"
folderA=""
#folderA=""
# 2019/4/4	10:16:26

StartT="20190404" * "103826"
StartTime= DtFormat(StartT)
TimePeriod=4*60 #s
ChannelID=1
StartFiberLength=16595-100 #meters
FLRange=200 #meters
TgrFolder="E:/a/"

################Para######################
ret = SnapPicedata2File(folderB,
                folderA,
                StartTime,
                TimePeriod,
                StartFiberLength,
                FLRange,
                ChannelID,
                TgrFolder  )
println(ret)

max_stars_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_issues_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_forks_repo_name: wangzhongtian/fiberAlarmAnalyseTool


text:
using Keras
using Flux
using Base.Test
using DataFlow:Call, constant, vertex, syntax

vcall(a...) = vertex(Call(), constant.(a)...)
Keras.new_type(a, b) = Keras.new_type(a, nothing, b)
ip = rand(32,32,3,1)

#BatchNorm Test
d = Dict{Any, Any}()
d["momentum"] = 0.9
d["epsilon"] = 0.1
a = Keras.new_type(:BatchNormalization, d)
@testset "Normalize" begin
    @test Keras.ops[:BatchNormalization](a)(ip) == BatchNorm(3, ϵ=0.1, momentum=0.9)(ip)
end


text:
# Recode: function signature -> Pattern creating AST
module Recode
using ..Meta, ..Patterns, ..Nodes
export recode, @qpat, @ipat
export refnode, egalpred, typepred, lengthnode

macro qpat(ex)
    recode(ex)[1]
end
macro ipat(ex)
    c = Context()
    recode(c, quot(argnode), ex)
    quote
        $(c.code...)
        intension($(c.preds...))
    end
end

type Context
    code::Vector
    preds::Vector
    bindings::Vector
    Context() = new([], [], [])
end

function recode(ex)
    c = Context()
    recode(c, quot(argnode), ex)
    p_ex = quote
        $(c.code...)
        Pattern(intension($(c.preds...)),
                $(Expr(:typed_dict, :(Symbol=>Node), c.bindings...)))
    end
    syms = Symbol[b.args[1].args[1] for b in c.bindings]
    if length(syms) != length(Set(syms))
        error("@pattern: Non-tree patterns not yet implemented. Use unique argument names.")
    end
    p_ex, syms
end

recode(c::Context, arg, ex) = push!(c.preds, :(egalpred($arg,$(quot(ex)))))
recode(c::Context, arg, ex::Symbol) = push!(c.bindings, :($(quot(ex))=>$arg))

const at_tilde_symbol = symbol("@~")

function recode(c::Context, arg, ex::Expr)
    head, args = ex.head, ex.args
    nargs = length(args)
    if head === :(::)
        @assert 1 <= nargs <= 2
        if nargs == 1
            push!(c.preds, :( typepred($arg, $(esc(args[1]))) ))
        else
            push!(c.preds, :( typepred($arg, $(esc(args[2]))) ))
            recode(c, arg, args[1])
        end
    elseif head === :tuple || head === :vect
        argtype = (head === :tuple) ? Tuple : Vector
        push!(c.preds, :( typepred($arg, $(quot(argtype))) ))
        push!(c.preds, :( egalpred(lengthnode($arg), $(quot(nargs))) ))
        for (k, p) in enumerate(args)
            node = gensym("e$k")
            push!(c.code, :( $node = refnode($arg, $k) ))
            recode(c, node, p)
        end
    elseif head === :cell1d
        error("@pattern: use [] for array patterns")
    elseif head === :macrocall && args[1] == at_tilde_symbol
        for p in args[2:end]; recode(c, arg, p); end
    elseif head === :$ && nargs == 1
        push!(c.preds, :(egalpred($arg, $(esc(args[1])))))
    elseif head === :... && nargs == 1
        error("@pattern: varargs are not yet implemented")
    else
        error("recode: unimplemented: ex = $ex")
    end
end

end # module