text:
# gausslaguerre from FastGaussQuadrature [FastGaussQuadrature.jl - Gauss quadrature nodes and weights in Julia, Townsend, A.", 2016, https://github.com/ajt60gaibb/FastGaussQuadrature.jl ] with some additional testing options

# (x,w) = gausslaguerre(n) returns n Gauss-Laguerre nodes and weights.
# (x,w) = gausslaguerre(n,alpha) allows generalized Gauss-Laguerre quadrature.
# (x,w) = gausslaguerre(n,alpha,method) allows the user to select which method to use.
# (x,w) = gausslaguerre(n,alpha,method,qm,m) allows the generalised weight function x^alpha*exp(-qm*x^m) where qm and m should be one when method = GW, rec(W) or exp(W)
# METHOD = "GW" will use the traditional Golub-Welsch eigenvalue method, which is best for when N is small.
# METHOD = "gen" can generate an arbitrary number of terms of the asymptotic expansion of Laguerre-type polynomials, orthogonal with respect to x^alpha*exp(-qm*x^m). "genW" does the same, but stops as the weights underflow: it is O(sqrt(n)) when m is one. gausslaguerre(round(Int64, (n/17)^2), alpha, "genW") returns about n nodes and weights above floatmin(Float64) for large n.
# METHOD = "exp" will use explicit expansions and METHOD = "expW" will only compute the weights above floatmin.
# METHOD = "rec" uses forward recurrence and METHOD = "recW" stops when the weights underflow.
# METHOD = "default" uses "gen" when m or qm are not one, explicit formulae when n <= 2, "rec" when 2 < n < 128 and else "exp".

function gausslaguerreTest( n::Int64, alpha::Float64=0.0, method::AbstractString="default",  qm::Float64=1.0, m::Int64=1, T::Int64=ceil(Int64, 34/log(n) ) )

    if ( imag(alpha) != 0 ) || ( alpha < -1 )
        error(string("alpha = ", alpha, " is not allowed.") )
    elseif ( (m != 1) || (qm != 1) ) && ( (method != "gen") && (method != "genW") && (method != "default") )
        error(string("Method ", method, " is not implemented for generalised weights.") )
    end

    if ( (method == "default") && ( (m != 1) || (qm != 1.0) ) ) || (method == "gen") || (method == "genW")
        (x,w) = asyRHgen(n, method == "genW", alpha, m, qm)
        w = w/sum(w)*gamma((alpha+1)/m)*qm^(-(alpha+1)/m)/m # We left out a constant factor while computing w in case we use finite differences
        (x,w)
    elseif (method == "default") && (n == 0)
        Float64[], Float64[]
    elseif (method == "default") && (n == 1)
        [1 + alpha], [1.0]
    elseif (method == "default") && (n == 2)
        [alpha+2 -sqrt(alpha+2); alpha+2+sqrt(alpha+2)], [((alpha-sqrt(alpha+2)+2)*gamma(alpha+2))/(2*(alpha+2)*(sqrt(alpha+2)-1)^2); ((alpha+sqrt(alpha+2)+2)*gamma(alpha+2))/(2*(alpha+2)*(sqrt(alpha+2)+1)^2)]
    elseif method == "GW"
        laguerreGW( n, alpha )         # Use Golub-Welsch
    elseif (method == "GLR")
	if (alpha != 0.0)
            error("GLR not implemented for associated Laguerre polynomials.")
	end
	laguerreGLR( n )               # Use GLR
    elseif ( (method == "default") && (n < 128 ) ) || (method == "rec") || (method == "recW")
        laguerreRec( n, method == "recW", alpha)   # Use forward recurrence
    elseif ( (method == "default") && (n >= 128) ) || (method == "exp") || (method == "expW")
        laguerreExp( n, method == "expW", alpha, T)   # Use explicit expansions and possibly only compute the representable weights
    else
        error(string("Wrong method string, got ", method) )
    end
end

function laguerreGW( n::Int64, alpha::Float64 )
# Calculate Gauss-Laguerre nodes and weights based on Golub-Welsch

    alph = 2*(1:n)-1+alpha           # 3-term recurrence coeffs
    beta = sqrt.( (1:n-1).*(alpha .+ (1:n-1) ) )
    T = SymTridiagonal(alph, beta)  # Jacobi matrix
    x, V = eig( T )                  # eigenvalue decomposition
    w = gamma(alpha+1)*V[1,:].^2     # Quadrature weights
    x, vec(w)
end

########################## Routines for the forward recurrence ##########################

function laguerreRec( n::Int64, compRepr::Bool, alpha::Float64)

    if compRepr
        # Get a heuristic for the indices where the weights are about above floatmin.
        mn = min(ceil(Int64, 17*sqrt(n)), n)
    else
        mn = n
    end
    itric = min(mn, 7)
    bes = besselroots(alpha, itric).^2/(4*n + 2*alpha+2) # [DLMF 18.16.10] says this is lower than the zero, so we do not risk skipping a zero (if none of the results coincide)

    factorw = (n^2 +alpha*n)^(-1/2) # Ratio of leading order coefficients
    w = zeros(mn)
    x = [ bes; zeros(mn - itric)]
    noUnderflow = true
    for k = 1:mn
        if ( x[k] == 0 ) # Use sextic extrapolation for the initial guesses.
            x[k] = 7*x[k-1] -21*x[k-2] +35*x[k-3] -35*x[k-4] +21*x[k-5] -7*x[k-6] +x[k-7]
        end
        step = x[k]
        l = 0 # Newton-Raphson iteration number
        ov = floatmax(Float64) # Previous/old value
        ox = x[k] # Old x
        while ( ( abs(step) > eps(Float64)*40*x[k] ) && ( l < 20) )
            l = l + 1
            # pe = Float64(lagpnRec(big(n), big(alpha), big(x[k]) ) ) #Calculations in BigFloat
            pe = lagpnRec(n, alpha, x[k] )
            if (abs(pe) >= abs(ov)*(1-35*eps(Float64)) )
                # The function values do not decrease enough any more due to roundoff errors.
                x[k] = ox # Set to the previous value and quit.
                break
            end
            step = pe/lagpnRecDer(n, alpha, x[k])
            ox = x[k]
            x[k] = x[k] -step
            ov = pe
    end
    if ( x[k] < 0 ) || ( x[k] > 4*n + 2*alpha + 2 ) || ( l == 20 ) || ( ( k != 1 ) && ( x[k - 1] >= x[k] ) )
        error("Newton method may not have converged.")
    end
    if noUnderflow
        w[k] = factorw/lagpnRec(n -1, alpha, x[k] )/lagpnRecDer(n, alpha, x[k] )
    end
    if noUnderflow && ( w[k] == 0 ) && ( k > 1 ) && ( w[k-1] > 0 ) # We could stop now as the weights underflow.
        if compRepr
		x = x[1:k-1]
		w = w[1:k-1]
		return (x,w)
	    else
                noUnderflow = false
	    end
        end
    end
    x, w
end

# Orthonormal associated Laguerre polynomial with positive leading coefficient, allows BigFloats
function lagpnRec(n,alpha,x)
    pnprev = 0*alpha
    pn= 1/sqrt(gamma(alpha+1) )
    for k=1:n
        pnold = pn
        pn = (x -2*k -alpha+1)/sqrt(k*(alpha+k))*pn-sqrt((k-1+alpha)*(k-1)/k/(k+alpha))*pnprev
        pnprev = pnold
    end
    return pn
end

# Derivative of the orthonormal associated Laguerre polynomial, allows BigFloats
function lagpnRecDer(n,alpha,x)
    pnprev = 0*alpha
    pn= 1/sqrt(gamma(alpha+1) )
    pndprev = 0*alpha
    pnd = 0*alpha
    for k=1:n
        pnold = pn
        pn = (x -2*k -alpha+1)/sqrt(k*(alpha+k))*pnold-sqrt((k-1+alpha)*(k-1)/k/(k+alpha))*pnprev
        pnprev = pnold
        pndold = pnd
        pnd = (pnold+(x-2*k-alpha+1)*pndold)/sqrt(k*(alpha+k)) -sqrt((k-1+alpha)*(k-1)/k/(alpha+k))*pndprev
        pndprev = pndold
    end
    return pnd
end

########################## Routines for the GLR algorithm ##########################

function laguerreGLR( n::Int64 )

    x, ders = ScaledLaguerreGLR( n )  # Nodes and L_n'(x)
    w = exp(-x)./(x.*ders.^2)         # Quadrature weights
    x, w
end

function  ScaledLaguerreGLR( n::Int64 )
    # Calculate the nodes and L_n'(x).

    ders = Array(Float64, n )
    x = Array(Float64,n)
    xs = 1/(2*n+1)
    n1 = 20
    n1 = min(n1, n)
    for k = 1:n1
        xs, dxs = BoundaryLaguerreGLR(n, xs)
        ders[k] = dxs
        x[k] = xs
        xs *= 1.1
    end
    x, ders = InteriorLaguerreGLR(n, x, ders, n1)
end

function InteriorLaguerreGLR(n::Int64, roots, ders, n1::Int64)

    m = 30
    hh1 = ones(m+1)
    zz = Array(Float64,m)
    u = Array(Float64,m+1)
    up = Array(Float64,m+1)
    x = roots[ n1 ]
    for j = n1 : (n - 1)
        # initial approx
        h = rk2_Lag(pi/2, -pi/2, x, n)
        h = h - x

        # scaling:
        M = 1/h
        M2 = M^2
        M3 = M^3
        M4 = M^4

        # recurrence relation for Laguerre polynomials
        r = x*(n + .5 - .25*x)
        p = x^2
        u[1] = 0.
        u[2] = ders[j]/M
        u[3] = -.5*u[2]/(M*x) - (n + .5 - .25*x)*u[1]/(x*M2)
        u[4] = -u[3]/(M*x) + ( -(1+r)*u[2]/6/M2 - (n+.5-.5*x)*u[1]/M3 ) / p
        up[1] = u[2]; up[2] = 2*u[3]*M; up[3] = 3*u[4]*M

        for k = 2:(m - 2)
              u[k+3] = ( -x*(2*k+1)*(k+1)*u[k+2]/M - (k*k+r)*u[k+1]/M2 -
                (n+.5-.5*x)*u[k]/M3 + .25*u[k-1]/M4 ) / (p*(k+2)*(k+1))
              up[k+2] = (k+2)*u[k+3]*M
        end
        up[m+1] = 0

        # Flip for more accuracy in inner product calculation.
        u = u[m+1:-1:1]
        up = up[m+1:-1:1]

        # Newton iteration
        hh = hh1
        hh[m+1] = M
        step = Inf
        l = 0
        if M == 1
            Mhzz = M*h + zz
            hh = [M ; cumprod( Mhzz )]
            hh = hh[end:-1:1]
        end
        while ( (abs(step) > eps(Float64)) && (l < 10) )
            l = l + 1
            step = dot(u,hh)/dot(up,hh)
            h = h - step
            Mhzz = (M*h) + zz
            # Powers of h (This is the fastest way!)
            hh = [M ; cumprod(Mhzz)]
            # Flip for more accuracy in inner product
            hh = hh[end:-1:1]
        end

        # Update
        x = x + h
        roots[j+1] = x
        ders[j+1] = dot(up,hh)
    end
    roots, ders
end

function BoundaryLaguerreGLR(n::Int64, xs)

    u, up = eval_Lag(xs, n)
    theta = atan(sqrt(xs/(n + .5 - .25*xs))*up/u)
    x1 = rk2_Lag(theta, -pi/2, xs, n)

    # Newton iteration
    step = Inf
    l = 0
    while ( (abs(step) > eps(Float64) || abs(u) > eps(Float64)) && (l < 200) )
        l = l + 1
        u, up = eval_Lag(x1, n)
        step = u/up
        x1 = x1 - step
    end

    ignored, d1 = eval_Lag(x1, n)

    x1, d1
end

function eval_Lag(x, n::Int64)
    # Evauate Laguerre polynomial via recurrence

    L = 0.
    Lp = 0.
    Lm2 = 0.
    Lm1 = exp(-x/2)
    Lpm2 = 0.
    Lpm1 = 0.
    for k = 0:n-1
        L = ( (2*k+1-x).*Lm1 - k*Lm2 ) / (k + 1)
        Lp = ( (2*k+1-x).*Lpm1 - Lm1 - k*Lpm2 ) / (k + 1)
        Lm2 = Lm1
        Lm1 = L
        Lpm2 = Lpm1
        Lpm1 = Lp
    end
    L, Lp
end

function rk2_Lag(t, tn, x, n::Int64)
    # Runge-Kutta for Laguerre Equation

    m = 10
    h = (tn - t)/m
    for j = 1:m
        f1 = (n + .5 - .25*x)
        k1 = -h/( sqrt(abs(f1/x)) + .25*(1/x-.25/f1)*sin(2*t) )
        t = t + h
        x = x + k1
        f1 = (n + .5 - .25*x)
        k2 = -h/( sqrt(abs(f1/x)) + .25*(1/x-.25/f1)*sin(2*t) )
        x = x + .5*(k2 - k1)
    end
    x
end

############## Routines for the "gen(W)" algorithm for computing an arbitrary number of terms with general w(x) = x^alpha*exp(-qm*x^m) ##########################

function asyRHgen(n, compRepr, alpha, m, qm)

    T = ceil(Int64, 34/log(n) ) # Heuristic for number of terms, should be scaled by the logarithm of eps(Float64) over the machine precision.
    UQ0 = getUQ(alpha, qm ,m, T)
    if compRepr
        mn = min(n,ceil(Int64, exp(exp(1/m)*1.05)*n^(1-1/2/m) ))
    else
        mn = n
    end
    itric = max(ceil(Int64, 3.6*n^0.188), 7)
    # Heuristics to switch between Bessel, extrapolation and Airy initial guesses.
    igatt = ceil(Int64, mn + 1.31*n^0.4 - n)

    A = zeros(m+1)
    for k =0:m
        A[k+1] = prod((2*(1:k)-1)/2/(1:k))
    end
    softEdge = (n*2/m/qm/A[m+1] )^(1/m)
    # Use finite differences for derivative of polynomial when not x^alpha*exp(-x) and use other initial approximations
    useFinDiff = (m != 1) || (qm != 1.0)
    bes = besselroots(alpha, itric).^2 # [Tricomi 1947 pg. 296]
    w = zeros(mn)
    if useFinDiff
        x = [bes*(2*m-1)^2/16/m^2/n^2*softEdge ; zeros(mn-itric) ]
    else
        ak = [-13.69148903521072; -12.828776752865757; -11.93601556323626;    -11.00852430373326; -10.04017434155809; -9.02265085340981; -7.944133587120853;    -6.786708090071759; -5.520559828095551; -4.08794944413097; -2.338107410459767]
        t = 3*pi/2*( (igatt:-1:12)-0.25) # [DLMF (9.9.6)]
        ak = [-t.^(2/3).*(1 + 5/48/t.^2 - 5/36/t.^4 + 77125/82944/t.^6     -10856875/6967296/t.^8); ak[max(1,12-igatt):11] ]
        nu = 4*n+2*alpha+2 # [Gatteshi 2002 (4.9)]
        air = (nu+ak*(4*nu)^(1/3)+ ak.^2*(nu/16)^(-1/3)/5 + (11/35-alpha^2-12/175*ak.^3)/nu + (16/1575*ak+92/7875*ak.^4)*2^(2/3)*nu^(-5/3) -(15152/3031875*ak.^5+1088/121275*ak.^2)*2^(1/3)*nu^(-7/3))
        x = [ bes/(4*n + 2*alpha+2).*(1 + (bes + 2*(alpha^2 - 1) )/(4*n + 2*alpha+2)^2/3 ) ; zeros(mn - itric -max(igatt,0) ) ; air]
    end

    if !useFinDiff
        UQ1 = getUQ(alpha+1, qm, m, T)
        factor0 = 1-4im*4^alpha*sum((UQ0[1,2,1:(T-1),1, 2] + UQ0[1,2,1:(T-1),1])./n.^reshape(1:(T-1), (1,1,T-1)) )
        factor1 = 1-4im*4^(alpha+1)*sum((UQ1[1,2,1:(T-1),1, 2] + UQ1[1,2,1:(T-1),1])./n.^reshape(1:(T-1), (1,1,T-1)) )
        factorx = real(sqrt(factor0/factor1 )/2/(1 - 1/n)^(1+alpha/2))
        factorw = real( -(1 - 1/(n + 1) )^(n + 1+ alpha/2)*(1 - 1/n)^(1 + alpha/2)*exp(1 + 2*log(2) )*4^(1+alpha)*pi*n^alpha*sqrt(factor0*factor1)*(1 + 1/n)^(alpha/2) )
    end

    if ( alpha^2/n > 1 )
        warn("A large alpha may lead to inaccurate results because the weight is low and R(z) is not close to identity.")
    end
    noUnderflow = true
    for k = 1:mn
        if ( x[k] == 0 ) && useFinDiff # Use linear extrapolation for the initial guesses for robustness in generalised weights.
            x[k] = 2*x[k-1] -x[k-2]
        elseif ( x[k] == 0 ) # Use sextic extrapolation for the initial guesses.
            x[k] = 7*x[k-1] -21*x[k-2] +35*x[k-3] -35*x[k-4] +21*x[k-5] -7*x[k-6] +x[k-7]
        end
        step = x[k]
        l = 0 # Newton-Raphson iteration number
        ov = floatmax(Float64) # Previous/old value
        ox = x[k] # Old x
        # Accuracy of the expansions up to machine precision would lower this bound.
        while ( ( abs(step) > eps(Float64)*40*x[k] ) && ( l < 20) )
            l = l + 1
