text:
"""
    Interpolation{dim, ref_shape, order}()

Return an `Interpolation` of given dimension `dim`, reference shape
(see see [`AbstractRefShape`](@ref)) `ref_shape` and order `order`.
`order` corresponds to the highest order term in the polynomial.
The interpolation is used to define shape functions to interpolate
a function between nodes.

The following interpolations are implemented:

* `Lagrange{1,RefCube,1}`
* `Lagrange{1,RefCube,2}`
* `Lagrange{2,RefCube,1}`
* `Lagrange{2,RefCube,2}`
* `Lagrange{2,RefTetrahedron,1}`
* `Lagrange{2,RefTetrahedron,2}`
* `Lagrange{3,RefCube,1}`
* `Serendipity{2,RefCube,2}`
* `Serendipity{3,RefCube,2}`
* `Lagrange{3,RefTetrahedron,1}`
* `Lagrange{3,RefTetrahedron,2}`

# Examples
```jldoctest
julia> ip = Lagrange{2,RefTetrahedron,2}()
Ferrite.Lagrange{2,Ferrite.RefTetrahedron,2}()

julia> getnbasefunctions(ip)
6
```
"""
abstract type Interpolation{dim,shape,order} end

"""
Return the dimension of an `Interpolation`
"""
@inline getdim(ip::Interpolation{dim}) where {dim} = dim

"""
Return the reference shape of an `Interpolation`
"""
@inline getrefshape(ip::Interpolation{dim,shape}) where {dim,shape} = shape

"""
Return the polynomial order of the `Interpolation`
"""
@inline getorder(ip::Interpolation{dim,shape,order}) where {dim,shape,order} = order

"""
Compute the value of the shape functions at a point ξ for a given interpolation
"""
function value(ip::Interpolation{dim}, ξ::Vec{dim,T}) where {dim,T}
    [value(ip, i, ξ) for i in 1:getnbasefunctions(ip)]
end

"""
Compute the gradients of the shape functions at a point ξ for a given interpolation
"""
function derivative(ip::Interpolation{dim}, ξ::Vec{dim,T}) where {dim,T}
    [gradient(ξ -> value(ip, i, ξ), ξ) for i in 1:getnbasefunctions(ip)]
end

#####################
# Utility functions #
#####################

"""
Return the number of base functions for an [`Interpolation`](@ref) or `Values` object.
"""
getnbasefunctions

# struct that gathers all the information needed to distribute
# dofs for a given interpolation.
struct InterpolationInfo
    # TODO: Can be smaller than `Int` if that matters...
    nvertexdofs::Int
    nedgedofs::Int
    nfacedofs::Int
    ncelldofs::Int
    dim::Int
    InterpolationInfo(interpolation::Interpolation{dim}) where {dim} =
        new(nvertexdofs(interpolation), nedgedofs(interpolation),
            nfacedofs(interpolation),   ncelldofs(interpolation), dim)
end

# The following functions are used to distribute the dofs. Definitions:
#   vertexdof: dof on a "corner" of the reference shape
#   facedof: dof in the dim-1 dimension (line in 2D, surface in 3D)
#   edgedof: dof on a line between 2 vertices (i.e. "corners") (3D only)
#   celldof: dof that is local to the element

# Fallbacks for the interpolations which are used to distribute the dofs correctly
nvertexdofs(::Interpolation) = 0
nedgedofs(::Interpolation)   = 0
nfacedofs(::Interpolation)   = 0
ncelldofs(::Interpolation)   = 0

# Needed for distrubuting dofs on shells correctly (face in 2d is edge in 3d)
edges(ip::Interpolation{2}) = faces(ip)
nedgedofs(ip::Interpolation{2}) = nfacedofs(ip)

# Fallbacks for vertices
vertices(::Interpolation{2,RefCube}) = (1,2,3,4)
vertices(::Interpolation{3,RefCube}) = (1,2,3,4,5,6,7,8)
vertices(::Interpolation{2,RefTetrahedron}) = (1,2,3)
vertices(::Interpolation{3,RefTetrahedron}) = (1,2,3,4)

############
# Lagrange #
############
struct Lagrange{dim,shape,order} <: Interpolation{dim,shape,order} end

getlowerdim(::Lagrange{dim,shape,order}) where {dim,shape,order} = Lagrange{dim-1,shape,order}()
getlowerorder(::Lagrange{dim,shape,order}) where {dim,shape,order} = Lagrange{dim,shape,order-1}()

##################################
# Lagrange dim 1 RefCube order 1 #
##################################
getnbasefunctions(::Lagrange{1,RefCube,1}) = 2
nvertexdofs(::Lagrange{1,RefCube,1}) = 1

faces(::Lagrange{1,RefCube,1}) = ((1,), (2,))

function reference_coordinates(::Lagrange{1,RefCube,1})
    return [Vec{1, Float64}((-1.0,)),
            Vec{1, Float64}(( 1.0,))]
end

function value(ip::Lagrange{1,RefCube,1}, i::Int, ξ::Vec{1})
    ξ_x = ξ[1]
    i == 1 && return (1 - ξ_x) * 0.5
    i == 2 && return (1 + ξ_x) * 0.5
    throw(ArgumentError("no shape function $i for interpolation $ip"))
end

##################################
# Lagrange dim 1 RefCube order 2 #
##################################
getnbasefunctions(::Lagrange{1,RefCube,2}) = 3
nvertexdofs(::Lagrange{1,RefCube,2}) = 1
ncelldofs(::Lagrange{1,RefCube,2}) = 1

faces(::Lagrange{1,RefCube,2}) = ((1,), (2,))

function reference_coordinates(::Lagrange{1,RefCube,2})
    return [Vec{1, Float64}((-1.0,)),
            Vec{1, Float64}(( 1.0,)),
            Vec{1, Float64}(( 0.0,))]
end

function value(ip::Lagrange{1,RefCube,2}, i::Int, ξ::Vec{1})
    ξ_x = ξ[1]
    i == 1 && return ξ_x * (ξ_x - 1) * 0.5
    i == 2 && return ξ_x * (ξ_x + 1) * 0.5
    i == 3 && return 1 - ξ_x^2
    throw(ArgumentError("no shape function $i for interpolation $ip"))
end

##################################
# Lagrange dim 2 RefCube order 1 #
##################################
getnbasefunctions(::Lagrange{2,RefCube,1}) = 4
nvertexdofs(::Lagrange{2,RefCube,1}) = 1

faces(::Lagrange{2,RefCube,1}) = ((1,2), (2,3), (3,4), (4,1))

function reference_coordinates(::Lagrange{2,RefCube,1})
    return [Vec{2, Float64}((-1.0, -1.0)),
            Vec{2, Float64}(( 1.0, -1.0)),
            Vec{2, Float64}(( 1.0,  1.0,)),
            Vec{2, Float64}((-1.0,  1.0,))]
end

function value(ip::Lagrange{2,RefCube,1}, i::Int, ξ::Vec{2})
    ξ_x = ξ[1]
    ξ_y = ξ[2]
    i == 1 && return (1 - ξ_x) * (1 - ξ_y) * 0.25
    i == 2 && return (1 + ξ_x) * (1 - ξ_y) * 0.25
    i == 3 && return (1 + ξ_x) * (1 + ξ_y) * 0.25
    i == 4 && return (1 - ξ_x) * (1 + ξ_y) * 0.25
    throw(ArgumentError("no shape function $i for interpolation $ip"))
end

##################################
# Lagrange dim 2 RefCube order 2 #
##################################
getnbasefunctions(::Lagrange{2,RefCube,2}) = 9
nvertexdofs(::Lagrange{2,RefCube,2}) = 1
nfacedofs(::Lagrange{2,RefCube,2}) = 1
ncelldofs(::Lagrange{2,RefCube,2}) = 1

faces(::Lagrange{2,RefCube,2}) = ((1,2,5), (2,3,6), (3,4,7), (4,1,8))

function reference_coordinates(::Lagrange{2,RefCube,2})
    return [Vec{2, Float64}((-1.0, -1.0)),
            Vec{2, Float64}(( 1.0, -1.0)),
            Vec{2, Float64}(( 1.0,  1.0)),
            Vec{2, Float64}((-1.0,  1.0)),
            Vec{2, Float64}(( 0.0, -1.0)),
            Vec{2, Float64}(( 1.0,  0.0)),
            Vec{2, Float64}(( 0.0,  1.0)),
            Vec{2, Float64}((-1.0,  0.0)),
            Vec{2, Float64}(( 0.0,  0.0))]
end

function value(ip::Lagrange{2,RefCube,2}, i::Int, ξ::Vec{2})
    ξ_x = ξ[1]
    ξ_y = ξ[2]
    i == 1 && return (ξ_x^2 - ξ_x) * (ξ_y^2 - ξ_y) * 0.25
    i == 2 && return (ξ_x^2 + ξ_x) * (ξ_y^2 - ξ_y) * 0.25
    i == 3 && return (ξ_x^2 + ξ_x) * (ξ_y^2 + ξ_y) * 0.25
    i == 4 && return (ξ_x^2 - ξ_x) * (ξ_y^2 + ξ_y) * 0.25
    i == 5 && return (1 - ξ_x^2) * (ξ_y^2 - ξ_y) * 0.5
    i == 6 && return (ξ_x^2 + ξ_x) * (1 - ξ_y^2) * 0.5
    i == 7 && return (1 - ξ_x^2) * (ξ_y^2 + ξ_y) * 0.5
    i == 8 && return (ξ_x^2 - ξ_x) * (1 - ξ_y^2) * 0.5
    i == 9 && return (1 - ξ_x^2) * (1 - ξ_y^2)
    throw(ArgumentError("no shape function $i for interpolation $ip"))
end

#########################################
# Lagrange dim 2 RefTetrahedron order 1 #
#########################################
getnbasefunctions(::Lagrange{2,RefTetrahedron,1}) = 3
getlowerdim(::Lagrange{2, RefTetrahedron, order}) where {order} = Lagrange{1, RefCube, order}()
nvertexdofs(::Lagrange{2,RefTetrahedron,1}) = 1

vertices(::Lagrange{2,RefTetrahedron,1}) = (1,2,3)
faces(::Lagrange{2,RefTetrahedron,1}) = ((1,2), (2,3), (3,1))

function reference_coordinates(::Lagrange{2,RefTetrahedron,1})
    return [Vec{2, Float64}((1.0, 0.0)),
            Vec{2, Float64}((0.0, 1.0)),
            Vec{2, Float64}((0.0, 0.0))]
end

function value(ip::Lagrange{2,RefTetrahedron,1}, i::Int, ξ::Vec{2})
    ξ_x = ξ[1]
    ξ_y = ξ[2]
    i == 1 && return ξ_x
    i == 2 && return ξ_y
    i == 3 && return 1. - ξ_x - ξ_y
    throw(ArgumentError("no shape function $i for interpolation $ip"))
end

#########################################
# Lagrange dim 2 RefTetrahedron order 2 #
#########################################
getnbasefunctions(::Lagrange{2,RefTetrahedron,2}) = 6
nvertexdofs(::Lagrange{2,RefTetrahedron,2}) = 1
nfacedofs(::Lagrange{2,RefTetrahedron,2}) = 1

vertices(::Lagrange{2,RefTetrahedron,2}) = (1,2,3)
faces(::Lagrange{2,RefTetrahedron,2}) = ((1,2,4), (2,3,5), (3,1,6))

function reference_coordinates(::Lagrange{2,RefTetrahedron,2})
    return [Vec{2, Float64}((1.0, 0.0)),
            Vec{2, Float64}((0.0, 1.0)),
            Vec{2, Float64}((0.0, 0.0)),
            Vec{2, Float64}((0.5, 0.5)),
            Vec{2, Float64}((0.0, 0.5)),
            Vec{2, Float64}((0.5, 0.0))]
end

function value(ip::Lagrange{2,RefTetrahedron,2}, i::Int, ξ::Vec{2})
    ξ_x = ξ[1]
    ξ_y = ξ[2]
    γ = 1. - ξ_x - ξ_y
    i == 1 && return ξ_x * (2ξ_x - 1)
    i == 2 && return ξ_y * (2ξ_y - 1)
    i == 3 && return γ * (2γ - 1)
    i == 4 && return 4ξ_x * ξ_y
    i == 5 && return 4ξ_y * γ
    i == 6 && return 4ξ_x * γ
    throw(ArgumentError("no shape function $i for interpolation $ip"))
end

#########################################
# Lagrange dim 3 RefTetrahedron order 1 #
#########################################
getnbasefunctions(::Lagrange{3,RefTetrahedron,1}) = 4
nvertexdofs(::Lagrange{3,RefTetrahedron,1}) = 1

faces(::Lagrange{3,RefTetrahedron,1}) = ((1,2,3), (1,2,4), (2,3,4), (1,4,3))
edges(::Lagrange{3,RefTetrahedron,1}) = ((1,2), (2,3), (3,1), (1,4), (2,4), (3,4))

function reference_coordinates(::Lagrange{3,RefTetrahedron,1})
    return [Vec{3, Float64}((0.0, 0.0, 0.0)),
            Vec{3, Float64}((1.0, 0.0, 0.0)),
            Vec{3, Float64}((0.0, 1.0, 0.0)),
            Vec{3, Float64}((0.0, 0.0, 1.0))]
