text:
"""
    MatrixUtilityModule

Module for general utility matrix product functions.
"""
module MatrixUtilityModule

__precompile__(true)

using DelimitedFiles
using SparseArrays
using ..FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
using LoopVectorization

"""
    loc!(loc::FFltMat, ecoords::FFltMat, N::FFltMat)

Compute the location of the quadrature point.

Arguments:
`loc` = matrix of coordinates, overwritten  inside the function
`ecoords` = matrix of the node coordinates for the element.
`N` = matrix of basis function values
"""
function loc!(loc::FFltMat, ecoords::FFltMat, N::FFltMat)
	return mulCAtB!(loc, N, ecoords)
end

"""
    jac!(J::FFltMat, ecoords::FFltMat,gradNparams::FFltMat)

Compute the Jacobian matrix at the quadrature point.

Arguments:
`J` = Jacobian matrix, overwritten  inside the function
`ecoords` = matrix of the node coordinates for the element.
`gradNparams` = matrix of basis function gradients
"""
function jac!(J::FFltMat, ecoords::FFltMat, gradNparams::FFltMat)
    return mulCAtB!(J, ecoords, gradNparams)
end

"""
    locjac!(loc::FFltMat, J::FFltMat, ecoords::FFltMat, N::FFltMat, gradNparams::FFltMat)

Compute location and Jacobian matrix at the quadrature point.

Arguments:
`loc` = matrix of coordinates, overwritten  inside the function
`J` = Jacobian matrix, overwritten  inside the function
`ecoords` = matrix of the node coordinates for the element.
`N` = matrix of basis function values
`gradNparams` = matrix of basis function gradients
"""
function locjac!(loc::FFltMat, J::FFltMat, ecoords::FFltMat, N::FFltMat, gradNparams::FFltMat)
    return loc!(loc, ecoords, N), jac!(J, ecoords, gradNparams)
end

"""
    add_mggt_ut_only!(Ke::FFltMat, gradN::FFltMat, mult::FFlt)

Add the product `gradN*mult*gradNT` to the elementwise matrix `Ke`.
The argument `mult` is a scalar.
*Only upper triangle* is computed; the lower triangle is not touched.
(Use `complete_lt!` to complete the lower triangle, if needed.)

The matrix `Ke` is assumed to be suitably initialized.

The matrix `Ke` is modified.  The matrix `gradN` is not modified
inside this function.
"""
function add_mggt_ut_only!(Ke::FFltMat, gradN::FFltMat, mult::FFlt)
    Kedim = size(Ke, 1)
    @assert Kedim == size(Ke, 2) # Square matrix?
    nne, mdim = size(gradN)
    @assert nne == Kedim # compatible matrices?
    @inbounds for nx = 1:Kedim # Do: Ce  =  Ce + gradN*((Jac*w[j]))*gradN' ;
        @inbounds for px = 1:mdim
            a = (mult)*gradN[nx, px]
            @inbounds for mx = 1:nx # only the upper triangle
                Ke[mx, nx] +=  gradN[mx, px] * a
            end
        end
    end
    return true
end

"""
    add_gkgt_ut_only!(Ke::FFltMat, gradN::FFltMat, Jac_w::FFlt,
      kappa_bar::FFltMat, kappa_bargradNT::FFltMat)

Add the product `gradN*kappa_bar*gradNT*(Jac*w[j])` to the elementwise matrix `Ke`.
*Only upper triangle* is computed; the lower triangle is not touched.
(Use `complete_lt!` to complete the lower triangle, if needed.)

The matrix `Ke` is assumed to be suitably initialized.

Upon return,  the matrix `Ke` is updated.  The scratch buffer `kappa_bargradNT`
is overwritten during each call of this function. The matrices `gradN` and
`kappa_bar` are not modified inside this function.
"""
function add_gkgt_ut_only!(Ke::FFltMat, gradN::FFltMat, Jac_w::FFlt,
        kappa_bar::FFltMat, kappa_bargradNT::FFltMat)
    @assert size(Ke, 1) == size(Ke, 2)
    Kedim = size(Ke, 1)
    nne, mdim = size(gradN)
    @assert size(kappa_bar) == (mdim, mdim)
    @assert size(kappa_bargradNT) == (mdim, nne)
    # A_mul_Bt!(kappa_bargradNT, kappa_bar, gradN); # intermediate result
    @inbounds for nx = 1:nne
        @inbounds for mx = 1:mdim
	        accum::FFlt  = 0.0
	        @inbounds @simd for px = 1:mdim
	            accum += kappa_bar[mx, px]*gradN[nx, px]
	        end
	        kappa_bargradNT[mx, nx] = Jac_w*accum
        end
    end
    # Ke = Ke + gradN*(kappa_bar*(Jac*w[j]))*gradN'; only upper triangle
    @inbounds for nx = 1:Kedim
        @inbounds for mx = 1:nx # only the upper triangle
            accum::FFlt  = 0.0
            @inbounds @simd for px = 1:mdim
                accum += gradN[mx, px]*kappa_bargradNT[px, nx]
            end
            Ke[mx, nx] += accum
        end
    end
    return true
end

"""
    complete_lt!(Ke::FFltMat)

Complete the lower triangle of the elementwise matrix `Ke`.

The matrix `Ke` is modified  inside this function. The
upper-triangle  entries  are copied  across the diagonal
to the lower triangle.
"""
function complete_lt!(Ke::FFltMat)
    Kedim = size(Ke, 1)
    @assert Kedim == size(Ke, 2)
    @inbounds for nx = 1:Kedim # complete the lower triangle
        @inbounds for mx = nx+1:Kedim
            Ke[mx, nx] = Ke[nx, mx]
        end
    end
    return true
end

"""
    add_btdb_ut_only!(Ke::FFltMat, B::FFltMat, Jac_w::FFlt,
                  D::FFltMat, DB::FFltMat)

Add the product  `(B'*(D*(Jac*w[j]))*B)`, to the elementwise matrix Ke.
*Only upper triangle* is computed; the lower triangle is not touched.
(Use `complete_lt!` to complete the lower triangle, if needed.)

The matrix Ke is assumed to be suitably initialized.

The matrix Ke is modified.  The matrices B and D are not modified
inside this function. The scratch buffer DB is overwritten
during each call of this function.
"""
function add_btdb_ut_only!(Ke::FFltMat, B::FFltMat, Jac_w::FFlt, D::FFltMat, DB::FFltMat)
    @assert size(Ke, 1) == size(Ke, 2)
    @assert size(B, 1) == size(D, 1)
    nstr, Kedim = size(B)
    @assert size(D) == (nstr, nstr)
    @assert size(DB) == (nstr, Kedim)
    # A_mul_B!(DB, D, B)    # intermediate product
    @inbounds for nx = 1:Kedim
        @inbounds for mx = 1:nstr
            accum::FFlt  = 0.0
            @inbounds for px = 1:nstr
                accum += D[mx, px]*B[px, nx]
            end
            DB[mx, nx] = Jac_w*accum
        end
    end
    #  Ke = Ke + (B'*(D*(Jac*w[j]))*B); only the upper triangle
    @inbounds for nx = 1:Kedim
        @inbounds for mx = 1:nx # only the upper triangle
	        accum::FFlt  = 0.0
	        @inbounds for px = 1:nstr
	            accum += B[px, mx]*DB[px, nx]
	        end
	        Ke[mx, nx] += accum
        end
    end
    return true
end

"""
    add_btsigma!(Fe::FFltVec, B::FFltMat, coefficient::FFlt, sigma::FFltVec)

Add the product  `B'*(sigma*coefficient)`, to the elementwise vector `Fe`.

The vector `Fe` is assumed to be suitably initialized.

The vector `Fe` is modified.  The vector `sigma` is not modified
inside this function.
"""
function add_btsigma!(Fe::FFltVec, B::FFltMat, coefficient::FFlt, sigma::FFltVec)
    @assert size(B, 1) == length(sigma)
    nstr, Kedim = size(B)
    @inbounds for nx = 1:Kedim
	    accum::FFlt  = 0.0
	    @inbounds for px = 1:nstr
	    	accum += sigma[px]*B[px, nx]
		end
		Fe[nx] += coefficient*accum
	end
    return true
end

"""
    add_nnt_ut_only!(Ke::FFltMat, Nn::FFltMat, Jac_w_coeff::FFlt)

Add the product  `Nn*(Nn'*(coeff*(Jac*w(j)))`, to the elementwise matrix `Ke`.
*Only upper triangle* is computed; the lower triangle is not touched.

The matrix `Ke` is assumed to be suitably initialized. The
matrix `Nn` has a single column.

The matrix `Ke` is modified.  The matrix `Nn` is not modified
inside this function.
"""
function add_nnt_ut_only!(Ke::FMat{T}, N::FFltMat, Jac_w_coeff::T) where {T<:Number}
    Kedim = size(N, 1)
    @assert Kedim == size(Ke, 2)
    @inbounds for nx = 1:Kedim
        a = (Jac_w_coeff)*N[nx]
        @inbounds for mx = 1:nx # only the upper triangle
            Ke[mx,nx] += N[mx] * a
        end
    end
    return true
end

function add_nnt_ut_only!(Ke::FFltMat, N::FFltVec, Jac_w_coeff::FFlt)
    return add_nnt_ut_only!(Ke, reshape(N, length(N), 1), Jac_w_coeff)
end

"""
    adjugate3!(B, A)

Compute the adjugate matrix of `A`.
"""
function adjugate3!(B, A)
    # % A = rand(3)
    # % det(A)*inv(A)-B
    # B =0*A;
    B[1,1] =(A[5]*A[9]-A[8]*A[6]);
    B[1,2] =-(A[4]*A[9]-A[7]*A[6]);
    B[1,3] =(A[4]*A[8]-A[7]*A[5]);

    B[2,1] =-(A[2]*A[9]-A[8]*A[3]);
    B[2,2] =(A[1]*A[9]-A[7]*A[3]);
    B[2,3] =-(A[1]*A[8]-A[7]*A[2]);

    B[3,1] =(A[2]*A[6]-A[5]*A[3]);
    B[3,2] =-(A[1]*A[6]-A[4]*A[3]);
    B[3,3] =(A[1]*A[5]-A[4]*A[2]);
    return B
end

"""
    symmetrize!(a)

Make the matrix on input symmetric.

The operation is in-place.
"""
function symmetrize!(a)
	@assert size(a, 1) == size(a, 2)
    @inbounds for c in 1:size(a, 2)
    	@inbounds for r in c:size(a, 1)
    		a[c, r] =  a[r, c] += a[c, r]
    	end
    end
    a .= 0.5 * a
end

"""
    mulCAB!(::Val{3}, C, A, B)

Compute the product of 3X3 matrices `C = A * B`

"""
function mulCAB!(::Val{3}, C, A, B)
	#define MAT3DUTILMULTAB3X3(C, A, B)
	C[1, 1]  = (A[1, 1]*B[1, 1] + A[1, 2]*B[2, 1] + A[1, 3]*B[3, 1]);
	C[1, 2]  = (A[1, 1]*B[1, 2] + A[1, 2]*B[2, 2] + A[1, 3]*B[3, 2]);
	C[1, 3]  = (A[1, 1]*B[1, 3] + A[1, 2]*B[2, 3] + A[1, 3]*B[3, 3]);
	C[2, 1]  = (A[2, 1]*B[1, 1] + A[2, 2]*B[2, 1] + A[2, 3]*B[3, 1]);
	C[2, 2]  = (A[2, 1]*B[1, 2] + A[2, 2]*B[2, 2] + A[2, 3]*B[3, 2]);
	C[2, 3]  = (A[2, 1]*B[1, 3] + A[2, 2]*B[2, 3] + A[2, 3]*B[3, 3]);
	C[3, 1]  = (A[3, 1]*B[1, 1] + A[3, 2]*B[2, 1] + A[3, 3]*B[3, 1]);
	C[3, 2]  = (A[3, 1]*B[1, 2] + A[3, 2]*B[2, 2] + A[3, 3]*B[3, 2]);
	C[3, 3]  = (A[3, 1]*B[1, 3] + A[3, 2]*B[2, 3] + A[3, 3]*B[3, 3]);
	return C
end

"""
    mulCAtB!(::Val{3}, C, A, B)

Compute the product of 3X3 matrices `C = Transpose(A) * B`

"""
function mulCAtB!(::Val{3}, C, A, B)
	#define MAT3DUTILMULTATB3X3(C, A, B)
	C[1, 1]  = (A[1, 1]*B[1, 1] + A[2, 1]*B[2, 1] + A[3, 1]*B[3, 1]);
	C[1, 2]  = (A[1, 1]*B[1, 2] + A[2, 1]*B[2, 2] + A[3, 1]*B[3, 2]);
	C[1, 3]  = (A[1, 1]*B[1, 3] + A[2, 1]*B[2, 3] + A[3, 1]*B[3, 3]);
	C[2, 1]  = (A[1, 2]*B[1, 1] + A[2, 2]*B[2, 1] + A[3, 2]*B[3, 1]);
	C[2, 2]  = (A[1, 2]*B[1, 2] + A[2, 2]*B[2, 2] + A[3, 2]*B[3, 2]);
	C[2, 3]  = (A[1, 2]*B[1, 3] + A[2, 2]*B[2, 3] + A[3, 2]*B[3, 3]);
	C[3, 1]  = (A[1, 3]*B[1, 1] + A[2, 3]*B[2, 1] + A[3, 3]*B[3, 1]);
	C[3, 2]  = (A[1, 3]*B[1, 2] + A[2, 3]*B[2, 2] + A[3, 3]*B[3, 2]);
	C[3, 3]  = (A[1, 3]*B[1, 3] + A[2, 3]*B[2, 3] + A[3, 3]*B[3, 3]);
	return C
end