#            pe = polyAsyRHgen(n, x[k], alpha, T, qm, m, UQ0)
            pe = polyAsyRHgen(n, x[k], alpha, T, qm, m, UQ0,k)
            if (abs(pe) >= abs(ov)*(1-35*eps(Float64)) )
                # The function values do not decrease enough any more due to roundoff errors.
                x[k] = ox # Set to the previous value and quit.
                break
            end
            if useFinDiff
                hh = max(sqrt(eps(Float64))*x[k], sqrt(eps(Float64)) )
                step = pe*hh/(polyAsyRHgen(n, x[k]+hh, alpha, T, qm, m, UQ0,k) - pe)
            else
                # poly' = (p*exp(-Q/2) )' = exp(-Q/2)*(p' -p/2) with orthonormal p.
                step = pe/(polyAsyRHgen(n-1, x[k], alpha+1, T, qm, m, UQ1,k)*factorx - pe/2)
            end
            ox = x[k]
            x[k] = x[k] -step
            ov = pe
        end
        if ( x[k] < 0 ) || ( l == 20 ) || ( ( k != 1 ) && ( x[k - 1] >= x[k] ) ) || isnan(x[k])
            print(x[k], "=x[k], k=", k, ", l=", l, ", x[k-1]=", x[k-1], ", x[k-2]=", x[k-1], ", step=", step, ", ox=", ox, ", ov=", ov, ".\n") # Print some debugging information and throw an error.
            error("Newton method may not have converged.")
        elseif ( x[k] > softEdge)
            warn("Node is outside the support of the measure: inaccuracy is expected.")
        end
        if noUnderflow&& useFinDiff
            hh = max(sqrt(eps(Float64))*x[k], sqrt(eps(Float64)) )
            w[k] = hh/(polyAsyRHgen(n,x[k]+hh,alpha,T,qm,m,UQ0,k) -polyAsyRHgen(n,x[k],alpha,T,qm,m,UQ0,k))/polyAsyRHgen(n-1,x[k],alpha,T,qm,m,UQ0,k)/exp(qm*x[k]^m) # This leaves out a constant factor, given by a ratio of leading order coefficients and normalising constants
        elseif noUnderflow
            w[k] = factorw/polyAsyRHgen(n-1, x[k], alpha+1, T, qm, m, UQ1,k)/polyAsyRHgen(n+1, x[k], alpha, T,qm,m,UQ0,k)/exp( x[k] )
        end
        if noUnderflow && ( w[k] == 0 ) && ( k > 1 ) && ( w[k-1] > 0 ) # We could stop now as the weights underflow.
            if compRepr
	        x = x[1:k-1]
	        w = w[1:k-1]
	        return (x,w)
	    else
                noUnderflow = false
	    end
        end
    end
    x, w
end

# Compute the expansion of the orthonormal polynomial without e^(qm*x^m/2) nor a constant factor based on some heuristics.
function polyAsyRHgen(np, y, alpha, T::Int64, qm, m::Int64, UQ, k)
    if (qm == 1) && (m == 1)
        z = y/4/np
        mnxi = 2*np*( sqrt(z).*sqrt(1 - z) - acos(sqrt(z) ) ) # = -n*xin/i
    else
        A = zeros(m+1)
        for k =0:m
            A[k+1] = prod((2*(1:k)-1)/2/(1:k))
        end
        z = y/(np*2/m/qm/A[m+1] )^(1/m)
        # Also correct but much slower: Hn = 4*m/(2*m-1)*double(hypergeom([1, 1-m], 3/2-m, z))/m
        Hn = 2/A[m+1]*sum(z.^(0:m-1).*A[m-(0:m-1)])/m
        mnxi = np*(sqrt(z+0im).*sqrt(1-z+0im).*Hn/2 -2*acos(sqrt(z+0im)))
    end
    # We could avoid these tests by splitting the loop k=1:mn into three parts with heuristics for the bounding indices.
#    if y < sqrt(np)
    if k < sqrt(np)
        # The fixed delta in the Riemann-Hilbert Problem would mean this bound has to be proportional to n, but x(1:k) are O(1/n) so choose the bound in between them to make more use of the (cheap) expansion in the bulk.
        return asyBesselgen(np, z, alpha, T, qm, m, UQ, mnxi + pi*np, true)
#    elseif y > 3.7*np
    elseif k > 0.9*np
        # Use the expansion in terms of the (expensive) Airy function, although the corresponding weights will start underflowing for n >= 186 for standard associated Laguerre polynomials.
        return asyAirygen(np, z, alpha, T, qm, m, UQ, (mnxi*3im/2)^(2/3), true)
    end
    asyBulkgen(np, z, alpha, T, qm, m, UQ, mnxi)
end

function asyBulkgen(np, z, alpha, T::Int64, qm, m::Int64, UQ, mnxi)
    if T == 1
        return real( 2/(z+0im)^(1/4 + alpha/2)/(1 - z+0im)^(1/4)*cos(acos(2*z - 1+0im)*(1/2 + alpha/2) - mnxi - pi/4) )
    end
    R = [1 0]
    for k = 1:T-1
        for i = 1:ceil(Int64, 3*T/2)
            R = R + reshape((UQ[1,:,k,i,1]/(z-1)^i + UQ[1,:,k,i,2]/z^i)/np^k,(1,2))
        end
    end
    p = real( 2/(z+0im)^(1/4 + alpha/2)*(cos(acos(2*z-1+0im)*(1/2+alpha/2) - mnxi-pi/4)*R[1]-cos(acos(2*z-1+0im)*(-1/2+alpha/2)-mnxi-pi/4)*R[2]*1im*4^alpha)/(1 - z+0im)^(1/4) )
end

function asyBesselgen(np, z, alpha, T::Int64, qm, m::Int64, UQ, npb, useQ::Bool)
    if T == 1
        return real( sqrt(2*pi)*(-1)^np*sqrt(npb)/(z+0im)^(1/4+alpha/2)/(1 - z+0im)^(1/4)*(sin( (alpha + 1)/2*acos(2*z - 1+0im) - pi*alpha/2)*besselj(alpha,npb) + cos( (alpha + 1)/2*acos(2*z - 1+0im) - pi*alpha/2)*(besselj(alpha-1,npb) - alpha/(npb)*besselj(alpha, npb) ) ) )
    end
    R = (1+0im)*[1 0]
    # Use the series expansion of R because it is faster and we use asyBessel only very close to zero to have less calls to besselj.
    if useQ
        for k = 1:T-1
            for i = 1:min(size(UQ,4),9-k)
                R = R + reshape(UQ[1, :, k, i, 4]*z^(i-1)/np^k,(1,2))
            end
        end
    else
        d = z-1+0im
        phi = 2*z-1+2*sqrt(z)*sqrt(d)
        Rko = (1+0im)zeros(T-1,2)
        sL = (1+0im)zeros(2,2,T-1)
	for m = 1:T-1
            for i = 1:ceil(Int64,3*m/2)
                Rko[m,:] += UQ[1, :, m, i, 1]/d^i+UQ[1, :, m, i, 2]/z^i
            end
            sL[:,:,m] = brac(m-1,alpha)/2^(1+2*m)/(npb/2im/np)^m*( [2^(-alpha)  0 ; 0   2^(alpha)]*[sqrt(phi)    1im/sqrt(phi)  ;  -1im/sqrt(phi)   sqrt(phi)]/2/z^(1/2)/d^(1/2)*[(-phi)^(alpha/2)  0 ; 0   (-phi)^(-alpha/2) ]*[((-1)^m)/m*(alpha^2+m/2-1/4)     (m-1/2)*1im  ;  -((-1)^m)*(m-1/2)*1im    (alpha^2+m/2-1/4)/m]*[(-phi)^(-alpha/2)  0 ; 0   (-phi)^(alpha/2) ]*[sqrt(phi)     -1im/sqrt(phi) ;  1im/sqrt(phi)    sqrt(phi)]*[2^(alpha)   0 ; 0 2^(-alpha)] -mod(m+1,2)*(4*alpha^2+2*m-1)/m*eye(2,2) )
        end
        for k = 1:T-1
            R += reshape((Rko[k,:] - sL[1,:,k])/np^k,(1,2))
            for m = 1:k-1
                R -= reshape(reshape(Rko[k-m,:],(1,2))*sL[:,:,m]/np^k,(1,2))
            end
        end
    end
    p = real( sqrt(2*pi)*(-1)^np*sqrt(npb)/(z+0im)^(1/4+alpha/2)/(1 - z+0im)^(1/4)*( (sin( (alpha + 1)/2*acos(2*z - 1+0im) - pi*alpha/2)*R[1] -sin( (alpha - 1)/2*acos(2*z - 1+0im) - pi*alpha/2)*R[2]*1im*4^alpha)*besselj(alpha, npb) + (cos( (alpha + 1)/2*acos(2*z - 1+0im)- pi*alpha/2)*R[1] - cos( (alpha - 1)/2*acos(2*z - 1+0im) - pi*alpha/2)*R[2]*1im*4^alpha)*(besselj(alpha-1, npb) - alpha/npb*besselj(alpha, npb) ) ) )
end

function asyAirygen(np, z, alpha, T::Int64, qm, m::Int64, UQ, fn, useQ::Bool, xin=NaN+NaN*1im)
    d = z - 1.0 +0im
    if T == 1
        return real( 4*sqrt(pi)/(z+0im)^(1/4+alpha/2)/d^(1/4)*(cos( (alpha + 1)/2*acos(2*z - 1+0im) )*fn^(1/4)*airyai(fn) -1im*sin( (alpha + 1)/2*acos(2*z - 1+0im) )*ifelse(angle(z-1) <= 0, -one(z), one(z) )*fn^(-1/4)*airyaiprime(fn) ) )
    end
    R = (1+0im)*[1 0]
    if useQ
        for k = 1:T-1
            for i = 1:min(size(UQ,4),9-k)
                R = R + reshape(UQ[1, :, k, i, 3]*d^(i-1)/np^k,(1,2))
            end
        end
    else
        phi = 2*z-1+2*sqrt(z)*sqrt(d)
        Rko = (1+0im)zeros(T-1,2)
        sR = (1+0im)zeros(2,2,T-1)
	for m = 1:T-1
            for i = 1:ceil(Int64,3*m/2)
                Rko[m,:] += UQ[1, :, m, i, 1]/d^i+UQ[1, :, m, i, 2]/z^i
            end
            sR[:,:,m] = nuk(m)/xin^m*( [2^(-alpha)  0 ; 0   2^(alpha)]*[sqrt(phi)    1im/sqrt(phi)  ;  -1im/sqrt(phi)   sqrt(phi)]/8/z^(1/2)/d^(1/2)*[phi^(alpha/2)  0 ; 0   phi^(-alpha/2) ]*[(-1.0)^m  -m*6im ; 6im*m*(-1)^m    1.0]*[phi^(-alpha/2)  0 ; 0   phi^(alpha/2) ]*[sqrt(phi)     -1im/sqrt(phi) ;  1im/sqrt(phi)    sqrt(phi)]*[2^(alpha)   0 ; 0 2^(-alpha)] - mod(m+1,2)*eye(2,2) )
        end
        for k = 1:T-1
            R += reshape((Rko[k,:] -sR[1,:,k])/np^k,(1,2))
            for m = 1:k-1
                R -= reshape(reshape(Rko[k-m,:],(1,2))*sR[:,:,m]/np^k,(1,2))
            end
        end
    end
    p = real( 4*sqrt(pi)/(z+0im)^(1/4+alpha/2)/d^(1/4)*( (R[1]*cos( (alpha + 1)/2*acos(2*z - 1+0im) ) -cos( (alpha - 1)/2*acos(2*z - 1+0im) )*R[2]*1im*4^alpha)*fn^(1/4)*airyai(fn) + 1im*(-sin( (alpha + 1)/2*acos(2*z - 1+0im) )*R[1] +sin( (alpha - 1)/2*acos(2*z - 1+0im) )*R[2]*1im*4^alpha)*ifelse(angle(z-1) <= 0, -one(z), one(z) )*fn^(-1/4)*airyaiprime(fn) ) )
end

# Additional short functions
function poch(x,n) # pochhammer does not seem to exist yet in Julia
    p = prod(x+(0:(n-1)) )
end
function binom(x,n) # binomial only works for integer x
    b = 1.0
    for i = 1:n
        b *= (x-(n-i))/i
    end
    b
end
function nuk(n)
    nu = -gamma(3*n-1/2)*2^n/27^n/2/n/sqrt(pi)/gamma(n*2)
end
function brac(n,alpha)
    b = prod(4*alpha^2-(2*(1:n)-1).^2 )/(2^(2*n)*gamma(1.0+n))
end