end

function value(ip::Lagrange{3,RefTetrahedron,1}, i::Int, ξ::Vec{3})
    ξ_x = ξ[1]
    ξ_y = ξ[2]
    ξ_z = ξ[3]
    i == 1 && return 1.0 - ξ_x - ξ_y - ξ_z
    i == 2 && return ξ_x
    i == 3 && return ξ_y
    i == 4 && return ξ_z
    throw(ArgumentError("no shape function $i for interpolation $ip"))
end

#########################################
# Lagrange dim 3 RefTetrahedron order 2 #
#########################################
getnbasefunctions(::Lagrange{3,RefTetrahedron,2}) = 10
nvertexdofs(::Lagrange{3,RefTetrahedron,2}) = 1
nedgedofs(::Lagrange{3,RefTetrahedron,2}) = 1

faces(::Lagrange{3,RefTetrahedron,2}) = ((1,2,3,5,6,7), (1,2,4,5,9,8), (2,3,4,6,10,9), (1,4,3,8,10,7))
edges(::Lagrange{3,RefTetrahedron,2}) = ((1,5,2), (2,6,3), (3,7,1), (1,8,4), (2,9,4), (3,10,4))

function reference_coordinates(::Lagrange{3,RefTetrahedron,2})
    return [Vec{3, Float64}((0.0, 0.0, 0.0)),
            Vec{3, Float64}((1.0, 0.0, 0.0)),
            Vec{3, Float64}((0.0, 1.0, 0.0)),
            Vec{3, Float64}((0.0, 0.0, 1.0)),
            Vec{3, Float64}((0.5, 0.0, 0.0)),
            Vec{3, Float64}((0.5, 0.5, 0.0)),
            Vec{3, Float64}((0.0, 0.5, 0.0)),
            Vec{3, Float64}((0.0, 0.0, 0.5)),
            Vec{3, Float64}((0.5, 0.0, 0.5)),
            Vec{3, Float64}((0.0, 0.5, 0.5))]
end

# http://www.colorado.edu/engineering/CAS/courses.d/AFEM.d/AFEM.Ch09.d/AFEM.Ch09.pdf
# http://www.colorado.edu/engineering/CAS/courses.d/AFEM.d/AFEM.Ch10.d/AFEM.Ch10.pdf
function value(ip::Lagrange{3,RefTetrahedron,2}, i::Int, ξ::Vec{3})
    ξ_x = ξ[1]
    ξ_y = ξ[2]
    ξ_z = ξ[3]
    i == 1  && return (-2 * ξ_x - 2 * ξ_y - 2 * ξ_z + 1) * (-ξ_x - ξ_y - ξ_z + 1)
    i == 2  && return ξ_x * (2 * ξ_x - 1)
    i == 3  && return ξ_y * (2 * ξ_y - 1)
    i == 4  && return ξ_z * (2 * ξ_z - 1)
    i == 5  && return ξ_x * (-4 * ξ_x - 4 * ξ_y - 4 * ξ_z + 4)
    i == 6  && return 4 * ξ_x * ξ_y
    i == 7  && return 4 * ξ_y * (-ξ_x - ξ_y - ξ_z + 1)
    i == 8  && return ξ_z * (-4 * ξ_x - 4 * ξ_y - 4 * ξ_z + 4)
    i == 9  && return 4 * ξ_x * ξ_z
    i == 10 && return 4 * ξ_y * ξ_z
    throw(ArgumentError("no shape function $i for interpolation $ip"))
end

##################################
# Lagrange dim 3 RefCube order 1 #
##################################
getnbasefunctions(::Lagrange{3,RefCube,1}) = 8
nvertexdofs(::Lagrange{3,RefCube,1}) = 1

faces(::Lagrange{3,RefCube,1}) = ((1,4,3,2), (1,2,6,5), (2,3,7,6), (3,4,8,7), (1,5,8,4), (5,6,7,8))
edges(::Lagrange{3,RefCube,1}) = ((1,2), (2,3), (3,4), (4,1), (1,5), (2,6), (3,7), (4,8), (5,6), (6,7), (7,8), (8,5))

function reference_coordinates(::Lagrange{3,RefCube,1})
    return [Vec{3, Float64}((-1.0, -1.0, -1.0)),
            Vec{3, Float64}(( 1.0, -1.0, -1.0)),
            Vec{3, Float64}(( 1.0,  1.0, -1.0)),
            Vec{3, Float64}((-1.0,  1.0, -1.0)),
            Vec{3, Float64}((-1.0, -1.0,  1.0)),
            Vec{3, Float64}(( 1.0, -1.0,  1.0)),
            Vec{3, Float64}(( 1.0,  1.0,  1.0)),
            Vec{3, Float64}((-1.0,  1.0,  1.0))]
end

function value(ip::Lagrange{3,RefCube,1}, i::Int, ξ::Vec{3})
    ξ_x = ξ[1]
    ξ_y = ξ[2]
    ξ_z = ξ[3]
    i == 1 && return 0.125(1 - ξ_x) * (1 - ξ_y) * (1 - ξ_z)
    i == 2 && return 0.125(1 + ξ_x) * (1 - ξ_y) * (1 - ξ_z)
    i == 3 && return 0.125(1 + ξ_x) * (1 + ξ_y) * (1 - ξ_z)
    i == 4 && return 0.125(1 - ξ_x) * (1 + ξ_y) * (1 - ξ_z)
    i == 5 && return 0.125(1 - ξ_x) * (1 - ξ_y) * (1 + ξ_z)
    i == 6 && return 0.125(1 + ξ_x) * (1 - ξ_y) * (1 + ξ_z)
    i == 7 && return 0.125(1 + ξ_x) * (1 + ξ_y) * (1 + ξ_z)
    i == 8 && return 0.125(1 - ξ_x) * (1 + ξ_y) * (1 + ξ_z)
    throw(ArgumentError("no shape function $i for interpolation $ip"))
end

###############
# Serendipity #
###############
struct Serendipity{dim,shape,order} <: Interpolation{dim,shape,order} end

#####################################
# Serendipity dim 2 RefCube order 2 #
#####################################
getnbasefunctions(::Serendipity{2,RefCube,2}) = 8
getlowerdim(::Serendipity{2,RefCube,2}) = Lagrange{1,RefCube,2}()
getlowerorder(::Serendipity{2,RefCube,2}) = Lagrange{2,RefCube,1}()
nvertexdofs(::Serendipity{2,RefCube,2}) = 1
nfacedofs(::Serendipity{2,RefCube,2}) = 1

faces(::Serendipity{2,RefCube,2}) = ((1,2,5), (2,3,6), (3,4,7), (4,1,8))

function reference_coordinates(::Serendipity{2,RefCube,2})
    return [Vec{2, Float64}((-1.0, -1.0)),
            Vec{2, Float64}(( 1.0, -1.0)),
            Vec{2, Float64}(( 1.0,  1.0)),
            Vec{2, Float64}((-1.0,  1.0)),
            Vec{2, Float64}(( 0.0, -1.0)),
            Vec{2, Float64}(( 1.0,  0.0)),
            Vec{2, Float64}(( 0.0,  1.0)),
            Vec{2, Float64}((-1.0,  0.0))]
end

function value(ip::Serendipity{2,RefCube,2}, i::Int, ξ::Vec{2})
    ξ_x = ξ[1]
    ξ_y = ξ[2]
    i == 1 && return (1 - ξ_x) * (1 - ξ_y) * 0.25(-ξ_x - ξ_y - 1)
    i == 2 && return (1 + ξ_x) * (1 - ξ_y) * 0.25( ξ_x - ξ_y - 1)
    i == 3 && return (1 + ξ_x) * (1 + ξ_y) * 0.25( ξ_x + ξ_y - 1)
    i == 4 && return (1 - ξ_x) * (1 + ξ_y) * 0.25(-ξ_x + ξ_y - 1)
    i == 5 && return 0.5(1 - ξ_x * ξ_x) * (1 - ξ_y)
    i == 6 && return 0.5(1 + ξ_x) * (1 - ξ_y * ξ_y)
    i == 7 && return 0.5(1 - ξ_x * ξ_x) * (1 + ξ_y)
    i == 8 && return 0.5(1 - ξ_x) * (1 - ξ_y * ξ_y)
    throw(ArgumentError("no shape function $i for interpolation $ip"))
end

##################################
# Serendipity dim 3 RefCube order 2 #
##################################
getnbasefunctions(::Serendipity{3,RefCube,2}) = 20
getlowerdim(::Serendipity{3,RefCube,2}) = Serendipity{2,RefCube,2}()
getlowerorder(::Serendipity{3,RefCube,2}) = Lagrange{3,RefCube,1}()
nvertexdofs(::Serendipity{3,RefCube,2}) = 1
nedgedofs(::Serendipity{3,RefCube,2}) = 1

faces(::Serendipity{3,RefCube,2}) = ((1,4,3,2,12,11,10,9), (1,2,6,5,9,18,13,17), (2,3,7,6,10,19,14,18), (3,4,8,7,11,20,15,19), (1,5,8,4,17,16,20,12), (5,6,7,8,13,14,15,16))

function reference_coordinates(::Serendipity{3,RefCube,2})
    return [Vec{3, Float64}((-1.0, -1.0, -1.0)),
            Vec{3, Float64}(( 1.0, -1.0, -1.0)),
            Vec{3, Float64}(( 1.0,  1.0, -1.0)),
            Vec{3, Float64}((-1.0,  1.0, -1.0)),
            Vec{3, Float64}((-1.0, -1.0,  1.0)),
            Vec{3, Float64}(( 1.0, -1.0,  1.0)),
            Vec{3, Float64}(( 1.0,  1.0,  1.0)),
            Vec{3, Float64}((-1.0,  1.0,  1.0)),
            Vec{3, Float64}((0.0, -1.0, -1.0)),
            Vec{3, Float64}((1.0, 0.0, -1.0)),
            Vec{3, Float64}((0.0, 1.0, -1.0)),
            Vec{3, Float64}((-1.0, 0.0, -1.0)),
            Vec{3, Float64}((0.0, -1.0, 1.0)),
            Vec{3, Float64}((1.0, 0.0, 1.0)),
            Vec{3, Float64}((0.0, 1.0, 1.0)),
            Vec{3, Float64}((-1.0, 0.0, 1.0)),
            Vec{3, Float64}((-1.0, -1.0, 0.0)),
            Vec{3, Float64}((1.0, -1.0, 0.0)),
            Vec{3, Float64}((1.0, 1.0, 0.0)),
            Vec{3, Float64}((-1.0, 1.0, 0.0)),]
end

function value(ip::Serendipity{3,RefCube,2}, i::Int, ξ::Vec{3})
    ξ_x = ξ[1]
    ξ_y = ξ[2]
    ξ_z = ξ[3]
    i == 1 && return 0.125(1 - ξ_x) * (1 - ξ_y) * (1 - ξ_z) - 0.5(value(ip,12,ξ) + value(ip,9,ξ) + value(ip,17,ξ))
    i == 2 && return 0.125(1 + ξ_x) * (1 - ξ_y) * (1 - ξ_z) - 0.5(value(ip,9,ξ) + value(ip,10,ξ) + value(ip,18,ξ))
    i == 3 && return 0.125(1 + ξ_x) * (1 + ξ_y) * (1 - ξ_z) - 0.5(value(ip,10,ξ) + value(ip,11,ξ) + value(ip,19,ξ))
    i == 4 && return 0.125(1 - ξ_x) * (1 + ξ_y) * (1 - ξ_z) - 0.5(value(ip,11,ξ) + value(ip,12,ξ) + value(ip,20,ξ))
    i == 5 && return 0.125(1 - ξ_x) * (1 - ξ_y) * (1 + ξ_z) - 0.5(value(ip,16,ξ) + value(ip,13,ξ) + value(ip,17,ξ))
    i == 6 && return 0.125(1 + ξ_x) * (1 - ξ_y) * (1 + ξ_z) - 0.5(value(ip,13,ξ) + value(ip,14,ξ) + value(ip,18,ξ))
    i == 7 && return 0.125(1 + ξ_x) * (1 + ξ_y) * (1 + ξ_z) - 0.5(value(ip,14,ξ) + value(ip,15,ξ) + value(ip,19,ξ))
    i == 8 && return 0.125(1 - ξ_x) * (1 + ξ_y) * (1 + ξ_z) - 0.5(value(ip,15,ξ) + value(ip,16,ξ) + value(ip,20,ξ))
    i == 9 && return 0.25(1 - ξ_x^2) * (1 - ξ_y) * (1 - ξ_z)
    i == 10 && return 0.25(1 + ξ_x) * (1 - ξ_y^2) * (1 - ξ_z)
    i == 11 && return 0.25(1 - ξ_x^2) * (1 + ξ_y) * (1 - ξ_z)
    i == 12 && return 0.25(1 - ξ_x) * (1 - ξ_y^2) * (1 - ξ_z)
    i == 13 && return 0.25(1 - ξ_x^2) * (1 - ξ_y) * (1 + ξ_z)
    i == 14 && return 0.25(1 + ξ_x) * (1 - ξ_y^2) * (1 + ξ_z)
    i == 15 && return 0.25(1 - ξ_x^2) * (1 + ξ_y) * (1 + ξ_z)
    i == 16 && return 0.25(1 - ξ_x) * (1 - ξ_y^2) * (1 + ξ_z)
    i == 17 && return 0.25(1 - ξ_x) * (1 - ξ_y) * (1 - ξ_z^2)
    i == 18 && return 0.25(1 + ξ_x) * (1 - ξ_y) * (1 - ξ_z^2)
    i == 19 && return 0.25(1 + ξ_x) * (1 + ξ_y) * (1 - ξ_z^2)
    i == 20 && return 0.25(1 - ξ_x) * (1 + ξ_y) * (1 - ξ_z^2)
    throw(ArgumentError("no shape function $i for interpolation $ip"))