"""
    mulCABt!(::Val{3}, C, A, B)

Compute the product of 3X3 matrices `C = A * Transpose(B)`

"""
function mulCABt!(::Val{3}, C, A, B)
	#define MAT3DUTILMULTABT3X3(C, A, B)
	C[1, 1]  = (A[1, 1]*B[1, 1] + A[1, 2]*B[1, 2] + A[1, 3]*B[1, 3]);
	C[1, 2]  = (A[1, 1]*B[2, 1] + A[1, 2]*B[2, 2] + A[1, 3]*B[2, 3]);
	C[1, 3]  = (A[1, 1]*B[3, 1] + A[1, 2]*B[3, 2] + A[1, 3]*B[3, 3]);
	C[2, 1]  = (A[2, 1]*B[1, 1] + A[2, 2]*B[1, 2] + A[2, 3]*B[1, 3]);
	C[2, 2]  = (A[2, 1]*B[2, 1] + A[2, 2]*B[2, 2] + A[2, 3]*B[2, 3]);
	C[2, 3]  = (A[2, 1]*B[3, 1] + A[2, 2]*B[3, 2] + A[2, 3]*B[3, 3]);
	C[3, 1]  = (A[3, 1]*B[1, 1] + A[3, 2]*B[1, 2] + A[3, 3]*B[1, 3]);
	C[3, 2]  = (A[3, 1]*B[2, 1] + A[3, 2]*B[2, 2] + A[3, 3]*B[2, 3]);
	C[3, 3]  = (A[3, 1]*B[3, 1] + A[3, 2]*B[3, 2] + A[3, 3]*B[3, 3]);
	return C
end

"""
    mulCAtB!(C, A, B)

Compute the matrix `C = A' * B`

The use of BLAS is purposefully avoided in order to eliminate contentions of
multi-threaded execution of the library code with the user-level threads.
"""
function mulCAtB!(C, A, B)
    M, N = size(C); K = size(B,1)
    @assert size(C, 1) == size(A, 2)
    @assert size(C, 2) == size(B, 2)
    @assert size(A, 1) == size(B, 1)
    # When the @avx macro is available, this code is faster:
    # z = zero(eltype(C))
    @avx for n in 1:size(C,2), m in 1:size(C,1)
        Cmn = zero(eltype(C)) # Cmn = z # For some reason this bombed in LoopVectorization 0.4.0
        for k in 1:size(A,1)
            Cmn += A[k,m] * B[k,n]
        end
        C[m,n] = Cmn
    end
    # Otherwise, this code works just fine:
 	# if mod(M, 2) == 0 && mod(N, 2) == 0
	#     @inbounds for m in 1:2:M
	#     	m1 = m + 1
	#     	@inbounds for n in 1:2:N
	#     		n1 = n + 1
	# 	    	C11, C21, C12, C22 = 0.0, 0.0, 0.0, 0.0
	# 	    	@inbounds for k in 1:K
	# 	    		C11 += A[k,m] * B[k,n]
	# 	    		C21 += A[k,m1] * B[k,n]
	# 	    		C12 += A[k,m] * B[k,n1]
	# 	    		C22 += A[k,m1] * B[k,n1]
	# 	    	end
	# 	    	C[m,n] = C11
	# 	    	C[m1,n] = C21
	# 	    	C[m,n1] = C12
	# 	    	C[m1,n1] = C22
	# 	    end
	#     end
	# else
	# 	@inbounds for n in 1:N, m in 1:M
	#     	Cmn = 0.0
	#     	@inbounds for k in 1:K
	#     		Cmn += A[k,m] * B[k,n]
	#     	end
	#     	C[m,n] = Cmn
	#     end
	# end
    return C
end

"""
    mulCAB!(C, A, B)

Compute the matrix `C = A * B`

The use of BLAS is purposefully avoided in order to eliminate contentions of
multi-threaded execution of the library code with the user-level threads.

Note: See the thread https://discourse.julialang.org/t/ann-loopvectorization/32843/36
"""
function mulCAB!(C, A, B)
	M, N = size(C); K = size(B,1)
	@assert M == size(A, 1)
    @assert N == size(B, 2)
    @assert size(A, 2) == K
    # Surprisingly the code below (nkm) is slower than the one in the order mnk.
	# C .= 0
	# @avx  for n in 1:N, k in 1:K
	# 	Bkn = B[k,n]
	# 	for m in 1:M
	# 		C[m,n] += A[m,k] * Bkn
	# 	end
	# end
	# z = zero(eltype(C))
	@avx  for m in 1:M, n in 1:N
		Cmn = zero(eltype(C))
	    for k in 1:K
	        Cmn += A[m,k] * B[k,n]
	    end
	    C[m,n] = Cmn
	end
	return C
end

"""
    mulCABt!(C, A, B)

Compute the matrix `C = A * B'`

The use of BLAS is purposefully avoided in order to eliminate contentions of
multi-threaded execution of the library code with the user-level threads.
"""
function mulCABt!(C, A, B)
    M, N = size(C); K = size(B,2)
	@assert size(C, 1) == size(A, 1)
	@assert size(C, 2) == size(B, 1)
	@assert size(A, 2) == size(B, 2)
	C .= 0
   	@avx for n ∈ 1:N, k ∈ 1:K # for k ∈ 1:K, n ∈ 1:N #
   		Bnk = B[n,k]
   		for m ∈ 1:M
   			C[m,n] += A[m,k] * Bnk
   		end
   	end
    return C
end

"""
    detC(::Val{3}, C::FFltMat)

Compute determinant of 3X3 `C`.
"""
function detC(::Val{3}, C::FFltMat)
	#define MAT3DUTIL_DET_3X3(T3X3)
	return (
		(C[1, 1] * C[2, 2] * C[3, 3]) +
		(C[1, 2] * C[2, 3] * C[3, 1]) +
		(C[1, 3] * C[2, 1] * C[3, 2]) -
		(C[1, 3] * C[2, 2] * C[3, 1]) -
		(C[1, 2] * C[2, 1] * C[3, 3]) -
		(C[1, 1] * C[2, 3] * C[3, 2]))
end

"""
    mulCAB!(C::Vector{T}, A, B::Vector{T})  where {T}

Compute the product `C = A * B`, where `C` and `B` are "vectors".

The use of BLAS is purposefully avoided in order to eliminate contentions of
multi-threaded execution of the library code with the user-level threads.
"""
function mulCAB!(C::Vector{T}, A, B::Vector{T})  where {T}
	M, N = size(A);
	@assert length(C) == M
    @assert length(B) == N
	C .= 0.0
	@avx  for n in 1:N
		for m in 1:M
			C[m] += A[m,n] * B[n]
		end
	end
	return C
end

"""
    export_sparse(filnam, M)

Export sparse matrix to a text file.
"""
function export_sparse(filnam, M)
#readdlm(source, ',')
    open(filnam, "w") do io
        out = [(size(M)..., 0.0)]
        I, J, V = findnz(M)
        out = vcat(out, [(I[i], J[i], V[i]) for i in 1:length(I)])
        writedlm(io, out, ',')
    end;
end

"""
    import_sparse(filnam)

Import sparse matrix from a text file.
"""
function import_sparse(filnam)
    B = open(filnam, "r") do io
        rows = readlines(io)
        s = split(replace(rows[1], "," => " "))
        m, n = parse(Int64, s[1]), parse(Int64, s[2])
        I = fill(0, length(rows)-1)
        J = fill(0, length(rows)-1)
        V = fill(0.0, length(rows)-1)
        for i in 1:length(I)
            s = split(replace(rows[i+1], "," => " "))
            I[i] = parse(Int64, s[1])
            J[i] = parse(Int64, s[2])
            V[i] = parse(Float64, s[3])
        end
        B = sparse(I, J, V, m, n)
    end
    return B
end

end


text:
# # Bar and stack plots with [PowerGraphics.jl](github.com/nrel-siip/PowerGraphics.jl)
# PowerGraphics also provides some basic specifications for plotting `SimulationResults`.
# This example demonstrates some simple plotting capabilities using different Plots.julia
# backends.
#
# The plotting capabilities use the Julia Plots package which can generate plots using
# several different graphics packages. We'll use GR.jl and PlotlyJS.jl.
#
# ## Dependencies
using SIIPExamples #for path locations
pkgpath = dirname(dirname(pathof(SIIPExamples)))
using PowerSystems #to load results
using PowerSimulations #to load results
using PowerGraphics

# ### Results file
# If you have already run some of the other examples, you should have generated some results
# (If you haven't run some of the other simulaitons, you can run
# `include(joinpath(pkgpath, "test", "3_PowerSimulations_examples", "2_sequential_simulations.jl"))`).
# You can load the results into memory with:
simulation_folder = joinpath(pkgpath, "RTS-GMLC-master", "rts-test")
simulation_folder =
    joinpath(simulation_folder, "$(maximum(parse.(Int64,readdir(simulation_folder))))")
res = load_simulation_results(simulation_folder, "UC")

# ## Plots
# By default, PowerGraphics uses the GR graphics package as the backend for Plots.jl to
# generate figures. This creates static plots and should execute without any extra steps.
# For example, we can create a stacked bar_plot:
gr() # loads the GR backend
bar_plot(res)

# However, interactive plotting can generate much more insightful figures, especially when
# creating somewhat complex stacked figures. So, we can use the PlotlyJS backend for Plots,
# but it requires that PlotlyJS.jl, and ORCA.jl (if in a notebook, WebIO.jl is required too)
# are installed in your Project.toml. To startup the PlotlyJS backend, run:
plotlyjs()

# Now we can create stacked bar plots that can be inspected interactively.
bar_plot(res)

# Similarly, we can create a stack plot for any combination of variable to see the time
# series values. *Note: the `load = true` kwarg populates a line for the total system load.*

stack_plot(res, [Symbol("P__ThermalStandard"), Symbol("P__RenewableDispatch")], load = true)

# Or, we can create a series of stack plots for every variable in the dictionary:
# ```julia
# stack_plot(res)
# ```

# Generator fuel types are not stored in the model, or the associated results files. So,
# to make aggregated fuel plots, we need to load the `System` as well. The simulation routine
# automatically serializes the `System` data into the results directory, so we just need to
# load it.
uc_sys = System(joinpath(
    simulation_folder,
    "models_json",
    "stage_UC_model",
    "Stage1_sys_data.json",
))

# Now we can make a set of aggregated plots by fuel type.
fuel_plot(res, uc_sys, load = true)

ext: jl


text:
const LEFT_PAREN  = Val{:LEFT_PAREN}
const RIGHT_PAREN = Val{:RIGHT_PAREN}
const PLUS        = Val{:PLUS}
const MINUS       = Val{:MINUS}
const ASTERISK    = Val{:ASTERISK}
const SLASH       = Val{:SLASH}
const CARET       = Val{:CARET}
const LITERAL     = Val{:LITERAL}
const INT_VAL     = Val{:INT_VAL}

const TokenType = Union{
	Type{LEFT_PAREN},
	Type{RIGHT_PAREN},
	Type{PLUS},
	Type{MINUS},
	Type{ASTERISK},
	Type{SLASH},
	Type{CARET},
	Type{LITERAL},
	Type{INT_VAL}
}

@enum Precedences begin
	NONE        = 0
	SUM         = 10
	PRODUCT     = 20
	EXPONENT_R  = 29
	EXPONENT    = 30
	PREFIX      = 40
end

struct Token{TokenType}
	text::String
end

# left binding powers

lbp(token::Union{Token{PLUS}, Token{MINUS}})             = SUM::Precedences
lbp(token::Union{Token{ASTERISK}, Token{SLASH}})         = PRODUCT::Precedences
lbp(token::Token{CARET})                                 = EXPONENT::Precedences
lbp(token::Union{Token{LEFT_PAREN}, Token{RIGHT_PAREN}}) = NONE::Precedences

# null denotations

nud(token::Token,          itr) = error("unexpected token: ", token)
nud(token::Token{PLUS},    itr) = expression(itr, PREFIX::Precedences)
nud(token::Token{MINUS},   itr) = Expr(:call, :-, expression(itr, PREFIX::Precedences))
nud(token::Token{INT_VAL}, itr) = parse(Int, token.text)
nud(token::Token{LITERAL}, itr) = Symbol(token.text)

function nud(token::Token{LEFT_PAREN}, itr)
	expr = expression(itr, NONE::Precedences)
	lookahead = popfirst!(itr)

	if typeof(lookahead) != Token{RIGHT_PAREN}
		error(Token(RIGHT_PAREN, ")"), " expected, got: ", lookahead)
	end

	expr
end

# left denotations

led(token::Token,           left, itr) = error("unexpected token: ", token)
led(token::Token{PLUS},     left, itr) = Expr(:call, :+, left, expression(itr, SUM::Precedences))
led(token::Token{MINUS},    left, itr) = Expr(:call, :-, left, expression(itr, SUM::Precedences))
led(token::Token{ASTERISK}, left, itr) = Expr(:call, :*, left, expression(itr, PRODUCT::Precedences))
led(token::Token{SLASH},    left, itr) = Expr(:call, :/, left, expression(itr, PRODUCT::Precedences))
led(token::Token{CARET},    left, itr) = Expr(:call, :^, left, expression(itr, EXPONENT_R::Precedences))

# main parser

function expression(itr, rbp = NONE::Precedences)
	token = popfirst!(itr)
	left = nud(token, itr)

	while Base.peek(itr) != nothing && rbp < lbp(Base.peek(itr))
		token = popfirst!(itr)
		left = led(token, left, itr)
	end

	left
end

# dummy token stream

struct TokenStream
	token_stream::Vector{Token}
end

Base.iterate(stream::TokenStream, count = 1) = count > size(stream.token_stream, 1) ? nothing : (stream.token_stream[count], count + 1)

stream = TokenStream(
		[
			Token{LITERAL}("A"),
			Token{ASTERISK}("*"),
			Token{LEFT_PAREN}("("),
			Token{INT_VAL}("1"),
			Token{PLUS}("+"),
			Token{INT_VAL}("2"),
			Token{CARET}("^"),
			Token{INT_VAL}("3"),
			Token{SLASH}("/"),
			Token{INT_VAL}("3"),
			Token{RIGHT_PAREN}(")")
		]
	)

itr = Iterators.Stateful(stream)

expr = expression(itr)

println("result: ", expr)
dump(expr)


text:
# RUN: julia -e "import Brutus; Brutus.lit(:emit_lowered)" --startup-file=no %s 2>&1 | FileCheck %s

add(x, y) = x + y
emit(add, Int64, Int64)