# Compute the W or V-matrices to construct the asymptotic expansion of R.
# Input
#   alpha, qm, m - Factors in the weight function w(x) = x^alpha*exp(-qm*x^m)
#   maxOrder     - The maximum order of the error
#   r            - 1 when computing Wright, -1 when computing Wleft
#   isW          - Whether we compute W(iso V)-matrices
# Output
#   WV           - Coefficient matrices for (z + 1/2 \pm 1/2)^m of Delta_k(z) or s_k(z)
function getV(alpha,qm,m::Int64,maxOrder::Int64,r)
    mo = ceil(Int64, 3*maxOrder/2) + 4
    ns = 0:mo
    f = NaN*zeros(mo+1) # Coefficients in the expansion of \bar{phi}_n(z) or \xi_n(z)
    g = NaN*zeros(maxOrder-1,mo+1)

    A = zeros(m+1)
    for k =0:m
        A[k+1] = prod((2*(1:k)-1.0)/2/(1:k))
    end
    if (r == 1) # Right disk: near z=1
        f = NaN*zeros(mo+2)
        ns = [ns; ns[mo+1]+1] # Extend by one because f(1) = 0 while not for left
        u = zeros(ns[mo+2]+1,ns[mo+2]+2)
        v = zeros(ns[mo+2]+1,ns[mo+2]+2)
        u[1,1] = 1.0
        v[1,1] = 1.0
        for n = [ns; ns[mo+2]+1]
            u[2,n+1] = binom(1/2,n+1)
            v[2,n+1] = binom(1/2,n+2)
        end
        for kt = 2:ns[mo+2]
            for n = ns
                u[kt+1,n+1] = sum(u[kt,(0:n)+1].*u[2,n-(0:n)+1])
                v[kt+1,n+1] = sum(v[kt,(0:n)+1].*v[2,n-(0:n)+1])
            end
        end
        q = zeros(ns[mo+2]+1)
        rr = zeros(ns[mo+2]+1) # Coeffs in the expansion of sqrt(2-2*sqrt(1-w))
        for kt = ns
            for l = 0:kt
                q[kt+1] = q[kt+1] + poch(1/2,kt-l)*u[kt-l+1,l+1]/(-2)^(kt-l)/gamma(1.0+kt-l)/(1+2*(kt-l) )
                rr[kt+1] = rr[kt+1] + binom(1/2,kt-l)*v[kt-l+1,l+1]*2^(kt-l)
            end
        end
        if (m == 1)
            for n = ns
                f[n+1,1] = -2*binom(1/2,n)
                for l = 0:n
                    f[n+1,1] = f[n+1,1] + 2*q[l+1]*rr[n-l+1]
                end
            end
        else
            for j = ns
                f[j+1,1] = 0
                for i=0:min(j,m-1)
                    f[j+1,1] = f[j+1,1] + binom(1/2,j-i)*(-1)^(m-i-1)*gamma(-1/2-i)/gamma(1/2-m)/gamma(m-i)
                end
                f[j+1,1] = -f[j+1,1]/m/A[m+1]
                for l = 0:j
                    f[j+1,1] = f[j+1,1] + 2*q[l+1]*rr[j-l+1]
                end
            end
        end
        if(abs(f[1]) > 10*eps(Float64) )
            error("xi_n should be O( (z-1)^(3/2) ): Expected f[1] to be zero")
        end
        ns = ns[1:mo+1] # Reset ns to its value before computing f's
        g[1,1,1] = -1/f[2,1]
        for n = 1:mo
            g[1,n+1,1] = -sum(reshape(g[1,1:n,1],(n,1)).*reshape(f[(n+2):-1:3,1],(n,1) ) )/f[2,1]
        end
    else # Left disk: near z=0
        if (m == 1)
            for n = ns
                f[n+1,1] = -(binom(1/2,n)*(-1)^n + poch(1/2,n)./(1+2*n)./gamma(1.0+n))
            end
        else
            for n = ns
                f[n+1,1] = 0.0
                for k = 0:min(m-1,n)
                    f[n+1,1] = f[n+1,1] + binom(1/2,n-k)*(-1)^(n-k)*A[m-k]
                end
                f[n+1,1] = -f[n+1,1]/2/m/A[m+1]-poch(1/2,n)./(1+2*n)./gamma(1.0+n)
            end
        end
        g[1,1,1] = 1/f[1,1]
        for n = 1:mo
            g[1,n+1,1] = -sum(reshape(g[1,1:n,1],(n,1)).*reshape(f[(n+1):-1:2,1],(n,1)) )/f[1,1]
        end
    end
    rho = (1+1im)*zeros(2*mo+3,mo+1)
    for n = ns
        rho[2,n+1] = poch(1/2,n)/gamma(1.0+n)/(1+2*n)*(-r)^n
    end
    rho[1,1] = 1
    for i = 2:(maxOrder-1)
        for n = ns
            g[i,n+1] = sum(g[i-1,1:(n+1) ].*g[1,(n+1):-1:1] )
        end
    end
    for i = 2:(mo*2+2)
        for n = ns
            rho[i+1,n+1] = sum(rho[i,1:(n+1) ].*rho[2,(n+1):-1:1] )
        end
    end
    OmOdd = (1+1im)*zeros(mo+1); OmEven = (1+1im)*zeros(mo+1)
    XiOdd = (1+1im)*zeros(mo+1); XiEven = (1+1im)*zeros(mo+1)
    ThOdd = (1+1im)*zeros(mo+1); ThEven = (1+1im)*zeros(mo+1)
    OmO = (1+1im)*zeros(mo+1); OmE = (1+1im)*zeros(mo+1)
    XiO = (1+1im)*zeros(mo+1); XiE = (1+1im)*zeros(mo+1)
    ThO = (1+1im)*zeros(mo+1); ThE = (1+1im)*zeros(mo+1)
    for n = ns
        js = 0:n
        for j = js
            OmOdd[n+1] = OmOdd[n+1] + (-1)^j/gamma(1.0+2.0*j)*(-2*alpha/sqrt(-r+0.0im))^(2*j)*rho[2*j+1,n-j+1]
            XiOdd[n+1] = XiOdd[n+1] + (-1)^j/gamma(1.0+2.0*j)*(-2*(alpha+1)/sqrt(-r+0.0im))^(2*j)*rho[2*j+1,n-j+1]
            ThOdd[n+1] = ThOdd[n+1] + (-1)^j/gamma(1.0+2.0*j)*(-2*(alpha-1)/sqrt(-r+0.0im))^(2*j)*rho[2*j+1,n-j+1]
            OmEven[n+1] = OmEven[n+1] + (-1)^j/gamma(1.0+2*j+1.0)*(-2*alpha/sqrt(-r+0.0im))^(2*j+1)*rho[2*j+2,n-j+1]
            XiEven[n+1] = XiEven[n+1] + (-1)^j/gamma(1.0+2*j+1.0)*(-2*(alpha+1)/sqrt(-r+0.0im))^(2*j+1)*rho[2*j+2,n-j+1]
            ThEven[n+1] = ThEven[n+1] + (-1)^j/gamma(1.0+2*j+1.0)*(-2*(alpha-1)/sqrt(-r+0.0im))^(2*j+1)*rho[2*j+2,n-j+1]
        end
        for j = js
            OmO[n+1] = OmO[n+1] + binom(-1/2,j)*(r)^j*OmOdd[n-j+1]
            XiO[n+1] = XiO[n+1] + binom(-1/2,j)*(r)^j*XiOdd[n-j+1]
            ThO[n+1] = ThO[n+1] + binom(-1/2,j)*(r)^j*ThOdd[n-j+1]
            OmE[n+1] = OmE[n+1] + binom(-1/2,j)*(r)^j*OmEven[n-j+1]
            XiE[n+1] = XiE[n+1] + binom(-1/2,j)*(r)^j*XiEven[n-j+1]
            ThE[n+1] = ThE[n+1] + binom(-1/2,j)*(r)^j*ThEven[n-j+1]
        end
    end
    Ts = (1+1im)*zeros(2,2,mo+1) # = G_{k,n}^{odd/even} depending on k, overwritten on each new k
    WV = (1+1im)*zeros(2,2,maxOrder-1,mo+1)
    for k = 1:(maxOrder-1)
        Ts[:,:,:] = 0
        if r == 1
            if mod(k,2) == 1
                for n = 0:mo
                    Ts[:,:,n+1] = nuk(k)*[-2*(2*binom(-1/2,n-1)*(n>0)+binom(-1/2,n))     2im*4^(-alpha)*binom(-1/2,n)    ;    2im*4^(alpha)*binom(-1/2,n)    (2*(2*binom(-1/2,n-1)*(n>0) +binom(-1/2,n)))] -6*k*nuk(k)*[-2*OmO[n+1]   4^(-alpha)*2im*XiO[n+1]  ;   4^(alpha)*2im*ThO[n+1]    2*OmO[n+1]]
                    WV[:,:,k,n+1] = sum(repeat(reshape(g[k,1:(n+1) ], (1,1,n+1) ), outer=[2,2,1]).*Ts[:,:,(n+1):-1:1],3)/8
                end
            else
                for n = 0:mo
                     Ts[:,:,n+1] = nuk(k)*4*(n==0)*eye(2) +6*k*nuk(k)*[-2im*OmE[n+1]    -2*4^(-alpha)*XiE[n+1]  ;   -2*4^alpha*ThE[n+1]   2im*OmE[n+1]]
                     WV[:,:,k,n+1] = sum(repeat(reshape(g[k,1:(n+1) ], (1,1,n+1) ), outer=[2,2,1]).*Ts[:,:,(n+1):-1:1],3)/8
                end
            end
        else
            if mod(k,2) == 1
                for n = 0:mo
                    Ts[:,:,n+1] = -(alpha^2+k/2-1/4)/k*[-(-1)^n*(2*binom(-1/2,n-1)*(n>0)+binom(-1/2,n))*2    -1im*4^(-alpha)*2*(-1)^n*binom(-1/2,n)  ;  -1im*4^(alpha)*2*(-1)^n*binom(-1/2,n)     ( (-1)^n*(2*binom(-1/2,n-1)*(n>0) +binom(-1/2,n))*2)] - (k-1/2)*[2*OmO[n+1]   4^(-alpha)*2im*XiO[n+1]  ;   4^(alpha)*2im*ThO[n+1]   -2*OmO[n+1]] # binom(-1/2,-1) should be zero
                    WV[:,:,k,n+1] = -(-1)^(ceil(Int64, k/2)+1)*(1im*sqrt(2))^k*(-2+0im)^(-k/2)/4^(k+1)*brac(k-1,alpha)*sum(repeat(reshape(g[k,1:(n+1) ], (1,1,n+1) ), outer=[2,2,1]).*Ts[:,:,(n+1):-1:1],3)
                end
            else
                for n = 0:mo
                    Ts[:,:,n+1] = (alpha^2+k/2-1/4)/k*4*(n==0)*eye(2)  -2*(k-1/2)*[ OmE[n+1]   4^(-alpha)*1im*XiE[n+1]  ;   4^alpha*1im*ThE[n+1]   -OmE[n+1] ]
                    WV[:,:,k,n+1] = -(-1)^(ceil(Int64, k/2)+1)*(1im*sqrt(2))^k*(-2)^(-k/2)/4^(k+1)*brac(k-1,alpha)*sum(repeat(reshape(g[k,1:(n+1) ], (1,1,n+1) ), outer=[2,2,1]).*Ts[:,:,(n+1):-1:1],3)
                end
            end
        end
    end
    WV
end

# Get the U-matrices to construct the asymptotic expansion of R using the procedure with the convolutions with a specified method.
# Input
#   alpha, qm, m - Parts of the weight function
#   maxOrder     - The maximal order of the error
# Output
#   UQ           - Coefficient matrices of R_k(z) for (z-1)^(-m) [Uright], or z^(-m) [Uleft] of R_k^{right}(z) for (z-1)^n [Qright] and of R_k^{left}(z) for z^n [Qleft]
function getUQ(alpha, qm, m::Int64, maxOrder::Int64)
    Vr = getV(alpha, qm, m, maxOrder, 1)
    Vl = getV(alpha, qm, m, maxOrder, -1)
    UQ = (1+1im)*zeros(2,2,maxOrder-1,ceil(Int64,3*maxOrder/2)+2, 4)
    for kt = 0:(maxOrder-2)
        # Uright(:,:,(maxOrder-1)+1,:) will not be used later on because first term in expansions is without U's
        for mt = 0:(ceil(Int64,3*(kt+1)/2)-1)
            UQ[:,:,kt+1,mt+1,1] = Vr[:,:,kt+1,ceil(Int64,3*(kt+1)/2)-mt]
            for j = 0:(kt-1)
                for l = 0:(ceil(Int64,3*(j+1)/2)-mt-1)
                    UQ[:,:,kt+1,mt+1,1] = UQ[:,:,kt+1,mt+1,1] + UQ[:,:,kt-j,l+1,3]*Vr[:,:,j+1,ceil(Int64,3*(j+1)/2)-l-mt]
                end
            end
        end
        for mt = 0:(ceil(Int64,(kt+1)/2)-1)
            UQ[:,:,kt+1,mt+1,2] = Vl[:,:,kt+1,ceil(Int64,(kt+1)/2)-mt]
            for j= 0:(kt-1)
                for l = 0:(ceil(Int64,(j+1)/2)-mt-1)
                    UQ[:,:,kt+1,mt+1,2] = UQ[:,:,kt+1,mt+1,2] + UQ[:,:,kt-j,l+1,4]*Vl[:,:,j+1,ceil(Int64,(j+1)/2)-l-mt]
                end
            end
        end
        for n = 0:(ceil(Int64,3*(maxOrder-kt+1)/2)-1)
            UQ[:,:,kt+1,n+1,3] = -Vr[:,:,kt+1,ceil(Int64,3*(kt+1)/2)+1+n]
            UQ[:,:,kt+1,n+1,4] = -Vl[:,:,kt+1,ceil(Int64,(kt+1)/2)+1+n]
            for i = 0:(ceil(Int64,(kt+1)/2)-1)
                UQ[:,:,kt+1,n+1,3] = UQ[:,:,kt+1,n+1,3] + binom(-i-1,n)*UQ[:,:,kt+1,i+1,2]
            end
            for i = 0:(ceil(Int64,3*(kt+1)/2)-1)
                UQ[:,:,kt+1,n+1,4] = UQ[:,:,kt+1,n+1,4] + binom(-i-1,n)*(-1.0)^(-i-1-n)*UQ[:,:,kt+1,i+1,1]
            end
            for j = 0:(kt-1)
                for l = 0:(ceil(Int64, (j+1)/2)+n)
                    UQ[:,:,kt+1,n+1,4] = UQ[:,:,kt+1,n+1,4] -UQ[:,:,kt-j,l+1,4]*Vl[:,:,j+1,n-l+1+ceil(Int64,(j+1)/2) ]
                end
                for l = 0:(ceil(Int64, 3*(j+1)/2)+n)
                    UQ[:,:,kt+1,n+1,3] = UQ[:,:,kt+1,n+1,3] -UQ[:,:,kt-j,l+1,3]*Vr[:,:,j+1,n-l+1+ceil(Int64, 3*(j+1)/2) ]
                end
            end
        end
    end
    UQ
end

########################## Routines for the explicit expansions ##########################

function laguerreExp( n::Int64, compRepr::Bool, alpha::Float64, T::Int64)

    if compRepr
        # Get a heuristic for the indices where the weights are about above floatmin.
        mn = min(ceil(Int64, 17*sqrt(n)), n)
    else
        mn = n
    end
    if ( alpha^2/n > 1 )
        warn("A large alpha may lead to inaccurate results because the weight is low and R(z) is not close to identity.")
    end
    d = 1/(4*n+2*alpha+2)

    itric = max(ceil(Int64, sqrt(n)), 7)
    # Heuristics to switch between Bessel, extrapolation and Airy initial guesses.

    iair = floor(Int64,0.9*n)
    pt = (4*n .- 4*((itric+1):min(mn, iair-1)) .+ 3)*d

    t = pi^2/16 * (pt .- 1).^2
    # This t is not a very good initial approximation of the inverse function of f(y) = (4*n -4*k +3)*d*pi +2*sqrt(y)*sqrt(1-y) -acos(2*y-1)
    for it = 1:6
		t = @. t - (pt*pi + 2*sqrt(t-t^2) - acos(2t - 1) ) * sqrt(t/(1-t))/2
    end
    # A second option is to use the Roots package and t = fzero(f, 0, 1); for each k
    # A third option is to improve the initial approximations for k starting from ichebSer = Int64(round(n/10-0.45*alpha+0.3)) with barycentric interpolation on a chebyshev grid as