end


text:
__precompile__() # this module is safe to precompile
module IMinuit

using PyCall: PyObject, pycall, PyNULL, PyAny, PyVector, pyimport_conda, pyimport, pytype_mapping
import PyCall: PyObject, pycall
# import PyCall: hasproperty # Base.hasproperty in Julia 1.2
import Base: convert, ==, isequal, hash, hasproperty,  haskey

using ForwardDiff: gradient

export Minuit, migrad, minos, hesse, matrix, iminuit, args, model_fit, @model_fit
export AbstractFit, Fit, ArrayFit, func_argnames, Data, chisq, @plt_data, @plt_data!, @plt_best, @plt_best!
export gradient, LazyHelp, contour, mncontour, mnprofile, draw_mncontour, profile,
    draw_contour, draw_profile
export get_contours, get_contours_all, contour_df, get_contours_given_parameter
export contour_df_given_parameter, get_contours_samples, contour_df_samples

# copied from PyPlot.jl
# that lazily looks up help from a PyObject via zero or more keys.
# This saves us time when loading iminuit, since we don't have
# to load up all of the documentation strings right away.
struct LazyHelp
    o # a PyObject or similar object supporting getindex with a __doc__ property
    keys::Tuple{Vararg{String}}
    LazyHelp(o) = new(o, ())
    LazyHelp(o, k::AbstractString) = new(o, (k,))
    LazyHelp(o, k1::AbstractString, k2::AbstractString) = new(o, (k1,k2))
    LazyHelp(o, k::Tuple{Vararg{AbstractString}}) = new(o, k)
end
function show(io::IO, ::MIME"text/plain", h::LazyHelp)
    o = h.o
    for k in h.keys
        o = getproperty(o, k)
    end
    if hasproperty(o, "__doc__")
        print(io, "Docstring pulled from the Python `iminuit`:\n\n")
        print(io, convert(AbstractString, o."__doc__"))
    else
        print(io, "no Python docstring found for ", h.k)
    end
end

Base.show(io::IO, h::LazyHelp) = show(io, "text/plain", h)
function Base.Docs.catdoc(hs::LazyHelp...)
    Base.Docs.Text() do io
        for h in hs
            show(io, MIME"text/plain"(), h)
        end
    end
end

@doc raw"""
    method_argnames(m::Method)

Extracting the argument names of a method as an array.
Modified from [`methodshow.jl`](https://github.com/JuliaLang/julia/blob/master/base/methodshow.jl)
(`Vector{Any}` changed to `Vector{Symbol}`).
"""
function method_argnames(m::Method)
    argnames = ccall(:jl_uncompress_argnames, Vector{Symbol}, (Any,), m.slot_syms)
    isempty(argnames) && return argnames
    return argnames[1:m.nargs]
end

"""
    func_argnames(f::Function)

Extracting the argument names of a function as an array.
"""
function func_argnames(f::Function)
    ms = collect(methods(f))
    return method_argnames(last(ms))[2:end]
end

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

include("init.jl")
include("FitStructs.jl")

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

# Wrappers of the iminuit functions
@doc raw"""
    Minuit(fcn; kwds...)
    Minuit(fcn, start; kwds...)
    Minuit(fcn, m::AbatractFit; kwds...)

Wrapper of the `iminuit` function `Minuit`.
* `fcn` is the function to be optimized.
* `start`: an array/tuple of the starting values of the parameters.
* `kwds` is the list of keyword arguments of `Minuit`. For more information, refer to the `iminuit` manual.
* `m`: a fit that was previously defined; its parameters at the latest stage can be passed to the new fit.

Example:
```
fit = Minuit(fcn, [1, 0]; name = ["a", "b"], error = 0.1*ones(2),
            fix_a = true, limit_b = (0, 50) )
migrad(fit)
```
where the parameters are collected in an array `par` which is the argument of `fcn(par)`. In this case,
one can use external code (e.g., `using ForwardDiff: gradient`) to compute the gradient as
`gradfun(par) = gradient(fcn, par)`, and include `grad = gradfun` as a keyword argument.

If `fcn` is defined as `fcn(a, b)`, then the starting values need to be
set as `Minuit(fcn, a = 1, b = 0)`.

From `iminuit`:

`Minuit(fcn, throw_nan=False, pedantic=True, forced_parameters=None, print_level=0,
        errordef=None, grad=None, use_array_call=False, **kwds)`

"""
function Minuit(fcn; kwds...)::Fit
    forced_parameters = Tuple(func_argnames(fcn)) # get the argument lists of fcn
    return iminuit.Minuit(fcn; forced_parameters= forced_parameters, pedantic = false, kwds...)
end

Minuit(fcn, start::AbstractVector; kwds...)::ArrayFit =  iminuit.Minuit.from_array_func(fcn, start; pedantic = false, kwds...)
Minuit(fcn, start::Tuple; kwds...)::ArrayFit =  iminuit.Minuit.from_array_func(fcn, start; pedantic = false, kwds...)

Minuit(fcn, m::ArrayFit; kwds...) = Minuit(fcn, args(m); (Symbol(k) => v for (k,v) in m.fitarg)...,  kwds...)
Minuit(fcn, m::Fit; kwds...) = Minuit(fcn; (Symbol(k) => v for (k,v) in m.fitarg)...,  kwds...)

# nsplit = 1, nsplit option removed since iminuit v1.5.0
function migrad(f::AbstractFit; ncall = 1000, resume = true, precision = nothing)
    return pycall(f.migrad, PyObject, ncall, resume, precision)
end

hesse(f::AbstractFit; maxcall = 0) = pycall(f.hesse, PyObject, maxcall)

function minos(f::AbstractFit; var = nothing, sigma = 1, maxcall = 0)
    return pycall(f.minos, PyObject, var, sigma, maxcall)
end

matrix(f::AbstractFit; kws...) = pycall(PyObject(f).matrix, PyObject, kws...)

function args(o::AbstractFit)::Vector{Float64}
    _a::PyObject = o.args;  _n::Int = length(_a)
    _res = zeros(_n)
    @views for i = 1:_n
        _res[i] = get(_a, i-1)
    end
    return _res
end

for fun in [:contour, :mncontour, :draw_contour, :draw_mncontour]
    :( ($fun)(f::AbstractFit, par1, par2; kws...) = f.$fun(par1, par2; kws...) ) |> eval
end

for fun in [:profile, :draw_profile, :mnprofile, :draw_mnprofile]
    :( ($fun)(f::AbstractFit, par1; kws...) = f.$fun(par1; kws...) ) |> eval
end

for f in [:migrad, :minos, :hesse, :matrix, :args, :contour, :mncontour, :profile,
          :mnprofile, :draw_mncontour, :draw_contour, :draw_profile, :draw_mnprofile]
    sf = string(f)
    @eval @doc LazyHelp(mMinuit, $sf)  function $f(ars...; kws...)
        if !hasproperty(mMinuit, $sf)
            error("iminuit ", version, " does not have iminuit.Minuit", $sf)
        end
        return pycall(mMinuit.$sf, PyAny, ars...; kws...)
    end
end

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

include("Data.jl")
include("contour.jl")

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

"""
    model_fit(model::Function, data::Data, start_values; kws...)

convenient wrapper for fitting a `model` to `data`; the returning stype is `ArrayFit`, which can be passed
to `migrad`, `minos` etc.
* `model` is the function to be fitted to `data`; it should be of the form `model(x, params)` with `params` given either as an array or a tuple.
* `start_values` can be either the initial values of parameters (`<: AbstractArray` or `<:Tuple`) or a previous fit of type `AbstractFit`.
"""
function model_fit(model::Function, data::Data, start_values; kws...)
    # _model(x, par) = length(func_argnames(model)) > 2 ? model(x, par...) : model(x, par)
    _chisq(par) = chisq(model, data, par)
    _fit = Minuit(_chisq, start_values; kws...)
    return _fit
end

"""
    @model_fit model data start_values kws...

convenient wrapper for fitting a `model` to `data`; see [`model_fit`](@ref).
"""
macro model_fit(model, data, start_values, kws...)
    _expr = quote
        _chisq(par) = chisq($model, $data, par)
        _fit = isempty($kws) ? Minuit(_chisq, $start_values) : Minuit(_chisq, $start_values; $(kws...))
    end
    esc(_expr)
end

end

ext: jl


text:
"Explicit Euler integrator."
struct IntegratorExplicitEuler{DT, TT, AT, D, ET <: NamedTuple} <: DeterministicIntegrator{DT,TT}
    equs::ET
    Δt::TT
    v::AT

    function IntegratorExplicitEuler{DT,AT,D}(equations::ET, Δt::TT) where {DT, TT, AT, D, ET}
        # create cache for vector field
        v = AT(zeros(DT, D))

        # create integrator
        new{DT,TT,AT,D,ET}(equations, Δt, v)
    end

    function IntegratorExplicitEuler{DT,AT,D}(v::Function, Δt::TT; kwargs...) where {DT,TT,AT,D}
        IntegratorExplicitEuler{DT,AT,D}(NamedTuple{(:v,)}((v,)), Δt; kwargs...)
    end

    function IntegratorExplicitEuler(equation::ODE{DT,TT,AT}, Δt; kwargs...) where {DT,TT,AT}
        IntegratorExplicitEuler{DT, AT, axes(equation)}(get_functions(equation), Δt; kwargs...)
    end
end

equation(int::IntegratorExplicitEuler, i::Symbol) = int.equs[i]
equations(int::IntegratorExplicitEuler) = int.equs
timestep(int::IntegratorExplicitEuler) = int.Δt

function integrate_step!(int::IntegratorExplicitEuler{DT,TT,AT}, sol::AtomicSolutionODE{DT,TT,AT}) where {DT,TT,AT}
    # reset atomic solution
    reset!(sol, timestep(int))

    # compute vector field
    equations(int)[:v](sol.t̄, sol.q, int.v)

    # compute update
    sol.q .+= timestep(int) .* int.v
end


text:
const ValueMap = Dict{Symbol,Union{Int,Rational}}
const RExprMap = Dict{Symbol,RExpr}
const ExprVec = Vector{Expr}

function extract(branches::Vector{Expr}, init::ValueMap)
    [__extract_single(ExprVec(b.args), init) for b in branches]
end

# Extract recurrence system from single-path loop
function __extract_single(exprs::Vector{Expr}, init::ValueMap, lc::Symbol = Recurrences.gensym_unhashed(:n))
    lhss, rhss = _parallel(split_assign(exprs)...)
    @debug "Splitted and parallel assignments" rhss lhss

    lines = Expr[]
    mvars = Symbol[]
    loop_counters = Dict{Symbol,RExpr}()
    for (l, r) in zip(lhss, rhss)
        _mvars = find_multiplications(r)
        filter!(x->x!=l, _mvars)
        append!(mvars, _mvars)
        cform = __is_loop_counter(l, r, init, lc)
        if !isnothing(cform)
            loop_counters[l] = cform
        end
    end
    # if !isempty(setdiff(mvars, keys(loop_counters)))
    #     @info "" mvars loop_counters
    #     error("Unsupported multiplication between variables")
    # end
    replacelc = !isempty(mvars)
    for (l, r) in zip(lhss, rhss)
        # skip loop counters
        # haskey(loop_counters, l) && continue
        _r = r
        for (v, cform) in loop_counters
            _r = Recurrences.symbol_walk(x->:($x($(lc))), r)
        end
        _l = Recurrences.symbol_walk(x->:($x($(lc)+1)), l)
        _r = Recurrences.symbol_walk(r) do x
            if x != l && haskey(loop_counters, x) && replacelc
                loop_counters[x]
            else
                :($x($(lc)))
            end
        end
        push!(lines, :($(_l) = $(_r)))
    end
    @debug "Lines for rec system creation" lines
    lrs = Recurrences.lrs(lines)
    Recurrences.invertible_system!(lrs)
    cfs = Recurrences.solve(lrs)
    return cfs
end

function __is_loop_counter(l::Symbol, r::RExpr, init::ValueMap, lc::Symbol)
    inc = Basic(r) - Basic(l)
    if isempty(SymEngine.free_symbols(inc)) && haskey(init, l)
        return :($(init[l]) + $(convert(Expr, inc))*($lc))
    end
    return nothing
end


text:
using PyPlot
using Random, Images, JLD2, LinearAlgebra
using JOLI, Statistics, FFTW
using Printf
using JUDI
using MECurvelets