# CHECK: Core.MethodMatch(Tuple{typeof(Main.Main.add), Int64, Int64}, svec(), add(x, y) in Main.Main at /{{.*}}/test/Codegen/lower/add_int.jl:3, true)after translating to MLIR in JLIR dialect:module  {
# CHECK-NEXT:   func nested @"Tuple{typeof(Main.add), Int64, Int64}"(%arg0: !jlir<"typeof(Main.add)">, %arg1: !jlir.Int64, %arg2: !jlir.Int64) -> !jlir.Int64 attributes {llvm.emit_c_interface} {
# CHECK-NEXT:     "jlir.goto"()[^bb1] : () -> ()
# CHECK-NEXT:   ^bb1:  // pred: ^bb0
# CHECK-NEXT:     %0 = "jlir.constant"() {value = #jlir<"#<intrinsic #2 add_int>">} : () -> !jlir<"typeof(Core.IntrinsicFunction)">
# CHECK-NEXT:     %1 = "jlir.call"(%0, %arg1, %arg2) : (!jlir<"typeof(Core.IntrinsicFunction)">, !jlir.Int64, !jlir.Int64) -> !jlir.Int64
# CHECK-NEXT:     "jlir.return"(%1) : (!jlir.Int64) -> ()
# CHECK-NEXT:   }
# CHECK-NEXT: }

# CHECK: module  {
# CHECK-NEXT:   func nested @"Tuple{typeof(Main.add), Int64, Int64}"(%arg0: !jlir<"typeof(Main.add)">, %arg1: i64, %arg2: i64) -> i64 attributes {llvm.emit_c_interface} {
# CHECK-NEXT:     %0 = addi %arg1, %arg2 : i64
# CHECK-NEXT:     return %0 : i64
# CHECK-NEXT:   }
# CHECK-NEXT: }

# CHECK:   llvm.func @"Tuple{typeof(Main.add), Int64, Int64}"(%arg0: !llvm.ptr<struct<"struct_jl_value_type", opaque>>, %arg1: i64, %arg2: i64) -> i64 attributes {llvm.emit_c_interface, sym_visibility = "nested"} {
# CHECK-NEXT:     %0 = llvm.add %arg1, %arg2  : i64
# CHECK-NEXT:     llvm.return %0 : i64
# CHECK-NEXT:   }
# CHECK-NEXT:   llvm.func @"_mlir_ciface_Tuple{typeof(Main.add), Int64, Int64}"(%arg0: !llvm.ptr<struct<"struct_jl_value_type", opaque>>, %arg1: i64, %arg2: i64) -> i64 attributes {llvm.emit_c_interface, sym_visibility = "nested"} {
# CHECK-NEXT:     %0 = llvm.call @"Tuple{typeof(Main.add), Int64, Int64}"(%arg0, %arg1, %arg2) : (!llvm.ptr<struct<"struct_jl_value_type", opaque>>, i64, i64) -> i64
# CHECK-NEXT:     llvm.return %0 : i64
# CHECK-NEXT:   }
# CHECK-NEXT: }

max_stars_repo_name: stephenneuendorffer/brutus

max_issues_repo_name: stephenneuendorffer/brutus

max_forks_repo_name: stephenneuendorffer/brutus


text:
facts("Leadfield") do

    L = 1

    context("Create") do

        L = Leadfield(rand(2000, 3, 6), vec(rand(2000, 1)), vec(rand(2000, 1)), vec(rand(2000, 1)), vec(["Cz","80Hz_SWN_70dB_R","20Hz_SWN_70dB_R","40Hz_SWN_70dB_R","10Hz_SWN_70dB_R","_4Hz_SWN_70dB_R"]))

        @fact size(L.L) --> (2000, 3, 6)
    end

    context("Show") do

        show(L)
    end

    context("Match") do

        s = read_SSR(joinpath(dirname(@__FILE__), "../../data", "test_Hz19.5-testing.bdf"))

        L1 = match_leadfield(L, s)
        @fact size(L1.L) --> (2000, 3, 6)

        keep_channel!(s, ["Cz","80Hz_SWN_70dB_R","20Hz_SWN_70dB_R","40Hz_SWN_70dB_R"])
        L2 = match_leadfield(deepcopy(L), s)
        @fact size(L2.L) --> (2000, 3, 4)
        @fact L2.L[:, :, 1] --> L1.L[:, :, 1]
        @fact L2.L[:, :, 2] --> L1.L[:, :, 3]
        @fact L2.L[:, :, 4] --> L1.L[:, :, 6]
        @fact L2.sensors --> channelnames(s)

        s = merge_channels(s, "Cz", "garbage")
        @fact_throws BoundsError match_leadfield(L, s)

        L = Leadfield(rand(2000, 3, 5), vec(rand(2000, 1)), vec(rand(2000, 1)), vec(rand(2000, 1)), vec(["Cz","80Hz_SWN_70dB_R","10Hz_SWN_70dB_R","_4Hz_SWN_70dB_R"]))
        @fact_throws ErrorException  match_leadfield(L, s)

    end

    context("Operations") do
        ldf = Leadfield(rand(5, 3, 2), collect(1:5.0), collect(1:5.0), collect(1:5.0), ["ch1", "ch2"])

        context("Find location") do
            for n = 1:size(ldf.L, 1)
                @fact n --> find_location(ldf, Talairach(ldf.x[n], ldf.y[n], ldf.z[n]))
            end

            for n = 1:size(ldf.L, 1)
                @fact n --> find_location(ldf, Talairach(ldf.x[n]+0.1, ldf.y[n]-0.1, ldf.z[n]))
            end
        end

    end
end


text:
## Basic functions

RealVector = Vector{<:Real}
RealMatrix = Matrix{<:Real}
# Base.@irrational ° 0.0174532925199432957692369076848861271344 (big(pi)/big(180))

function DeleteDuplicates(v::Array{T,1})::Array{T,1} where T <: Any
    w = Array{eltype(v),1}()
    if length(v) > 0
        for e ∈ v
            if e ∉ w
                push!(w,e)
            end
        end
    end
    return w
end

function OrthogonalVector(v::RealVector) where T <: Real
    if norm(v) == 0
        error("zero vector")
    else
        v₁ = normalize(v)
        i = argmin(abs.(v))
        v₂ = zeros(length(v))
        v₂[i] = 1.0
        v₂ = normalize(v₂-dot(v₁,v₂)*v₁)
        return v₂
    end
end

function isorthogonal(M::AbstractMatrix{T}) where T<:Real
    size(M)
    return M*M'≈oneunit(M)
end

function rotatematrix(v::RealVector,θ::Real)
    n = normalize(v)
    N = [0 -n[3] n[2];n[3] 0 -n[1];-n[2] n[1] 0]
    return oneunit(N)+N^2+(-cos(θ)*N^2+sin(θ)*N)
end

function Circumcenter(p₁::RealVector,p₂)
    c₁ = 1
    c₂ = 1
    center = (c₁*p₁+c₂*p₂)/(c₁+c₂)
    return center
end

function Circumcenter(p₁::RealVector,p₂::RealVector,p₃::RealVector)
    l₁² = (norm(p₂-p₃))^2
    l₂² = (norm(p₃-p₁))^2
    l₃² = (norm(p₁-p₂))^2
    c₁ = l₁²*(l₂²+l₃²-l₁²)
    c₂ = l₂²*(l₃²+l₁²-l₂²)
    c₃ = l₃²*(l₁²+l₂²-l₃²)
    center = (c₁*p₁+c₂*p₂+c₃*p₃)/(c₁+c₂+c₃)
    return center
end

function Circumcenter(p₁::RealVector,p₂::RealVector,p₃::RealVector,p₄::RealVector)
    n = 3
    𝒑 = hcat(p₁,p₂,p₃,p₄)
    p² = [dot(𝒑[1:n,i],𝒑[1:n,i]) for i in 1:(n+1)]
    P = vcat(𝒑,reshape([1. for i in 1:(n+1)],1,n+1))
    Q = inv(P)
    𝒄 = (Q'*p²)[1:n]/2
    r = sqrt(dot(𝒄,𝒄)+dot(p²,Q[:,end]))
    return 𝒄
end

function NormalVector(p₁::RealVector,p₂::RealVector,p₃::RealVector)
    return normalize(cross(p₁,p₂)+cross(p₂,p₃)+cross(p₃,p₁))
end

function Angles(p₁::RealVector,p₂::RealVector,p₃::RealVector)
    l₁ = norm(p₂-p₃)
    l₂ = norm(p₃-p₁)
    l₃ = norm(p₁-p₂)
    cosα₁ = dot(p₂-p₁,p₃-p₁)/(l₂*l₃)
    cosα₂ = dot(p₃-p₂,p₁-p₂)/(l₃*l₁)
    cosα₃ = dot(p₁-p₃,p₂-p₃)/(l₁*l₂)
    α₁ = acos(clamp(cosα₁,-1,1))
    α₂ = acos(clamp(cosα₂,-1,1))
    α₃ = acos(clamp(cosα₃,-1,1))
    return α₁,α₂,α₃
end

max_stars_repo_name: hyrodium/JuliRay

max_issues_repo_name: hyrodium/JuliRay

max_forks_repo_name: hyrodium/JuliRay


text:
using VUMPS
using KrylovKit
using CUDA
using Test
using OMEinsum
using Random
CUDA.allowscalar(false)

@testset "OMEinsum with $atype{$dtype} " for atype in [Array, CuArray], dtype in [Float64, ComplexF64]
    Random.seed!(100)
    d = 9
    D = 20
    AL = atype(rand(dtype, D, d, D))
    M = atype(rand(dtype, d, d, d, d))
    FL = atype(rand(dtype, D, d, D))
    @time ein"((γcη,ηpβ),csap),γsα -> αaβ"(FL,AL,M,conj(AL))
end

@testset "KrylovKit with $atype{$dtype}" for atype in [Array, CuArray], dtype in [Float64, ComplexF64]
    Random.seed!(100)
    d = 4
    D = 10
    FL = atype(rand(dtype, D, d, D))
    M = atype(rand(dtype, d, d, d, d))
    AL = atype(rand(dtype, D, d, D))
    @time λs, FLs, info = eigsolve(FL -> ein"((γcη,ηpβ),csap),γsα -> αaβ"(FL,AL,M,conj(AL)), FL, 1, :LM; ishermitian = false)

    λl,FL = real(λs[1]),real(FLs[1])
    dFL = atype(rand(dtype, D, d, D))
    @time ξl, info = linsolve(FR -> ein"((ηpβ,βaα),csap),γsα -> ηcγ"(AL, FR, M, conj(AL)), permutedims(dFL, (3, 2, 1)), -λl, 1)
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "Pythia8Builder"
version = v"0.1.0"

# Collection of sources required to build Pythia8Builder
sources = [
    "http://home.thep.lu.se/~torbjorn/pythia8/pythia8243.tgz" =>
    "f8ec27437d9c75302e192ab68929131a6fd642966fe66178dbe87da6da2b1c79",

]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/pythia8243/
./configure --prefix=$prefix --enable-shared
make -j $nproc
make install

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:x86_64, libc=:glibc),
    Linux(:x86_64, libc=:musl),
    MacOS(:x86_64)
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libpythia8", :libpythia8)
]

# Dependencies that must be installed before this package can be built
dependencies = [

]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
using Dice
using Dice: num_flips, num_nodes, to_dice_ir

#TODO: the code below is not correct
code = @dice begin
    # triangle distribution

    function triangle(b::Int)
        s = false
        m = 2^(2b)
        n = 2^b
        x = Vector(undef, b)
        y = Vector(undef, b)
        for i = b:-1:1
            x[i] = if s flip(1/2) else flip((m/2-n/2)/(m-n)) end
            y[i] = x[i] || flip(1/(3/2-n/m))
            s = s || (!x[i] && y[i])
            m = m/4
            n = n/2
        end
        return DistInt(x)
    end

    for b = 1:16
        t = triangle(b)
        println("$b-bit triangle uses $(num_flips(t)) flips and $(num_nodes(t)) BDD nodes")
    end
    triangle(2)

end
bdd = compile(code)
num_flips(bdd)
num_nodes(bdd)
infer(code, :bdd)


text:
abstract type NaturalEvolutionStrategyOpt <: PopulationOptimizer end

"""
Separable Natural Evolution Strategy (sNES) optimizer.
"""
mutable struct SeparableNESOpt{F,E<:EmbeddingOperator} <: NaturalEvolutionStrategyOpt
    embed::E                        # operator embedding into the search space
    lambda::Int                     # number of samples to take per iteration
    mu::Vector{Float64}             # center of the population
    sigma::Vector{Float64}          # average std deviation for sampling in position
    distr::Distribution             # distribution to sample the step sizes from
    mu_learnrate::Float64
    sigma_learnrate::Float64
    max_sigma::Float64

    last_s::Matrix{Float64}         # `s` values sampled in the last call to `ask()`
    candidates::Vector{Candidate{F}}# the last sampled values, now being evaluated
    sortedUtilities::Vector{Float64}# the fitness shaping utility vector
    tmp_Utilities::Vector{Float64}   # the fitness shaping utility vector sorted by current population fitness

    function SeparableNESOpt{F}(
            embed::E;
            lambda::Int = 0,
            mu_learnrate::Float64 = 1.0,
            sigma_learnrate::Float64 = 0.0,
            ini_x = nothing, max_sigma::Float64 = 1.0E+10
    ) where {F, E<:EmbeddingOperator}
        d = numdims(search_space(embed))

        if lambda == 0
            lambda = 4 + ceil(Int, log(3*d)) # default lambda
        end
        if sigma_learnrate == 0.0
            sigma_learnrate = calc_sigma_learnrate_for_snes(d) # default sigma learn rate
        end
        if ini_x === nothing
            ini_x = rand_individual(search_space(embed))
        else
            ini_x = copy(ini_x::Individual)
        end

        new{F,E}(embed, lambda,
                 ini_x, fill(1.0, d),
                 Normal(0, 1),
                 mu_learnrate, sigma_learnrate, max_sigma,
                 fill(0.0, d, lambda),
                 [Candidate{F}(fill!(Individual(undef, d), NaN), i) for i in 1:lambda],
                 # Most modern NES papers use log rather than linear fitness shaping.
                 fitness_shaping_utilities_log(lambda),
                 Vector{Float64}(undef, lambda))
    end
end

population(o::NaturalEvolutionStrategyOpt) = o.candidates
numdims(o::NaturalEvolutionStrategyOpt) = numdims(search_space(o.embed))

set_candidate!(o::SeparableNESOpt, x0) = (o.mu = x0)
candidate(o::SeparableNESOpt) = o.mu