#    xk = 0.5+ 0.49*cos((2*(1:16) -1)/2/16*pi)'
#    baryc = [8.392374553903038e+02, -1.573904756067641e+04,  1.217801459223059e+05,  -5.241187241822941e+05,   1.561671400871400e+06, -3.632334022840342e+06,   7.056947602530619e+06,  -1.192343707574069e+07,   1.795854669340786e+07,  -2.446733743951402e+07,   3.036493740995492e+07,  -3.430589960928102e+07,   3.489836980231517e+07,  -3.097800484853177e+07,   2.192179297293409e+07,  -8.079248247153265e+06]
#    for mix = 1:16
#	t[ichebSer+1:end] = t[ichebSer+1:end].*(pt[ichebSer+1:end] -xk[mix])
#	asum = asum + baryc[mix]./(pt[ichebSer+1:end] -xk[end])
#    end
#    t[ichebSer+1:end] = t[ichebSer+1:end].*asum
    # The third option ensures that only 3 Newton iterations are required. However, both are slower in total than the uncommented version, which still takes more than half of the total execution time of laguerreExp. Maybe a planned Clenshaw algorithm can improve this ...

    jak = besselroots(alpha, itric)
    bes = jak*0
    wbes = 0*bes

    bulk = t*0
    wbulk = t*0

    ak = [-13.69148903521072; -12.828776752865757; -11.93601556323626;    -11.00852430373326; -10.04017434155809; -9.02265085340981; -7.944133587120853;    -6.786708090071759; -5.520559828095551; -4.08794944413097; -2.338107410459767]

    tair = 3*pi/2*( ((mn-iair+1):-1:12).-0.25)
    ak = [-tair.^(2/3).*(1 .+ 5/48 ./tair.^2 - 5/36 ./tair.^4 .+ 77125/82944 ./tair.^6 -10856875/6967296 ./tair.^8); ak[max(1,12-mn+iair-1):11] ]
    air = 0*ak

    if (T >= 9)
		bes = bes + (10644*jak^8 + 60*(887*alpha^2 - 2879)*jak^6 + (125671*alpha^4 -729422*alpha^2 + 1456807)*jak^4 + 3*(63299*alpha^6 - 507801*alpha^4 + 1678761*alpha^2 - 2201939)*jak^2 + 2*(107959*alpha^8 - 1146220*alpha^6 + 5095482*alpha^4 -10087180*alpha^2 + 6029959) )*d^8/42525
		wbes = wbes + (215918*alpha^8 + 53220*jak.^8 + 240*(887*alpha^2 - 2879)*jak.^6 -2292440*alpha^6 + 3*(125671*alpha^4 - 729422*alpha^2 + 1456807)*jak.^4 + 10190964*alpha^4 + 6*(63299*alpha^6 - 507801*alpha^4 + 1678761*alpha^2 -2201939)*jak.^2 - 20174360*alpha^2 + 12059918)/42525*d^8

		bulk = bulk + d^7/10886400*(1-t).^3/t.^3*(43222750000*(1-t).^(-14) - 241928673000*(1-t).^(-13) + 566519158800*(1-t).^(-12) -714465642135*(1-t).^(-11) + 518401904799*(1-t).^(-10) + 672*(12000*alpha^4 - 24000*alpha^2 +64957561)*(1-t).^(-8) - 212307298152*(1-t).^(-9) + 24883200*alpha^8 - 192*(103425*alpha^4 -206850*alpha^2 + 15948182)*(1-t).^(-7) + 3360*(4521*alpha^4 - 9042*alpha^2 - 7823)*(1-t).^(-6) -232243200*alpha^6 - 1792*(3375*alpha^6 - 13905*alpha^4 + 17685*alpha^2 - 1598)*(1-t).^(-5)+ 16128*(450*alpha^6 - 2155*alpha^4 + 2960*alpha^2 - 641)*(1-t).^(-4) + 812851200*alpha^4 -768*(70875*alpha^8 - 631260*alpha^6 + 2163630*alpha^4 - 2716980*alpha^2 +555239)*(1-t).^(-3) + 768*(143325*alpha^8 - 1324260*alpha^6 + 4613070*alpha^4 -5826660*alpha^2 + 1193053)*(1-t).^(-2) - 1028505600*alpha^2 - 5806080*(15*alpha^8 -140*alpha^6 + 490*alpha^4 - 620*alpha^2 + 127)*(1-t).^(-1) + 210677760)
	#Computing this term for wbulk takes too much memory
    end

    if (T >= 7)
        # These higher order terms in the left and bulk region are derived in [Opsomer 2017, in preparation]
		bes = bes + (657*jak.^6 +36*jak.^4*(73*alpha^2-181) +2*jak.^2*(2459*alpha^4 -10750*alpha^2 +14051) + 4*(1493*alpha^6 -9303*alpha^4 +19887*alpha^2 - 12077) )*d^6/2835

		wbes = wbes + (1493*alpha^6 + 657*jak.^6 + 27*(73*alpha^2 - 181)*jak.^4 - 9303*alpha^4 + (2459*alpha^4 -10750*alpha^2 + 14051)*jak.^2 + 19887*alpha^2 - 12077)*4/2835*d^6

		bulk = bulk -d^5/181440*(1-t).^2/t.^2*(10797500*(1-t).^(-10) - 43122800*(1-t).^(-9) + 66424575*(1-t).^(-8) -48469876*(1-t).^(-7) + 193536*alpha^6 + 16131880*(1-t).^(-6) + 80*(315*alpha^4 - 630*alpha^2 -221)*(1-t).^(-4) - 1727136*(1-t).^(-5) - 967680*alpha^4 - 320*(63*alpha^4 - 126*alpha^2 +43)*(1-t).^(-3) + 384*(945*alpha^6 - 4620*alpha^4 + 6405*alpha^2 - 1346)*(1-t).^(-2) +1354752*alpha^2 - 23040*(21*alpha^6 - 105*alpha^4 + 147*alpha^2 - 31)*(1-t).^(-1) -285696)

		wbulk = wbulk - (1-t).^3/90720/t.^3*d^6*(43190000*(1-t).^(-12) -204917300*(1-t).^(-11) + 393326325*(1-t).^(-10) - 386872990*(1-t).^(-9) + 201908326*(1-t).^(-8) +80*(315*alpha^4 - 630*alpha^2 + 53752)*(1-t).^(-6) - 50986344*(1-t).^(-7) - 320*(189*alpha^4 -378*alpha^2 - 89)*(1-t).^(-5) + 480*(63*alpha^4 - 126*alpha^2 + 43)*(1-t).^(-4) -384*(315*alpha^6 - 1470*alpha^4 + 1995*alpha^2 - 416)*(1-t).^(-3) + 2304*(21*alpha^6 -105*alpha^4 + 147*alpha^2 - 31)*(1-t).^(-2) )
    end

    if (T >= 5)
		bes = bes + (11*jak.^4 +3*jak.^2*(11*alpha^2-19) +46*alpha^4 -140*alpha^2 +94)*d^4/45
		wbes = wbes + (46*alpha^4 + 33*jak.^4 +6*jak.^2*(11*alpha^2 -19) -140*alpha^2 +94)/45*d^4

		air = air -(15152/3031875*ak.^5+1088/121275*ak.^2)*2^(1/3)*d^(7/3) # [Gatteshi 2002 (4.9)], Gives an O(n^{-4}) relative error

		bulk = bulk +d^3*(1-t)/t/720*(1600*(1-t).^(-6) - 3815*(1-t).^(-5) + 480*alpha^4 +2814*(1-t).^(-4) - 576*(1-t).^(-3) - 960*alpha^2 - 48*(15*alpha^4 - 30*alpha^2 + 7)*(1-t).^(-1) -16*(1-t).^(-2) + 224)
		wbulk = wbulk +(1-t).^2/720/t.^2*d^4*(8000*(1-t).^(-8) - 24860*(1-t).^(-7) + 27517*(1-t).^(-6) - 12408*(1-t).^(-5) + 1712*(1-t).^(-4) +16*(15*alpha^4 - 30*alpha^2 + 7)*(1-t).^(-2) + 32*(1-t).^(-3))
    end

    if (T >= 3)
        # From here, the terms are also in [Tricomi 1947 pg. 296] and
		bes = bes + (jak.^2 + 2*alpha^2 - 2)*d^2/3
		wbes = wbes + (alpha^2 + jak.^2 -1)*2/3*d^2

		air = air +  ak.^2*(d*16)^(1/3)/5 + (11/35-alpha^2-12/175*ak.^3)*d +  (16/1575*ak+92/7875*ak.^4)*2^(2/3)*d^(5/3)
		bulk = bulk -d/12*(12*alpha^2 + 5*(1-t).^(-2) - 4*(1-t).^(-1) - 4)
		wbulk = wbulk - d^2/6*(5*(1-t).^(-3) - 2*(1-t).^(-2) )
    end
    bes = jak.^2*d.*(1 + bes )
    air = 1/d +ak*(d/4)^(-1/3) + air
    bulk = bulk + t/d

    w = [ 4*d*bes.^alpha.*exp(-bes)./(besselj(alpha-1, jak)).^2 .*(1+ wbes); bulk.^alpha.*exp(-bulk)*2*pi.*sqrt(t./(1-t)).*(1 +wbulk) ; 2^(2*alpha+1/3)*n^(alpha+1/3)*exp.(-air)./(airyaiprime.(ak)).^2]
    # For the Airy region, only O(n^{-4}) relative error for x and O(n^{-2/3}) for w as the latter are extremely small or even underflow
    x = [ bes; bulk ; air]

    if ( minimum(x) < 0.0 ) || ( maximum(x) > 4*n + 2*alpha + 2 ) ||  ( (minimum(diff(x)) <= 0.0) && (T == ceil(Int64, 34/log(n) )) && (n <= 1e5) ) || (minimum(w) < 0.0) # Higher errors when too few terms: when switching from bulk to air and avoid error for extremely large n: previously (n <= 5e6)
		print(minimum(x), " =min, max= ", maximum(x), "=minmax x, ", 4*n + 2*alpha + 2, ", ",  minimum(diff(x)), "=mindiff, minw=",  minimum(w), '\n' )
		print(bulk[end-100:end], "=end bulk, start air = ", air[1:100], ", ", findmin(diff(x)))
		error("Wrong node or weight.")
    end

    if compRepr
		k = findlast(w)
		x = x[1:k]
		w = w[1:k]
    end
    x, w
end


text:
struct NaNChecker{T, F} <: AbstractDiagnostic
    schedule :: T
      fields :: F
end

"""
    NaNChecker(; schedule, fields)

Returns a `NaNChecker` that checks for a `NaN` anywhere in `fields` when `schedule` actuates.
`fields` should be a named tuple. The simulation is aborted if a `NaN` is found.
"""
NaNChecker(model=nothing; schedule, fields) = NaNChecker(schedule, fields)

function run_diagnostic!(nc::NaNChecker, model)
    for (name, field) in pairs(nc.fields)
        if any(isnan, field.data.parent)
            t, i = model.clock.time, model.clock.iteration
            error("time = $t, iteration = $i: NaN found in $name. Aborting simulation.")
        end
    end
end

ext: jl


text:
using jInvSeismic.FWI
using jInvSeismic.Utils
using Helmholtz
using jInv.Mesh
using jInv.Utils
using jInv.ForwardShare
using Multigrid.ParallelJuliaSolver
using Distributed

using DelimitedFiles

if nworkers() == 1
	addprocs(16);
elseif nworkers() < 16
	addprocs(16 - nworkers());
end

@everywhere begin
using jInv.InverseSolve
using jInv.LinearSolvers
using jInvSeismic.FWI
using jInv.Mesh
using jInv.Utils
using DelimitedFiles
using jInv.ForwardShare
using KrylovMethods
end

dim     = 2;
pad     = 30;
newSize = [600,300];

(m,Mr,mref,boundsHigh,boundsLow) = readModelAndGenerateMeshMref("examples","SEGmodel2Dsalt.dat",dim,pad,[0.0,13.5,0.0,4.2],newSize,1.752,2.7);
m = 1 ./ (m.^2);
mref = 1 ./ (mref.^2);

# attenuation for BC
padx = 34; padz = 34
a    = 2.0;
xc = getCellCenteredGrid(Mr);
gamma = getHelmholtzABL(Mr,true,[padx;padz],a);
# gamma = Helmholtz.getABL(Mr,true,[padx;padz],a);
println("Size gamma: ", size(gamma))
println("Size Mr: " , Mr.n)
println("Size m: " , size(m))
# parameters for the Helmholtz (units in km)
h = Mr.h;
n = Mr.n;
# omega = 0.3*2*pi*[2.0;5.0; 10.0;20.0;25.0]
omega = 2*pi*[2.0;2.5;3.5;4.5;6.0;]
nfreq = length(omega)
# generate sources
q = zeros(tuple(n.+1...)); q[padx+1:4:end-padx-1,1] .= 1e4
print(size(q))
Q = sdiag(vec(q))
print(size(Q))
Q = Q[:,(LinearIndices(sum(Q,dims=2) .!= 0))[findall(sum(Q,dims=2) .!= 0)]]
print(size(Q))
nsrc = size(Q,2)
# receivers
p = zeros(tuple(n.+1...)); p[padx+1:end-padx-1,1] .= 1
P = sdiag(vec(p))
P = P[:,(LinearIndices(sum(P,dims=2) .!= 0))[findall(sum(P,dims=2) .!= 0)]]
nrec = size(P,2)
Ainv  = getParallelJuliaSolver(ComplexF64,Int64,numCores=4,backend=1);

(pForp,contDivFWI,SourcesSubIndFWI)  = getFWIparam(omega,ones(ComplexF64,size(omega)),vec(gamma),Q,P,Mr,Ainv,workers())
# pForp = getBasicFWIparam(omega,gamma,Q,P,Mr,true)

# P = nothing
# Q = nothing
# inversion mesh and forward mesh are the same here
M2Mp = ones(length(pForp))


text:
export TotalDegreeSystem

"""
    TotalDegreeSystem(polynomials, vars) <: AbstractSystem

Create a tottal degree system
"""
struct TotalDegreeSystem{T} <: AbstractSystem
    degrees::Vector{Int}
    degree_idxs::Vector{Int}
    hom_idx::Int # 0 if affine
    scaling_factors::Vector{T}

    function TotalDegreeSystem(
        degrees::Vector{Int},
        degree_idxs::Vector{Int},
        hom_idx::Int = 0,
        scaling_factors::Vector{T} = ones(length(degrees));
        affine = false,
    ) where {T}
        @assert length(degrees) == length(degree_idxs) == length(scaling_factors)
        if affine
            hom_idx = 0
        else
            @assert 1 ≤ hom_idx ≤ length(degrees) + 1
        end
        new{T}(degrees, degree_idxs, hom_idx, scaling_factors)
    end
end
function TotalDegreeSystem(F::Vector{<:MP.AbstractPolynomialLike}, args...; kwargs...)
    TotalDegreeSystem(MP.maxdegree.(F), args...; kwargs...)
end
function TotalDegreeSystem(degrees::Vector{Int}, args...; kwargs...)
    TotalDegreeSystem(
        degrees,
        collect(1:length(degrees)),
        length(degrees) + 1,
        args...;
        kwargs...,
    )
end

function Base.size(F::TotalDegreeSystem)
    n = length(F.degrees)
    F.hom_idx == 0 ? (n, n) : (n, n + 1)
end

cache(::TotalDegreeSystem, x) = SystemNullCache()

function evaluate!(u, F::TotalDegreeSystem, x::ProjectiveVectors.PVector, ::SystemNullCache)
    λ = F.scaling_factors
    @inbounds for i = 1:length(F.degrees)
        d = F.degrees[i]
        u[i] = λ[i] * (x[i]^d - x[end]^d)
    end
    u