using ArgParse
include("../utils/parse_cmd.jl")
parsed_args = parse_commandline()
nv = 16
vm = "Standard_E16s_v3"
nth = 2
batchsize = 8
nsrc = 16

creds=joinpath(pwd(),"..","credentials.json")
Random.seed!(1234);

import MECurvelets.soft_thresholding

function soft_thresholding(c::Vector{Complex{vDt}}, λ::vDt) where {vDt}
    out = exp.(angle.(c)im) .* max.(abs.(c) .- convert(vDt, λ), 0)
    return out
end

function soft_thresholding(c::Vector{Complex{vDt}}, λ::Vector{vDt}) where {vDt}
    out = exp.(angle.(c)im) .* max.(abs.(c) .- convert.(vDt, λ), 0)
    return out
end

Random.seed!(1234);

JLD2.@load "../models/Compass_tti_625m.jld2"
JLD2.@load "../data/dobs$(nv)vint$(nsrc)nsrcnoisefree.jld2" dobs_stack q_stack

C = joMECurvelet2D(n; DDT=Float32, RDT=Complex{Float32})
support = BitVector(undef, size(C,1))   # initially zero

idx_wb = find_water_bottom(rho.-rho[1,1])

v = 1f0./sqrt.(m)
v0 = deepcopy(v)
v0[:,maximum(idx_wb)+1:end] = 1f0./convert(Array{Float32,2},imfilter(1f0./v[:,maximum(idx_wb)+1:end], Kernel.gaussian(10)))
m0 = 1f0./v0.^2f0

rho0 = deepcopy(rho)
rho0[:,maximum(idx_wb)+1:end] = 1f0./convert(Array{Float32,2},imfilter(1f0./rho[:,maximum(idx_wb)+1:end], Kernel.gaussian(10)))

model0 = Model(n,d,o,m0; rho=rho0, nb=80)

opt = JUDI.Options(isic=true)

# Preconditioner

Tm = judiTopmute(model0.n, maximum(idx_wb), 3)  # Mute water column
S = judiDepthScaling(model0)
Mr = Tm*S

niter = 8

for L = 2:nv

    #init_culsterless(L; credentials=creds, vm_size=vm, pool_name="JRM$(L)vintages", verbose=1, nthreads=nth, auto_scale=false, n_julia_per_instance=batchsize)
    @printf("Simultaneous source JRM Iteration at vintage %d \n", L)
    flush(Base.stdout)

    w = ones(Float32, length(support)) - 0.5f0 .* support # support is 0.5, non-support is 1
    W = joDiag(1f0./w; DDT=Complex{Float32}, RDT=Complex{Float32})
    x = [zeros(Complex{Float32}, size(C,1)) for i = 1:L+1]
    z = [zeros(Complex{Float32}, size(C,1)) for i = 1:L+1]
    γ   = Float32.(sqrt(L)) # hyperparameter to tune

    # sim src acquisition
    xsrc_stack = [[[q_stack[i].geometry.xloc[s][1] for s = 1:dobs_stack[i].nsrc] for k = 1:batchsize] for i = 1:L]
    ysrc_stack = [[[0.0f0] for k = 1:batchsize] for i = 1:L]
    zsrc_stack = [[[q_stack[i].geometry.zloc[s][1] for s = 1:dobs_stack[i].nsrc] for k = 1:batchsize] for i = 1:L]

    src_geometry_stack = [Geometry(xsrc_stack[i], ysrc_stack[i], zsrc_stack[i]; dt=q_stack[i].geometry.dt[1], t=q_stack[i].geometry.t[1]) for i = 1:L]

    # Main loop

    for  j=1:niter

        # Set up weights for current iteration
        weights = [randn(Float32, batchsize, nsrc) for i = 1:L]

        # Create wavelet
        wavelet = [[hcat(q_stack[i].data...) .* weights[i][k:k,:] for k = 1:batchsize] for i = 1:L]

        # Create sim src
        q_j = [judiVector(src_geometry_stack[i], wavelet[i]) for i = 1:L]  # super shot for current iteration

        # Create sim data
        data_j = [[sum(weights[i][k,:].*dobs_stack[i].data) for k = 1:batchsize] for i = 1:L]
        dobs_j = [judiVector(Geometry(dobs_stack[i].geometry.xloc[1],dobs_stack[i].geometry.yloc[1],dobs_stack[i].geometry.zloc[1]; dt=dobs_stack[i].geometry.dt[1], t=dobs_stack[i].geometry.t[1], nsrc=batchsize),data_j[i]) for i = 1:L]

        dmx = [Mr*C'*W*(1f0/γ*x[1]+x[i+1]) for i = 1:L]

        phi, g1 = lsrtm_objective([model0 for i = 1:L], q_j, dobs_j, dmx; options=opt, nlind=true)
        g = [1f0/γ*W*C*Mr'*vec(sum(g1)), [W*C*Mr'*vec(g1[i]) for i=1:L]...]

	    @printf("At iteration %d function value is %2.2e \n", j, phi)
	    flush(Base.stdout)

        # Step size and update variable
        t = 2f-3*phi/norm(g)^2f0; # fixed step
        z -= t .* g;
        if L == 2 && j == 1
	        global threshold = 30*[Float32.(quantile(abs.(z[i]), .98)) for i = 1:L+1]
            global threshold1 = threshold[1]
            global threshold2 = maximum(threshold[2:end])
        end
        if j == 1
            z -= 5 * t .* g;
        end
        x[1] = soft_thresholding(z[1], threshold1)
        x[2:L+1] = [soft_thresholding(z[i], threshold2) for i = 2:L+1]
    end
    JLD2.@save "../results/JRMsimnv$(nv)no$(L)IncreJRMwl1.jld2" x z threshold support
    for i = 1:L+1
        global support = support .| (abs.(x[i]).>0f0)
    end
end


text:
using ITensors
using ITensorParallel

ITensors.BLAS.set_num_threads(1)
ITensors.Strided.disable_threads()

# Turn block sparse multithreading on or off
ITensors.disable_threaded_blocksparse()

@show Threads.nthreads()

threaded_projmpo = true

Nx, Ny = 10, 5
N = Nx * Ny

sites = siteinds("S=1/2", N; conserve_qns=true)

lattice = square_lattice(Nx, Ny; yperiodic=false)

function heisenberg_hamiltonian(lattice)
  ℋ = OpSum()
  for b in lattice
    ℋ += 0.5, "S+", b.s1, "S-", b.s2
    ℋ += 0.5, "S-", b.s1, "S+", b.s2
    ℋ += "Sz", b.s1, "Sz", b.s2
  end
  return ℋ
end
ℋ = heisenberg_hamiltonian(lattice)

# Define how the Hamiltonian will be partitioned
# into a sum of sub-Hamiltonians.
function in_partition(sites::Tuple{Int,Int}, p, nparts)
  i, j = sites
  return p == mod1(i, nparts)
end

ℋs = opsum_sum(ℋ, Threads.nthreads(); in_partition=in_partition)

H = Vector{MPO}(undef, Threads.nthreads())
Threads.@threads for n in 1:Threads.nthreads()
  H[Threads.threadid()] = splitblocks(linkinds, MPO(ℋs[n], sites))
end

PH = if threaded_projmpo
  ThreadedProjMPOSum(H)
else
  ProjMPOSum(H)
end

state = [isodd(n) ? "Up" : "Dn" for n in 1:N]
psi0 = randomMPS(sites, state, 10)

sweeps = Sweeps(10)
setmaxdim!(sweeps, 20, 60, 100, 100, 200, 400, 800)
setcutoff!(sweeps, 1E-10)
@show sweeps

energy, psi = dmrg(PH, psi0, sweeps)
@show energy

max_stars_repo_name: mtfishman/ITensorParallel

max_issues_repo_name: mtfishman/ITensorParallel

max_forks_repo_name: mtfishman/ITensorParallel


text:
# TODO : QuBase.jl - WIP : complex and expectation.

@doc """
Quantum Monte-Carlo Wave Function Method

Step Propagation using the exponential solver Expokit.expmv.

### Fields :

* eps  :: Float64

  Random number used to realize (quantum) jumps.

* options :: Dict

  Dictionary to set the solver used to proagate the state generated using `draw`.
""" ->
type QuMCWF <: QuPropagatorMethod
    eps::Float64
    options::Dict{Symbol, Any}
end

QuMCWF(options::Dict=Dict()) = QuMCWF(rand(), options)

@doc """
Ensemble of state, number of trajectories, decomposition based on the state.

### Fields :

* rho <: QuBase.AbstractQuArray

  State of the system
* ntraj :: Int

  Number of trajectories
* decomp

  Decomposition is of state `rho` if `rho` is a `QuBase.AbstractQuMatrix`.
""" ->
type QuMCWFEnsemble{QA<:QuBase.AbstractQuArray}
    rho::QA
    ntraj::Int
    decomp
    QuMCWFEnsemble(rho, ntraj, decomp) = new(rho, ntraj, decomp)
end

QuMCWFEnsemble{QM<:QuBase.AbstractQuMatrix}(rho::QM, ntraj=500) = QuMCWFEnsemble{QM}(rho, ntraj, eigfact(full(coeffs(rho))))
QuMCWFEnsemble{QV<:QuBase.AbstractQuVector}(psi::QV, ntraj=500) = QuMCWFEnsemble{QV}(psi, ntraj, nothing)

Base.length(en::QuMCWFEnsemble) = en.ntraj

@doc """
Iterator for QuMCWFEnsemble

Iterator `start` method of the Ensemble

### Arguments

Inputs :
* mcwfensemble :: QuMCWFEnsemble

Output :
* 1
""" ->
Base.start(mcwfensemble::QuMCWFEnsemble) = 1

@doc """
Iterator for QuMCWFEnsemble

Iterator `next` method of the Ensemble

### Arguments

Inputs :
* mcwfensemble :: QuMCWFEnsemble

* state :: Int

Output :
* State using `draw` and next state.
""" ->
function Base.next(mcwfensemble::QuMCWFEnsemble, i::Int)
    state = draw(mcwfensemble)
    return (state, i+1)
end

@doc """
Iterator for QuMCWFEnsemble

Iterator `done` method of the Ensemble

### Arguments

Inputs :
* mcwfensemble :: QuMCWFEnsemble

* state ::  Int

Output :
* Compares the iterator state with the number of trajectories and
  returns true or false depending on the comparison.
""" ->
function Base.done(mcwfensemble::QuMCWFEnsemble, i::Int)
    i > mcwfensemble.ntraj
end

@doc """
Method to draw a vector from the decomposition, using random number for selection.

### Arguments

Inputs :
* mcwfensemble :: QuMCWFEnsemble{QuBase.AbstractQuMatrix}

  Ensemble on which decomposition is performed

Output :
* QuArray with eigen vector from decomposition passed as an argument.
""" ->
function draw{QM<:QuBase.AbstractQuMatrix}(mcwfensemble::QuMCWFEnsemble{QM})
    r = rand() # draw random number from [0,1)
    pacc = 0.
    for i=1:length(mcwfensemble.decomp[:values])
        pacc = pacc + mcwfensemble.decomp[:values][i]
        if pacc >= r
            return QuArray(mcwfensemble.decomp[:vectors][:,i])
        end
    end
end

@doc """
Method to draw a vector, which returns a copy of the vector instead of decomposing.

### Arguments

Inputs :
* mcwfensemble :: QuMCWFEnsemble{QuBase.AbstractQuVector}

  Ensemble which is constructed using a QuVector

Output :
* Copy of the QuVector.
""" ->
function draw{QV<:QuBase.AbstractQuVector}(mcwfensemble::QuMCWFEnsemble{QV})
    return copy(mcwfensemble.rho)
end