const NES_DefaultOptions = ParamsDict(
    :lambda => 0,              # If 0 it will be set based on the number of dimensions
    :ini_x => nothing,         # starting point, "nothing" generates random point in a search space
    :mu_learnrate => 1.0,
    :sigma_learnrate => 0.0,   # If 0.0 it will be set based on the number of dimensions
    :max_sigma => 1.0E+10      # Maximal sigma
)

function separable_nes(problem::OptimizationProblem, parameters)
    params = chain(NES_DefaultOptions, parameters)
    embed = RandomBound(search_space(problem))
    SeparableNESOpt{fitness_type(problem)}(embed,
        lambda = params[:lambda],
        mu_learnrate = params[:mu_learnrate],
        sigma_learnrate = params[:sigma_learnrate],
        ini_x = params[:ini_x],
        max_sigma = params[:max_sigma])
end

calc_sigma_learnrate_for_snes(d) = (3 + log(d)) / (5 * sqrt(d))
calc_sigma_learnrate_for_nes(d) = (9 + 3 * log(d)) / (5 * d * sqrt(d))

function ask(snes::SeparableNESOpt)
    # Sample from N(0, 1)
    randn!(snes.last_s)

    for i in eachindex(snes.candidates)
        candi = snes.candidates[i]
        candi.index = i # reset ordering
        @inbounds for j in eachindex(candi.params)
            candi.params[j] = snes.mu[j] + snes.sigma[j] * snes.last_s[j, i]
        end
        apply!(snes.embed, candi.params, snes.mu)
        candi.fitness = NaN # FIXME: use nafitness()
    end

    return snes.candidates
end

function tell!(snes::SeparableNESOpt{F}, rankedCandidates::AbstractVector{<:Candidate{F}}) where F
    u = assign_weights!(snes.tmp_Utilities, rankedCandidates, snes.sortedUtilities)

    # Calc gradient
    gradient_mu = snes.last_s * u
    gradient_sigma = (abs2.(snes.last_s) .- 1) * u

    # Update the mean and sigma vectors based on the gradient
    old_mu = copy(snes.mu)
    snes.mu += snes.mu_learnrate * (snes.sigma .* gradient_mu)
    @inbounds for i in eachindex(snes.sigma)
        snes.sigma[i] = clamp(snes.sigma[i] * exp(0.5 * snes.sigma_learnrate * gradient_sigma[i]),
                              0.0, snes.max_sigma)
    end
    apply!(snes.embed, snes.mu, old_mu)

    # There is no notion of how many was better in NES so return 0
    return 0
end

"""
    assign_weights!(weights, rankedCandidates, sortedWeights)

Assigns the candidate `weights` according to the candidate index.
`rankedCandidates` are ranked by their fitness, `sortedWeights` are
the corresponding weights.

Returns candidate weights sorted by the individual's index in the population.
"""
function assign_weights!(weights::Vector{Float64},
                         rankedCandidates::AbstractVector{<:Candidate},
                         sortedWeights::Vector{Float64})
    @assert length(weights) == length(sortedWeights) && length(rankedCandidates) == length(weights)
    for i in eachindex(rankedCandidates)
        weights[rankedCandidates[i].index] = sortedWeights[i]
    end
    return weights
end

#  Traces the current `sNES` optimization state.
#  Called by `OptRunController` `trace_progress()`.
function trace_state(io::IO, snes::SeparableNESOpt, mode::Symbol)
    println(io, "|σ|=", norm(snes.sigma))
end

"""
Abstract type for a family of NES methods that represent population as
```
x = μ + σ B⋅Z,
```
where `B` is an exponential of some symmetric matrix `lnB`, `tr(lnB)==0.0`
"""
abstract type ExponentialNaturalEvolutionStrategyOpt <: NaturalEvolutionStrategyOpt end

function update_candidates!(exnes::ExponentialNaturalEvolutionStrategyOpt, Z::Matrix)
    B = exp(exnes.ln_B)
    sBZ = mul!(exnes.tmp_sBZ, B, Z)
    sBZ .*= exnes.sigma
    for i in eachindex(exnes.candidates)
        candi = exnes.candidates[i]
        candi.index = i # reset ordering
        @inbounds for j in 1:length(candi.params)
            candi.params[j] = exnes.x[j] + sBZ[j, i]
        end
        apply!(exnes.embed, candi.params, exnes.x)
        candi.fitness = NaN # FIXME: use nafitness()
    end

    return exnes.candidates
end

function update_parameters!(exnes::ExponentialNaturalEvolutionStrategyOpt, u::AbstractVector)
    # TODO use syrk(Z,dA) to speed-up multiplication
    Zu = broadcast!(*, exnes.tmp_Zu, exnes.Z, transpose(u))
    ln_dB = mul!(exnes.tmp_lndB, Zu, transpose(exnes.Z))
    sumU = sum(u)
    dSigma = 0.0
    @inbounds for i in 1:size(ln_dB, 1)
        ln_dB[i, i] -= sumU
        dSigma += ln_dB[i, i]
    end
    dSigma /= size(ln_dB, 1)
    @inbounds for i in 1:size(ln_dB, 1)
        ln_dB[i, i] -= dSigma
    end
    ln_dB .*= 0.5 * exnes.B_learnrate

    prev_x = copyto!(exnes.tmp_x, exnes.x)
    dx = mul!(exnes.tmp_dx, exnes.tmp_sBZ, u)
    dx .*= exnes.x_learnrate
    exnes.x += dx
    apply!(exnes.embed, exnes.x, prev_x)

    exnes.ln_B += ln_dB

    new_sigma = exnes.sigma * exp(0.5 * exnes.sigma_learnrate * dSigma)
    if isfinite(new_sigma) && new_sigma > 0
        exnes.sigma = min(new_sigma, exnes.max_sigma)
    end
    return exnes
end

" Identity for generic search space "
ini_xnes_B(ss::SearchSpace) = Matrix{Float64}(I, numdims(ss), numdims(ss))

"""
Calculates the initial ``log B`` matrix for `xNES` based on the deltas of each dimension.
"""
function ini_xnes_B(ss::RectSearchSpace)
    diag = log.(dimdelta(ss))
    diag .-= mean(diag)
    return Matrix{Float64}(Diagonal(diag))
end

"""
`xNES` method.

Nice but scales badly with increasing dimensions.
"""
mutable struct XNESOpt{F,E<:EmbeddingOperator} <: ExponentialNaturalEvolutionStrategyOpt
    embed::E                        # operator embedding into the search space
    lambda::Int                     # number of samples to take per iteration
    sortedUtilities::Vector{Float64}# the fitness shaping utility vector
    tmp_Utilities::Vector{Float64}   # the fitness shaping utility vector sorted by current population fitness
    x_learnrate::Float64
    sigma_learnrate::Float64
    B_learnrate::Float64
    max_sigma::Float64

    # TODO use Symmetric{Float64} to improve exponent etc calculation
    ln_B::Matrix{Float64}           # `log` of "covariation" matrix
    sigma::Float64                  # step size
    x::Individual                   # center of the population (aka most likely value, `mu` etc)
    Z::Matrix{Float64}              # current `N(0,I)` samples
    candidates::Vector{Candidate{F}}# the last sampled values, now being evaluated

    # temporary variables to minimize GC overhead
    tmp_x::Individual
    tmp_dz::Individual
    tmp_dx::Individual
    tmp_lndB::Matrix{Float64}
    tmp_Zu::Matrix{Float64}
    tmp_sBZ::Matrix{Float64}

    function XNESOpt{F}(embed::E; lambda::Int = 0, mu_learnrate::Float64 = 1.0,
            sigma_learnrate = 0.0, B_learnrate::Float64 = 0.0,
            ini_x = nothing, ini_sigma::Float64 = 1.0, ini_lnB = nothing,
            max_sigma::Float64 = 1.0E+10) where {F, E<:EmbeddingOperator}
        d = numdims(search_space(embed))
        if lambda == 0
            lambda = 4 + 3*floor(Int, log(d))
        end
        if sigma_learnrate == 0.0
            sigma_learnrate = 0.5 * min(1.0 / d, 0.25) # 0.6 * (3+log(d))/(d*sqrt(d)) in other literature
        end
        if B_learnrate == 0.0
            B_learnrate = 0.6 * (3+log(d))/(d*sqrt(d)) # 0.5 * min(1.0 / d, 0.25)
        end
        if ini_x === nothing
            ini_x = rand_individual(search_space(embed))
        else
            ini_x = copy(ini_x::Individual)
            apply!(embed, ini_x, rand_individual(search_space(embed)))
        end

        new{F,E}(embed, lambda, fitness_shaping_utilities_log(lambda),
                 fill!(Vector{Float64}(undef, lambda), NaN),
                 mu_learnrate, sigma_learnrate, B_learnrate, max_sigma,
                 ini_lnB === nothing ? ini_xnes_B(search_space(embed)) : ini_lnB, ini_sigma, ini_x, fill(0.0, d, lambda),
                 [Candidate{F}(fill!(Individual(undef, d), NaN), i) for i in 1:lambda],
                  # temporaries
                  Vector{Float64}(undef, d), Vector{Float64}(undef, d), Vector{Float64}(undef, d),
                  Matrix{Float64}(undef, d, d),
                  Matrix{Float64}(undef, d, lambda), Matrix{Float64}(undef, d, lambda),
        )
    end
end

set_candidate!(o::XNESOpt, x0) = (o.x = x0)
candidate(o::XNESOpt) = o.x

const XNES_DefaultOptions = chain(NES_DefaultOptions, ParamsDict(
    :B_learnrate => 0.0,   # If 0.0 it will be set based on the number of dimensions
    :ini_sigma => 1.0,     # Initial sigma (step size)
    :ini_lnB => nothing    # Initial log(B)
))

function xnes(problem::OptimizationProblem, parameters)
    params = chain(XNES_DefaultOptions, parameters)
    embed = RandomBound(search_space(problem))
    F = fitness_type(problem)
    XNESOpt{F}(embed; lambda = params[:lambda],
               mu_learnrate = params[:mu_learnrate],
               sigma_learnrate = params[:sigma_learnrate],
               B_learnrate = params[:B_learnrate],
               ini_x = params[:ini_x],
               ini_sigma = params[:ini_sigma],
               ini_lnB = params[:ini_lnB],
               max_sigma = params[:max_sigma])
end

function ask(xnes::XNESOpt)
    randn!(xnes.Z)
    update_candidates!(xnes, xnes.Z)
end

function tell!(xnes::XNESOpt{F}, rankedCandidates::Vector{Candidate{F}}) where F
    u = assign_weights!(xnes.tmp_Utilities, rankedCandidates, xnes.sortedUtilities)

    update_parameters!(xnes, u)
    return 0
end

function trace_state(io::IO, xnes::XNESOpt, mode::Symbol)
    println(io, "sigma=", xnes.sigma,
                " |trace(ln_B)|=", tr(xnes.ln_B))
end

"""
    fitness_shaping_utilities_log(n)

Calculate the `n`-dimensional fitness shaping utilities vector using the "log" method.
"""
function fitness_shaping_utilities_log(n::Int)
    u = max.(0.0, log(n / 2 + 1.0) .- log.(1:n))
    return u./sum(u) .- 1/n
end

"""
    fitness_shaping_utilities_linear(n)

Calculate the `n`-dimensional fitness shaping utilities vector
using the "steps" method.
"""
function fitness_shaping_utilities_linear(n::Int)
    # Second half has zero utility.
    threshold = n÷2
    second_half = fill(0.0, n - threshold)

    # While first half's utility decreases in linear steps
    first_half = range(1.0, stop=1/threshold, length=threshold)

    # But the utilities should sum to 0, so we normalize and return
    u = vcat(first_half, second_half)
    return u./sum(u) .- 1/n
end


text:
function test_shape()
    cell = load_mesh(joinpath(data_folder,"cell.obj"))
    new_cell = mesh_centering(cell)
    @test all(faces(cell) .== faces(new_cell))
    new_centroid = convert(Vector{Float64},mean(coordinates(new_cell)))
    @test isapprox(new_centroid, [0., 0., 0.],atol =0.0001)
    new_cell2 = volume_normalizing(cell)
    @test all(faces(cell) .== faces(new_cell2))
    @test isapprox(mesh_volume(new_cell2), 1.0,atol =0.0001)
end


text:
#=

# Setting preferences

To set the [preferences](@ref Preferences) `LOG_DESTRUCTION` to "true" for debugging, you can do it either by having a `LocalPreferences.toml` in your environment or by calling `Vulkan.set_preferences!`:

=#

using Vulkan

Vulkan.set_preferences!("LOG_DESTRUCTION" => "true")


text:
module SVM

using StatsBase
import StatsBase: predict

export svm, cddual, pegasos, predict

type SVMFit
    w::Vector{Float64}
    pass::Int
    converged::Bool
end

function Base.show(io::IO, fit::SVMFit)
    @printf(io, "Fitted linear SVM\n")
    @printf(io, " * Non-zero weights: %d\n", countnz(fit.w))
    @printf(io, " * Iterations: %d\n", fit.pass)
    @printf(io, " * Converged: %s\n", string(fit.converged))
end

function predict{T<:Real}(fit::SVMFit, X::AbstractMatrix{T})
    n, l = size(X)
    preds = Array{Float64,1}(l)
    for i in 1:l
        tmp = 0.0
        for j in 1:n
            tmp += fit.w[j] * X[j, i]
        end
        preds[i] = sign(tmp)
    end
    return preds
end

include("pegasos.jl")
include("cddual.jl")

function svm{T<:Real}(X::AbstractMatrix{T},
                      Y::AbstractVector{T};
                      k::Integer = 5,
                      lambda::Real = 0.1,
                      maxpasses::Integer = 100)
    pegasos(X, Y, k = k, lambda = lambda, maxpasses = maxpasses)
end

end # module


text:
import Base.Operators: +, -, ^, /, //, \, *, ==

## equality
function ==(b1::SymbolicType, b2::SymbolicType)
    b1,b2 = map(Basic, (b1, b2))
    ccall((:basic_eq, libsymengine), Int, (Ptr{Basic}, Ptr{Basic}), &b1, &b2) == 1
end