end

function evaluate!(u, F::TotalDegreeSystem, x::AbstractVector, ::SystemNullCache)
    λ = F.scaling_factors
    @inbounds for i = 1:length(F.degrees)
        d = F.degrees[i]
        dix = F.degree_idxs[i]
        u[i] = λ[i] * x[dix]^d - λ[i]
    end
    u
end
function evaluate(F::TotalDegreeSystem, x, cache::SystemNullCache)
    u = similar(x, size(F, 1))
    evaluate!(u, F, x, cache)
    u
end

function jacobian!(U, F::TotalDegreeSystem, x::PVector{<:Number,1}, ::SystemNullCache)
    U .= zero(eltype(x))
    λ = F.scaling_factors
    N = length(x)
    @inbounds for i = 1:length(F.degrees)
        d = F.degrees[i]
        didx = F.degree_idxs[i]
        if d == 1
            U[i, didx] = λ[i]
            U[i, N] = -λ[i]
        elseif d > 1
            U[i, didx] = λ[i] * d * x[didx]^(d - 1)
            U[i, N] = -λ[i] * d * x[N]^(d - 1)
        end
    end
    U
end
function jacobian!(U, F::TotalDegreeSystem, x::AbstractVector, ::SystemNullCache)
    U .= zero(eltype(x))
    λ = F.scaling_factors
    @inbounds for i = 1:length(F.degrees)
        d = F.degrees[i]
        didx = F.degree_idxs[i]
        if d == 1
            U[i, didx] = λ[i]
        elseif d > 1
            U[i, didx] = λ[i] * d * x[didx]^(d - 1)
        end
    end
    U
end

function jacobian(F::TotalDegreeSystem, x, cache::SystemNullCache)
    U = similar(x, size(F))
    jacobian!(U, F, x, cache)
    U
end

function evaluate_and_jacobian!(
    u,
    U,
    F::TotalDegreeSystem,
    x::PVector{<:Number,1},
    ::SystemNullCache,
)
    U .= zero(eltype(x))
    N = length(x)
    λ = F.scaling_factors
    @inbounds for i = 1:length(F.degrees)
        d = F.degrees[i]
        didx = F.degree_idxs[i]
        if d == 1
            u[i] = λ[i] * (x[didx] - x[N])
            U[i, didx] = λ[i]
            U[i, N] = -λ[i]
        elseif d > 1
            xd = x[didx]^(d - 1)
            xh = x[N]^(d - 1)
            u[i] = λ[i] * (xd * x[didx] - xh * x[N])
            U[i, didx] = λ[i] * d * xd
            U[i, N] = -λ[i] * d * xh
        end
    end
    nothing
end
function evaluate_and_jacobian!(
    u,
    U,
    F::TotalDegreeSystem,
    x::AbstractVector,
    ::SystemNullCache,
)
    U .= zero(eltype(x))
    λ = F.scaling_factors
    @inbounds for i = 1:length(F.degrees)
        d = F.degrees[i]
        didx = F.degree_idxs[i]
        if d == 1
            u[i] = λ[i] * x[didx] - λ[i]
            U[i, didx] = λ[i]
        elseif d > 1
            xd = x[didx]^(d - 1)
            u[i] = λ[i] * xd * x[didx] - λ[i]
            U[i, didx] = λ[i] * d * xd
        end
    end
    nothing
end

function evaluate_and_jacobian(F::TotalDegreeSystem, x, cache::SystemNullCache)
    u = similar(x, size(F, 1))
    U = similar(x, size(F))
    evaluate_and_jacobian!(u, U, F, x, cache)
    u, U
end


text:
using FeedbackParticleFilters, Statistics

println("Testing PoissonEquation.jl:")

@testset "PoissonEquation.jl" begin

    h(x)      = 2x[1:2]
    pos       = randn(3, 10)
    H         = randn(2, 10)
    mean_H    = Statistics.mean(H, dims=2)
    potential = zeros(2, 10)
    gain      = zeros(3, 10, 2);

    print("  inner constructor")
    eq = PoissonEquation(h, pos, H, mean_H, potential, gain)
    print("-")
    @test eq isa PoissonEquation{Float64,typeof(h)}
    print("-")
    @test_throws DimensionMismatch PoissonEquation(h, randn(2, 10), H, mean_H, potential, gain) # wrong size of pos
    print("-")
    @test_throws DimensionMismatch PoissonEquation(h, randn(3, 11), H, mean_H, potential, gain) # wrong size of pos
    print("-")
    @test_throws DimensionMismatch PoissonEquation(h, pos, randn(3,10), mean_H, potential, gain) # wrong size of H
    print("-")
    @test_throws DimensionMismatch PoissonEquation(h, pos, randn(2,11), mean_H, potential, gain) # wrong size of H
    print("-")
    @test_throws DimensionMismatch PoissonEquation(h, pos, H, randn(3,1), potential, gain) # wrong size of mean_H
    print("-")
    @test_throws DimensionMismatch PoissonEquation(h, pos, H, randn(2,2), potential, gain) # wrong size of mean_H
    print("-")
    @test_throws DimensionMismatch PoissonEquation(h, pos, H, mean_H, randn(3,10), gain) # wrong size of potential
    print("-")
    @test_throws DimensionMismatch PoissonEquation(h, pos, H, mean_H, randn(2,11), gain) # wrong size of potential
    print("-")
    @test_throws DimensionMismatch PoissonEquation(h, pos, H, mean_H, potential, randn(2,10,2)) # wrong size of gain
    print("-")
    @test_throws DimensionMismatch PoissonEquation(h, pos, H, mean_H, potential, randn(3,11,2)) # wrong size of gain
    print("-")
    @test_throws DimensionMismatch PoissonEquation(h, pos, H, mean_H, potential, randn(3,10,3)) # wrong size of gain
    println("DONE")

    print("  outer constructor")
    eq2 = PoissonEquation(h, pos)
    print("-")
    @test eq2 isa PoissonEquation{Float64,typeof(h)}
    print("-")
    @test eq2.H == 2*pos[1:2, :]
    print("-")
    @test eq2.mean_H == Statistics.mean(eq2.H, dims=2)
    ens = UnweightedParticleEnsemble(pos)
    eq3 = PoissonEquation(h, ens)
    print("-")
    @test eq3.positions == pos
    print("-")
    @test eq3.H == eq2.H
    print("-")
    @test eq3.mean_H == eq2.mean_H
    println("DONE")

    print("  method state_dim")
    print("-")
    @test state_dim(eq) == 3
    println("DONE")

    print("  method obs_dim")
    print("-")
    @test obs_dim(eq) == 2
    println("DONE")

    print("  method Htilde")
    print("-")
    @test Htilde(eq2) == eq2.H .- eq2.mean_H
    println("DONE")

    print("  method update!")
    update!(eq)
    print("-")
    @test eq.positions == pos
    print("-")
    @test eq.H == 2*pos[1:2, :]
    print("-")
    @test eq.mean_H == Statistics.mean(eq.H, dims=2)
    pos = randn(3, 10)
    update!(eq, pos)
    print("-")
    @test eq.positions == pos
    print("-")
    @test eq.H == 2*pos[1:2, :]
    print("-")
    @test eq.mean_H == Statistics.mean(eq.H, dims=2)
    pos = randn(3, 10)
    ens = UnweightedParticleEnsemble(pos)
    update!(eq, ens)
    print("-")
    @test eq.positions == pos
    print("-")
    @test eq.H == 2*pos[1:2, :]
    print("-")
    @test eq.mean_H == Statistics.mean(eq.H, dims=2)
    println("DONE")

end; #PoissonEquation.jl

max_stars_repo_name: simsurace/FeedbackParticleFilters

max_issues_repo_name: simsurace/FeedbackParticleFilters

max_forks_repo_name: simsurace/FeedbackParticleFilters


text:
# run bandit algorithms on piecewise stationary bandit instance (and possibly save results)
# the choice of parameters given in this example correspond to the Bayesian experiment (average of many random problem instances) studied in Table 1 of [Besson et al. 20] with alpha0=0.5

using HDF5
using Statistics
using PyPlot

include("Algos.jl")
include("Benchmarks.jl")

# SAVING RESULTS
Save = "on" # should we save results ?

# TYPE OF EXPERIMENT
Type = "Bayesian" # Frequentist or Bayesian experiment?
# NUMBER OF REPETITIONS
N = 100
# TIME HORIZON
T = 20000

# NAME YOUR INSTANCE
# default choices : PB1, PB2, PB3, PB4 [see Besson et al. 2020], Random
# or use an other name for Bayesian experiments or for a custom problem
PB = "Bayesian"

# PARAMETERS OF THE RANDOM INSTANCES SAMPLER
K = 5
nbBreaks = 6
p=0.5 # probability that a breakpoint is a CP for each arm
Deltamin = 0.05 # minimal magnitude of a change
Deltamax = 0.3 # maximal magnitude of a change
Spacemin=1000 # minimal space between two breakpoints

# CUSTOM NON-STATIONARY PROBLEM (only necessary if you didn't choose one of the default names)
# (see different ways to specify them in Benchmarks.jl)
Problem,BreakPoints = ProblemUnif(T,MeansMatrix1)

if (PB == "PB1")
	# Problem 1 in [Besson et al. 2020]
	Problem,BreakPoints = ProblemUnif(T,MeansMatrix1)
elseif (PB == "PB2")
	# Problem 2 in [Besson et al. 2020]
	Problem, BreakPoints = ProblemUnif(T,MeansMatrix2)
elseif (PB =="PB3")
	# Problem 3 in [Besson et al. 2020]
	Problem = h5read("benchmarks/WorseBecomeGood","Problem")
	BreakPoints = h5read("benchmarks/WorseBecomeGood","BreakPoints")
elseif (PB =="PB4")
	# Problem 4 in [Besson et al. 2020]
	Problem = h5read("benchmarks/OneOptChange","Problem")
	BreakPoints = h5read("benchmarks/OneOptChange","BreakPoints")
elseif (PB[1:6] == "Random")
	# select a problem at random from the above sampler
	Problem, BreakPoints = RandomProblem(K,T,nbBreaks,p,Deltamin,Deltamax,Spacemin)
end

# display the arms means
if Type=="Frequentist"
	nbBreaks=length(BreakPoints)
	(K,T)=size(Problem)
	clf()
	PlotProblem(Problem)
end

## parameters for the CPD-based algorithms

# TWEAK
tweak = 0.5 # multiplicative constant in front of the exploration parameter alpha for all CDP-based algorithms

alpha0 = tweak*sqrt(K*log(T)/T)
alpha_opt = alpha0*sqrt(nbBreaks) # for algorithms using the knwoledge of nbBreaks
delta_GLR = 1/sqrt(T) # GLR-klUCB
M_CUSUM = 200 # CUMSUM (two variants tried in the paper)
epsilon_CUSUM = 0.1
M_CUSUM2 = 400
epsilon_CUSUM2 = Deltamin
h_CUSUM = log(T/nbBreaks)
w_MUCB = 200 # MUCB (two variants tried in the paper)
b_MUCB = sqrt(w_MUCB*log(2*K*T^2)/2)
w_MUCB2 = 800
b_MUCB2 = sqrt(w_MUCB2*log(2*K*T^2)/2)

# print a summary of the tuning
print("Exploration parameters are alpha0=$(alpha0), alpha_opt=$(alpha_opt)\n")
print("Parameters for GLR-klUCB: delta = $(delta_GLR)\n")
print("Parameters for CUSUM-klUCB: M = $(M_CUSUM), epsilon = $(epsilon_CUSUM), h = $(h_CUSUM)\n")
print("Parameters for M-klUCB: windows w = $(w_MUCB), threshold b = $(b_MUCB)\n")
print("Parameters for CUSUM-klUCB 2: M = $(M_CUSUM2), epsilon = $(epsilon_CUSUM2), h = $(h_CUSUM)\n")
print("Parameters for M-klUCB 2: windows w = $(w_MUCB2), threshold b = $(b_MUCB2)\n")

# name of the expes
series = "$(PB)_tweak_$(tweak)_"

## defining algorithms

# GLR-klUCB
beta(n,delta)=log(n^(3/2)/delta)
# Bernoulli GLR
GLRKLUCB(Table) = GLRklUCBGlobal(Table,beta,kl,alpha0,delta_GLR)
GLRKLUCBLoc(Table) = GLRklUCBLocal(Table,beta,kl,alpha0,delta_GLR)
# Gaussian GLR
GLRKLUCBGauss(Table) = GLRklUCBGlobal(Table,beta,klSG,alpha0,delta_GLR)
GLRKLUCBLocGauss(Table) = GLRklUCBLocal(Table,beta,klSG,alpha0,delta_GLR)
# Bernoulli GLR with constant, optimized, exploration
GLRKLUCBConstant(Table) = GLRklUCBGlobal(Table,beta,kl,alpha_opt,delta_GLR,"constant")
GLRKLUCBLocConstant(Table) = GLRklUCBLocal(Table,beta,kl,alpha_opt,delta_GLR,"constant")

# Oracle tuned as a function of the positions of breakpoints
Oracle(Table) = OracleklUCB(Table,BreakPoints)

# Other algorithms
SWKLUCB(Table) = SWklUCB(Table,nbBreaks)
DKLUCB(Table) = DklUCB(Table,nbBreaks)
EXP3SOpt(Table) = EXP3S(Table,nbBreaks)
MKLUCB(Table) = MklUCB(Table,alpha_opt,w_MUCB,b_MUCB)
MKLUCB2(Table) = MklUCB(Table,alpha_opt,w_MUCB2,b_MUCB2)
CUSUMKLUCB(Table) = CUSUMklUCB(Table,alpha_opt,M_CUSUM,epsilon_CUSUM,h_CUSUM)
CUSUMKLUCB2(Table) = CUSUMklUCB(Table,alpha_opt,M_CUSUM2,epsilon_CUSUM2,h_CUSUM)

print("Parameters for SW-klUCB: windows size tau = $(floor(Int,2*sqrt(T*log(T)/nbBreaks)))\n")
print("Parameters for D-klUCB: discount gamma = $(1- 0.25*sqrt(nbBreaks/T))\n\n")

## SELECT THE POLICIES TO RUN
#policies=[Oracle,klUCB,EXP3SOpt,SWKLUCB,DKLUCB,MKLUCB,MKLUCB2,CUSUMKLUCB,CUSUMKLUCB2,GLRKLUCB,GLRKLUCBConstant,GLRKLUCBLoc,GLRKLUCBLocConstant]
#policies=[Oracle,klUCB,EXP3SOpt,SWKLUCB,DKLUCB,MKLUCB2,CUSUMKLUCB2,GLRKLUCB,GLRKLUCBLoc]
policies=[klUCB,EXP3SOpt,SWKLUCB,DKLUCB,MKLUCB2,CUSUMKLUCB2,GLRKLUCB,GLRKLUCBLoc]