function propagate(prob::QuMCWF, eq::QuLindbladMasterEq, t_end, t_start, psi)
    const jtol = 1.e-6  # jump tolerance
    # get information of QME
    heff = eff_hamiltonian(eq)
    c_ops = eq.collapse_ops
    dt = t_end - t_start
    accdt = 0.
    while accdt < t_end - t_start
        # propagate one time-step
        # the following only works for time independent cases (WIP for time dependent)
        psi1 = propagate(get(prob.options, :solver, QuExpmV()), QuSchrodingerEq(heff), t_start + dt, t_start, psi)
        if norm(psi1)^2 > prob.eps + jtol/2
            psi = copy(psi1)  # no jump
            accdt = accdt + dt
            dt = (t_end - t_start) - accdt
        elseif norm(psi1)^2 < prob.eps - jtol/2
            dt = dt/2		# jump in the current interval
        else
            if isempty(c_ops)
                psi = copy(psi1)
            else
                PnS = 0.
                for cind in c_ops
                    # PnS = PnS + real(expectation(psi1, cind'*cind))
                    PnS =  PnS + real(dot(psi1, cind'*cind*psi1))
                end
                Pn = 0.
                for cind in c_ops
                    # Pn = Pn + real(expectation(psi1, cind'*cind))/PnS
                    Pn = Pn + real(dot(psi1, cind'*cind*psi1))/PnS
                    if Pn >= prob.eps
                        psi = cind*psi1
                        break
                    end
                end
            end
            normalize!(psi)
            prob.eps =  rand()
            accdt = accdt + dt
            dt =  (t_end - t_start) - accdt
        end
    end
    return psi
end

export QuMCWFEnsemble,
      QuMCWF


text:
#=

Pkg.add("JLD2")
Pkg.add("Random")
Pkg.add("LinearAlgebra")
Pkg.add("Statistics")
Pkg.add("CSV")
Pkg.add("DataFrames")
Pkg.add("FreqTables")

Pkg.add("Kronecker")

=#

import JLD2
import Random
import LinearAlgebra
import Statistics
import CSV
import DataFrames
import FreqTables

import Kronecker

###################################################################
### Question 1

#function q1()

## (a)

Random.seed!(1234)

# i matrix 10* 7 with uniform values between -5 and 10

A=rand(-5:10,(10,7)) # -5:10 is the range, 10,7 is the size of the matrix

# ii matrix 10*7 with normal values (-2,15) => not re-centered, do not know how

B=randn(10,7)*15 # 15 is the standard deviation

# iii put 5 rows * 5 columns of the first together with 5 rows * 2 columns (last) of the second

A_small=A[1:5,:1:5]
A_small
B_small=B[1:5,6:7]
B_small

C=[A_small B_small] # same on the row dimension only 5*5 and 5*2 matrix

size(C)

# iv 10*7 matrix where values are equal to As if positive, else 0s

A

D=A

A # it seems that this changes M1 too - Why?

for row=1:10
    for column=1:7  # goes through all arguments of the matrix
        if D[row,column]<0
            D[row,column]=0
        end
    end
end

D

## (b) get the number of elements of A

length(A)

## (c) list the number of unique elements in D

nr_unique=length(unique(D))

# (d)

println(reshape(B,(70,1))) # reshaped it into a vector 70*1

size(reshape(B,(70,1))) # 70,1

size(reshape(B,(1,70))) # 1,70

# (e)

E=rand(2,10,7) # a matrix with the right dimensions but random numbers

E[1,1:10,1:7]=A

E[2,1:10,1:7]=B

E

# (f)

E=permutedims(E, [2,3,1]) # what used to be the second dimension now is first ...

E

# (g)

F = A ⊗ B

#F = C ⊗ E # no method matching error - Why?

size(A) # 10,7
size(B) # 10,7
size(C) # 5,7
size(E) # 10,7,2 => dimensions do not match

# (h)

using JLD2
@save("C:/Users/tsankova/Dropbox/Structural Course Material/Week 1/matrixpractice.jld", "arr", A, B, C, D, E, F)

# (i)

using JLD2
@save("C:/Users/tsankova/Dropbox/Structural Course Material/Week 1/firstmatrix.jld", "arr", A, B, C, D)

# (j)

using JLD2
@save("C:/Users/tsankova/Dropbox/Structural Course Material/Week 1/Cmatrix.csv", "arr", C)

# (k)

using DataFrames
D=convert(DataFrame, D)

using JLD2
@save("C:/Users/tsankova/Dropbox/Structural Course Material/Week 1/Dmatrix.dat", D)

# (l)

#end

#q1()

###################################################################
### Question 2

#function q2()

# (a)

AB=zeros(10,7)

for row=1:10
    for column=1:7  # goes through all arguments of the matrix
        a=A[row,column]
        println(a)
        b=B[row,column]
        println(b)
        c=a*b
        println(c)
        AB[row,column]=c
    end
end

# (b)

C

Cprime=reshape(C,(35,1))

Cprime=Vector{Float64}(vec(Cprime))

Cprime

try
    for row=1:length(Cprime) # length changes as drop obs. so in the end there is an error message
        if Cprime[row]>5
            deleteat!(Cprime,row)
        elseif Cprime[row]<-5
            deleteat!(Cprime,row)
        end
    end
catch
finally
end

Cprime # 26 elements left

# Do it not in a loop => left for later

# (c)

# (d)

# (e)

# (f)

#end

#q2()

###################################################################
### Question 3

#function q3()

# (a)

using CSV
nlsw88 = CSV.read("nlsw88.csv")

names(nlsw88) # lists variable names, seem file

using Statistics
mean(nlsw88.grade) # mean of grade variable - missing?

# missing values?

using JLD2
@save("C:/Users/tsankova/Dropbox/Structural Course Material/Week 1/nlsw88.jld", "df", nlsw88)

# (b)

# share married - married variable

using Statistics
mean(nlsw88.married) # 0.642

# share with a college degree - collgrad variable

using Statistics
mean(nlsw88.collgrad) # 0.237

# (c) do not know how to show percentag ???

using FreqTables
freqtable(nlsw88.race)

# (d)

summarystats=describe(nlsw88) # min, max, median, mean values

nlsw88

using FreqTables
freqtable(nlsw88.grade) # 2 missing

# (e)

using FreqTables
freqtable(nlsw88.industry)
freqtable(nlsw88.occupation)
freqtable(nlsw88.industry,nlsw88.occupation)

# (f)

ind_wage=select(nlsw88,:industry,:wage) # select data wanted
ind_wage_grouped=groupby(ind_wage,1) # create a grouped data frame
ind_wage_selected=combine(ind_wage_grouped,:wage => mean) # calculate mean wage by group

occ_wage=select(nlsw88,:occupation,:wage)
occ_wage_grouped=groupby(occ_wage,1)
occ_wage_selected=combine(occ_wage_grouped,:wage => mean)

#end

#q3()

###################################################################
### Question 4

#function q4()

# (a)

@load "firstmatrix.jld"

# (b) and (c)

function matrixops(A_some,B_some)

    #=

    The function takes A and B matrices as agruments of the same
    dimension, then creates three matrices - one element-by-element
    multiplication, one matrix multiplication and one matrix
    summation

    =#

    size_A=size(A_some)

    size_A_r=size_A[1] # number of rows
    size_A_c=size_A[2] # number of columns

    e_by_e_product_A_B=zeros(size_A_r,size_A_c)

    for row in 1:size_A_r
        for column in 1:size_A_c
            a=A_some[row,column]
            b=B_some[row,column]
            c=a*b
            e_by_e_product_A_B[row,column]=c
        end
    end

    A_transpose_B=transpose(A_some)*B_some

    A_plus_B=A_some+B_some

    return e_by_e_product_A_B, A_transpose_B, A_plus_B

end

# (d)

matrixops(A,B) # returns three matrices of the correct size

# (e)

function matrixops(A_some,B_some)

    size_A=size(A_some)
    size_B=size(B_some)

    if size_A==size_B # like stata

        size_A_r=size_A[1] # number of rows
        size_A_c=size_A[2] # number of columns

        e_by_e_product_A_B=zeros(size_A_r,size_A_c)

        for row in 1:size_A_r
            for column in 1:size_A_c
                a=A_some[row,column]
                b=B_some[row,column]
                c=a*b
                e_by_e_product_A_B[row,column]=c
            end
        end

        A_transpose_B=transpose(A_some)*B_some

        A_plus_B=A_some+B_some

        return e_by_e_product_A_B, A_transpose_B, A_plus_B

    else

        println("Matrices of different sizes")

    end

end

# (e)

matrixops(C,D) # prints the message

# (f)

ttl_exp=convert(Array,nlsw88.ttl_exp)

wage=convert(Array,nlsw88.wage)

ttl_exp
wage # vectors

ttl_exp_array=reshape(ttl_exp, length(ttl_exp), 1)
wage_array=reshape(wage, length(wage), 1)

matrixops(ttl_exp_array,wage_array) # yes

#end

#q4()

#=

List of issues:

Do not know how to fix the normal distribution
to have a different mean and standard deviation

Have not done the second part of Question 2 -
need to see answers

Q: Do we write using ...
when the packages have sub-packaged / commands with the same name
only or we should always do it? Seems wordy

Q: Why do we need three dimensional arrays?

Q: If I want to replace the values of a matrix with something
one by one then I first need to create an empty / zeros matrix -
is this true?

Do not know how to declare the functions
without inputs and outpits
- using expressions need to come forward
- some other pieces of code do not run (kron product)
- unclear when to use return as well

Q: Is it good practice to wrap everything in a function?

Q: Should I have se a home directory?

max_stars_repo_path: PS1 Answers - TTsankova 31.08.2020.jl

max_issues_repo_path: PS1 Answers - TTsankova 31.08.2020.jl

max_forks_repo_path: PS1 Answers - TTsankova 31.08.2020.jl


text:
# =============================================================================
# graphics.jl -> Funciones útiles para visualizar gráficamente los resultados de una RNA:
#   -
# =============================================================================

using Plots;

function print_train_results(
    trainingLosses::Array{Float64,1},validationLosses::Array{Float64,1},testLosses::Array{Float64,1},
    trainingAccuracies::Array{Float64,1},validationAccuracies::Array{Float64,1},testAccuracies::Array{Float64,1})
    x=1:size(trainingLosses,1)
    losses=plot(x,[trainingLosses validationLosses testLosses],title = "Loss",label = ["Training" "Validation" "Test"],);
    xlabel!("Epoch");
    ylabel!("Error");
    accuracy=plot(x,[trainingAccuracies validationAccuracies testAccuracies],title = "Accuracies",label = ["Training" "Validation" "Test"],);
    xlabel!("Epoch");
    ylabel!("% Accuracies");
    display(plot(losses,accuracy));
end;

function printAccStd(mean_acc::Array{Any,1}, sdev::Array{Any,1}, mean_f1::Array{Any,1}, sdev_f1::Array{Any,1}, N::Int64, xlabel::String, rep::Symbol)
    if rep==:All
        m = plot([1:N],[mean_acc mean_f1],title = "Mean Values",label = ["Accurracies" "F1-Score"]);
        xlabel!(xlabel);
        ylabel!("Precision");
        stdd = plot([1:N],[sdev sdev_f1],title = "Standard Deviation",label = ["Accurracies" "F1-Score"]);
        xlabel!(xlabel);
        display(plot(m,stdd));
    elseif rep==:AccStd
        m = plot([1:N],mean_acc,title = "Mean Values",label = "Accurracies");
        xlabel!(xlabel);
        ylabel!("Precision");
        stdd = plot([1:N],sdev,title = "Standard Deviation",label = "Standard Deviation");
        xlabel!(xlabel);
        display(plot(m,stdd));
    elseif rep==:F1Std
        m = plot([1:N],mean_f1,title = "Mean Values",label = "F1-Score");
        xlabel!(xlabel);
        ylabel!("Precision");
        stdd = plot([1:N],sdev_f1,title = "Standard Deviation",label = "Standard Deviation");
        xlabel!(xlabel);
        display(plot(m,stdd));
    elseif rep==:AccF1
        m = plot([1:N],[mean_acc mean_f1],title = "Mean Values",label = ["Accurracies" "F1-Score"]);
        xlabel!(xlabel);
        ylabel!("Precision");
        display(m);
    end;
end;

function printAccStdRNA(mean_acc::Array{Any,1}, sdev::Array{Any,1}, mean_f1::Array{Any,1}, sdev_f1::Array{Any,1}, topologyarr::Array{Any,1}, rep::Symbol)
    if rep==:All
        m = plot(topologyarr,[mean_acc mean_f1],title = "Mean Values",label = ["Accurracies" "F1-Score"]);
        xlabel!("Topology");
        ylabel!("Precision");
        stdd = plot(topologyarr,[sdev sdev_f1],title = "Standard Deviation",label = ["Accurracies" "F1-Score"]);
        xlabel!("Topology");
        display(plot(m,stdd));
    elseif rep==:AccStd
        m = plot(topologyarr,mean_acc,title = "Mean Values",label = "Accurracies");
        xlabel!("Topology");
        ylabel!("Precision");
        stdd = plot(topologyarr,sdev,title = "Standard Deviation",label = "Standard Deviation");
        xlabel!("Topology");
        display(plot(m,stdd));
    elseif rep==:F1Std
        m = plot(topologyarr,mean_f1,title = "Mean Values",label = "F1-Score");
        xlabel!("Topology");
        ylabel!("Precision");
        stdd = plot(topologyarr,sdev_f1,title = "Standard Deviation",label = "Standard Deviation");
        xlabel!("Topology");
        display(plot(m,stdd));
    elseif rep==:AccF1
        m = plot(topologyarr,[mean_acc mean_f1],title = "Mean Values",label = ["Accurracies" "F1-Score"]);
        xlabel!("Topology");
        ylabel!("Precision");
        display(m);
    end;
end;

function print_train_results(num_ciclos, trainingLosses, testLosses, trainingAccuracies, testAccuracies)
    losses=plot(num_ciclos,[trainingLosses testLosses],title = "Loss",label = ["Training" "Test"],);
    xlabel!("Ciclo");
    ylabel!("Error(%)");
    accuracy=plot(num_ciclos,[trainingAccuracies testAccuracies],title = "Accuracies",label = ["Training" "Test"],);
    xlabel!("Ciclo");
    ylabel!("Accuracies(%)");
    display(plot(losses,accuracy));
end;

function oneMetric(x::Array{Any,1}, result::Array{Any,1})
    m = plot(x, result, title = "Convolutional Neural Network",label = "Accurracy");
    xlabel!("Cycle");
    ylabel!("Precision");
    display(m);
end;


text:
module CSDP

using LinearAlgebra # For Diagonal
using SparseArrays # For SparseMatrixCSC

if VERSION < v"1.3"
    if isfile(joinpath(dirname(@__FILE__),"..","deps","deps.jl"))
        include("../deps/deps.jl")
    else
        error("CSDP not properly installed. Please run Pkg.build(\"CSDP\")")
    end
else
    import CSDP_jll: libcsdp
end

# This is the size of int used by the LAPACK library used by CSDP.
# If libcsdp is patched to use a 64-bit integer LAPACK, this should be replaced by `Clong`.
const CSDP_INT = Cint

export Blockmatrix

include("blockdiag.jl")
include("blockmat.h.jl")
include("blockmat.jl")
include("declarations.h.jl")
include("declarations.jl")
include("debug-mat.jl")
include("options.jl")
include("MOI_wrapper.jl")

end # module


text:
"""
    sig = signature(sigsv::SimpleVector)

Compute a method signature from a suitable `SimpleVector`: `sigsv[1]` holds the signature
and `sigsv[2]` the `TypeVar`s.

# Example:

For `f(x::AbstractArray{T}) where T`, the corresponding `sigsv` is constructed as

    T = TypeVar(:T)
    sig1 = Core.svec(typeof(f), AbstractArray{T})
    sig2 = Core.svec(T)
    sigsv = Core.svec(sig1, sig2)
    sig = signature(sigsv)
"""
function signature(sigsv::SimpleVector)
    sigp, sigtv = sigsv
    sig = Tuple{sigp...}
    for i = length(sigtv):-1:1
        sig = UnionAll(sigtv[i], sig)
    end
    return sig
end

"""
    sigt, lastpc = signature(recurse, frame, pc)
    sigt, lastpc = signature(frame, pc)

Compute the signature-type `sigt` of a method whose definition in `frame` starts at `pc`.
Generally, `pc` should point to the `Expr(:method, methname)` statement, in which case
`lastpc` is the final statement number in `frame` that is part of the signature
(i.e, the line above the 3-argument `:method` expression).
Alternatively, `pc` can point to the 3-argument `:method` expression,
as long as all the relevant SSAValues have been assigned.
In this case, `lastpc == pc`.

If no 3-argument `:method` expression is found, `sigt` will be `nothing`.
"""
function signature(@nospecialize(recurse), frame::Frame, @nospecialize(stmt), pc)
    mod = moduleof(frame)
    lastpc = frame.pc = pc
    while !isexpr(stmt, :method, 3)  # wait for the 3-arg version
        if isanonymous_typedef(stmt)
            lastpc = pc = step_through_methoddef(recurse, frame, stmt)   # define an anonymous function
        elseif isexpr(stmt, :call) && is_quotenode(stmt.args[1], Core.Typeof) &&
               (sym = stmt.args[2]; isa(sym, Symbol) && !isdefined(mod, sym))
            return nothing, pc
        else
            lastpc = pc
            pc = step_expr!(recurse, frame, stmt, true)
            pc === nothing && return nothing, lastpc
        end
        stmt = pc_expr(frame, pc)
    end
    sigsv = @lookup(frame, stmt.args[2])::SimpleVector
    sigt = signature(sigsv)
    return sigt, lastpc
end
signature(@nospecialize(recurse), frame::Frame, pc) = signature(recurse, frame, pc_expr(frame, pc), pc)
signature(frame::Frame, pc) = signature(finish_and_return!, frame, pc)

function minid(@nospecialize(node), stmts, id)
    if isa(node, SSAValue)
        id = min(id, node.id)
        stmt = stmts[node.id]
        return minid(stmt, stmts, id)
    elseif isa(node, Expr)
        for a in node.args
            id = minid(a, stmts, id)
        end
    end
    return id
end

function signature_top(frame, stmt::Expr, pc)
    @assert ismethod3(stmt)
    return minid(stmt.args[2], frame.framecode.src.code, pc)
end

function step_through_methoddef(@nospecialize(recurse), frame, @nospecialize(stmt))
    while !isexpr(stmt, :method)
        pc = step_expr!(recurse, frame, stmt, true)
        stmt = pc_expr(frame, pc)
    end
    return step_expr!(recurse, frame, stmt, true)  # also define the method
end

"""
    MethodInfo(start, stop, refs)

Given a frame and its CodeInfo, `start` is the line of the first `Expr(:method, name)`,
whereas `stop` is the line of the last `Expr(:method, name, sig, src)` expression for `name`.
`refs` is a vector of line numbers of other references.
Some of these will be the location of the "declaration" of a method,
the `:thunk` expression containing a CodeInfo that just returns a 1-argument `:method` expression.
Others may be `:global` declarations.

In some cases there may be more than one method with the same name in the `start:stop` range.
"""
mutable struct MethodInfo
    start::Int
    stop::Int
    refs::Vector{Int}
end
MethodInfo(start) = MethodInfo(start, -1, Int[])

"""
    methodinfos, selfcalls = identify_framemethod_calls(frame)

Analyze the code in `frame` to locate method definitions and "self-calls," i.e., calls
to methods defined in `frame` that occur within `frame`.

`methodinfos` is a Dict of `name=>info` pairs, where `info` is a [`MethodInfo`](@ref).

`selfcalls` is a list of `(linetop, linebody, callee, caller)` tuples that holds the location of
calls the methods defined in `frame`. `linetop` is the line in `frame` (top meaning "top level"),
which will correspond to a 3-argument `:method` expression containing a `CodeInfo` body.
`linebody` is the line within the `CodeInfo` body from which the call is made.
`callee` is the Symbol of the called method.
"""
function identify_framemethod_calls(frame)
    refs = Pair{Symbol,Int}[]
    methodinfos = Dict{Symbol,MethodInfo}()
    selfcalls = NamedTuple{(:linetop, :linebody, :callee, :caller),Tuple{Int64,Int64,Symbol,Union{Symbol,Bool}}}[]
    for (i, stmt) in enumerate(frame.framecode.src.code)
        if isexpr(stmt, :global, 1)
            key = stmt.args[1]
            # We don't know for sure if this is a reference to a method, but let's
            # tentatively cue it
            push!(refs, key=>i)
        elseif isexpr(stmt, :thunk, 1) && stmt.args[1] isa CodeInfo
            tsrc = stmt.args[1]::CodeInfo
            if length(tsrc.code) == 1
                tstmt = tsrc.code[1]
                if isexpr(tstmt, :return, 1)
                    tex = tstmt.args[1]
                    if isexpr(tex, :method)
                        push!(refs, tex.args[1]=>i)
                    end
                end
            end
        elseif ismethod1(stmt)
            key = stmt.args[1]::Symbol
            mi = get(methodinfos, key, nothing)
            if mi === nothing
                methodinfos[key] = MethodInfo(i)
            else
                mi.stop == -1 && (mi.start = i) # advance the statement # unless we've seen the method3
            end
        elseif ismethod3(stmt)
            key = stmt.args[1]
            if key isa Symbol
                mi = methodinfos[key]
                mi.stop = i
            elseif key isa Expr   # this is a module-scoped call. We don't have to worry about these because they are named
                continue
            end
            msrc = stmt.args[3]
            if msrc isa CodeInfo
                for (j, mstmt) in enumerate(msrc.code)
                    if isexpr(mstmt, :call)
                        mkey = mstmt.args[1]
                        if isa(mkey, Symbol)
                            # Could be a GlobalRef but then it's outside frame
                            haskey(methodinfos, mkey) && push!(selfcalls, (linetop=i, linebody=j, callee=mkey, caller=key))
                        end
                    elseif isexpr(mstmt, :meta) && mstmt.args[1] == :generated
                        newex = mstmt.args[2]
                        if isexpr(newex, :new) && length(newex.args) >= 2 && is_global_ref(newex.args[1], Core, :GeneratedFunctionStub)
                            mkey = newex.args[2]::Symbol
                            haskey(methodinfos, mkey) && push!(selfcalls, (linetop=i, linebody=j, callee=mkey, caller=key))
                        end
                    end
                end
            end
        end
    end
    for r in refs
        mi = get(methodinfos, r.first, nothing)
        mi === nothing || push!(mi.refs, r.second)
    end
    return methodinfos, selfcalls
end

function callchain(selfcalls)
    calledby = Dict{Symbol,Union{Symbol,Bool}}()
    for sc in selfcalls
        startswith(String(sc.callee), '#') || continue
        caller = get(calledby, sc.callee, nothing)
        if caller === nothing
            calledby[sc.callee] = sc.caller
        elseif caller == sc.caller
        else
            error("unexpected multiple callers, ", caller, " and ", sc.caller)
        end
    end
    return calledby
end

function set_to_running_name!(@nospecialize(recurse), replacements, frame, methodinfos, calledby, callee, caller)
    if isa(caller, Symbol) && startswith(String(caller), '#')
        rep = get(replacements, caller, nothing)
        if rep === nothing
            parentcaller = get(calledby, caller, nothing)
            if parentcaller !== nothing
                set_to_running_name!(recurse, replacements, frame, methodinfos, calledby, caller, parentcaller)
            end
        else
            caller = rep
        end
    end
    if isa(caller, Symbol)
        mi = methodinfos[caller]
        cname, _pc, _ = get_running_name(recurse, frame, mi.start, callee, get(replacements, caller, caller))
    else
        # For generated constructors (which have no name), we just assume they immediately follow their callee
        mi = methodinfos[callee]
        cname, _pc, _ = get_running_name(recurse, frame, mi.stop+1, callee, get(replacements, caller, caller))
    end
    replacements[callee] = cname
    mi = methodinfos[cname] = methodinfos[callee]
    src = frame.framecode.src
    replacename!(src.code[mi.start:mi.stop], callee=>cname)        # the method itself
    for r in mi.refs                                               # the references
        replacename!((src.code[r])::Expr, callee=>cname)
    end
    return replacements
end

"""
    methranges = rename_framemethods!(frame)
    methranges = rename_framemethods!(recurse, frame)

Rename the gensymmed methods in `frame` to match those that are currently active.
The issues are described in https://github.com/JuliaLang/julia/issues/30908.
`frame` will be modified in-place as needed.

Returns a vector of `name=>start:stop` pairs specifying the range of lines in `frame`
at which method definitions occur. In some cases there may be more than one method with
the same name in the `start:stop` range.
"""
function rename_framemethods!(@nospecialize(recurse), frame::Frame, methodinfos, selfcalls, calledby)
    src = frame.framecode.src
    replacements = Dict{Symbol,Symbol}()
    for (callee, caller) in calledby
        (!startswith(String(callee), '#') || haskey(replacements, callee)) && continue
        try
            set_to_running_name!(recurse, replacements, frame, methodinfos, calledby, callee, caller)
        catch err
            @warn "skipping callee $callee (called by $caller) due to $err"
        end
    end
    for (linetop, linebody, callee, caller) in selfcalls
        cname = get(replacements, callee, nothing)
        if cname !== nothing && cname !== callee
            replacename!((src.code[linetop].args[3])::CodeInfo, callee=>cname)
        end
    end
    return methodinfos
end

function rename_framemethods!(@nospecialize(recurse), frame::Frame)
    pc0 = frame.pc
    methodinfos, selfcalls = identify_framemethod_calls(frame)
    calledby = callchain(selfcalls)
    rename_framemethods!(recurse, frame, methodinfos, selfcalls, calledby)
    frame.pc = pc0
    return methodinfos
end
rename_framemethods!(frame::Frame) = rename_framemethods!(finish_and_return!, frame)

"""
    pctop, isgen = find_corrected_name(recurse, frame, pc, name, parentname)

Scans forward from `pc` in `frame` until a method is found that calls `name`.
`pctop` points to the beginning of that method's signature.
`isgen` is true if `name` corresponds to sa GeneratedFunctionStub.

Alternatively, this returns `nothing` if `pc` does not appear to point to either
a keyword or generated method.
"""
function find_corrected_name(@nospecialize(recurse), frame, pc, name, parentname)
    stmt = pc_expr(frame, pc)
    while true
        pc0 = pc
        while !ismethod3(stmt)
            pc = next_or_nothing(frame, pc)
            pc === nothing && return nothing
            stmt = pc_expr(frame, pc)
        end
        body = stmt.args[3]
        if stmt.args[1] !== name && isa(body, SSAValue)
            # OK, we can't skip all the stuff that might define the body
            # See https://github.com/timholy/Revise.jl/issues/398
            pc = pc0
            stmt = pc_expr(frame, pc)
            while !ismethod3(stmt)
                pc = step_expr!(recurse, frame, stmt, true)
                pc === nothing && return nothing
                stmt = pc_expr(frame, pc)
            end
            body = @lookup(frame, stmt.args[3])
        end
        if stmt.args[1] !== name && isa(body, CodeInfo)
            # This might be the GeneratedFunctionStub for a @generated method
            for (i, bodystmt) in enumerate(body.code)
                if isexpr(bodystmt, :meta) && bodystmt.args[1] === :generated
                    return signature_top(frame, stmt, pc), true
                end
                i >= 5 && break  # we should find this early
            end
            if length(body.code) > 1
                bodystmt = body.code[end-1]  # the line before the final return
                iscallto(bodystmt, name) && return signature_top(frame, stmt, pc), false
            end
        end
        pc = next_or_nothing(frame, pc)
        pc === nothing && return nothing
        stmt = pc_expr(frame, pc)
    end
end

"""
    replacename!(stmts, oldname=>newname)

Replace a Symbol `oldname` with `newname` in `stmts`.
"""
function replacename!(ex::Expr, pr)
    replacename!(ex.args, pr)
    return ex
end

replacename!(src::CodeInfo, pr) = replacename!(src.code, pr)

function replacename!(args::Vector{Any}, pr)
    oldname, newname = pr
    for i = 1:length(args)
        a = args[i]
        if isa(a, Expr)
            replacename!(a, pr)
        elseif isa(a, CodeInfo)
            replacename!(a.code, pr)
        elseif isa(a, QuoteNode) && a.value == oldname
            args[i] = QuoteNode(newname)
        elseif isa(a, Vector{Any})
            replacename!(a, pr)
        elseif a == oldname
            args[i] = newname
        end
    end
    return args
end

function get_running_name(@nospecialize(recurse), frame, pc, name, parentname)
    # Get the correct name (the one that's actively running)
    nameinfo = find_corrected_name(recurse, frame, pc, name, parentname)
    if nameinfo === nothing
        pc = skip_until(stmt->isexpr(stmt, :method, 3), frame, pc)
        pc = next_or_nothing(frame, pc)
        return name, pc, nothing
    end
    pctop, isgen = nameinfo
    # Sometimes signature_top---which follows SSAValue links backwards to find the first
    # line needed to define the signature---misses out on a SlotNumber assignment.
    # Fix https://github.com/timholy/Revise.jl/issues/422
    stmt = pc_expr(frame, pctop)
    while pctop > 1 && isa(stmt, SlotNumber) && !isassigned(frame.framedata.locals, pctop)
        pctop -= 1
        stmt = pc_expr(frame, pctop)
    end   # end fix
    sigtparent, lastpcparent = signature(recurse, frame, pctop)
    sigtparent === nothing && return name, pc, lastpcparent
    methparent = whichtt(sigtparent)
    methparent === nothing && return name, pc, lastpcparent  # caller isn't defined, no correction is needed
    if isgen
        cname = nameof(methparent.generator.gen)
    else
        bodyparent = Base.uncompressed_ast(methparent)
        bodystmt = bodyparent.code[end-1]
        @assert isexpr(bodystmt, :call)
        ref = getcallee(bodystmt)
        isa(ref, GlobalRef) || @show ref typeof(ref)
        @assert isa(ref, GlobalRef)
        @assert ref.mod == moduleof(frame)
        cname = ref.name
    end
    return cname, pc, lastpcparent
end

"""
    ret = methoddef!(recurse, signatures, frame; define=true)
    ret = methoddef!(signatures, frame; define=true)

Compute the signature of a method definition. `frame.pc` should point to a
`:method` expression. Upon exit, the new signature will be added to `signatures`.

There are several possible return values:

    pc, pc3 = ret

is the typical return. `pc` will point to the next statement to be executed, or be `nothing`
if there are no further statements in `frame`. `pc3` will point to the 3-argument `:method`
expression.

Alternatively,

    pc = ret

occurs for "empty method" expressions, e.g., `:(function foo end)`. `pc` will be `nothing`.

By default the method will be defined (evaluated). You can prevent this by setting `define=false`.
This is recommended if you are simply extracting signatures from code that has already been evaluated.
"""
function methoddef!(@nospecialize(recurse), signatures, frame::Frame, @nospecialize(stmt), pc::Int; define::Bool=true)
    framecode, pcin = frame.framecode, pc
    if ismethod3(stmt)
        pc3 = pc
        sigt, pc = signature(recurse, frame, stmt, pc)
        if sigt === nothing && define
            step_expr!(recurse, frame, stmt, true)
        end
        sigt, pc = signature(recurse, frame, stmt, pc)
        meth = whichtt(sigt)
        if (meth === nothing || !(meth.sig <: sigt && sigt <: meth.sig)) && define
            step_expr!(recurse, frame, stmt, true)
            meth = whichtt(sigt)
        end
        if isa(meth, Method)
            push!(signatures, meth.sig)
        elseif stmt.args[1] === false
            # If it's anonymous and not defined, define it
            pc = step_expr!(recurse, frame, stmt, true)
            meth = whichtt(sigt)
            isa(meth, Method) && push!(signatures, meth.sig)
            return pc, pc3
        else
            # guard against busted lookup, e.g., https://github.com/JuliaLang/julia/issues/31112
            code = framecode.src
            codeloc = codelocation(code, pc)
            loc = code.linetable[codeloc]
            ft = Base.unwrap_unionall(Base.unwrap_unionall(sigt).parameters[1])
            if !startswith(String(ft.name.name), "##")
                @warn "file $(loc.file), line $(loc.line): no method found for $sigt"
            end
        end
        frame.pc = pc
        return ( define ? step_expr!(recurse, frame, stmt, true) : next_or_nothing!(frame) ), pc3
    end
    ismethod1(stmt) || error("expected method opening, got ", stmt)
    name = stmt.args[1]
    if isa(name, Bool)
        error("not valid for anonymous methods")
    end
    name = name::Symbol
    while true  # methods containing inner methods may need multiple trips through this loop
        sigt, pc = signature(recurse, frame, stmt, pc)
        stmt = pc_expr(frame, pc)
        while !isexpr(stmt, :method, 3)
            pc = next_or_nothing(frame, pc)  # this should not check define, we've probably already done this once
            pc === nothing && return nothing   # this was just `function foo end`, signal "no def"
            stmt = pc_expr(frame, pc)
        end
        pc3 = pc
        name3 = stmt.args[1]
        sigt === nothing && (error("expected a signature"); return next_or_nothing(frame, pc)), pc3
        # Methods like f(x::Ref{<:Real}) that use gensymmed typevars will not have the *exact*
        # signature of the active method. So let's get the active signature.
        frame.pc = pc
        pc = define ? step_expr!(recurse, frame, stmt, true) : next_or_nothing!(frame)
        meth = whichtt(sigt)
        isa(meth, Method) && push!(signatures, meth.sig) # inner methods are not visible
        name === name3 && return pc, pc3     # if this was an inner method we should keep going
        stmt = pc_expr(frame, pc)  # there *should* be more statements in this frame
    end
end
methoddef!(@nospecialize(recurse), signatures, frame::Frame, pc::Int; define=true) =
    methoddef!(recurse, signatures, frame, pc_expr(frame, pc), pc; define=define)
function methoddef!(@nospecialize(recurse), signatures, frame::Frame; define=true)
    pc = frame.pc
    stmt = pc_expr(frame, pc)
    if !ismethod(stmt)
        pc = next_until!(ismethod, recurse, frame, true)
    end
    pc === nothing && error("pc at end of frame without finding a method")
    methoddef!(recurse, signatures, frame, pc; define=define)
end
methoddef!(signatures, frame::Frame; define=true) =
    methoddef!(finish_and_return!, signatures, frame; define=define)

function methoddefs!(@nospecialize(recurse), signatures, frame::Frame, pc; define=true)
    ret = methoddef!(recurse, signatures, frame, pc; define=define)
    pc = ret === nothing ? ret : ret[1]
    return _methoddefs!(recurse, signatures, frame, pc; define=define)
end
function methoddefs!(@nospecialize(recurse), signatures, frame::Frame; define=true)
    ret = methoddef!(recurse, signatures, frame; define=define)
    pc = ret === nothing ? ret : ret[1]
    return _methoddefs!(recurse, signatures, frame, pc; define=define)
end
methoddefs!(signatures, frame::Frame; define=true) =
    methoddefs!(finish_and_return!, signatures, frame; define=define)

function _methoddefs!(@nospecialize(recurse), signatures, frame::Frame, pc; define=define)
    while pc !== nothing
        stmt = pc_expr(frame, pc)
        if !ismethod(stmt)
            pc = next_until!(ismethod, recurse, frame, true)
        end
        pc === nothing && break
        ret = methoddef!(recurse, signatures, frame, pc; define=define)
        pc = ret === nothing ? ret : ret[1]
    end
    return pc
end

function is_self_call(@nospecialize(stmt), slotnames, argno=1)
    if isa(stmt, Expr)
        if stmt.head == :call
            a = stmt.args[argno]
            if isa(a, SlotNumber) || isa(a, Core.SlotNumber)
                sn = slotnames[a.id]
                if sn == Symbol("#self#") || sn == Symbol("") # allow empty to fix https://github.com/timholy/CodeTracking.jl/pull/48
                    return true
                end
            end
        end
    end
    return false
end

"""
    mbody = bodymethod(m::Method)

Return the "body method" for a method `m`. `mbody` contains the code of the function body
when `m` was defined.
"""
function bodymethod(mkw::Method)
    m = mkw
    local src
    while true
        framecode = JuliaInterpreter.get_framecode(m)
        fakeargs = Any[nothing for i = 1:framecode.scope.nargs]
        frame = JuliaInterpreter.prepare_frame(framecode, fakeargs, isa(m.sig, UnionAll) ? sparam_ub(m) : Core.svec())
        src = framecode.src
        (length(src.code) > 1 && is_self_call(src.code[end-1], src.slotnames)) || break
        # Build the optional arg, so we can get its type
        pc = frame.pc
        while pc < length(src.code) - 1
            pc = step_expr!(frame)
        end
        val = pc > 1 ? frame.framedata.ssavalues[pc-1] : src.code[1].args[end]
        sig = Tuple{Base.unwrap_unionall(m.sig).parameters..., typeof(val)}
        m = whichtt(sig)
    end
    length(src.code) > 1 || return m
    stmt = src.code[end-1]
    if isexpr(stmt, :call) && (f = stmt.args[1]; isa(f, QuoteNode))
        # Check that it has a #self# call
        hasself = false
        for i = 2:length(stmt.args)
            hasself |= is_self_call(stmt, src.slotnames, i)
        end
        hasself || return m
        f = f.value
        mths = methods(f)
        if length(mths) == 1
            return first(mths)
        end
    end
    return m
end


text:
# This file was generated, do not modify it. # hide
using MLJ
import RDatasets: dataset
using PrettyPrinting
import DataFrames: DataFrame, select, Not
MLJ.color_off() # hide
DTC = @load DecisionTreeClassifier pkg=DecisionTree

carseats = dataset("ISLR", "Carseats")

first(carseats, 3) |> pretty

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
"""
    identity_vector(l::Integer)
"""
function identity_vector(l::Integer)
  IdentityVector(l)
end

struct IdentityVector{T<:Integer} <: AbstractVector{T}
  length::T
end

function getindex(c::IdentityVector{T},i::Integer) where T
  @assert i > 0
  @assert i <= c.length
  j::T = i
  j
end

size(c::IdentityVector) = (c.length,)

IndexStyle(::Type{<:IdentityVector}) = IndexLinear()

function reindex(values::AbstractArray, indices::IdentityVector)
  @assert length(values) == length(indices)
  values
end

function reindex(a::AppliedArray,b::IdentityVector)
  @assert length(a) == length(b)
  a
end

function reindex(a::Fill,b::IdentityVector)
  @assert length(a) == length(b)
  a
end

function reindex(a::CompressedArray,b::IdentityVector)
  @assert length(a) == length(b)
  a
end

function reindex(a::Reindexed,b::IdentityVector)
  @assert length(a) == length(b)
  a
end


text:
using TailRiskScenGen
using Distributions
## using Gadfly

mu = [0.0, 0.5]
cov = [[1 0.5], [0.5 1]]
num_scen = 1000

normal = MvNormal(mu, cov)
scenarios = rand(normal, num_scen)
probs = fill(1.0/num_scen, num_scen)

# print(scenarios)

new_scen, new_probs = aggregate_scenarios(scenarios, probs, mu, cov, eye(2), quantile(Normal(), 0.95))
# print(new_scenarios)

## x_min, y_min = minimum(scenarios,2)
## x_max, y_max = maximum(scenarios,2)
## x_scale = Scale.x_continuous(minvalue=x_min, maxvalue=x_max)
## y_scale = Scale.y_continuous(minvalue=y_min, maxvalue=y_max)
## p1 = plot(x=scenarios[1,:], y=scenarios[2,:], x_scale, y_scale, Guide.title("Full Scenario Set"))
## p2 = plot(x=new_scen[1,:], y=new_scen[2,:], x_scale, y_scale, Guide.title("Aggregated Scenario Set"))
## draw(PDF("data.pdf", 6inch, 6inch), vstack(p1,p2))


text:
R_trues = [
    # Case 1: τ = 0.10, (λ = 530 nm), ρ = 0.0
    [
        [
            [0.03916, 0.04235, 0.04378, 0.04593, 0.04940, 0.05438, 0.05816, 0.06541, 0.07639, 0.08773, 0.09549, 0.11842, 0.13619, 0.17836, 0.22393, 0.27069],
            [0.03916, 0.03922, 0.03930, 0.03952, 0.04021, 0.04208, 0.04409, 0.04894, 0.05748, 0.06698, 0.07367, 0.09388, 0.10977, 0.14790, 0.18964, 0.23441],
            [0.03916, 0.03634, 0.03533, 0.03414, 0.03316, 0.03381, 0.03557, 0.04073, 0.05062, 0.06186, 0.06981, 0.09387, 0.11280, 0.15826, 0.20833, 0.26423]
        ],
        [
            [0.04337, 0.05057, 0.05428, 0.06032, 0.07077, 0.08615, 0.09743, 0.11784, 0.14627, 0.17387, 0.19217, 0.24472, 0.28457, 0.37797, 0.47857, 0.58547],
            [0.04337, 0.04384, 0.04434, 0.04543, 0.04800, 0.05316, 0.05781, 0.06762, 0.08322, 0.09957, 0.11078, 0.14393, 0.16958, 0.23062, 0.29730, 0.37033],
            [0.04337, 0.03865, 0.03753, 0.03691, 0.03854, 0.04531, 0.05258, 0.06882, 0.09508, 0.12247, 0.14114, 0.19590, 0.23803, 0.33792, 0.44741, 0.57244]
        ],
        [
            [0.08599, 0.09709, 0.10459, 0.11842, 0.14567, 0.19035, 0.22510, 0.28998, 0.38196, 0.47165, 0.53118, 0.70246, 0.83299, 1.14328, 1.49145, 1.94280],
            [0.08599, 0.08786, 0.08979, 0.09388, 0.10311, 0.12041, 0.13518, 0.16489, 0.21010, 0.25623, 0.28751, 0.37921, 0.45007, 0.61998, 0.81151, 1.05733],
            [0.08599, 0.08458, 0.08701, 0.09387, 0.11184, 0.14747, 0.17799, 0.23845, 0.32807, 0.41745, 0.47732, 0.65075, 0.78354, 1.10034, 1.45755, 1.92791]
        ]
    ],
    # Case 2: τ = 0.10, (λ = 530 nm), ρ = 0.25
    [
        [
            [0.27073, 0.27049, 0.27026, 0.26981, 0.26902, 0.26814, 0.26785, 0.26811, 0.26976, 0.27228, 0.27427, 0.28078, 0.28616, 0.29929, 0.31335, 0.32473],
            [0.27073, 0.27049, 0.27026, 0.26981, 0.26902, 0.26814, 0.26785, 0.26811, 0.26976, 0.27228, 0.27427, 0.28078, 0.28616, 0.29929, 0.31335, 0.32473],
            [0.27073, 0.27049, 0.27026, 0.26981, 0.26902, 0.26814, 0.26785, 0.26811, 0.26976, 0.27228, 0.27427, 0.28078, 0.28616, 0.29929, 0.31335, 0.32473]
        ],
        [
            [0.26868, 0.27355, 0.27572, 0.27897, 0.28401, 0.29064, 0.29518, 0.30312, 0.31405, 0.32471, 0.33181, 0.35233, 0.36795, 0.40442, 0.44282, 0.47762],
            [0.26868, 0.26865, 0.26864, 0.26864, 0.26881, 0.26960, 0.27062, 0.27331, 0.27837, 0.28415, 0.28828, 0.30085, 0.31077, 0.33449, 0.35998, 0.38391],
            [0.26868, 0.26443, 0.26290, 0.26106, 0.25934, 0.25939, 0.26087, 0.26563, 0.27490, 0.28540, 0.29280, 0.31502, 0.33239, 0.37385, 0.41909, 0.46776]
        ],
        [
            [0.26977, 0.27849, 0.28337, 0.29166, 0.30666, 0.32947, 0.34640, 0.37706, 0.41944, 0.46019, 0.48705, 0.56379, 0.62178, 0.75766, 0.90507, 1.06945],
            [0.26977, 0.27046, 0.27119, 0.27275, 0.27632, 0.28318, 0.28915, 0.30137, 0.32028, 0.33976, 0.35302, 0.39193, 0.42193, 0.49318, 0.57123, 0.65836],
            [0.26977, 0.26536, 0.26492, 0.26589, 0.27116, 0.28448, 0.29699, 0.32305, 0.36302, 0.40351, 0.43078, 0.50990, 0.57036, 0.71331, 0.87044, 1.05478]
        ]
    ],
    # Case 3: τ = 0.25, (λ = 440 nm), ρ = 0.0
    [
        [
            [0.09775, 0.09765, 0.09758, 0.09756, 0.09802, 0.10032, 0.10331, 0.11101, 0.12470, 0.13925, 0.14893, 0.17503, 0.19223, 0.22240, 0.23988, 0.24541],
            [0.09775, 0.09765, 0.09758, 0.09756, 0.09802, 0.10032, 0.10331, 0.11101, 0.12470, 0.13925, 0.14893, 0.17503, 0.19223, 0.22240, 0.23988, 0.24541],
            [0.09775, 0.09765, 0.09758, 0.09756, 0.09802, 0.10032, 0.10331, 0.11101, 0.12470, 0.13925, 0.14893, 0.17503, 0.19223, 0.22240, 0.23988, 0.24541]
        ],
        [
            [0.09854, 0.11077, 0.11651, 0.12536, 0.13970, 0.15938, 0.17315, 0.19709, 0.22857, 0.25691, 0.27443, 0.31878, 0.34654, 0.39284, 0.41701, 0.41911],
            [0.09854, 0.09902, 0.09954, 0.10073, 0.10375, 0.11035, 0.11662, 0.13014, 0.15139, 0.17252, 0.18619, 0.22222, 0.24563, 0.28647, 0.31039, 0.31815],
            [0.09854, 0.08882, 0.08570, 0.08246, 0.08101, 0.08597, 0.09344, 0.11195, 0.14275, 0.17386, 0.19405, 0.24750, 0.28254, 0.34534, 0.38580, 0.40830]
        ],
        [
            [0.12483, 0.14482, 0.15623, 0.17577, 0.21126, 0.26488, 0.30406, 0.37292, 0.46272, 0.54253, 0.59155, 0.71594, 0.79516, 0.93532, 1.02629, 1.07987],
            [0.12483, 0.12707, 0.12939, 0.13430, 0.14534, 0.16580, 0.18297, 0.21654, 0.26486, 0.31056, 0.33949, 0.41484, 0.46385, 0.55189, 0.60965, 0.64215],
            [0.12483, 0.11563, 0.11527, 0.11871, 0.13315, 0.16710, 0.19778, 0.25918, 0.34774, 0.43090, 0.48327, 0.61925, 0.70784, 0.86948, 0.98212, 1.06400]
        ]
    ],
    # Case 4: τ = 0.25, (λ = 440 nm), ρ = 0.25
    [
        [
            [0.30216, 0.30941, 0.31233, 0.31637, 0.32204, 0.32875, 0.33312, 0.34069, 0.35103, 0.36074, 0.36686, 0.38248, 0.39207, 0.40655, 0.41048, 0.40057],
            [0.30216, 0.30181, 0.30148, 0.30089, 0.30002, 0.29967, 0.30031, 0.30317, 0.30956, 0.31696, 0.32203, 0.33586, 0.34489, 0.35978, 0.36608, 0.36139],
            [0.30216, 0.29480, 0.29182, 0.28782, 0.28299, 0.27991, 0.28010, 0.28406, 0.29397, 0.30555, 0.31346, 0.33517, 0.34962, 0.37511, 0.38988, 0.39347]
        ],
        [
            [0.29778, 0.31426, 0.32250, 0.33561, 0.35756, 0.38810, 0.40924, 0.44504, 0.49028, 0.52970, 0.55360, 0.61322, 0.65012, 0.71160, 0.74439, 0.74824],
            [0.29778, 0.29855, 0.29936, 0.30113, 0.30528, 0.31349, 0.32075, 0.33560, 0.35786, 0.37938, 0.39310, 0.42884, 0.45181, 0.49139, 0.51373, 0.51745],
            [0.29778, 0.28635, 0.28331, 0.28104, 0.28288, 0.29467, 0.30776, 0.33658, 0.38085, 0.42367, 0.45092, 0.52194, 0.56801, 0.65033, 0.70385, 0.73402]
        ],
        [
            [0.32057, 0.34195, 0.35634, 0.38273, 0.43408, 0.51642, 0.57872, 0.69088, 0.84057, 0.97630, 1.06099, 1.28186, 1.42882, 1.71233, 1.94102, 2.17632],
            [0.32057, 0.32427, 0.32809, 0.33611, 0.35394, 0.38641, 0.41326, 0.46520, 0.53946, 0.60987, 0.65476, 0.77392, 0.85412, 1.00911, 1.13208, 1.24962],
            [0.32057, 0.31777, 0.32239, 0.33541, 0.36924, 0.43508, 0.49017, 0.59577, 0.74397, 0.88185, 0.96899, 1.19862, 1.35283, 1.65336, 1.89994, 2.16054]
        ]
    ],
    # Case 5: τ = 0.50, (λ = 360 nm), ρ = 0.0
    [
        [
            [0.18734, 0.20137, 0.20744, 0.21632, 0.22977, 0.24706, 0.25869, 0.27821, 0.30219, 0.32145, 0.33202, 0.35358, 0.36290, 0.37017, 0.36753, 0.35823],
            [0.18734, 0.18747, 0.18765, 0.18821, 0.19020, 0.19583, 0.20183, 0.21531, 0.23594, 0.25474, 0.26575, 0.29009, 0.30211, 0.31567, 0.31882, 0.31623],
            [0.18734, 0.17461, 0.16996, 0.16436, 0.15939, 0.16066, 0.16636, 0.18268, 0.21033, 0.23658, 0.25232, 0.28844, 0.30757, 0.33310, 0.34470, 0.35054]
        ],
        [
            [0.19653, 0.22598, 0.24104, 0.26527, 0.30603, 0.36223, 0.40016, 0.46146, 0.53206, 0.58626, 0.61560, 0.67672, 0.70575, 0.73877, 0.74768, 0.74280],
            [0.19653, 0.19878, 0.20112, 0.20610, 0.21732, 0.23808, 0.25527, 0.28776, 0.33113, 0.36789, 0.38887, 0.43500, 0.45838, 0.48763, 0.49832, 0.49827],
            [0.19653, 0.17753, 0.17316, 0.17113, 0.17840, 0.20530, 0.23221, 0.28715, 0.36368, 0.42992, 0.46830, 0.55516, 0.60175, 0.66795, 0.70337, 0.72761]
        ],
        [
            [0.26496, 0.29541, 0.31615, 0.35420, 0.42791, 0.54424, 0.63025, 0.77987, 0.96802, 1.12673, 1.22007, 1.44405, 1.57909, 1.81703, 2.00445, 2.22063],
            [0.26496, 0.27112, 0.27745, 0.29066, 0.31960, 0.37085, 0.41186, 0.48782, 0.58912, 0.67767, 0.73051, 0.85810, 0.93470, 1.06678, 1.16606, 1.27144],
            [0.26496, 0.26195, 0.26922, 0.28899, 0.33910, 0.43410, 0.51149, 0.65474, 0.84430, 1.00909, 1.10743, 1.34664, 1.49277, 1.75350, 1.96156, 2.20439]
        ]
    ],
    # Case 6: τ = 0.50, (λ = 360 nm), ρ = 0.25
    [
        [
            [0.36020, 0.35915, 0.35814, 0.35621, 0.35283, 0.34915, 0.34788, 0.34828, 0.35225, 0.35726, 0.36044, 0.36727, 0.36996, 0.37015, 0.36572, 0.35639],
            [0.36020, 0.35915, 0.35814, 0.35621, 0.35283, 0.34915, 0.34788, 0.34828, 0.35225, 0.35726, 0.36044, 0.36727, 0.36996, 0.37015, 0.36572, 0.35639],
            [0.36020, 0.35915, 0.35814, 0.35621, 0.35283, 0.34915, 0.34788, 0.34828, 0.35225, 0.35726, 0.36044, 0.36727, 0.36996, 0.37015, 0.36572, 0.35639]
        ],
        [
            [0.35140, 0.37273, 0.38227, 0.39640, 0.41781, 0.44412, 0.46053, 0.48561, 0.51301, 0.53301, 0.54331, 0.56260, 0.56960, 0.57122, 0.56274, 0.54467],
            [0.35140, 0.35158, 0.35180, 0.35241, 0.35431, 0.35923, 0.36427, 0.37527, 0.39171, 0.40636, 0.41478, 0.43275, 0.44098, 0.44817, 0.44638, 0.43685],
            [0.35140, 0.33315, 0.32682, 0.31953, 0.31365, 0.31621, 0.32377, 0.34399, 0.37666, 0.40685, 0.42473, 0.46542, 0.48688, 0.51554, 0.52829, 0.53301]
        ],
        [
            [0.35266, 0.38470, 0.40275, 0.43327, 0.48745, 0.56599, 0.62061, 0.71069, 0.81666, 0.90002, 0.94627, 1.04741, 1.10019, 1.17386, 1.21231, 1.23573],
            [0.35266, 0.35608, 0.35961, 0.36699, 0.38322, 0.41211, 0.43526, 0.47792, 0.53371, 0.58077, 0.60779, 0.66857, 0.70092, 0.74601, 0.76813, 0.77720],
            [0.35266, 0.33743, 0.33652, 0.34136, 0.36268, 0.41221, 0.45569, 0.53880, 0.64951, 0.74376, 0.79837, 0.92371, 0.99327, 1.09940, 1.16471, 1.21897]
        ]
    ]
];