## main ops
for (op, libnm) in ((:+, :add), (:-, :sub), (:*, :mul), (:/, :div), (://, :div), (:^, :pow))
    tup = (Base.Symbol("basic_$libnm"), libsymengine)
    @eval begin
        function ($op)(b1::Basic, b2::Basic)
            a = Basic()
            ccall($tup, Void, (Ptr{Basic}, Ptr{Basic}, Ptr{Basic}), &a, &b1, &b2)
            return a
        end
        ($op)(b1::BasicType, b2::BasicType) = ($op)(Basic(b1), Basic(b2))
    end
end

^{T<:SymbolicType, S <: Integer}(a::T, b::S) = Basic(a)^Basic(b)
^{T<:SymbolicType, S <: Rational}(a::T, b::S) = Basic(a)^Basic(b)
+(b::SymbolicType) = b
-(b::SymbolicType) = 0 - b
\(b1::SymbolicType, b2::SymbolicType) = b2 / b1

## ## constants
Base.zero(x::Basic) = Basic(0)
Base.zero{T<:Basic}(::Type{T}) = Basic(0)
Base.one(x::Basic) = Basic(1)
Base.one{T<:Basic}(::Type{T}) = Basic(1)

Base.zero(x::BasicType) = BasicType(Basic(0))
Base.zero{T<:BasicType}(::Type{T}) = BasicType(Basic(0))
Base.one(x::BasicType) = BasicType(Basic(1))
Base.one{T<:BasicType}(::Type{T}) = BasicType(Basic(1))

## Math constants
## no oo!

for op in [:IM, :PI, :E, :EulerGamma, :Catalan, :oo, :zoo, :NAN]
    @eval begin
        const $op = Basic(C_NULL)
    end
    eval(Expr(:export, op))
end

macro init_constant(op, libnm)
    tup = (Base.Symbol("basic_const_$libnm"), libsymengine)
    alloc_tup = (:basic_new_stack, libsymengine)
    :(
        begin
            ccall($alloc_tup, Void, (Ptr{Basic}, ), &($op))
            ccall($tup, Void, (Ptr{Basic}, ), &($op))
            finalizer($op, basic_free)
        end
    )
end

function init_constants()
    @init_constant IM I
    @init_constant PI pi
    @init_constant E E
    @init_constant EulerGamma EulerGamma
    @init_constant Catalan Catalan
    @init_constant oo infinity
    @init_constant zoo complex_infinity
    @init_constant NAN nan
end

## ## Conversions
Base.convert(::Type{Basic}, x::Irrational{:π}) = PI
Base.convert(::Type{Basic}, x::Irrational{:e}) = E
Base.convert(::Type{Basic}, x::Irrational{:γ}) = EulerGamma
Base.convert(::Type{Basic}, x::Irrational{:catalan}) = Catalan
Base.convert(::Type{Basic}, x::Irrational{:φ}) = (1 + Basic(5)^Basic(1//2))/2
Base.convert(::Type{BasicType}, x::Irrational) = BasicType(convert(Basic, x))

## Logical operators
@compat Base.:<(x::SymbolicType, y::SymbolicType) = N(x) < N(y)
@compat Base.:<(x::SymbolicType, y) = <(promote(x,y)...)
@compat Base.:<(x, y::SymbolicType) = <(promote(x,y)...)


text:
using Documenter, Univariates

makedocs(;
    modules=[Univariates],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/wakakusa/Univariates.jl/blob/{commit}{path}#L{line}",
    sitename="Univariates.jl",
    authors="wakakusa",
    assets=String[],
)

deploydocs(;
    repo="github.com/wakakusa/Univariates.jl",
)


text:
@generated function _index(::Val{name}) where {name}
    #:(Int((ccall("llvm.nvvm.read.ptx.sreg.ctaid.x", llvmcall, UInt32, ()))+UInt32(1)))

    s = "llvm.nvvm.read.ptx.sreg.ctaid.$name"
    :(Int((ccall($s, llvmcall, UInt32, ()))+UInt32(1)))
end

#@inline blockIdx_x() = Int((ccall("llvm.nvvm.read.ptx.sreg.ctaid.x", llvmcall, UInt32, ()))+UInt32(1))
@inline blockIdx_x() = _index(Val(:x))
@inline blockIdx_y() = _index(Val(:y))
@inline blockIdx_z() = _index(Val(:z))

@inline blockIdx() = (x=blockIdx_x(), y=blockIdx_y(), z=blockIdx_z())

function kernel_vadd()
    blockIdx().x-1
end

code_llvm(stdout, kernel_vadd, Tuple{})

using TinyCassette
struct Ctx end

code_llvm(stdout, TinyCassette.execute,
          Tuple{Ctx,
                typeof(kernel_vadd)})


text:
# LP programming using MATLAB-like construction
#
#   maximize 1000 x + 350 y
#
#       s.t. x >= 30, y >= 0
#            x - 1.5y >= 0  (i.e. -x + 1.5 y <= 0)
#            12 x + 8 y <= 1000
#            1000 x + 300 y <= 70000
#
#   solution: (59.0909, 36.3636)
#   objv: 71818.1818

using CPLEX, Base.Test

@testset "LP 02" begin
    env = CPLEX.Env()

    model = CPLEX.cplex_model(env;
        name="lp_02",
        sense=:Max,
        f = [1000., 350.],
        A = [-1. 1.5; 12. 8.; 1000. 300.],
        b = [0., 1000., 70000.],
        lb = [0., 30.])

    CPLEX.optimize!(model)

    sol = CPLEX.get_solution(model)
    @test sol[1] ≈ 650/11
    @test sol[2] ≈ 400/11

    @test CPLEX.get_objval(model) ≈ 790000/11
end


text:
# 2^15 = 32768 and the sum of its digits is 3 + 2 + 7 + 6 + 8 = 26.
#
# What is the sum of the digits of the number 2^1000?

using ProjectEulerSolutions

# Simple with Julia's built-in BigInt and the digits function.
function p016solution(mantissa::Integer=2, exponent::Integer=3)::Integer
    return sum(digits(BigInt(mantissa) ^ exponent))
end

p016 = Problems.Problem(p016solution)

Problems.benchmark(p016, 2, 1000)

max_stars_repo_name: gnujosh/julia-euler-project

max_issues_repo_name: gnujosh/julia-euler-project

max_forks_repo_name: gnujosh/julia-euler-project


text:
using Test

@testset "VSCodeDebugger" begin
end


text:
#Subject
"""
    setdosetime!(data::T, dosetime::DoseTime) where T <: PKSubject

Set dose time `dosetime` for subject `data`.
"""
function setdosetime!(data::T, dosetime::DoseTime) where T <: PKSubject
    data.dosetime = dosetime
    data
end
#DS ind Int
"""
    setdosetime!(data::DataSet{T}, dosetime::DoseTime, ind::Int) where T <: PKSubject

* `ind` - index in DataSet.
"""
function setdosetime!(data::DataSet{T}, dosetime::DoseTime, ind::Int) where T <: PKSubject
    setdosetime!(data[ind], dosetime)
    data
end
#DS iter Int
"""
    setdosetime!(data::DataSet{T}, dosetime::DoseTime, inds::Union{Vector{Int}, UnitRange{Int}, Tuple{Vararg{Int}}}) where T <: PKSubject

* `inds` - indexes in DataSet.
"""
function setdosetime!(data::DataSet{T}, dosetime::DoseTime, inds::Union{Vector{Int}, UnitRange{Int}, Tuple{Vararg{Int}}}) where T <: PKSubject
    for i in inds
        setdosetime!(data[i], dosetime)
    end
    data
end
#DS all
"""
    setdosetime!(data::DataSet{T}, dosetime::DoseTime) where T <: PKSubject

For all subjects in DataSet.
"""
function setdosetime!(data::DataSet{T}, dosetime::DoseTime) where T <: PKSubject
    for i = 1:length(data)
        setdosetime!(data[i], dosetime)
    end
    data
end
#DS Dict
"""
    setdosetime!(data::DataSet{T}, dosetime::DoseTime, sort::Dict) where T <: PKSubject

Set dose time `dosetime` for subjects if `sort` ⊆ subject's `id`.
"""
function setdosetime!(data::DataSet{T}, dosetime::DoseTime, sort::Dict) where T <: PKSubject
    for i = 1:length(data)
        if sort ⊆ data[i].id setdosetime!(data[i], dosetime) end
    end
    data
end
#GET subj
"""
    getdosetime(data::T) where T <: PKSubject

Return dosetime.
"""
function getdosetime(data::T) where T <: PKSubject
    data.dosetime
end

#Subject
#DS ind Int
#DS iter Int
#DS all
#DS Dict
#GET subj


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

using FEMBasis
using LinearAlgebra
using Test

# interpolate scalar field in unit square:
# T(X,t) = t*(1 + X[1] + 3*X[2] - 2*X[1]*X[2])
B = Quad4()
X = ((0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0))
T = (1.0, 2.0, 3.0, 4.0)
T_known(X) = 1 + X[1] + 3*X[2] - 2*X[1]*X[2]
T_interpolated = interpolate(B, T, (0.0, 0.0))
@test isapprox(T_interpolated, T_known((0.5, 0.5)))

# calculate gradient dB/dX
B = Quad4()
X = ([0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0])
dBdX = grad(B, X, (0.0, 0.0))
@test isapprox(dBdX, 1/2*[-1 1 1 -1; -1 -1 1 1])

# interpolate gradient of scalar field in unit square:
# grad(T)(X) = [1-2X[2], 3-2*X[1]]
B = Quad4()
X = ([0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0])
T = (1.0, 2.0, 3.0, 4.0)
dTdX = grad(B, T, X, (0.0, 0.0))
dTdX_expected(X) = [1-2*X[2] 3-2*X[1]]
@test isapprox(dTdX, dTdX_expected([0.5, 0.5]))

# interpolate vector field in unit square:
# u(X,t) = [1/4*t*X[1]*X[2], 0, 0]
B = Quad4()
u = Vector[[0.0, 0.0], [0.0, 0.0], [1/4, 0.0], [0.0, 0.0]]
u_known(X) = [1/4*X[1]*X[2], 0]
u_interpolated = interpolate(B, u, (0.0, 0.0))
@test isapprox(u_interpolated, u_known((0.5, 0.5)))

# interpolate gradient of vector field in unit square:
# u(X) = t*[X[1]*(X[2]+1), X[1]*(4*X[2]-1)]
# => u_i,j = t*[X[2]+1 X[1]; 4*X[2]-1 4*X[1]]
B = Quad4()
X = ([0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0])
u = ([0.0, 0.0], [1.0, -1.0], [2.0, 3.0], [0.0, 0.0])
dudX = grad(B, u, X, (0.0, 0.0))
dudX_expected(X) = [X[2]+1 X[1]; 4*X[2]-1 4*X[1]]
@test isapprox(dudX, dudX_expected([0.5, 0.5]))

# test BasisInfo
B = BasisInfo(Seg2)
@test length(B) == 2
@test size(B) == (1, 2)
X = ((0.0,), (1.1))
xi = (0.0,)
eval_basis!(B, X, xi)
@test isapprox(B.invJ, inv(B.J))
@test isapprox(B.detJ, det(B.J))
B = BasisInfo(Quad4)
X = ((0.0,0.0), (1.1,0.0), (1.0,1.0), (0.0,1.0))
xi = (0.0,0.0)
eval_basis!(B, X, xi)
@test isapprox(B.invJ, inv(B.J))
@test isapprox(B.detJ, det(B.J))
B = BasisInfo(Hex8)
X = ((0.0,0.0,0.0), (1.1,0.0,0.0), (1.0,1.0,0.0), (0.0,1.0,0.0),
     (0.0,0.0,1.0), (1.0,0.0,1.0), (1.0,1.0,1.0), (0.0,1.0,1.0))
xi = (0.0,0.0,0.0)
eval_basis!(B, X, xi)
@test isapprox(B.invJ, inv(B.J))
@test isapprox(B.detJ, det(B.J))

# evaluate gradient using BasisInfo
B = BasisInfo(Quad4)
X = ((0.0,0.0), (1.0,0.0), (1.0,1.0), (0.0,1.0))
xi = (0.0, 0.0)
eval_basis!(B, X, xi)
u = ((0.0, 0.0), (1.0, -1.0), (2.0, 3.0), (0.0, 0.0))
gradu = zeros(2, 2)
grad!(B, gradu, u)
gradu_expected = [1.5 0.5; 1.0  2.0]
@test isapprox(gradu, gradu_expected)

# test curves
bi = BasisInfo(Seg2)
X1 = ((0.0,0.0,0.0), (1.0,0.0,0.0))
X2 = ((0.0,0.0,0.0), (0.0,1.0,0.0))
X3 = ((0.0,0.0,0.0), (0.0,0.0,1.0))
xi = (0.0, 0.0)
eval_basis!(bi, X1, xi)
@test isapprox(bi.detJ, 0.5)
eval_basis!(bi, X2, xi)
@test isapprox(bi.detJ, 0.5)
eval_basis!(bi, X3, xi)
@test isapprox(bi.detJ, 0.5)
@test isapprox(jacobian(Seg2(), X1, xi), [0.5 0.0 0.0])
@test isapprox(jacobian(Seg2(), X2, xi), [0.0 0.5 0.0])
@test isapprox(jacobian(Seg2(), X3, xi), [0.0 0.0 0.5])

# test manifolds
bi = BasisInfo(Quad4)
X1 = ((0.0,0.0,0.0), (1.0,0.0,0.0), (1.0,1.0,0.0), (0.0,1.0,0.0))
X2 = ((0.0,0.0,0.0), (0.0,1.0,0.0), (0.0,1.0,1.0), (0.0,0.0,1.0))
X3 = ((0.0,0.0,0.0), (0.0,0.0,1.0), (1.0,0.0,1.0), (1.0,0.0,0.0))
xi = (0.0, 0.0)
eval_basis!(bi, X1, xi)
@test isapprox(bi.detJ, 0.25)
eval_basis!(bi, X2, xi)
@test isapprox(bi.detJ, 0.25)
eval_basis!(bi, X3, xi)
@test isapprox(bi.detJ, 0.25)
J1 = jacobian(Quad4(), X1, xi)
J2 = jacobian(Quad4(), X2, xi)
J3 = jacobian(Quad4(), X3, xi)
@test isapprox(J1, [0.5 0.0 0.0; 0.0 0.5 0.0])
@test isapprox(J2, [0.0 0.5 0.0; 0.0 0.0 0.5])
@test isapprox(J3, [0.0 0.0 0.5; 0.5 0.0 0.0])


text:
using TropicalGEMM
using Test

@testset "TropicalGEMM.jl" begin
    include("gemm.jl")
end


text:
function main(args)
    return args
end
# Load inline arguments for worker unit
num_workers, this = main(ARGS)
num_workers = parse(Int64, num_workers)
this = parse(Int64, this)
# num_workers, this = 10,1
println("Worker " * string(this) * " of " * string(num_workers))

include("../src/init.jl")
include("../src/fct.jl")
include("../src/fitting-fct.jl")
include("../src/worker_fct.jl")
data_type = "cal";
set = "sim";
high_cut = 4;

plots_path = "../plots/sim/"
base_path  = joinpath("../../waveforms/sim/", "wf")
base_path_AE  = joinpath("../../waveforms/sim/", "calib_AE")

calib_filepath = "../dicts/calib.json"
calib = JSON.parsefile(calib_filepath)
cut_lib_filepath = "../dicts/cut_lib.json"
cut_lib = JSON.parsefile(cut_lib_filepath)
AE_cal_filepath = "../dicts/AE_cal.json"
AE_cal = JSON.parsefile(AE_cal_filepath)
det_lib_filepath = "../dicts/det_lib.json"
det_lib = JSON.parsefile(det_lib_filepath)

channels = []
for ch in 0:1:36
    if ctb[ch] && ch != 6
        push!(channels, ch)
    end
end
for ch in channels
    ch_str = lpad(ch, 2, "0");
    if !channel_to_bege[ch]
        println("STOP here! This is a coax detector.")
    end
    if !haskey(calib, channel_to_name[ch])
        calib[channel_to_name[ch]] = Dict()
        calib[channel_to_name[ch]]["data"] = Dict()
        calib[channel_to_name[ch]]["sim"] = Dict()
    end

    filepath = joinpath(base_path_AE, ch_str * "-" * ctn[ch] * "-" * string(this) * ".h5")
    !isdir(dirname(filepath)) ? mkpath(dirname(filepath)) : ""
    if !isfile(filepath)
        numofele_A = 5
        BackDelta5 = div(numofele_A,2)
        ForwardDelta5 = isodd(numofele_A) ? div(numofele_A,2) : div(numofele_A,2) - 1
        numofele_E = 201
        BackDelta201 = div(numofele_E,2)
        ForwardDelta201 = isodd(numofele_E) ? div(numofele_E,2) : div(numofele_E,2) - 1

        files = glob(joinpath(base_path, ch_str * "-" * channel_to_name[ch] * "/*/*.h5"));
        files = get_share_for_worker(files, num_workers, this)
        data = Table(energy = [], A = [], E = [])
        @showprogress 1 "Read data and calculate A & E for Ch$ch_str" for file in files
            tmp = HDF5.h5open(file, "r") do h5f
                LegendHDF5IO.readdata(h5f, "data")
            end
            A = []
            E = []
            for wf in tmp.waveform
                push!(A, maximum(movingaverage(diff(wf.value),numofele_A,BackDelta5,ForwardDelta5,3)))
                push!(E, maximum(movingaverage(diff(wf.value),numofele_E,BackDelta201,ForwardDelta201,13)))
            end
            if size(tmp,1) > 1
                append!(data, Table(energy = sum.(tmp.energy)*1000, A = A, E = E))
            end
        end
        HDF5.h5open(filepath, "w") do h5f
            LegendHDF5IO.writedata(h5f, "data", Table(
                A = float.(data.A),
                E = float.(data.E),
                energy = float.(data.energy)
            ))
        end
        data = nothing
    end
end


text:
using DataFrames
using Test

m = page_model()
include("../src/components/AdaptationCosts.jl")

adaptationcosts_noneconomic = MimiPAGE2009.addadaptationcosts_noneconomic(m)
adaptationcosts_noneconomic[:y_year_0] = 2008.
adaptationcosts_noneconomic[:y_year] = Mimi.dim_keys(m.md, :time)
adaptationcosts_noneconomic[:gdp] = readpagedata(m, "test/validationdata/gdp.csv")

adaptationcosts_economic = MimiPAGE2009.addadaptationcosts_economic(m)
adaptationcosts_economic[:y_year_0] = 2008.
adaptationcosts_economic[:y_year] = Mimi.dim_keys(m.md, :time)
adaptationcosts_economic[:gdp] = readpagedata(m, "test/validationdata/gdp.csv")

adaptationcosts_sealevel = MimiPAGE2009.addadaptationcosts_sealevel(m)
adaptationcosts_sealevel[:y_year_0] = 2008.
adaptationcosts_sealevel[:y_year] = Mimi.dim_keys(m.md, :time)
adaptationcosts_sealevel[:gdp] = readpagedata(m, "test/validationdata/gdp.csv")

p = load_parameters(m)
set_leftover_params!(m, p)

run(m)

autofac = m[:AdaptiveCostsEconomic, :autofac_autonomouschangefraction]
atl_economic = m[:AdaptiveCostsEconomic, :atl_adjustedtolerablelevel]
imp_economic = m[:AdaptiveCostsEconomic, :imp_adaptedimpacts]
acp_economic = m[:AdaptiveCostsEconomic, :acp_adaptivecostplateau]
aci_economic = m[:AdaptiveCostsEconomic, :aci_adaptivecostimpact]

autofac_compare = readpagedata(m, "test/validationdata/autofac_autonomouschangefraction.csv")
atl_economic_compare = readpagedata(m, "test/validationdata/atl_1.csv")
imp_economic_compare = readpagedata(m, "test/validationdata/imp_1.csv")
acp_economic_compare = readpagedata(m, "test/validationdata/acp_adaptivecostplateau_economic.csv")
aci_economic_compare = readpagedata(m, "test/validationdata/aci_adaptivecostimpact_economic.csv")

@test autofac ≈ autofac_compare rtol=1e-8
@test atl_economic ≈ atl_economic_compare rtol=1e-8
@test imp_economic ≈ imp_economic_compare rtol=1e-8
@test acp_economic ≈ acp_economic_compare rtol=1e-3
@test aci_economic ≈ aci_economic_compare rtol=1e-3

ac_noneconomic = m[:AdaptiveCostsNonEconomic, :ac_adaptivecosts]
ac_economic = m[:AdaptiveCostsEconomic, :ac_adaptivecosts]
ac_sealevel = m[:AdaptiveCostsSeaLevel, :ac_adaptivecosts]

ac_noneconomic_compare = readpagedata(m, "test/validationdata/ac_adaptationcosts_noneconomic.csv")
ac_economic_compare = readpagedata(m, "test/validationdata/ac_adaptationcosts_economic.csv")
ac_sealevel_compare = readpagedata(m, "test/validationdata/ac_adaptationcosts_sealevelrise.csv")

@test ac_noneconomic ≈ ac_noneconomic_compare rtol=1e-2
@test ac_economic ≈ ac_economic_compare rtol=1e-3
@test ac_sealevel ≈ ac_sealevel_compare rtol=1e-2


text:
@testset "Hubbard Model" begin
    @test HubbardModel <: Model

    # constructors
    m = HubbardModel(8, 1)
    @test length(m.l) == 8
    @test typeof(m) == HubbardModel{Chain}

    l = SquareLattice(5)
    m = HubbardModel(l = l);
    @test length(m.l) == 25
    @test typeof(m) == HubbardModel{SquareLattice}

    m = HubbardModel(l, t = 0.5);
    @test length(m.l) == 25
    @test m.t == 0.5
    @test typeof(m) == HubbardModel{SquareLattice}

    d = Dict(:l => l, :U => HubbardModel == HubbardModelAttractive ? 0.5 : -0.5)
    m = HubbardModel(d)
    @test typeof(m) == HubbardModel{SquareLattice}
    @test m.U == (HubbardModel == HubbardModelAttractive ? 0.5 : -0.5)
end


text:
function compareField(Allc, Bllc, syms)::Bool
  (!isdefined(Allc, syms) && !isdefined(Bllc, syms)) && return true
  !isdefined(Allc, syms) && return false
  !isdefined(Bllc, syms) && return false
  return eval(:($Allc.$syms == $Bllc.$syms))
end

"""
    $(SIGNATURES)

Compare the all fields of T that are not in `skip` for objects `Al` and `Bl`.

TODO > add to func_ref.md
"""
function compareFields(Al::T,
                       Bl::T;
                       show::Bool=true,
                       skip::Vector{Symbol}=Symbol[]  )::Bool where {T}
  for field in fieldnames(T)
    (field in skip) && continue
    tp = compareField(Al, Bl, field)
    show && @debug("  $tp : $field")==nothing
    !tp && return false
  end
  return true
end

function compareFields(Al::T,
                       Bl::T;
                       show::Bool=true,
                       skip::Vector{Symbol}=Symbol[]  )::Bool where {T <: Union{Number, AbstractString}}
  #
  return Al == Bl
end

function compareAll(Al::T,
                    Bl::T;
                    show::Bool=true,
                    skip::Vector{Symbol}=Symbol[]  )::Bool where {T <: Union{AbstractString,Symbol}}
  #
  return Al == Bl
end

function compareAll(Al::T,
                    Bl::T;
                    show::Bool=true,
                    skip::Vector{Symbol}=Symbol[]  )::Bool where {T <: Union{Array{<:Number}, Number}}
  #
  (length(Al) != length(Bl)) && return false
  return norm(Al - Bl) < 1e-6
end

function compareAll(Al::T,
                    Bl::T;
                    show::Bool=true,
                    skip::Vector{Symbol}=Symbol[]  )::Bool where {T <: Array}
  #
  (length(Al) != length(Bl)) && return false
  for i in 1:length(Al)
    !compareAll(Al[i],Bl[i], show=false) && return false
  end
  return true
end

"""
    $(SIGNATURES)

Recursively compare the all fields of T that are not in `skip` for objects `Al` and `Bl`.

TODO > add to func_ref.md
"""
function compareAll(Al::T,
                    Bl::T;
                    show::Bool=true,
                    skip::Vector{Symbol}=Symbol[]  )::Bool where {T <: Tuple}
  #
  length(Al) != length(Bl) && return false
  for i in 1:length(Al)
    !compareAll(Al[i], Bl[i], show=show, skip=skip) && return false
  end
  return true
end

function compareAll(Al::T,
                    Bl::T;
                    show::Bool=true,
                    skip::Vector{Symbol}=Symbol[]  )::Bool where {T <: Dict}
  #
  (length(Al) != length(Bl)) && return false
  for (id, val) in Al
    (Symbol(id) in skip) && continue
    !compareAll(val, Bl[id], show=show, skip=skip) && return false
  end
  return true
end

function compareAll(Al::T, Bl::T; show::Bool=true, skip::Vector{Symbol}=Symbol[])::Bool where T
  @debug "Comparing types $T:"
  @debug "  Al = $Al"
  @debug "  Bl = $Bl"
  !compareFields(Al, Bl, show=show, skip=skip) && return false
  for field in fieldnames(T)
    field in skip && continue
    @debug("  Checking field: $field")
    (!isdefined(Al, field) && !isdefined(Al, field)) && return true
    !isdefined(Al, field) && return false
    !isdefined(Bl, field) && return false
    Ad = eval(:($Al.$field))
    Bd = eval(:($Bl.$field))
    !compareAll(Ad, Bd, show=show, skip=skip) && return false
  end
  return true
end

"""
    $SIGNATURES

Compare that all fields are the same in a `::FactorGraph` variable.
"""
function compareVariable(A::DFGVariable,
                         B::DFGVariable;
                         skip::Vector{Symbol}=Symbol[],
                         show::Bool=true,
                         skipsamples::Bool=true  )::Bool
  #
  skiplist = union([:attributes;:solverDataDict;:_internalId],skip)
  TP = compareAll(A, B, skip=skiplist, show=show)
  varskiplist = skipsamples ? [:val; :bw] : Symbol[]
  skiplist = union([:softtype;],varskiplist)
  union!(skiplist, skip)
  TP = TP && compareAll(A.solverDataDict, B.solverDataDict, skip=skiplist, show=show)

  Ad = solverData(A)
  Bd = solverData(B)

  # TP = TP && compareAll(A.attributes, B.attributes, skip=[:softtype;], show=show)
  varskiplist = union(varskiplist, [:softtype;:_internalId])
  union!(varskiplist, skip)
  TP = TP && compareAll(Ad, Bd, skip=varskiplist, show=show)
  TP = TP && typeof(Ad.softtype) == typeof(Bd.softtype)
  TP = TP && compareAll(Ad.softtype, Bd.softtype, show=show, skip=skip)
  return TP
end

function compareAllSpecial(A::T1,
                           B::T2;
                           skip=Symbol[],
                           show::Bool=true) where {T1 <: GenericFunctionNodeData, T2 <: GenericFunctionNodeData}
  if T1 != T2
    return false
  else
    return compareAll(A, B, skip=skip, show=show)
  end
end

"""
    $SIGNATURES

Compare that all fields are the same in a `::FactorGraph` factor.
"""
function compareFactor(A::DFGFactor,
                       B::DFGFactor;
                       show::Bool=true,
                       skip::Vector{Symbol}=Symbol[],
                       skipsamples::Bool=true,
                       skipcompute::Bool=true  )
  #
  TP =  compareAll(A, B, skip=union([:attributes;:data;:_variableOrderSymbols;:_internalId],skip), show=show)
  # TP = TP & compareAll(A.attributes, B.attributes, skip=[:data;], show=show)
  TP = TP & compareAllSpecial(solverData(A), solverData(B), skip=union([:fnc;:_internalId], skip), show=show)
  TP = TP & compareAllSpecial(solverData(A).fnc, solverData(B).fnc, skip=union([:cpt;:measurement;:params;:varidx;:threadmodel], skip), show=show)
  TP = TP & (skipsamples || compareAll(solverData(A).fnc.measurement, solverData(B).fnc.measurement, show=show, skip=skip))
  TP = TP & (skipcompute || compareAll(solverData(A).fnc.params, solverData(B).fnc.params, show=show, skip=skip))
  TP = TP & (skipcompute || compareAll(solverData(A).fnc.varidx, solverData(B).fnc.varidx, show=show, skip=skip))

  return TP
end
  # Ad = solverData(A)
  # Bd = solverData(B)
  # TP =  compareAll(A, B, skip=[:attributes;:data], show=show)
  # TP &= compareAll(A.attributes, B.attributes, skip=[:data;], show=show)
  # TP &= compareAllSpecial(solverData(A).fnc, solverData(B).fnc, skip=[:cpt;], show=show)
  # TP &= compareAll(solverData(A).fnc.cpt, solverData(B).fnc.cpt, show=show)

"""
    $SIGNATURES

Compare all variables in both `::FactorGraph`s A and B.

Notes
- A and B should all the same variables and factors.

Related:

`compareFactorGraphs`, `compareSimilarVariables`, `compareVariable`, `ls`
"""
function compareAllVariables(fgA::G1,
                             fgB::G2;
                             skip::Vector{Symbol}=Symbol[],
                             show::Bool=true,
                             skipsamples::Bool=true )::Bool where {G1 <: AbstractDFG, G2 <: AbstractDFG}
  # get all the variables in A or B
  xlA =  getVariableIds(fgA)
  xlB =  getVariableIds(fgB)
  vars = union(xlA, xlB)

  # compare all variables exist in both A and B
  TP = length(xlA) == length(xlB)
  for xla in xlA
    TP &= xla in xlB
  end
  # slightly redundant, but repeating opposite direction anyway
  for xlb in xlB
    TP &= xlb in xlA
  end

  # compare each variable is the same in both A and B
  for var in vars
    TP = TP && compareVariable(getVariable(fgA, var), getVariable(fgB, var), skipsamples=skipsamples, skip=skip)
  end

  # return comparison result
  return TP
end

"""
    $SIGNATURES

Compare similar labels between `::FactorGraph`s A and B.

Notes
- At least one variable label should exist in both A and B.

Related:

`compareFactorGraphs`, `compareAllVariables`, `compareSimilarFactors`, `compareVariable`, `ls`.
"""
function compareSimilarVariables(fgA::G1,
                                 fgB::G2;
                                 skip::Vector{Symbol}=Symbol[],
                                 show::Bool=true,
                                 skipsamples::Bool=true )::Bool where {G1 <: AbstractDFG, G2 <: AbstractDFG}
  #
  xlA = getVariableIds(fgA)
  xlB = getVariableIds(fgB)

  # find common variables
  xlAB = intersect(xlA, xlB)
  TP = length(xlAB) > 0

  # compare the common set
  for var in xlAB
    @info var
    TP &= compareVariable(getVariable(fgA, var), getVariable(fgB, var), skipsamples=skipsamples, skip=skip)
  end

  # return comparison result
  return TP
end

"""
    $SIGNATURES

Compare similar factors between `::FactorGraph`s A and B.

Related:

`compareFactorGraphs`, `compareSimilarVariables`, `compareAllVariables`, `ls`.
"""
function compareSimilarFactors(fgA::G1,
                               fgB::G2;
                               skipsamples::Bool=true,
                               skipcompute::Bool=true,
                               show::Bool=true  )::Bool where {G1 <: AbstractDFG, G2 <: AbstractDFG}
  #
  xlA = getFactorIds(fgA)
  xlB = getFactorIds(fgB)

  # find common variables
  xlAB = intersect(xlA, xlB)
  TP = length(xlAB) > 0

  # compare the common set
  for var in xlAB
    TP = TP && compareFactor(getFactor(fgA, var), getFactor(fgB, var), skipsamples=skipsamples, skipcompute=skipcompute, show=show)
  end

  # return comparison result
  return TP
end

"""
    $SIGNATURES

Compare and return if two factor graph objects are the same, by comparing similar variables and factors.

Notes:
- Default items to skip with `skipsamples`, `skipcompute`.
- User defined fields to skip can be specified with `skip::Vector{Symbol}`.

Related:

`compareSimilarVariables`, `compareSimilarFactors`, `compareAllVariables`, `ls`.
"""
function compareFactorGraphs(fgA::G1,
                             fgB::G2;
                             skipsamples::Bool=true,
                             skipcompute::Bool=true,
                             skip::Vector{Symbol}=Symbol[],
                             show::Bool=true  )::Bool where {G1 <: AbstractDFG, G2 <: AbstractDFG}
  #
  skiplist = Symbol[:g;:bn;:IDs;:fIDs;:id;:nodeIDs;:factorIDs;:fifo;:solverParams]
  skiplist = union(skiplist, skip)
  @warn "compareFactorGraphs will skip comparisons on: $skiplist"

  TP = compareAll(fgA, fgB, skip=skiplist, show=show)
  TP = TP && compareSimilarVariables(fgA, fgB, skipsamples=skipsamples, show=show, skip=skiplist )
  TP = TP && compareSimilarFactors(fgA, fgB, skipsamples=skipsamples, skipcompute=skipcompute, show=show )
  TP = TP && compareAll(fgA.solverParams, fgB.solverParams, skip=skiplist)

  return TP
end


text:
using Libdl # Why do we need ???
run(`make`)

max_issues_repo_name: AtelierArith/gomalish

max_forks_repo_name: AtelierArith/gomalish


text:
export vdoc

"""
    vdoc(array_path_bam::Array{String,1}, path_bed12::String;
					output_prefix::String = "")

Calculates Variability of Depth of Coverage (VDoC) score across samples.

A VDoC score is defined for a transcript as follows:
1. For each base on a transcript, a depth of coverage is calculated.
2. For each base, the coefficient of variation (CV) of the depths of coverage is calculated
   across samples.
3. The average of CVs across the transcript is calculated, which we designate as a VDoC score.

# Arguments
---------
- `array_path_bam::Array{String,1}`: Array of paths to BAM files.
- `path_bed12::String`: Path to a reference gene model file with BED12 format.
- `output_prefix::String`: Prefix for output files. If this keyword is not specified or set to
  "" (defalut), no output files are saved.
"""
function vdoc(array_path_bam::Array{String,1}, path_bed12::String;
					output_prefix::String = "")

	if output_prefix != ""
		path_out = output_prefix * ".vdoc.txt"
		path_out_plot = output_prefix * ".vdoc.pdf"
		println(@sprintf "Start calculating absolute gene body coverage...\n- bam: %s\n- bed12: %s\n- output: %s\n          %s\n" join(array_path_bam, ", ") path_bed12 path_out path_out_plot)
	end

	# File check
	for i in 1:length(array_path_bam)
		path_bam = array_path_bam[i]
		if !isfile(path_bam)
			error(@sprintf "No such file: %s\n" path_bam)
		end
	end
	if !isfile(path_bed12)
		error(@sprintf "No such file: %s\n" path_bed12)
	end
	if output_prefix != ""
		if !uv_access_writable(dirname(output_prefix))
			error(@sprintf "Output files are not writable to: %s\n" (dirname(output_prefix) != "" ? dirname(output_prefix) : "."))
		end
	end

	# Load transcripts
	transcripts = load_transcript(path_bed12)
	println(@sprintf "transcripts: %d" length(transcripts))

	# Define containers of results
	array_vdoc = zeros(length(transcripts))
	array_transcript_length = zeros(Int64, length(transcripts))
	array_transcript_name = Array{String,1}(undef, length(transcripts))

	# Caclutate the read coverage
	for (j, t) in enumerate(transcripts)
		# Define matrix
		transcript_length = sum(BED.blocksizes(t))
		coverages = zeros(Int, length(array_path_bam), transcript_length)

		for (i, path_bam) in enumerate(array_path_bam)
			open(BAM.Reader, path_bam, index=path_bam*".bai") do bam_reader
				# Load coverage
				coverages[i,:] = readcoverage_transcript_bam(bam_reader::BAM.Reader, t::BED.Record)
			end
		end
		# Calculate mean of CV
		mean_coverages = mean(coverages, dims=1)
		std_coverages = std(coverages, dims=1)
		cv_coverages = std_coverages ./ mean_coverages
		vdoc = mean(cv_coverages[mean_coverages .> 0])

		array_vdoc[j] = vdoc
		array_transcript_length[j] = transcript_length
		array_transcript_name[j] = BED.name(t)
	end

	if output_prefix != ""
		# Write output (TSV) (Sample_ID, Bin, Coverage)
		println("Write results...")
		open(path_out, "w") do io
			writedlm(io, [array_transcript_name array_transcript_length array_vdoc], '\t')
		end

		# Message
		println(@sprintf "Finished! Check output files:\n- %s" path_out)
	end

    return(array_transcript_name, array_transcript_length, array_vdoc)
end


text:
using IterativeSolvers
using Test

A = Float64[
  16 4 8;
  4 5 -4;
  8 -4 22;
]
a = Float64[16; 18; -22]
x = Float64[1; 1; 1]
jacobi!(x, A, a, maxiter=100)

@test round.(A * x, digits=5) == a

x = Float64[1; 1; 1]
gauss_seidel!(x, A, a, maxiter=100)

@test round.(A * x, digits=5) == a

x = Float64[1; 1; 1]
ω = 1.0
ssor!(x, A, a, ω, maxiter=100)

@test round.(A * x, digits=5) == a


text:
# Computes I_v(k) and I_v(k) vectors from Luger et al. (2018), along
# the derivatives with respect to k using recursion.
include("cel_bulirsch.jl")

#using GSL

function Iv_series(k2::T,v::Int64) where {T <: Real}
# Use series expansion to compute I_v:
nmax = 100
n = 2; error = Inf; if k2 < 1; tol = eps(k2); else; tol = eps(inv(k2)); end
# Computing leading coefficient (n=0):
coeff = 2/(2v+1)
# Add leading term to I_v:
Iv = one(k2)*coeff
# Now, compute higher order terms until desired precision is reached:
while n < nmax && abs(error) > tol
  coeff *= (n-1)*(n+2v-1)
  coeff /= n*(n+2v+1)
  coeff *= k2
  Iv += coeff
#  error = coeff/Iv
  error = coeff
  n += 2
end
return Iv*k2^v*sqrt(k2)
end

# Compute I_v with hypergeometric function (this requires GSL library,
# which can't use BigFloat or ForwardDiff types):
function Iv_hyp(k2::T,v::Int64) where {T <: Real}
a = 0.5*one(k2); b=v+0.5*one(k2); c=v+1.5*one(k2);  fac = 2/(1+2v)
return sqrt(k2)*k2^v*fac*hypergeom([a,b],c,k2)
end

# Compute J_v with hypergeometric function:
function Jv_hyp(k2::T,v::Int64) where {T <: Real}
if k2 < 1
  a = 0.5; b=v+0.5; c=v+3.0;  fac = 3pi/(4*(v+1)*(v+2))
  for i=2:2:2v
    fac *= (i-1)/i
  end
  return sqrt(k2)*k2^v*fac*hypergeom([a,b],c,k2)
else # k^2 >=1
  k2inv = inv(k2)
  # Found a simpler expression than the one in paper (and perhaps more stable for large k^2):
  return  sqrt(pi)*gamma(v+.5)*(sf_hyperg_2F1_renorm(-.5,v+.5,v+1.,k2inv)-(.5+v)*k2inv*sf_hyperg_2F1_renorm(-.5,v+1.5,v+2.,k2inv))
end
end

function Jv_series(k2::T,v::Int64) where {T <: Real}
# Use series expansion to compute J_v:
nmax = 100
n = 2; error = Inf; if k2 < 1; tol = eps(k2); else; tol = eps(inv(k2)); end
# Computing leading coefficient (n=0):
#coeff = 3pi/(2^(2+v)*factorial(v+2))
if k2 < 1
  coeff = 0.75*pi/exp(lfact(v+2))
# multiply by (2v-1)!!
  for i=2:2:2v
    coeff *= (i-1)/2
  end
# Add leading term to J_v:
  Jv = convert(typeof(k2),coeff)
# Now, compute higher order terms until desired precision is reached:
  while n < nmax && abs(error) > tol
    coeff *= (n-1)*(n+2v-1)
    coeff /= n*(n+2v+4)
    coeff *= k2
    Jv += coeff
    error = coeff/Jv
    n += 2
  end
  return Jv*k2^v*sqrt(k2)
else # k^2 >= 1
  coeff = convert(typeof(k2),pi)
  # Compute (2v-1)!!/(2^v v!):
  for i=2:2:2v
    coeff *= (i-1)/i
  end
  Jv = convert(typeof(k2),coeff)
  k2inv = inv(k2)
  while n < nmax && abs(error) > tol
#    coeff *= (1.-2.5/n)*(1.-.5/(n+v))/k2
#    coeff *= (1-5/(n))*(1-1/(n+2v))/k2
    coeff *= (n-5)*(n+2v-1)
    coeff /= (n*(n+2v))
    coeff *= k2inv
    Jv += coeff
    error = coeff/Jv
    n += 2
  end
  return Jv
end
end

function dJv_seriesdk(k2::T,v::Int64) where {T <: Real}
# Use series expansion to compute J_v:
nmax = 100
n = 2; error = Inf; if k2 < 1; tol = eps(k2); else; tol = eps(inv(k2)); end
# Computing leading coefficient (n=0):
#coeff = 3pi/(2^(2+v)*factorial(v+2))
coeff = zero(k2)
if k2 < 1
  coeff = 3pi/(2^(2+v)*exp(lfact(v+2)))
#  println("coefficient: ",coeff)
# multiply by (2v-1)!!
  for i=2:2:2v
    coeff *= i-1
  end
# Add leading term to J_v:
  Jv = one(k2)*coeff
  dJvdk = one(k2)*coeff*(2v+1)
# Now, compute higher order terms until desired precision is reached:
  while n < nmax && abs(error) > tol
    coeff *= (n-1)*(n+2v-1)
    coeff /= n*(n+2v+4)
    coeff *= k2
    Jv += coeff
    dJvdk += coeff*(n+2v+1)
#    error = coeff/Jv
    error = coeff
    n += 2
  end
  dJvdk *= k2^v
  Jv *= k2^v*sqrt(k2)
#  println("Jv: ",Jv," dJv/dk: ",dJvdk)
  return Jv,dJvdk
else # k^2 >= 1
  coeff = convert(typeof(k2),pi)
  # Compute (2v-1)!!/(2^v v!):
  for i=2:2:2v
    coeff *= (i-1)/i
  end
  Jv = one(k2)*coeff
  dJvdk = zero(k2)
  k2inv = inv(k2)
  while n < nmax && abs(error) > tol
#    coeff *= (1-5/n)*(1-1/(n+2v))*k2inv
    coeff *= (n-5)*(n+2v-1)
    coeff /= (n*(n+2v))
    coeff *= k2inv

    Jv += coeff
    dJvdk -= n*coeff
#    error = coeff/Jv
    error = coeff
    n += 2
  end
  dJvdk /= sqrt(k2)
  return Jv,dJvdk
end
end

function IJv_raise!(v_max::Int64,k2::T,kck::T,kc::T,kap::T,Iv::Array{T,1},Jv::Array{T,1})  where {T <: Real}
# This function needs debugging. [ ]
# Compute I_v, J_v for 0 <= v <= v_max = l_max+2
# Define k:
k = sqrt(k2)
# Iterate upwards in v:
v = v_max
# Compute I_v via upward iteration on v:
if k2 < 1
# First, compute value for v=0:
  if k2 < 0.5
    Iv[1] = 2*asin(sqrt(k2))
#  Iv[1] = acos(1.-2k2)
  else
    Iv[1] = 2*acos(kc)
  end
# Try something else:
#  Iv[1] = asin(2*kck)
  Iv[1] = kap
# Next, iterate upwards in v:
#  f0 = kc/k
  f0 = kck
  v = 1
# Loop over v, computing I_v and J_v from higher v:
  while v <= v_max
    Iv[v+1]=((2v-1)*Iv[v]/2-f0)/v
    f0 *= k2
    v += 1
  end
else # k^2 >= 1
  # Compute v=0
  Iv[1] = pi
  for v=1:v_max
    Iv[v+1]=Iv[v]*(1-1/(2v))
  end
end
# Need to compute J_v for v=0, 1:
v= 0
if k2 < 1
  # Use cel_bulirsch:
  if k2 > 0
    Jv[v+1]=2/(3k2*k)*cel_bulirsch(k2,kc,one(k2),k2*(3k2-1),k2*(1-k2))
    Jv[v+2]= 2/(15k2*k)*cel_bulirsch(k2,kc,one(k2),2k2*(3k2-2),k2*(4-7k2+3k2*k2))
  else
    Jv[v+1]= 0.0
    Jv[v+2]= 0.0
  end
else # k^2 >=1
  k2inv = inv(k2)
  Jv[v+1]=2/3*cel_bulirsch(k2inv,kc,one(k2),3-k2inv,3-5k2inv+2k2inv^2)
  Jv[v+2]=cel_bulirsch(k2inv,kc,one(k2),12-8*k2inv,2*(9-8k2inv)*(1-k2inv))/15
end
v=2
while v <= v_max
  Jv[v+1] = (2*(v+1+(v-1)*k2)*Jv[v]-k2*(2v-3)*Jv[v-1])/(2v+3)
  v += 1
end
return
end

function dIJv_raise_dk!(v_max::Int64,k2::T,kck::T,kc::T,kap::T,Iv::Array{T,1},Jv::Array{T,1},dIvdk::Array{T,1},dJvdk::Array{T,1})  where {T <: Real}
# Compute I_v, J_v for 0 <= v <= v_max = l_max+2
# Define k:
k = sqrt(k2)
# Iterate upwards in v:
v = v_max
# Compute I_v via upward iteration on v:
if k2 < 1
# First, compute value for v=0:
  if k2 < 0.5
    Iv[1] = 2*asin(sqrt(k2))
  else
    Iv[1] = 2*acos(kc)
#    Iv[1] = 2*convert(Float64,acos(big(kc)))
#    Iv[1] = convert(Float64,acos((big(b)^2+big(r)^2-big(1))/(big(2)*big(b)*big(r))))
  end
  Iv[1] = kap
# Try something else:
# Next, iterate upwards in v:
  f0 = kck
  v = 1
# Loop over v, computing I_v and J_v from higher v:
  while v <= v_max
    Iv[v+1]=((2v-1)*Iv[v]/2-f0)/v
    f0 *= k2
    v += 1
  end
  # Now compute compute derivatives:
  dIvdk[1] = 2/kc
#  dIvdk[1] = 2*k2/kck
  for v=1:v_max
    dIvdk[v+1] = k2*dIvdk[v]
  end
else # k^2 >= 1
  # Compute v=0
  Iv[1] = pi
  for v=1:v_max
    Iv[v+1]=Iv[v]*(1-1/(2v))
  end
  # Derivatives of I_v are zero:
  fill!(dIvdk,zero(k2))
end
# Need to compute J_v for v=0, 1:
v= 0
if k2 < 1
  # Use cel_bulirsch:
  if k2 > 0
    Jv[v+1]=2/(3k2*k)*cel_bulirsch(k2,kc,one(k2),k2*(3k2-1),k2*(1-k2))
    Jv[v+2]= 2/(15k2*k)*cel_bulirsch(k2,kc,one(k2),2k2*(3k2-2),k2*(4-7k2+3k2*k2))
    dJvdk[v+1] = 2*cel_bulirsch(k2,kc,one(k2),inv(k2),(1-inv(k2)))
    dJvdk[v+2] = -3*Jv[v+2]/k+k2*dJvdk[v+1]
  else
    Jv[v+1]= 0.0
    Jv[v+2]= 0.0
    dJvdk[v+1] = 0.0
    dJvdk[v+2] = 0.0
  end
else # k^2 >=1
  k2inv = inv(k2)
  Jv[v+1] = 2/3*cel_bulirsch(k2inv,kc,one(k2),3-k2inv,3-5k2inv+2k2inv^2)
  Jv[v+2] = cel_bulirsch(k2inv,kc,one(k2),12-8*k2inv,2*(9-8k2inv)*(1-k2inv))/15
  dJvdk[v+1] = 2/(k2*k)*cel_bulirsch(k2inv,kc,one(k2),one(k2),2*(1-k2inv))
  dJvdk[v+2] = -3*Jv[v+2]/k+k2*dJvdk[v+1]
end
v=2
while v <= v_max
  Jv[v+1] = (2*(v+1+(v-1)*k2)*Jv[v]-k2*(2v-3)*Jv[v-1])/(2v+3)
  dJvdk[v+1] = -3*Jv[v+1]/k+k2*dJvdk[v]
  v += 1
end
return
end

function IJv_lower!(v_max::Int64,k2::T,kck::T,kc::T,kap::T,Iv::Array{T,1},Jv::Array{T,1})  where {T <: Real}
# Compute I_v, J_v for 0 <= v <= v_max = l_max+2
# Define k:
k = sqrt(k2)
# Iterate downwards in v:
v = v_max
# Add in k2 > 1 cases [ ]
# First, compute approximation for large v:
#Iv[v+1]=Iv_hyp(k2,v)
if k2 < 1
  Iv[v+1]=Iv_series(k2,v)
# Next, iterate downwards in v:
  f0 = k2^(v-1)*kck
# Loop over v, computing I_v and J_v from higher v:
  while v >= 2
    Iv[v] = 2/(2v-1)*(v*Iv[v+1]+f0)
    f0 /= k2
    v -= 1
  end
  Iv[1] = kap
else # k^2 >= 1
  # Compute v=0 (no need to iterate downwards in this case):
  Iv[1] = pi
  for v=1:v_max
    Iv[v+1]=Iv[v]*(1-.5/v)
  end
end
v= v_max
# Need to compute top two for J_v:
#if typeof(k2) == BigFloat
  Jv[v]=Jv_series(k2,v-1); Jv[v+1]=Jv_series(k2,v)
#else
#  Jv[v]=Jv_hyp(k2,v-1); Jv[v+1]=Jv_hyp(k2,v)
#end
#if typeof(k2) == Float64
#  println("v ",v," k2 ",k2," Jv_ser ",convert(Float64,Jv_series(big(k2),v-1))," ",convert(Float64,Jv_series(big(k2),v))," Jv_hyp ",Jv_hyp(k2,v-1)," ",Jv_hyp(k2,v))
#end
# Iterate downwards in v (lower):
while v >= 2
  f2 = k2*(2v-3); f1 = 2*(v+1+(v-1)*k2)/f2; f3 = (2v+3)/f2
#  k2inv = inv(k2)
#  f2 = (2v-3); f1 = 2*((v+1)*k2inv+(v-1))/f2; f3 = (2v+3)/f2*k2inv
  Jv[v-1] = f1*Jv[v]-f3*Jv[v+1]
  v -= 1
end
# Compute first two exactly:
v= 0
if k2 < 1
#  # Use cel_bulirsch:
#  if k2 > 0
#    Jv[v+1]=2/(3*k2*k)*cel_bulirsch(k2,kc,one(k2),k2*(3k2-1),k2*(1-k2))
#    Jv[v+2]=2/(15*k2*k)*cel_bulirsch(k2,kc,one(k2),2*k2*(3k2-2),k2*(4-7k2+3k2*k2))
##    Jv[v+3] =-2/(35*k)*cel_bulirsch(k2,kc,one(k2),8-11k2+k2*k2,(-8+5k2+2k2*k2)*(1-k2))
#  else
#    Jv[v+1]= 0.0
#    Jv[v+2]= 0.0
#  end
else # k^2 >=1
  k2inv = inv(k2)
  Jv[v+1]=2/3*cel_bulirsch(k2inv,kc,one(k2),3-k2inv,3-5k2inv+2k2inv^2)
  Jv[v+2]=cel_bulirsch(k2inv,kc,one(k2),12-8*k2inv,2*(9-8k2inv)*(1-k2inv))/15
#  Jv[v+3]=2/(35*k2)*cel_bulirsch(k2inv,kc,one(k2),-k2^2+11k2-8,(k2^2+16*k2-16)*(1-k2inv))
end
return
end

function dIJv_lower_dk!(v_max::Int64,k2::T,kck::T,kc::T,kap::T,Iv::Array{T,1},Jv::Array{T,1},dIvdk::Array{T,1},dJvdk::Array{T,1})  where {T <: Real}
# Compute I_v, J_v for 0 <= v <= v_max = l_max+2
# Define k:
k = sqrt(k2)
# Iterate downwards in v:
v = v_max
# Add in k2 > 1 cases [ ]
# First, compute approximation for large v:
#Iv[v+1]=Iv_hyp(k2,v)
if k2 < 1
  Iv[v+1]=Iv_series(k2,v)
# Next, iterate downwards in v:
  f0 = k2^(v-1)*kck
# Loop over v, computing I_v and J_v from higher v:
  while v >= 2
    Iv[v] = 2/(2v-1)*(v*Iv[v+1]+f0)
    f0 /= k2
    v -= 1
  end
  Iv[1] = kap
  # Now compute compute derivatives:
  dIvdk[1] = 2/kc
  for v=1:v_max
    dIvdk[v+1] = k2*dIvdk[v]
  end
else # k^2 >= 1
  # Compute v=0 (no need to iterate downwards in this case):
  Iv[1] = pi
  for v=1:v_max
    Iv[v+1]=Iv[v]*(1-.5/v)
  end
  # Derivatives of I_v are zero:
  fill!(dIvdk,zero(k2))
end
v= v_max
# Need to compute top two for J_v:
#if typeof(k2) == BigFloat
if Jv[v+1] == 0.0
  while Jv[v+1] == 0.0  # Loop downward in v until we get a non-zero Jv[v]
    dJvdk0 = zero(typeof(k2)); dJvdk1 = zero(typeof(k2))
    Jv[v],dJvdk0 = dJv_seriesdk(k2,v-1); Jv[v+1],dJvdk1=dJv_seriesdk(k2,v)
    dJvdk[v] = dJvdk0; dJvdk[v+1] = dJvdk1
#  println("Jv: ",Jv[v]," Jv+1: ",Jv[v+1])
    v -=1
  end
  v +=1
end
#else
#  Jv[v]=Jv_hyp(k2,v-1); Jv[v+1]=Jv_hyp(k2,v)
#end
#if typeof(k2) == Float64
#  println("v ",v," k2 ",k2," Jv_ser ",convert(Float64,Jv_series(big(k2),v-1))," ",convert(Float64,Jv_series(big(k2),v))," Jv_hyp ",Jv_hyp(k2,v-1)," ",Jv_hyp(k2,v))
#end
# Iterate downwards in v (lower):
while v >= 2
  f2 = k2*(2v-3); f1 = 2*(v+1+(v-1)*k2)/f2; f3 = (2v+3)/f2
#  k2inv = inv(k2)
#  f2 = (2v-3); f1 = 2*((v+1)*k2inv+(v-1))/f2; f3 = (2v+3)/f2*k2inv
  Jv[v-1] = f1*Jv[v]-f3*Jv[v+1]
#  dJvdk[v-1] = (dJvdk[v]+3/k*Jv[v])*k2inv
  dJvdk[v-1] = (dJvdk[v]+3/k*Jv[v])/k2
  v -= 1
end
# Compute first two exactly:
v= 0
if k2 < 1
#  # Use cel_bulirsch:
#  if k2 > 0
#    Jv[v+1]=2/(3*k2*k)*cel_bulirsch(k2,kc,one(k2),k2*(3k2-1),k2*(1-k2))
#    Jv[v+2]=2/(15*k2*k)*cel_bulirsch(k2,kc,one(k2),2*k2*(3k2-2),k2*(4-7k2+3k2*k2))
##    Jv[v+3] =-2/(35*k)*cel_bulirsch(k2,kc,one(k2),8-11k2+k2*k2,(-8+5k2+2k2*k2)*(1-k2))
#  else
#    Jv[v+1]= 0.0
#    Jv[v+2]= 0.0
#  end
else # k^2 >=1
  k2inv = inv(k2)
  Jv[v+1]=2/3*cel_bulirsch(k2inv,kc,one(k2),3-k2inv,3-5k2inv+2k2inv^2)
  Jv[v+2]=cel_bulirsch(k2inv,kc,one(k2),12-8*k2inv,2*(9-8k2inv)*(1-k2inv))/15
#  Jv[v+3]=2/(35*k2)*cel_bulirsch(k2inv,kc,one(k2),-k2^2+11k2-8,(k2^2+16*k2-16)*(1-k2inv))
end
return
end

max_stars_repo_name: tagordon/limbdark

max_issues_repo_name: tagordon/limbdark

max_forks_repo_name: tagordon/limbdark