function MultiExpes(Problem,N,policies,series="test",display="on",saveRes="off")
	# run N repetitions of the chosen policies on the bandit with means given in Problem
	# optional: regret plot / saving the data
	(K,T)=size(Problem)
	lP = length(policies)
	# sub-sample the time
	tsave=floor.(Int,collect(range(1,T,length=500)))'
	ts = length(tsave)
	# draw N Tables of reward ahead of time
	Tables = [DrawTable(Problem) for n in 1:N]
	if display=="on"
		figure()
	end
	for imeth in 1:lP
		starttime = time()
		policy = policies[imeth]
		name = "results/$(series)_$(policy)_T_$(T)_N_$(N)"
		# storing regret and restarts
		Regret = zeros(N,ts)
		Restarts = zeros(N)
		for n in 1:N
			ChosenArms,ReceivedRewards,episode,ChangePoints= policy(Tables[n])
			reg = ComputeCumRegret(Problem,ChosenArms)
			Regret[n,:]=vec(reg[tsave])
			Restarts[n]=length(ChangePoints)
		end
		print("Results for $(policy)\n")
		print("Mean final regret is $(mean(Regret[:,end]))\n")
		print("Standard deviation of the final regret is $(std(Regret[:,end]))\n")
		print("Mean number of restarts is $(mean(Restarts))\n")
		print("Elapsed time is $(time()-starttime)\n\n")
		if (display == "on")
			# plot regret
			plot(tsave,mean(Regret,dims=1)',label="$(policy)")
		end
		if (saveRes == "on")
			# save date
			h5write(name,"Problem",Problem)
        	h5write(name,"Regret",Regret)
			h5write(name,"Restarts",Restarts)
		end
	end
	if (display == "on")
		legend()
	end
	return
end

function BayesianExpes(N,policies,K,T,nbBreaks,p,Deltamin,Deltamax,Spacemin,series="test",saveres="on")
	# run each policy on N random instances drawn from the sampler
	# optional : save the results
	lP = length(policies)
	# sub-sample the time
	tsave=floor.(Int,collect(range(1,T,length=500)))'
	ts = length(tsave)
	# draw the N Tables for N different random instances ahead of time
	Problems = [RandomProblem(K,T,nbBreaks,p,Deltamin,Deltamax,Spacemin)[1] for n in 1:N]
	Tables = [DrawTable(Pb) for Pb in Problems]
	for imeth in 1:lP
		starttime = time()
		policy = policies[imeth]
		name = "results/$(series)_$(policy)_T_$(T)_N_$(N)_K_$(K)_U_$(nbBreaks)"
		Regret=zeros(N,ts)
		Restarts=zeros(N)
		for n in 1:N
			Problem = Problems[n]
			ChosenArms,ReceivedRewards,episode,cpd= policy(Tables[n])
			reg = ComputeCumRegret(Problem,ChosenArms)
			Regret[n,:] = reg[tsave]
			Restarts[n] = length(cpd)
		end
		print("Results for $(policy)\n")
		print("Mean final regret is $(mean(Regret[:,end]))\n")
		print("Mean number of restarts is $(mean(Restarts))\n")
		print("Elapsed time is $(time()-starttime)\n\n")
		# save data
		if (saveres=="on")
	    	h5write(name,"Regret",Regret)
			h5write(name,"Restarts",Restarts)
			h5write(name,"Deltamin",Deltamin)
			h5write(name,"Deltamax",Deltamax)
			h5write(name,"Spacemin",Spacemin)
			h5write(name,"probaCP",p)
		end
	end
	return
end

start=time()

if (Type=="Frequentist")
	if Save == "on"
		# saving the data (no visualization)
		MultiExpes(Problem,N,policies,series,"off","on")
	else
		# run expes and display a regret plot
		MultiExpes(Problem,N,policies,series,"on","off")
	end
else
	if Save == "on"
		BayesianExpes(N,policies,K,T,nbBreaks,p,Deltamin,Deltamax,Spacemin,series)
	else
		BayesianExpes(N,policies,K,T,nbBreaks,p,Deltamin,Deltamax,Spacemin,series,"off")
	end
end

print("total elapsed time is $(time()-start)")


text:
@recipe function f(network::T) where {T<:AbstractEcologicalNetwork}
    if plotattributes[:seriestype] == :heatmap
        adjacency(network)
    end
end

@recipe function f(
    layout::Dict{K,NodePosition},
    network::T;
    nodesize=nothing,
    nodefill=nothing,
    bipartite=false,
) where {T<:AbstractEcologicalNetwork} where {K}

    # Node positions
    X = [layout[s].x for s in species(network)]
    Y = [layout[s].y for s in species(network)]

    # Default values
    framestyle --> :none
    legend --> false

    if typeof(network) <: QuantitativeNetwork
        int_range = extrema(network.edges.nzval)
    end

    if get(plotattributes, :seriestype, :plot) == :plot
        for interaction in network
            y = [layout[interaction.from].y, layout[interaction.to].y]
            x = [layout[interaction.from].x, layout[interaction.to].x]
            @series begin
                seriestype := :line
                linecolor --> :darkgrey
                if typeof(network) <: QuantitativeNetwork
                    linewidth --> EcologicalNetworksPlots._scale_value(
                        interaction.strength, int_range, (0.5, 3.5)
                    )
                end
                if typeof(network) <: ProbabilisticNetwork
                    seriesalpha --> interaction.probability
                end
                x, y
            end
        end
    end

    if get(plotattributes, :seriestype, :plot) == :scatter
        @series begin
            if nodesize !== nothing
                nsi_range = (minimum(values(nodesize)), maximum(values(nodesize)))
                markersize := [
                    EcologicalNetworksPlots._scale_value(nodesize[s], nsi_range, (2, 8)) for
                    s in species(network)
                ]
            end

            if nodefill !== nothing
                nfi_range = (minimum(values(nodefill)), maximum(values(nodefill)))
                marker_z := [
                    EcologicalNetworksPlots._scale_value(nodefill[s], nfi_range, (0, 1)) for
                    s in species(network)
                ]
            end

            if bipartite
                m_shape = Symbol[]
                for (i, s) in enumerate(species(network))
                    this_mshape = s ∈ species(network; dims=1) ? :dtriangle : :circle
                    push!(m_shape, this_mshape)
                end
                marker := m_shape
            end

            seriestype := :scatter
            seriescolor --> :white
            X, Y
        end
    end
end


text:
CD = ComparisonDataset

@testset "Rect" begin
  @testset "Basic" begin
    a = CD.Rect(10, 20, 30, 40)
    @test CD.top(a) == 10
    @test CD.left(a) == 20
    @test CD.bottom(a) == 30
    @test CD.right(a) == 40
    @test CD.width(a) == 20
    @test CD.height(a) == 20
  end

  @testset "Reordering" begin
    a = CD.Rect(30, 40, 10, 20)
    @test CD.top(a) == 10
    @test CD.left(a) == 20
    @test CD.bottom(a) == 30
    @test CD.right(a) == 40
    @test CD.width(a) == 20
    @test CD.height(a) == 20
  end
end
# Test overlapping
o1 = CD.Rect(0,0,10,10)
o2 = CD.Rect(-1,-1,4,4)
o3 = CD.Rect(5,-1,15,4)
o4 = CD.Rect(-1,8,4,15)
o5 = CD.Rect(6,8,15,15)

@test CD.isoverlapping(o1,o2) == true
@test CD.isoverlapping(o2,o1) == true
@test CD.isoverlapping(o1,o3) == true
@test CD.isoverlapping(o3,o1) == true
@test CD.isoverlapping(o1,o4) == true
@test CD.isoverlapping(o4,o1) == true
@test CD.isoverlapping(o1,o5) == true
@test CD.isoverlapping(o5,o1) == true
@test CD.isoverlapping(o2,o3) == false
@test CD.isoverlapping(o2,o4) == false
@test CD.isoverlapping(o2,o5) == false
@test CD.isoverlapping(o3,o4) == false
@test CD.isoverlapping(o3,o5) == false
@test CD.isoverlapping(o4,o5) == false

# Test random Rect generation
r = CD.randrect(20, 20, 10, 15)
@test CD.width(r) == 10
@test CD.height(r) == 15

rects = (CD.randrect(20, 20, 10, 15) for i in 1:100)
borders = (CD.left(r) > 0 &&
           CD.right(r) <= 20 &&
           CD.top(r) > 0 &&
           CD.bottom(r) <= 20 for r in rects)
@test false ∉ borders

# Test random non-overlapping Rect generation
rects = CD.nonoverlappingrandrect(100,100,10,10,10)
@test length(rects) == 10
borders = (CD.left(r) > 0 &&
           CD.right(r) <= 100 &&
           CD.top(r) > 0 &&
           CD.bottom(r) <= 100 for r in rects)
@test false ∉ borders
# Ensure none of the rects are overlapping
overlapping = false
for (a,b) in ((a,b) for a in rects, b in rects)
  if a ≠ b && CD.isoverlapping(a,b)
    overlapping = true
    break
  end
end
@test overlapping == false
# Test for exhaustion during search
@test_throws CD.ExhaustionError CD.nonoverlappingrandrect(20,20,10,10,100)


text:
module MiCRM

    using Reexport
    @reexport using ModelingToolkit
    @reexport using DifferentialEquations

    using Distributions, IfElse

    #community generation
    include("crm_parameters.jl")

    #ModelingToolkit wrappers
    include("MTK/crm_systems.jl")

    #Differential equations wrappers
    include("DiffEq/crm_problem.jl")

    include("./coalescence.jl")

    export random_micrm_params, micrm_system, micrm_params

end


text:
@doc raw"""
    SeriesAccelerations

This modules provides an algorithm for calculating infinite series efficiently.
Currently, two algorithms are provided. `accel_cohen_villegas_zagier()` for an
alternating series, and `accel_wynn_eps()` for more general series.
"""
module SeriesAccelerations

export accel_cohen_villegas_zagier
export accel_wynn_eps

################# Cohen, Villegas, Zagier 2000 algorithm #####################

@doc raw"""
    accel_cohen_villegas_zagier(ak)

Estimate the sum over the array `ak`. It is assumed without check that the
elements in the array have alternating sign and the series converges.

This algorithm is
[presented](https://people.mpim-bonn.mpg.de/zagier/files/exp-math-9/fulltext.pdf)
in [Cohen et al 2000]
(https://www.tandfonline.com/doi/abs/10.1080/10586458.2000.10504632).
"""
function accel_cohen_villegas_zagier(ak)
    T = eltype(ak)
    n = length(ak)
    d = (3 + √T(8))^n
    d = (d + 1/d) / 2
    b = -T(1)
    c = -d
    s = T(0)
    for k=0:n-1
        c = b - c
        s += (-1)^k * c * ak[k+1]
        b = (k + n) * (k - n) * b / ((k + 1/T(2)) * (k + 1))
    end
    return s / d
end

################ Wynn 1956 epsilon-algorithm #################################

function next_j!(c2, ϵ, n)
    c1 = (c2 & 1) + 1  # c1 is the other column (row, actually)
    for i=1:n
        ϵ[c1,i] = ϵ[c1,i+1] + 1 / (ϵ[c2,i+1] - ϵ[c2,i])
    end
    return c1
end

function accel_wynn_eps!(ϵ)
    N = size(ϵ, 2)
    c = 2
    for j=1:N-1
        c = next_j!(c, ϵ, N - j)
    end
    return ϵ[c,1]
end

@doc raw"""
    accel_wynn_eps(ak)

Estimate the sum over the array `ak`. It is assumed without check that the
series converges.

The algorithm used is [Wynn's
epsilon-algorithm](https://www.ams.org/journals/mcom/1956-10-054/S0025-5718-1956-0084056-6/S0025-5718-1956-0084056-6.pdf)
presented in [Wynn 1956](https://www.jstor.org/stable/2002183).
"""
function accel_wynn_eps(ak)
    ΔN = 1 - (length(ak) & 1)  # 1 for even length, 0 for odd
    N = length(ak) - ΔN
    ϵ = Array{eltype(ak)}(undef, 2, N)
    for n=1:N
        ϵ[1,n] = zero(eltype(ak))
        ϵ[2,n] = sum(ak[1:n+ΔN])
    end
    return accel_wynn_eps!(ϵ)
end

end

# vim: set sw=4 et sts=4 :


text:
function Update(P, N)

    if P == 0
        Res = N
    elseif N == 0
        Res = P
    else
        Res = [P; N]
    end

    return Res
end


text:
import Pkg

Pkg.activate(Pkg.dir("SearchBall"))
Pkg.add(Pkg.PackageSpec(url="https://github.com/control13/Chipmunk.jl"))
println("Chipmunk wurde hinzugefügt")


text:
# ---
# title: 867. Transpose Matrix
# id: problem867
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Array
# link: <https://leetcode.com/problems/transpose-matrix/description/>
# hidden: true
# ---
#
# Given a matrix `A`, return the transpose of `A`.
#
# The transpose of a matrix is the matrix flipped over it's main diagonal,
# switching the row and column indices of the matrix.
#
#
# ![](https://assets.leetcode.com/uploads/2019/10/20/hint_transpose.png)
#
#
#
# **Example 1:**
#
#
#
#     Input: [[1,2,3],[4,5,6],[7,8,9]]
#     Output: [[1,4,7],[2,5,8],[3,6,9]]
#
#
# **Example 2:**
#
#
#
#     Input: [[1,2,3],[4,5,6]]
#     Output: [[1,4],[2,5],[3,6]]
#
#
#
#
# **Note:**
#
#   1. `1 <= A.length <= 1000`
#   2. `1 <= A[0].length <= 1000`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
# # Nearly Incompressible Hyperelasticity
#
# ![](quasi_incompressible_hyperelasticity.gif)
#-
#md # !!! tip
#md #     This example is also available as a Jupyter notebook:
#md #     [`quasi_incompressible_hyperelasticity.ipynb`](@__NBVIEWER_ROOT_URL__/examples/quasi_incompressible_hyperelasticity.ipynb)
#-
# ## Introduction
#
# In this example we study quasi- or nearly-incompressible hyperelasticity using the stable Taylor-Hood approximation. In spirit, this example is the nonlinear analogue of
# [`incompressible_elasticity`](@__NBVIEWER_ROOT_URL__/examples/incompressible_elasticity.ipynb) and the incompressible analogue of
# [`hyperelasticity`](@__NBVIEWER_ROOT_URL__/examples/hyperelasticity.ipynb). Much of the code therefore follows from the above two examples.
# The problem is formulated in the undeformed or reference configuration with the displacement `u` and pressure `p` being the unknown fields. We now briefly outline
# the formulation. Consider the standard hyperelasticity problem
#
# ```math
#   u = \argmin_{v\in\mathcal{K}}\Pi(v),\quad \text{where}\quad \Pi(v)  = \int_\Omega \Psi(v) \ \mathrm{d}\Omega\ .
# ```
# For clarity of presentation we ignore any non-zero surface tractions and body forces and instead consider only
# applied displacements (i.e. non-homogeneous dirichlet boundary conditions). Moreover we stick our attention to the standard Neo-Hookean stored energy density
#
# ```math
#     \Psi(u) = \frac{\mu}{2}\left(I_1 - 3 \right) - \mu \log(J) + \frac{\lambda}{2}\left( J - 1\right){}^2,
# ```
# where $I_1 = F:F = F_{ij}F_{ij}$ and J = det(F) denote the standard invariants of the deformation gradient tensor $F = \mathbb{I}+\nabla u$.
# The above problem is ill-posed in the limit of incompressibility (or near-incompressibility), namely when
# ```math
#     \lambda/\mu \rightarrow +\infty.
# ```
# In order to alleviate the problem, we consider the partial legendre transform of the strain energy density Ψ with respect to J = det(F), namely
# ```math
#   \widehat{\Psi}(u, p) = \sup_{J} \left[ p(J - 1) - \frac{\mu}{2}\left(I_1 - 3 \right) + \mu \log(J) - \frac{\lambda}{2}\left( J - 1\right){}^2 \right].
# ```
# The supremum, say $J^\star$, can be calculated in closed form by the first order optimailty condition $\partial\widehat{\Psi}/\partial J = 0$. This gives
# ```math
#   J^\star(p) = \frac{\lambda + p + \sqrt{(\lambda + p){}^2 + 4 \lambda \mu }}{(2 \lambda)}.
# ```
# Furthermore, taking the partial legendre transform of $\widehat{\Psi}$ once again, gives us back the original problem, i.e.
# ```math
#     \Psi(u) = \Psi^\star(u, p) = \sup_{p} \left[ p(J - 1) - p(J^\star - 1) + \frac{\mu}{2}\left(I_1 - 3 \right) - \mu \log(J^\star) + \frac{\lambda}{2}\left( J^\star - 1\right){}^2 \right].
# ```
# Therefore our original hyperelasticity problem can now be reformulated as
# ```math
#   \inf_{u\in\mathcal{K}}\sup_{p} \int_\Omega\Psi^{\star} (u, p) \ \mathrm{d}\Omega.
# ```
# The total (modified) energy $\Pi^\star$ can then be written as
# ```math
#   \Pi^\star(u, p) = \int_\Omega p (J - J^\star) \ \mathrm{d}\Omega + \int_\Omega \frac{\mu}{2} \left( I_1 - 3\right) \ \mathrm{d}\Omega - \int_\Omega \mu\log(J^\star)\ \mathrm{d}\Omega + \int_\Omega \frac{\lambda}{2}\left( J^\star - 1 \right){}^2\ \mathrm{d}\Omega
# ```
# The euler-lagrange equations corresponding to the above energy give us our governing PDEs in the weak form, namely
# ```math
#   \int_\Omega \frac{\partial\Psi^\star}{\partial F}:\delta F \ \mathrm{d}\Omega = 0
# ```
# and
# ```math
#   \int_\Omega \frac{\partial \Psi^\star}{\partial p}\delta p \ \mathrm{d}\Omega = 0,
# ```
# where δF = δ∇u = ∇(δu) and δu and δp denote arbitrary variations with respect to displacement and pressure (or the test functions). See the references
# below for a more detailed explanation of the above mathematical trick. Now, in order to apply Newton's method to the
# above problem, we further need to linearize the above equations and calculate the respective hessians (or tangents), namely, $\partial^2\Psi^\star/\partial F^2$, $\partial^2\Psi^\star/\partial p^2$ and $\partial^2\Psi^\star/\partial F\partial p$
# which, using `Tensors.jl`, can be determined conveniently using automatic differentiation (see the code below). Hence we only need to define the above potential.
# The remaineder of the example follows similarly.
# ## References
# 1. [A paradigm for higher-order polygonal elements in finite elasticity using a gradient correction scheme, CMAME 2016, 306, 216–251](http://pamies.cee.illinois.edu/Publications_files/CMAME_2016.pdf)
# 2. [Approximation of incompressible large deformation elastic problems: some unresolved issues, Computational Mechanics, 2013](https://link.springer.com/content/pdf/10.1007/s00466-013-0869-0.pdf)
#
# We now get to the actual code. First, we import the respective packages

using Ferrite, Tensors, ProgressMeter
using BlockArrays, SparseArrays, LinearAlgebra

# and the corresponding `struct` to store our material properties.

struct NeoHooke
    μ::Float64
    λ::Float64
end

# We then create a function to generate a simple test mesh on which to compute FE solution. We also mark the boundaries
# to later assign Dirichlet boundary conditions
function importTestGrid()
    grid = generate_grid(Tetrahedron, (5, 5, 5), zero(Vec{3}), ones(Vec{3}));
    addfaceset!(grid, "myBottom", x -> norm(x[2]) ≈ 0.0);
    addfaceset!(grid, "myBack", x -> norm(x[3]) ≈ 0.0);
    addfaceset!(grid, "myRight", x -> norm(x[1]) ≈ 1.0);
    addfaceset!(grid, "myLeft", x -> norm(x[1]) ≈ 0.0);
    return grid
end;

# The function to create corresponding cellvalues for the displacement field `u` and pressure `p`
# follows in a similar fashion from the `incompressible_elasticity` example
function create_values(interpolation_u, interpolation_p)
    ## quadrature rules
    qr      = QuadratureRule{3,RefTetrahedron}(4)
    face_qr = QuadratureRule{2,RefTetrahedron}(4)

    ## geometric interpolation
    interpolation_geom = Lagrange{3,RefTetrahedron,1}()

    ## cell and facevalues for u
    cellvalues_u = CellVectorValues(qr, interpolation_u, interpolation_geom)
    facevalues_u = FaceVectorValues(face_qr, interpolation_u, interpolation_geom)

    ## cellvalues for p
    cellvalues_p = CellScalarValues(qr, interpolation_p, interpolation_geom)

    return cellvalues_u, cellvalues_p, facevalues_u
end;

# We now create the function for Ψ*
function Ψ(F, p, mp::NeoHooke)
    μ = mp.μ
    λ = mp.λ
    Ic = tr(tdot(F))
    J = det(F)
    Js = (λ + p + sqrt((λ + p)^2. + 4. * λ * μ ))/(2. * λ)
    return p * (Js - J) + μ / 2 * (Ic - 3) - μ * log(Js) + λ / 2 * (Js - 1)^2
end;

# and it's derivatives (required in computing the jacobian and hessian respectively)
function constitutive_driver(F, p, mp::NeoHooke)
    ## Compute all derivatives in one function call
    ∂²Ψ∂F², ∂Ψ∂F = Tensors.hessian(y -> Ψ(y, p, mp), F, :all)
    ∂²Ψ∂p², ∂Ψ∂p = Tensors.hessian(y -> Ψ(F, y, mp), p, :all)
    ∂²Ψ∂F∂p = Tensors.gradient(q -> Tensors.gradient(y -> Ψ(y, q, mp), F), p)
    return ∂Ψ∂F, ∂²Ψ∂F², ∂Ψ∂p, ∂²Ψ∂p², ∂²Ψ∂F∂p
end;

# The functions to create the `DofHandler` and `ConstraintHandler` (to assign corresponding boundary conditions) follow
# likewise from the incompressible elasticity example, namely

function create_dofhandler(grid, ipu, ipp)
    dh = DofHandler(grid)
    push!(dh, :u, 3, ipu) # displacement dim = 3
    push!(dh, :p, 1, ipp) # pressure dim = 1
    close!(dh)
    return dh
end;

# We are simulating a uniaxial tensile loading of a unit cube. Hence we apply a displacement field (`:u`) in `x` direction on the right face.
# The left, bottom and back faces are fixed in the `x`, `y` and `z` components of the displacement so as to emulate the uniaxial nature
# of the loading.
function create_bc(dh)
    dbc = ConstraintHandler(dh)
    add!(dbc, Dirichlet(:u, getfaceset(dh.grid, "myLeft"), (x,t) -> zero(Vec{1}), [1]))
    add!(dbc, Dirichlet(:u, getfaceset(dh.grid, "myBottom"), (x,t) -> zero(Vec{1}), [2]))
    add!(dbc, Dirichlet(:u, getfaceset(dh.grid, "myBack"), (x,t) -> zero(Vec{1}), [3]))
    add!(dbc, Dirichlet(:u, getfaceset(dh.grid, "myRight"), (x,t) -> t*ones(Vec{1}), [1]))
    close!(dbc)
    Ferrite.update!(dbc, 0.0)
    return dbc
end;

# Also, since we are considering incompressible hyperelasticity, an interesting quantity that we can compute is the deformed volume of the solid.
# It is easy to show that this is equal to ∫J*dΩ where J=det(F). This can be done at the level of each element (cell)
function calculate_element_volume(cell, cellvalues_u, ue)
    reinit!(cellvalues_u, cell)
    evol::Float64=0.0;
    @inbounds for qp in 1:getnquadpoints(cellvalues_u)
        dΩ = getdetJdV(cellvalues_u, qp)
        ∇u = function_gradient(cellvalues_u, qp, ue)
        F = one(∇u) + ∇u
        J = det(F)
        evol += J * dΩ
    end
    return evol
end;

# and then assembled over all the cells (elements)
function calculate_volume_deformed_mesh(w, dh::DofHandler, cellvalues_u)
    evol::Float64 = 0.0;
    @inbounds for cell in CellIterator(dh)
        global_dofs = celldofs(cell)
        nu = getnbasefunctions(cellvalues_u)
        global_dofs_u = global_dofs[1:nu]
        ue = w[global_dofs_u]
        δevol = calculate_element_volume(cell, cellvalues_u, ue)
        evol += δevol;
    end
    return evol
end;

# The function to assemble the element stiffness matrix for each element in the mesh now has a block structure like in
# `incompressible_elasticity`.
function assemble_element!(Ke, fe, cell, cellvalues_u, cellvalues_p, mp, ue, pe)
    ## Reinitialize cell values, and reset output arrays
    ublock, pblock = 1, 2
    reinit!(cellvalues_u, cell)
    reinit!(cellvalues_p, cell)
    fill!(Ke, 0.0)
    fill!(fe, 0.0)

    n_basefuncs_u = getnbasefunctions(cellvalues_u)
    n_basefuncs_p = getnbasefunctions(cellvalues_p)

    @inbounds for qp in 1:getnquadpoints(cellvalues_u)
        dΩ = getdetJdV(cellvalues_u, qp)
        ## Compute deformation gradient F
        ∇u = function_gradient(cellvalues_u, qp, ue)
        p = function_value(cellvalues_p, qp, pe)
        F = one(∇u) + ∇u

        ## Compute first Piola-Kirchhoff stress and tangent modulus
        ∂Ψ∂F, ∂²Ψ∂F², ∂Ψ∂p, ∂²Ψ∂p², ∂²Ψ∂F∂p = constitutive_driver(F, p, mp)

        ## Loop over the `u`-test functions to calculate the `u`-`u` and `u`-`p` blocks
        for i in 1:n_basefuncs_u
            ## gradient of the test function
            ∇δui = shape_gradient(cellvalues_u, qp, i)
            ## Add contribution to the residual from this test function
            fe[BlockIndex((ublock), (i))] += ( ∇δui ⊡ ∂Ψ∂F) * dΩ

            ∇δui∂S∂F = ∇δui ⊡ ∂²Ψ∂F²
            for j in 1:n_basefuncs_u
                ∇δuj = shape_gradient(cellvalues_u, qp, j)

                ## Add contribution to the tangent
                Ke[BlockIndex((ublock, ublock), (i, j))] += ( ∇δui∂S∂F ⊡ ∇δuj ) * dΩ
            end
            ## Loop over the `p`-test functions
            for j in 1:n_basefuncs_p
                δp = shape_value(cellvalues_p, qp, j)
                ## Add contribution to the tangent
                Ke[BlockIndex((ublock, pblock), (i, j))] += ( ∂²Ψ∂F∂p ⊡ ∇δui ) * δp * dΩ
            end
        end
        ## Loop over the `p`-test functions to calculate the `p-`u` and `p`-`p` blocks
        for i in 1:n_basefuncs_p
            δp = shape_value(cellvalues_p, qp, i)
            fe[BlockIndex((pblock), (i))] += ( δp * ∂Ψ∂p) * dΩ

            for j in 1:n_basefuncs_u
                ∇δuj = shape_gradient(cellvalues_u, qp, j)
                Ke[BlockIndex((pblock, ublock), (i, j))] += ∇δuj ⊡ ∂²Ψ∂F∂p * δp * dΩ
            end
            for j in 1:n_basefuncs_p
                δp = shape_value(cellvalues_p, qp, j)
                Ke[BlockIndex((pblock, pblock), (i, j))] += δp * ∂²Ψ∂p² * δp * dΩ
            end
        end
    end
end;

# The only thing that changes in the assembly of the global stiffness matrix is slicing the corresponding element
# dofs for the displacement (see `global_dofsu`) and pressure (`global_dofsp`).
function assemble_global!(K::SparseMatrixCSC, f, cellvalues_u::CellVectorValues{dim},
                         cellvalues_p::CellScalarValues{dim}, dh::DofHandler, mp::NeoHooke, w) where {dim}
    nu = getnbasefunctions(cellvalues_u)
    np = getnbasefunctions(cellvalues_p)

    ## start_assemble resets K and f
    fe = PseudoBlockArray(zeros(nu + np), [nu, np]) # local force vector
    ke = PseudoBlockArray(zeros(nu + np, nu + np), [nu, np], [nu, np]) # local stiffness matrix

    assembler = start_assemble(K, f)
    ## Loop over all cells in the grid
    for cell in CellIterator(dh)
        global_dofs = celldofs(cell)
        global_dofsu = global_dofs[1:nu]; # first nu dofs are displacement
        global_dofsp = global_dofs[nu + 1:end]; # last np dofs are pressure
        @assert size(global_dofs, 1) == nu + np # sanity check
        ue = w[global_dofsu] # displacement dofs for the current cell
        pe = w[global_dofsp] # pressure dofs for the current cell
        assemble_element!(ke, fe, cell, cellvalues_u, cellvalues_p, mp, ue, pe)
        assemble!(assembler, global_dofs, fe, ke)
    end
end;

# We now define a main function `solve`. For nonlinear quasistatic problems we often like to parameterize the
# solution in terms of a pseudo time like parameter, which in this case is used to gradually apply the boundary
# displacement on the right face. Also for definitenessm we consider λ/μ = 10⁴
function solve(interpolation_u, interpolation_p)

    ## import the mesh
    grid = importTestGrid()

    ## Material parameters
    μ = 1.
    λ = 1.E4 * μ
    mp = NeoHooke(μ, λ)

    ## Create the DofHandler and CellValues
    dh = create_dofhandler(grid, interpolation_u, interpolation_p)
    cellvalues_u, cellvalues_p, facevalues_u = create_values(interpolation_u, interpolation_p)

    ## Create the DirichletBCs
    dbc = create_bc(dh)

    ## Pre-allocation of vectors for the solution and Newton increments
    _ndofs = ndofs(dh)
    w  = zeros(_ndofs)
    ΔΔw = zeros(_ndofs)
    apply!(w, dbc)

    ## Create the sparse matrix and residual vector
    K = create_sparsity_pattern(dh)
    f = zeros(_ndofs)

    ## We run the simulation parameterized by a time like parameter. `Tf` denotes the final value
    ## of this parameter, and Δt denotes its increment in each step
    Tf = 2.0;
    Δt = 0.1;
    NEWTON_TOL = 1e-8

    pvd = paraview_collection("hyperelasticity_incomp_mixed.pvd");
    for t ∈ 0.0:Δt:Tf
        ## Perform Newton iterations
        Ferrite.update!(dbc, t)
        apply!(w, dbc)
        newton_itr = -1
        prog = ProgressMeter.ProgressThresh(NEWTON_TOL, "Solving @ time $t of $Tf;")
        fill!(ΔΔw, 0.0);
        while true; newton_itr += 1
            assemble_global!(K, f, cellvalues_u, cellvalues_p, dh, mp, w)
            norm_res = norm(f[Ferrite.free_dofs(dbc)])
            apply_zero!(K, f, dbc)
            ## Only display output at specific load steps
            if t%(5*Δt) == 0
                ProgressMeter.update!(prog, norm_res; showvalues = [(:iter, newton_itr)])
            end
            if norm_res < NEWTON_TOL
                break
            elseif newton_itr > 30
                error("Reached maximum Newton iterations, aborting")
            end
            ## Compute the incremental `dof`-vector (both displacement and pressure)
            ΔΔw .= K\f;

            apply_zero!(ΔΔw, dbc)
            w .-= ΔΔw
        end;

        ## Save the solution fields
        vtk_grid("hyperelasticity_incomp_mixed_$t.vtu", dh, compress=false) do vtkfile
            vtk_point_data(vtkfile, dh, w)
            vtk_save(vtkfile)
            pvd[t] = vtkfile
        end
    end;
    vtk_save(pvd);
    vol_def = calculate_volume_deformed_mesh(w, dh, cellvalues_u);
    print("Deformed volume is $vol_def")
    return vol_def;
end;

# We can now test the solution using the Taylor-Hood approximation
quadratic = Lagrange{3, RefTetrahedron, 2}()
linear = Lagrange{3, RefTetrahedron, 1}()
vol_def = solve(quadratic, linear);

# We can also check that the deformed volume is indeed close to 1 (as should be for a nearly incompressible material)
using Test                #src
@test isapprox(vol_def, 1.0, atol=1E-3) #src

#md # ## Plain Program
#md #
#md # Below follows a version of the program without any comments.
#md # The file is also available here: [`quasi_incompressible_hyperelasticity.jl`](quasi_incompressible_hyperelasticity.jl)
#md #
#md # ```julia
#md # @__CODE__
#md # ```


text:
parseresponse(response) = map(
    x -> parse(Float64, split(x, ' ')[2]),
    split(String(response.body), '\n')[4:end-1],
)

bodykey(::Type{PhotoelectricAbsorption}) = "photoelectric"
bodykey(::Type{Coherent}) = "coherent"
bodykey(::Type{Incoherent}) = "incoherent"
bodykey(::Type{InNuclearField}) = "nuclear"
bodykey(::Type{InElectronField}) = "electron"
bodykey(::Type{WithCoherent}) = "with"
bodykey(::Type{WithoutCoherent}) = "without"

formatenergies(energies::AbstractArray{<:Unitful.Energy}) =
    join(["$(Float64(uconvert(MeV, e).val))" for e in energies], ';')

const basebody = Dict{String,String}(
    "NumAdd" => "1",
    "OutOpt" => "PIC",
    "Output" => "",
    "character" => "space",
    "WindowXmin" => "0.001",
    "WindowXmax" => "100000",
)

const XCOM_URL = "https://www.physics.nist.gov/cgi-bin/Xcom/data.pl"

function XCOM(body)
    merge!(body, basebody)
    r = HTTP.request(
        "POST",
        XCOM_URL,
        [],
        join(["$k=$v" for (k, v) in body], '&'),
    )
    if r.status != 200
        desc = STATUS_CODES[r.status]
        error("XCOM error - try again")
    end
    parseresponse(r)
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMPlates.jl/blob/master/LICENSE

"""
    DKT - Discrete Kirchhoff Triangle element

# References

- Batoz, Jean‐Louis, Klaus‐JÜRgen Bathe, and Lee‐Wing Ho.
  "A study of three‐node triangular plate bending elements."
  International Journal for Numerical Methods in Engineering 15.12 (1980): 1771-1812.
- Lucena Neto, Eliseu, et al.
  "An Explicit Consistent Geometric Stiffness Matrix for the DKT Element."
  Latin American Journal of Solids and Structures 14.4 (2017): 613-628.
"""
struct DKT <: FieldProblem
    geometric_stiffness :: Bool
end

function DKT()
    geometric_stiffness = true
    return DKT(geometric_stiffness)
end

function FEMBase.get_unknown_field_name(::Problem{DKT})
    return "displacement"
end

"""
        basis_factory(::Problem{DKT}, element::Element{Tri3})

Return basis functions and their partial derivatives for DKT.
"""
function basis_factory(::Problem{DKT}, element::Element{Tri3}, time)
    X = element("geometry", time)
    x1, x2, x3 = (X[i][1] for i in 1:3)
    y1, y2, y3 = (X[i][2] for i in 1:3)

    # for ij = 23, 31, 12:
    # x_ij = x_i - x_j
    # y_ij = y_i - y_j
    # L_ij^2 = x_ij^2 - y_ij^2
    x23 = x2-x3
    x31 = x3-x1
    x12 = x1-x2
    y23 = y2-y3
    y31 = y3-y1
    y12 = y1-y2
    L23 = (x23^2 + y23^2)
    L31 = (x31^2 + y31^2)
    L12 = (x12^2 + y12^2)

    # for k = 4,5,6 and ij=23,31,12:

    # a_k = -x_ij/L_ij^2
    # b_k = 3*x_ij*y_ij / (4*L_ij^2)
    # c_k = (x_ij^2 - 2*y_ij^2) / 4*L_ij^2
    # d_k = -y_ij / L_ij^2
    # e_k = (y_ij^2 - 2*x_ij^2) / 4*L_ij^2
    a4 = -x23/L23
    a5 = -x31/L31
    a6 = -x12/L12
    b4 = 3*x23*y23/(4*L23)
    b5 = 3*x31*y31/(4*L31)
    b6 = 3*x12*y12/(4*L12)
    c4 = (x23^2 - 2*y23^2)/(4*L23)
    c5 = (x31^2 - 2*y31^2)/(4*L31)
    c6 = (x12^2 - 2*y12^2)/(4*L12)
    d4 = -y23/L23
    d5 = -y31/L31
    d6 = -y12/L12
    e4 = (y23^2 - 2*x23^2)/(4*L23)
    e5 = (y31^2 - 2*x31^2)/(4*L31)
    e6 = (y12^2 - 2*x12^2)/(4*L12)

    # p_k = -6*x_ij/l_ij^2
    # q_k = 3*x_ij*y_ij/L_ij^2
    # r_k = 3*y_ij^2/L_ij^2
    # t_k = -6*y_ij/L_ij^2
    p4 = -6*x23/L23
    p5 = -6*x31/L31
    p6 = -6*x12/L12
    q4 =  3*x23*y23/L23
    q5 =  3*x31*y31/L31
    q6 =  3*x12*y12/L12
    r4 =  3*y23^2/L23
    r5 =  3*y31^2/L31
    r6 =  3*y12^2/L12
    t4 = -6*y23/L23
    t5 = -6*y31/L31
    t6 = -6*y12/L12

    # println("(a4, a5, a6) = ($a4, $a5, $a6)")
    # println("(b4, b5, b6) = ($b4, $b5, $b6)")
    # println("(c4, c5, c6) = ($c4, $c5, $c6)")
    # println("(d4, d5, d6) = ($d4, $d5, $d6)")
    # println("(e4, e5, e6) = ($e4, $e5, $e6)")
    # println("(p4, p5, p6) = ($p4, $p5, $p6)")
    # println("(q4, q5, q6) = ($q4, $q5, $q6)")
    # println("(r4, r5, r6) = ($r4, $r5, $r6)")
    # println("(t4, t5, t6) = ($t4, $t5, $t6)")

    function N(ip)
        xi, eta = ip
        N1 = 2.0*(1.0-xi-eta)*(0.5-xi-eta)
        N2 = xi*(2.0*xi-1.0)
        N3 = eta*(2.0*eta-1.0)
        N4 = 4.0*xi*eta
        N5 = 4.0*eta*(1.0-xi-eta)
        N6 = 4.0*xi*(1.0-xi-eta)
        return (N1, N2, N3, N4, N5, N6)
    end

    function Hx(xi)
        N1, N2, N3, N4, N5, N6 = N(xi)
        Hx1 = 1.5*(a6*N6 - a5*N5)
        Hx2 = b5*N5 + b6*N6
        Hx3 = N1 - c5*N5 - c6*N6
        Hx4 = 1.5*(a4*N4 - a6*N6)
        Hx5 = b4*N4 + b6*N6
        Hx6 = N2 - c4*N4 - c6*N6
        Hx7 = 1.5*(a5*N5 - a4*N4)
        Hx8 = b4*N4 + b5*N5
        Hx9 = N3 - c4*N4 - c5*N5
        return [Hx1 Hx2 Hx3 Hx4 Hx5 Hx6 Hx7 Hx8 Hx9]
    end

    function Hy(xi)
        N1, N2, N3, N4, N5, N6 = N(xi)
        Hy1 = 1.5*(d6*N6 - d5*N5)
        Hy2 = -N1 + e5*N5 + e6*N6
        Hy3 = -b5*N5 - b6*N6
        Hy4 = 1.5*(d4*N4 - d6*N6)
        Hy5 = -N2 + e4*N4 + e6*N6
        Hy6 = -b4*N4 - b6*N6
        Hy7 = 1.5*(d5*N5 - d4*N4)
        Hy8 = -N3 + e4*N4 + e5*N5
        Hy9 = -b4*N4 - b5*N5
        return [Hy1 Hy2 Hy3 Hy4 Hy5 Hy6 Hy7 Hy8 Hy9]
    end

    function dHxdxi(ip)
        xi, eta = ip
        d1 =  p6*(1.0-2.0*xi) + (p5-p6)*eta
        d2 =  q6*(1.0-2.0*xi) - (q5+q6)*eta
        d3 = -4.0 + 6.0*(xi+eta) + r6*(1.0-2.0*xi) - eta*(r5+r6)
        d4 = -p6*(1.0-2.0*xi) + eta*(p4+p6)
        d5 =  q6*(1.0-2.0*xi) - eta*(q6-q4)
        d6 = -2.0 + 6.0*xi + r6*(1.0-2.0*xi) + eta*(r4-r6)
        d7 = -eta*(p5+p4)
        d8 =  eta*(q4-q5)
        d9 = -eta*(r5-r4)
        return [d1 d2 d3 d4 d5 d6 d7 d8 d9]
    end

    function dHydxi(ip)
        xi, eta = ip
        d1 =  t6*(1.0-2.0*xi) + eta*(t5-t6)
        d2 =  1.0 + r6*(1.0-2.0*xi) - eta*(r5+r6)
        d3 = -q6*(1.0-2.0*xi) + eta*(q5+q6)
        d4 = -t6*(1.0-2.0*xi) + eta*(t4+t6)
        d5 = -1.0 + r6*(1.0-2.0*xi) + eta*(r4-r6)
        d6 = -q6*(1.0-2*xi) - eta*(q4-q6)
        d7 = -eta*(t4+t5)
        d8 =  eta*(r4-r5)
        d9 = -eta*(q4-q5)
        return [d1 d2 d3 d4 d5 d6 d7 d8 d9]
    end

    function dHxdeta(ip)
        xi, eta = ip
        d1 = -p5*(1.0-2.0*eta) - xi*(p6-p5)
        d2 =  q5*(1.0-2.0*eta) - xi*(q5+q6)
        d3 = -4.0 + 6.0*(xi+eta) + r5*(1.0-2.0*eta) - xi*(r5+r6)
        d4 =  xi*(p4+p6)
        d5 =  xi*(q4-q6)
        d6 = -xi*(r6-r4)
        d7 =  p5*(1.0-2.0*eta) - xi*(p4+p5)
        d8 =  q5*(1.0-2.0*eta) + xi*(q4-q5)
        d9 = -2.0 + 6.0*eta + r5*(1.0-2.0*eta+xi*(r4-r5))
        return [d1 d2 d3 d4 d5 d6 d7 d8 d9]
    end

    function dHydeta(ip)
        xi, eta = ip
        d1 = -t5*(1.0-2.0*eta) - xi*(t6-t5)
        d2 =  1.0 + r5*(1.0-2.0*eta) - xi*(r5+r6)
        d3 = -q5*(1.0-2.0*eta) + xi*(q5+q6)
        d4 =  xi*(t4+t6)
        d5 =  xi*(r4-r6)
        d6 = -xi*(q4-q6)
        d7 =  t5*(1.0-2.0*eta) - xi*(t4+t5)
        d8 = -1.0 + r5*(1.0-2.0*eta) + xi*(r4-r5)
        d9 = -q5*(1.0-2.0*eta) - xi*(q4-q5)
        return [d1 d2 d3 d4 d5 d6 d7 d8 d9]
    end

    return Hx, Hy, dHxdxi, dHxdeta, dHydxi, dHydeta

end

function FEMBase.assemble_elements!(problem::Problem{DKT},
                                    assembly::Assembly,
                                    elements::Vector{Element{Tri3}},
                                    time::Float64)

    ndofs = 3*3
    Ke = zeros(ndofs, ndofs)
    fe = zeros(ndofs)

    for element in elements

        fill!(Ke, 0.0)
        fill!(fe, 0.0)

        X = element("geometry", time)
        x1, x2, x3 = (X[i][1] for i in 1:3)
        y1, y2, y3 = (X[i][2] for i in 1:3)
        x23 = x2-x3
        x31 = x3-x1
        x12 = x1-x2
        y23 = y2-y3
        y31 = y3-y1
        y12 = y1-y2
        Hx, Hy, dHxdxi, dHxdeta, dHydxi, dHydeta = basis_factory(problem, element, time)

        for ip in get_integration_points(element)

            detJ = element(ip, time, Val{:detJ})
            E = element("youngs modulus", ip, time)
            t = element("thickness", ip, time)
            nu = element("poissons ratio", ip, time)

            D = E*t^3/(12*(1-nu^2)) * [
                1.0  nu  0.0
                nu  1.0  0.0
                0.0 0.0  0.5*(1-nu)]

            B = 1.0/(x31*y12 - x12*y31) * [
                 y31*dHxdxi(ip) + y12*dHxdeta(ip)
                -x31*dHydeta(ip) - x12*dHydeta(ip)
                -x31*dHxdxi(ip) - x12*dHxdeta(ip) + y31*dHydxi(ip) + y12*dHydeta(ip)]

            Ke += ip.weight * B'*D*B * detJ

            # if problem.properties.geometric_stiffness
            #     # Assemble Kg
            # end

            if haskey(element, "distributed load")
                N = element(ip, time)
                p = element("distributed load", ip, time)
                fe[1:3:end] += ip.weight * p * N' * detJ
            end

        end

        gdofs = get_gdofs(problem, element)
        add!(assembly.K, gdofs, gdofs, Ke)
        add!(assembly.f, gdofs, fe)

    end

    return nothing

end


text:
function DoubleIntegrator(;N=21)
    opts = SolverOptions(
        penalty_scaling=1000.,
        penalty_initial=1.,
    )

    model = RobotZoo.DoubleIntegrator()
    n,m = RD.dims(model)

    # Task
    x0 = @SVector [0., 0.]
    xf = @SVector [1., 0]
    tf = 2.0

    # Discretization info
    dt = tf/(N-1)

    # Costs
    Q = 1.0*Diagonal(@SVector ones(n))
    Qf = 1.0*Diagonal(@SVector ones(n))
    R = 1.0e-1*Diagonal(@SVector ones(m))
    obj = LQRObjective(Q*dt,R*dt,Qf,xf,N)

    # Constraints
    u_bnd = 3.0
    x_bnd = [Inf,0.6]
    conSet = ConstraintList(n,m,N)
    bnd = BoundConstraint(n,m, u_min=-u_bnd, u_max=u_bnd, x_min=-x_bnd, x_max=x_bnd)
    goal = GoalConstraint(xf)
    add_constraint!(conSet, bnd, 1:N-1)
    add_constraint!(conSet, goal, N:N)

    doubleintegrator_static = Problem(model, obj, x0, tf, constraints=conSet, xf=xf)
    TO.rollout!(doubleintegrator_static)
    return doubleintegrator_static, opts
end


text:
### Simulate polyploid cells with XA repression hypothesis in stochastic Gillespie model
include("sim_functions_gil_Oct2018.jl")

println(ARGS[1])
println(ARGS[2])

pars=convert(Array{Float64},readcsv(ARGS[1]))

# Simulate all parameters from stochastic simulation for initial cond XaXa
nr_cells=100
#Simulation time
t = 100.0

ma_scan=Array{Float64}(size(pars,1)*nr_cells,33+4*(Int64(t)+1))
@time for n=1:size(pars,1)
  println(n)
  p=pars[n,:]
  println(p)
  ma_scan[(n-1)*nr_cells+1:n*nr_cells,1:33] = p'.*ones(nr_cells,33)
  ma_scan[(n-1)*nr_cells+1:n*nr_cells,34:33+4*(Int64(t)+1)] = sim_ba_gil_polypl_repXA(p,nr_cells,t)
end

writecsv(ARGS[2],ma_scan)