text:
"""
    FESetModule

Module for finite element sets.
"""
module FESetModule

__precompile__(true)

import Base.count
import Base.cat
import LinearAlgebra: mul!, Transpose
At_mul_B!(C, A, B) = mul!(C, Transpose(A), B)
import LinearAlgebra: norm, cross
using ..FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict

"""
    AbstractFESet{NODESPERELEM}

Abstract type of a finite element set. Parameterized with the number of of the nodes per element.
"""
abstract type AbstractFESet{NODESPERELEM} end
"""
    AbstractFESet0Manifold{NODESPERELEM} <: FESet{NODESPERELEM}

Abstract type of a finite element set for 0-dimensional manifolds (points).
Parameterized with the number of of the nodes per element.
"""
abstract type AbstractFESet0Manifold{NODESPERELEM} <: AbstractFESet{NODESPERELEM} end
"""
    AbstractFESet1Manifold{NODESPERELEM} <: FESet{NODESPERELEM}

Abstract type of a finite element set for 1-dimensional manifolds (curves).
Parameterized with the number of of the nodes per element.
"""
abstract type AbstractFESet1Manifold{NODESPERELEM} <: AbstractFESet{NODESPERELEM} end
"""
    AbstractFESet2Manifold{NODESPERELEM} <: FESet{NODESPERELEM}

Abstract type of a finite element set for 2-dimensional manifolds (surfaces).
Parameterized with the number of of the nodes per element.
"""
abstract type AbstractFESet2Manifold{NODESPERELEM} <: AbstractFESet{NODESPERELEM} end
"""
    AbstractFESet3Manifold{NODESPERELEM} <: FESet{NODESPERELEM}

Abstract type of a finite element set for 3-dimensional manifolds (solids).
Parameterized with the number of of the nodes per element.
"""
abstract type AbstractFESet3Manifold{NODESPERELEM} <: AbstractFESet{NODESPERELEM} end

"""
    add_FESet_fields(NODESPERELEM)

Generate standard fields for the finite element set.
"""
macro add_FESet_fields(NODESPERELEM)
    return esc(:(
    conn::Array{NTuple{$NODESPERELEM, FInt}, 1};
    label::FIntVec; )
    )
end
# show(macroexpand(:(@add_FESet_fields 6)))

"""
    define_FESet(NAME, MANIFOLD, NODESPERELEM)

Define the concrete type for a finite element set.
"""
macro define_FESet(NAME, MANIFOLD, NODESPERELEM)
    return esc(:(
        mutable struct $NAME <: $MANIFOLD{$NODESPERELEM}
            @add_FESet_fields $NODESPERELEM
            function $NAME(conn::FIntMat)
                self = new(NTuple{$NODESPERELEM, FInt}[], FInt[])
                self = fromarray!(self, conn)
                setlabel!(self, 0)
                return self
            end
        end
    ))
end
# show(macroexpand(:(@define_FESet FESetT10 AbstractFESet3Manifold 10)))

"""
    nodesperelem(fes::AbstractFESet{NODESPERELEM}) where {NODESPERELEM}

Provide the number of nodes per element.
"""
nodesperelem(fes::AbstractFESet{NODESPERELEM}) where {NODESPERELEM} = NODESPERELEM

"""
    nodesperelem(::Type{T}) where {NODESPERELEM, T<:AbstractFESet{NODESPERELEM}}

Provide the number of nodes per element for a given type.
"""
nodesperelem(::Type{T}) where {NODESPERELEM, T<:AbstractFESet{NODESPERELEM}} = NODESPERELEM

"""
    manifdim(me)

Get the manifold dimension.
"""
manifdim(me::AbstractFESet0Manifold{NODESPERELEM}) where {NODESPERELEM} = 0
manifdim(me::AbstractFESet1Manifold{NODESPERELEM}) where {NODESPERELEM} = 1
manifdim(me::AbstractFESet2Manifold{NODESPERELEM}) where {NODESPERELEM} = 2
manifdim(me::AbstractFESet3Manifold{NODESPERELEM}) where {NODESPERELEM} = 3

"""
    count(self::T)::FInt where {T<:AbstractFESet}

Get the number of individual connectivities in the FE set.
"""
count(self::T) where {T<:AbstractFESet} = length(self.conn)

"""
    fromarray!(self::AbstractFESet{NODESPERELEM}, conn::FIntMat) where {NODESPERELEM}

Set  the connectivity from an integer array.
"""
function fromarray!(self::AbstractFESet{NODESPERELEM}, conn::FIntMat) where {NODESPERELEM}
    @assert size(conn, 2) == NODESPERELEM
    self.conn = fill(tuple(fill(0, NODESPERELEM)...), size(conn, 1));
    for i = 1:length(self.conn)
        self.conn[i] = ntuple(y -> conn[i, y], NODESPERELEM);
    end
    return self
end

"""
    connasarray(self::AbstractFESet{NODESPERELEM}) where {NODESPERELEM}

Return the connectivity as an array.

Return the connectivity as an integer array (matrix), where the number of rows
matches the number of connectivities in the set.
"""
function connasarray(self::AbstractFESet{NODESPERELEM}) where {NODESPERELEM}
    conn = fill(zero(FInt), length(self.conn), NODESPERELEM)
    for i = 1:size(conn, 1)
        for j = 1:NODESPERELEM
            conn[i, j] = self.conn[i][j]
        end
    end
    return conn
end

"""
    boundaryconn(self::T) where {T<:AbstractFESet}

Get boundary connectivity.
"""
function boundaryconn(self::T) where {T<:AbstractFESet}
    return privboundaryconn(self)
end

"""
    boundaryfe(self::T) where {T<:AbstractFESet}

Return the constructor of the type of the boundary finite element.
"""
function boundaryfe(self::T) where {T<:AbstractFESet}
    return privboundaryfe(self)
end

"""
    bfun(self::T,  param_coords::FFltVec)::FFltMat where {T<:AbstractFESet}

Compute the values of the basis functions.

Compute the values of the basis functions at a given parametric coordinate.
One basis function per row.
"""
function bfun(self::T,  param_coords::FFltVec)::FFltMat where {T<:AbstractFESet}
    return privbfun(self, param_coords)
end

"""
    bfundpar(self::T,  param_coords::FFltVec)::FFltMat where {T<:AbstractFESet}

Compute the values of the basis function gradients.

Compute the values of the basis function gradients with respect to the
parametric coordinates at a given parametric coordinate. One basis function
gradients per row.
"""
function bfundpar(self::T,  param_coords::FFltVec)::FFltMat where {T<:AbstractFESet}
    return privbfundpar(self, param_coords)
end

"""
    setlabel!(self::T, val::FInt) where {T<:AbstractFESet}

Set the label of the entire finite elements set.

All elements are labeled with this number.
"""
function setlabel!(self::T, val::FInt) where {T<:AbstractFESet}
    self.label = zeros(FInt, size(self.conn, 1));
    fill!(self.label, val)
    return self
end

"""
    setlabel!(self::T, val::FIntVec) where {T<:AbstractFESet}

Set the labels of individual elements.
"""
function setlabel!(self::T, val::FIntVec) where {T<:AbstractFESet}
    #    Set the label of this set.
    @assert size(self.conn, 1)==length(val) "Must get one  label per finite element connectivity"
    self.label=zeros(FInt, size(self.conn, 1));
    copyto!(self.label, val);
    return self
end

"""
    subset(self::T, L::FIntVec) where {T<:AbstractFESet}

Extract a subset of the finite elements from the given finite element set.
"""
function subset(self::T, L::FIntVec) where {T<:AbstractFESet}
    result = deepcopy(self)
    result.conn = deepcopy(self.conn[L])
    result.label = deepcopy(self.label[L])
    return  result
end

"""
    cat(self::T,  other::T) where {T<:AbstractFESet}

Concatenate the connectivities of two FE sets.
"""
function cat(self::T,  other::T) where {T<:AbstractFESet}
    @assert nodesperelem(self) == nodesperelem(other)
    result =deepcopy(self)
    result.conn = vcat(self.conn, other.conn);
    setlabel!(result, vcat(self.label, other.label))
    return result
end

"""
    updateconn!(self::T, newids::FIntVec) where {T<:AbstractFESet}

Update the connectivity after the IDs of nodes changed.

`newids`= new node IDs. Note that indexes in the conn array "point"
_into_ the  `newids` array. After the connectivity was updated
this will no longer be true!
"""
function updateconn!(self::T, newids::FIntVec) where {T<:AbstractFESet}
    conn = connasarray(self)
    for i=1:size(conn, 1)
        for j=1:size(conn, 2)
            conn[i, j]=newids[conn[i, j]];
        end
    end
    return fromarray!(self, conn)
end

"""
    inparametric(self::AbstractFESet, param_coords::FFltVec)

Are given parametric coordinates inside the element parametric domain?

Return a Boolean: is the point inside, true or false?
"""
function inparametric(self::T, param_coords::FFltVec; tolerance = 0.0) where {T<:AbstractFESet}
    return privinparametric(self, param_coords, tolerance)
end

"""
    map2parametric(self::T, x::FFltMat, pt::FFltVec;
        tolerance = 0.001, maxiter =5) where {T<:AbstractFESet}

Map a spatial location to parametric coordinates.

- `x`=array of spatial coordinates of the nodes, size(x) = nbfuns x dim,
- `c`= spatial location
- `tolerance` = tolerance in parametric coordinates; default is 0.001.

# Return
- `success` = Boolean flag, true if successful, false otherwise.
- `pc` = Returns a row array of parametric coordinates if the solution was
  successful, otherwise NaN are returned.
"""
function map2parametric(self::T, x::FFltMat, pt::FFltVec; tolerance = 0.001, maxiter =5) where {T<:AbstractFESet}
    sdim = size(x, 2); # number of space dimensions
    mdim = manifdim(self); # manifold dimension of the element
    ppc = zeros(mdim);
    pc = deepcopy(ppc);
    J = [i==j ? one(FFlt) : zero(FFlt) for i=1:sdim, j=1:mdim] # Jac. matx: "identity"
    loc = fill(zero(FFlt), (1, sdim)) # point location buffer
    success = true
    for i = 1:maxiter
        gradNparams = privbfundpar(self, pc);
        At_mul_B!(J, x, gradNparams); # calculate the Jacobian matrix
        N = privbfun(self, pc);
        At_mul_B!(loc, N, x);# Iterated point location
        pc[:] = ppc .- (J\(vec(loc) - pt))
        if (norm(pc-ppc) < tolerance)
            return pc, success;
        end
        ppc[:] = pc[:];
    end
    success = false
    return pc, success;
end

"""
    centroidparametric(self::T) where {T<:AbstractFESet}

Return the parametric coordinates  of the centroid of the element.
"""
function centroidparametric(self::T) where {T<:AbstractFESet}
    return privcentroidparametric(self)
end

"""
    Jacobian(self::T, J::FFltMat)::FFlt where {T<:AbstractFESet0Manifold}

Evaluate the point Jacobian.

- `J` = Jacobian matrix, columns are tangent to parametric coordinates curves.
"""
function Jacobian(self::T, J::FFltMat)::FFlt where {T<:AbstractFESet0Manifold}
    return 1.0::FFlt;
end

"""
    Jacobian(self::T, J::FFltMat)::FFlt where {T<:AbstractFESet1Manifold}

Evaluate the curve Jacobian.

- `J` = Jacobian matrix, columns are tangent to parametric coordinates curves.
"""
function Jacobian(self::T, J::FFltMat)::FFlt where {T<:AbstractFESet1Manifold}
    sdim,  ntan = size(J);
    @assert ntan == 1 "Expected number of tangent vectors: 1"
    return norm(J)::FFlt;
end

"""
    gradN!(self::AbstractFESet1Manifold, gradN::FFltMat, gradNparams::FFltMat,
      redJ::FFltMat)

Compute the gradient of the basis functions with the respect to
the "reduced" spatial coordinates.

- `gradN`= output,  matrix of gradients,  one per row
- `gradNparams`= matrix of gradients with respect to parametric coordinates, one per row
- `redJ`= reduced Jacobian matrix `redJ=transpose(Rm)*J`
"""
function gradN!(self::AbstractFESet1Manifold, gradN::FFltMat, gradNparams::FFltMat, redJ::FFltMat)
    for r=1:size(gradN, 1)
        gradN[r, 1]=  gradNparams[r, 1]/redJ[1, 1]
    end
end

"""
    Jacobian(self::T, J::FFltMat)::FFlt where {T<:AbstractFESet2Manifold}

Evaluate the curve Jacobian.

- `J` = Jacobian matrix, columns are tangent to parametric coordinates curves.
"""
function Jacobian(self::T, J::FFltMat)::FFlt where {T<:AbstractFESet2Manifold}
    sdim,  ntan = size(J);
    @assert ntan == 2 "Expected number of tangent vectors: 2"
    if sdim == ntan
        @inbounds Jac = (J[1, 1]*J[2, 2] - J[2, 1]*J[1, 2])
        return Jac::FFlt;# is det(J);% Compute the Jacobian
    else
        return norm(cross(J[:, 1], J[:, 2]))::FFlt;
    end
end

"""
    gradN!(self::AbstractFESet2Manifold, gradN::FFltMat, gradNparams::FFltMat,
      redJ::FFltMat)

Compute the gradient of the basis functions with the respect to
the "reduced" spatial coordinates.

- `gradN`= output,  matrix of gradients,  one per row
- `gradNparams`= matrix of gradients with respect to parametric coordinates,
  one per row
- `redJ`= reduced Jacobian matrix `redJ=transpose(Rm)*J`
"""
function gradN!(self::AbstractFESet2Manifold, gradN::FFltMat, gradNparams::FFltMat, redJ::FFltMat)
    # This is the unrolled version that avoids allocation of a 3 x 3 matrix
    invdet=1.0/(redJ[1, 1]*redJ[2, 2] - redJ[1, 2]*redJ[2, 1]);
    invredJ11 =  (redJ[2, 2])*invdet;
    invredJ12 = -(redJ[1, 2])*invdet;
    invredJ21 = -(redJ[2, 1])*invdet;
    invredJ22 =  (redJ[1, 1])*invdet;
    @assert size(gradN, 1)==size(gradNparams, 1)
    @inbounds for r=1:size(gradN, 1)
        gradN[r, 1]= gradNparams[r, 1]*invredJ11 +gradNparams[r, 2]*invredJ21;
        gradN[r, 2]= gradNparams[r, 1]*invredJ12 +gradNparams[r, 2]*invredJ22;
    end
end

"""
    Jacobian(self::T, J::FFltMat)::FFlt where {T<:AbstractFESet3Manifold}

Evaluate the volume Jacobian.

`J` = Jacobian matrix, columns are tangent to parametric coordinates curves.
"""
function Jacobian(self::T, J::FFltMat)::FFlt where {T<:AbstractFESet3Manifold}
    sdim,  ntan = size(J);
    @assert (ntan == 3) && (sdim == 3) "Expected number of tangent vectors: 3"
    #Jac = det(J);# Compute the Jacobian
    # The unrolled version
    return (+J[1, 1]*(J[2, 2]*J[3, 3]-J[3, 2]*J[2, 3])
    -J[1, 2]*(J[2, 1]*J[3, 3]-J[2, 3]*J[3, 1])
    +J[1, 3]*(J[2, 1]*J[3, 2]-J[2, 2]*J[3, 1]) )::FFlt;
end

"""
    gradN!(self::AbstractFESet3Manifold, gradN::FFltMat, gradNparams::FFltMat,
      redJ::FFltMat)

Compute the gradient of the basis functions with the respect to
the "reduced" spatial coordinates.

- `gradN`= output,  matrix of gradients,  one per row
- `gradNparams`= matrix of gradients with respect to parametric coordinates,
  one per row
- `redJ`= reduced Jacobian matrix `redJ=transpose(Rm)*J`
"""
function gradN!(self::AbstractFESet3Manifold, gradN::FFltMat, gradNparams::FFltMat, redJ::FFltMat)
    invdet = 1.0 / ( +redJ[1, 1]*(redJ[2, 2]*redJ[3, 3]-redJ[3, 2]*redJ[2, 3])
                    -redJ[1, 2]*(redJ[2, 1]*redJ[3, 3]-redJ[2, 3]*redJ[3, 1])
                    +redJ[1, 3]*(redJ[2, 1]*redJ[3, 2]-redJ[2, 2]*redJ[3, 1]) );
    # This is the unrolled version that avoids allocation of a 3 x 3 matrix
    invredJ11 =  (redJ[2, 2]*redJ[3, 3]-redJ[3, 2]*redJ[2, 3])*invdet;
    invredJ12 = -(redJ[1, 2]*redJ[3, 3]-redJ[1, 3]*redJ[3, 2])*invdet;
    invredJ13 =  (redJ[1, 2]*redJ[2, 3]-redJ[1, 3]*redJ[2, 2])*invdet;
    invredJ21 = -(redJ[2, 1]*redJ[3, 3]-redJ[2, 3]*redJ[3, 1])*invdet;
    invredJ22 =  (redJ[1, 1]*redJ[3, 3]-redJ[1, 3]*redJ[3, 1])*invdet;
    invredJ23 = -(redJ[1, 1]*redJ[2, 3]-redJ[2, 1]*redJ[1, 3])*invdet;
    invredJ31 =  (redJ[2, 1]*redJ[3, 2]-redJ[3, 1]*redJ[2, 2])*invdet;
    invredJ32 = -(redJ[1, 1]*redJ[3, 2]-redJ[3, 1]*redJ[1, 2])*invdet;
    invredJ33 =  (redJ[1, 1]*redJ[2, 2]-redJ[2, 1]*redJ[1, 2])*invdet;
    @assert size(gradN, 1)==size(gradNparams, 1)
    @inbounds for r=1:size(gradN, 1)
        gradN[r, 1]= gradNparams[r, 1]*invredJ11 +gradNparams[r, 2]*invredJ21 +gradNparams[r, 3]*invredJ31;
        gradN[r, 2]= gradNparams[r, 1]*invredJ12 +gradNparams[r, 2]*invredJ22 +gradNparams[r, 3]*invredJ32;
        gradN[r, 3]= gradNparams[r, 1]*invredJ13 +gradNparams[r, 2]*invredJ23 +gradNparams[r, 3]*invredJ33;
    end
end

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

"""
    FESetP1

Type for sets of point-like of finite elements.
"""
@define_FESet FESetP1 AbstractFESet0Manifold 1

privbfun(self::FESetP1,  param_coords::FFltVec) = reshape([1.0], 1, 1) # make sure this is a matrix
privbfundpar(self::FESetP1,  param_coords::FFltVec) = zeros(1, 0)

function privinparametric(self::FESetP1, param_coords::FFltVec, tolerance::FFlt)
    return  param_coords[1] == 0.0;
end

function privcentroidparametric(self::FESetP1)
    return vec([0.0])
end

function map2parametric(self::FESetP1, x::FFltMat, pt::FFltVec;
    tolerance = 0.001, maxiter =5)
    success = false; pc = [0.0]
    if norm(vec(x) - pt) < tolerance
        success = true
    end
    return pc, success;
end

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

"""
    FESetL2

Type for sets of curve-like finite elements with two nodes.
"""
@define_FESet FESetL2 AbstractFESet1Manifold 2

privbfun(self::FESetL2,  param_coords::FFltVec) = reshape([(1. - param_coords[1]); (1. + param_coords[1])] / 2.0, 2, 1) # make sure this is a matrix
privbfundpar(self::FESetL2,  param_coords::FFltVec) = reshape([-1.0; +1.0]/2.0, 2, 1)

function privboundaryconn(self::FESetL2)
    conn = connasarray(self)
    return [conn[:, 1]; conn[:, 2]];
end

function privboundaryfe(self::FESetL2)
    return FESetP1;
end

function privinparametric(self::FESetL2, param_coords::FFltVec, tolerance::FFlt)
    s = findall(v -> (v >= (-1.0 - tolerance)) && (v <= (+1.0 + tolerance)), param_coords);
    return  (length(s) == length(param_coords));
end

function privcentroidparametric(self::FESetL2)
    return vec([0.0])
end

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

"""
    FESetL3

Type for sets of curve-like of finite elements with three nodes.
"""
@define_FESet FESetL3 AbstractFESet1Manifold 3

function privbfun(self::FESetL3,  param_coords::FFltVec)
    xi=param_coords[1];
    val = [(xi-1)*xi/2;  (xi+1)*xi/2;  -(xi+1)*(xi-1)];
    return reshape(val, 3, 1) # make sure this is a matrix
end

function privbfundpar(self::FESetL3,  param_coords::FFltVec)
    xi=param_coords[1];
    val = [(xi-1/2); (xi+1/2); -2*xi];
    return reshape(val, 3, 1) # make sure this is a matrix
end

function privboundaryconn(self::FESetL3)
    conn = connasarray(self)
    return [conn[:, 1]; conn[:, 2]];
end

function privboundaryfe(self::FESetL3)
    return FESetP1;
end

function privinparametric(self::FESetL3, param_coords::FFltVec, tolerance::FFlt)
    s = findall(v -> (v >= (-1.0 - tolerance)) && (v <= (+1.0 + tolerance)), param_coords);
    return  (length(s) == length(param_coords));
end

function privcentroidparametric(self::FESetL3)
    return vec([0.0])
end

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

"""
    FESetT3

Type for sets of surface-like triangular finite elements with three nodes.
"""
@define_FESet FESetT3 AbstractFESet2Manifold 3

function privbfun(self::FESetT3,  param_coords::FFltVec)
    # Evaluate the basis function matrix for an 3-node triangle.
    return reshape([(1 - param_coords[1] - param_coords[2]); param_coords[1]; param_coords[2]], 3, 1); # Make sure this is a matrix
end

function privbfundpar(self::FESetT3,  param_coords::FFltVec)
    # Evaluate the derivatives of the basis function matrix.
    return  [-1. -1.;  +1.  0.;  0. +1.];
end

function privboundaryconn(self::FESetT3)
    conn = connasarray(self)
    return [conn[:, 1:2]; conn[:, 2:3]; conn[:, [3, 1]]];
end

function privboundaryfe(self::FESetT3)
    # Get  the constructor of the class of the  boundary finite element.
    return FESetL2;
end

function privinparametric(self::FESetT3, param_coords::FFltVec, tolerance::FFlt)
    s = findall(v -> (v >= (-tolerance)) && (v <= (+1.0 + tolerance)), param_coords);
    return  (length(s) == length(param_coords)) &&
        (-3.0 * tolerance <= sum(param_coords) <= 1.0 + 3.0 * tolerance);
end

function privcentroidparametric(self::FESetT3)
    return vec([1.0/3 1.0/3])
end

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

"""
    FESetQ4

Type for sets of surface-like quadrilateral finite elements with four nodes.
"""
@define_FESet FESetQ4 AbstractFESet2Manifold 4

function privbfun(self::FESetQ4,  param_coords::FFltVec)
    # Evaluate the basis function matrix for an 4-node quadrilateral.
    val = [0.25 * (1. - param_coords[1]) * (1. - param_coords[2]);
           0.25 * (1. + param_coords[1]) * (1. - param_coords[2]);
           0.25 * (1. + param_coords[1]) * (1. + param_coords[2]);
           0.25 * (1. - param_coords[1]) * (1. + param_coords[2])];
    return reshape(val, 4, 1); # Make sure this is a matrix
end

function privbfundpar(self::FESetQ4,  param_coords::FFltVec)
    # Evaluate the derivatives of the basis function matrix.
    val = [-(1. - param_coords[2])*0.25 -(1. - param_coords[1])*0.25;
            (1. - param_coords[2])*0.25 -(1. + param_coords[1])*0.25;
            (1. + param_coords[2])*0.25 (1. + param_coords[1])*0.25;
           -(1. + param_coords[2])*0.25 (1. - param_coords[1])*0.25];
    return val
end

function privboundaryconn(self::FESetQ4)
    conn = connasarray(self)
    return [conn[:, 1:2]; conn[:, 2:3]; conn[:, 3:4]; conn[:, [4, 1]]];
end

function privboundaryfe(self::FESetQ4)
    # Get  the constructor of the class of the  boundary finite element.
    return FESetL2;
end

function privinparametric(self::FESetQ4, param_coords::FFltVec, tolerance::FFlt)
    s = findall(v -> (v >= (-1.0 - tolerance)) && (v <= (+1.0 + tolerance)), param_coords);
    return  (length(s) == length(param_coords));
end

function privcentroidparametric(self::FESetQ4)
    return vec([0.0 0.0])
end

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

"""
    FESetQ9

Type for sets of surface-like quadrilateral finite elements with nine nodes.
"""
@define_FESet FESetQ9 AbstractFESet2Manifold 9

function privbfun(self::FESetQ9,  param_coords::FFltVec)
    # Evaluate the basis function matrix for an 9-node quadrilateral.
    xi=param_coords[1];
    xis = [(xi-1)*xi/2;  (xi+1)*xi/2;  -(xi+1)*(xi-1)];
    eta=param_coords[2];
    etas = [(eta-1)*eta/2;  (eta+1)*eta/2;  -(eta+1)*(eta-1)];
    xisetas=(xis*etas');
    val = xisetas[[1     2     5     4     3     8     6     7     9]];
    return reshape(val, 9, 1); # Make sure this is a matrix
end

function privbfundpar(self::FESetQ9,  param_coords::FFltVec)
    # Evaluate the derivatives of the basis function matrix.
    xi=param_coords[1];
    xis = [(xi-1)*xi/2;  (xi+1)*xi/2;  -(xi+1)*(xi-1)];
    dxis = [(xi-1/2); (xi+1/2); -2*xi];
    eta=param_coords[2];
    etas = [(eta-1)*eta/2;  (eta+1)*eta/2;  -(eta+1)*(eta-1)];
    detas=[(eta-1/2); (eta+1/2); -2*eta];
    dxisetas=(dxis*etas');
    xisdetas=(xis*detas');
    ix=[1     2     5     4     3     8     6     7     9]
    val =  [dxisetas[ix]'             xisdetas[ix]'];
    return reshape(val, 9, 2); # Make sure this is a matrix
end

function privboundaryconn(self::FESetQ9)
    conn = connasarray(self)
    return [conn[:, [1, 2, 5]]; conn[:, [2, 3, 6]]; conn[:, [3, 4, 7]]; conn[:, [4, 1, 8]]];
end

function privboundaryfe(self::FESetQ9)
    return FESetL3;
end

function privinparametric(self::FESetQ9, param_coords::FFltVec, tolerance::FFlt)
    s = findall(v -> (v >= (-1.0 - tolerance)) && (v <= (+1.0 + tolerance)), param_coords);
    return  (length(s) == length(param_coords));
end

function privcentroidparametric(self::FESetQ9)
    return vec([0.0 0.0])
end

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

"""
    FESetQ8

Type for sets of surface-like quadrilateral finite elements with eight nodes.
"""
@define_FESet FESetQ8 AbstractFESet2Manifold 8

function privbfun(self::FESetQ8,  param_coords::FFltVec)
    # Evaluate the basis function matrix for an 4-node quadrilateral.
    xim    = (-1 + param_coords[1]);
    etam   = (-1 + param_coords[2]);
    xip    = (1 + param_coords[1]);
    etap   = (1 + param_coords[2]);
    val = [ -1.0/4*xim*etam*(1+param_coords[1]+param_coords[2]);
           1.0/4*xip*etam*(1-param_coords[1]+param_coords[2]);
           -1.0/4*xip*etap*(1-param_coords[1]-param_coords[2]);
           1.0/4*xim*etap*(1+param_coords[1]-param_coords[2]);
           1.0/2*xim*xip*etam;
           -1.0/2*etam*etap*xip;
           -1.0/2*xim*xip*etap;
           1.0/2*etam*etap*xim];
    return reshape(val, 8, 1); # Make sure this is a matrix
end

function privbfundpar(self::FESetQ8,  param_coords::FFltVec)
    # Evaluate the derivatives of the basis function matrix.
    xi =param_coords[1];
    eta =param_coords[2];
    val=zeros(FFlt, 8, 2)
    val[:, 1]= [1.0/4*(1-eta)*(1+xi+eta)-1.0/4*(1-xi)*(1-eta);
               -1.0/4*(1-eta)*(1-xi+eta)+1.0/4*(1+xi)*(1-eta);
               -1.0/4*(1+eta)*(1-xi-eta)+1.0/4*(1+xi)*(1+eta);
               1.0/4*(1+eta)*(1+xi-eta)-1.0/4*(1-xi)*(1+eta);
               -1.0/2*(1+xi)*(1-eta)+1.0/2*(1-xi)*(1-eta);
               1.0/2*(1-eta)*(1+eta);
               -1.0/2*(1+xi)*(1+eta)+1.0/2*(1-xi)*(1+eta);
               -1.0/2*(1-eta)*(1+eta)];
    val[:, 2] = [1.0/4*(1-xi)*(1+xi+eta)-1.0/4*(1-xi)*(1-eta);
                1.0/4*(1+xi)*(1-xi+eta)-1.0/4*(1+xi)*(1-eta);
                -1.0/4*(1+xi)*(1-xi-eta)+1.0/4*(1+xi)*(1+eta);
                -1.0/4*(1-xi)*(1+xi-eta)+1.0/4*(1-xi)*(1+eta);
                -1.0/2*(1-xi)*(1+xi);
                -1.0/2*(1+xi)*(1+eta)+1.0/2*(1+xi)*(1-eta);
                1.0/2*(1-xi)*(1+xi);
                -1.0/2*(1-xi)*(1+eta)+1.0/2*(1-xi)*(1-eta)];
    return  reshape(val, 8, 2); # Make sure this is a matrix
end

function privboundaryconn(self::FESetQ8)
    conn = connasarray(self)
    return [conn[:, [1, 2, 5]]; conn[:, [2, 3, 6]]; conn[:, [3, 4, 7]]; conn[:, [4, 1, 8]]];
end

function privboundaryfe(self::FESetQ8)
    return FESetL3;
end

function privinparametric(self::FESetQ8, param_coords::FFltVec, tolerance::FFlt)
    s = findall(v -> (v >= (-1.0 - tolerance)) && (v <= (+1.0 + tolerance)), param_coords);
    return  (length(s) == length(param_coords));
end

function privcentroidparametric(self::FESetQ8)
    return vec([0.0 0.0])
end

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

"""
    FESetT6

Type for sets of surface-like triangular finite elements with six nodes.
"""
@define_FESet FESetT6 AbstractFESet2Manifold 6

function privbfun(self::FESetT6,  param_coords::FFltVec)
    # Evaluate the basis function matrix for an 4-node quadrilateral.
    r=param_coords[1];
    s=param_coords[2];
    t = 1. - r - s;
    val = [t * (t + t - 1);
           r * (r + r - 1);
           s * (s + s - 1);
           4 * r * t;
           4 * r * s;
           4 * s * t];
    return reshape(val, 6, 1); # Make sure this is a matrix
end

function privbfundpar(self::FESetT6,  param_coords::FFltVec)
    # Evaluate the derivatives of the basis function matrix.
    r =param_coords[1];
    s =param_coords[2];
    t = 1. - r - s;
    val = [-3+4*r+4*s  -3+4*r+4*s;
           4*r-1  0.0;
           0.0 4*s-1;
           4-8*r-4*s  -4*r;
           4*s  4*r;
           -4*s  4-4*r-8*s];
    return  reshape(val, 6, 2); # Make sure this is a matrix
end

function privboundaryconn(self::FESetT6)
    conn = connasarray(self)
    return [conn[:, [1,  2,  4]]; conn[:, [2,  3,  5]]; conn[:, [3,  1,  6]]];
end

function privboundaryfe(self::FESetT6)
    # Get  the constructor of the class of the  boundary finite element.
    return FESetL3;
end

function privinparametric(self::FESetT6, param_coords::FFltVec, tolerance::FFlt)
    s = findall(v -> (v >= (-tolerance)) && (v <= (+1.0 + tolerance)), param_coords);
    return  (length(s) == length(param_coords)) &&
      (-3.0 * tolerance <= sum(param_coords) <= 1.0 + 3.0 * tolerance);
end

function privcentroidparametric(self::FESetT6)
    return vec([1.0/3 1.0/3])
end

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

"""
    FESetH8

Type for sets of volume-like hexahedral finite elements with eight nodes.
"""
@define_FESet FESetH8 AbstractFESet3Manifold 8

function privbfun(self::FESetH8,  param_coords::FFltVec)
    # Evaluate the basis function matrix for an 8-node hexahedron.
    one_minus_xi    = (1.0 - param_coords[1]);
    one_minus_eta   = (1.0 - param_coords[2]);
    one_minus_theta = (1.0 - param_coords[3]);
    one_plus_xi     = (1.0 + param_coords[1]);
    one_plus_eta    = (1.0 + param_coords[2]);
    one_plus_theta  = (1.0 + param_coords[3]);
    val = [one_minus_xi*one_minus_eta*one_minus_theta;
           one_plus_xi*one_minus_eta*one_minus_theta;
           one_plus_xi*one_plus_eta*one_minus_theta;
           one_minus_xi*one_plus_eta*one_minus_theta;
           one_minus_xi*one_minus_eta*one_plus_theta;
           one_plus_xi*one_minus_eta*one_plus_theta;
           one_plus_xi*one_plus_eta*one_plus_theta;
           one_minus_xi*one_plus_eta*one_plus_theta] / 8.0;
    return reshape(val, 8, 1); # Make sure this is a matrix
end

function privbfundpar(self::FESetH8,  param_coords::FFltVec)
    # Evaluate the derivatives of the basis function matrix.
     omxi    = (1.0 - param_coords[1]);
    ometa   = (1.0 - param_coords[2]);
    omtheta = (1.0 - param_coords[3]);
    opxi     = (1.0 + param_coords[1]);
    opeta    = (1.0 + param_coords[2]);
    optheta  = (1.0 + param_coords[3]);
   val = [-ometa*omtheta  ometa*omtheta  opeta*omtheta  -opeta*omtheta  -ometa*optheta  ometa*optheta  opeta*optheta   -opeta*optheta;
          -omxi*omtheta  -opxi*omtheta   opxi*omtheta   omxi*omtheta    -omxi*optheta   -opxi*optheta  opxi*optheta   omxi*optheta;
          -omxi*ometa    -opxi*ometa     -opxi*opeta    -omxi*opeta      omxi*ometa     opxi*ometa     opxi*opeta    omxi*opeta]'/8.0;
    return val
end

function privboundaryconn(self::FESetH8)
    conn = connasarray(self)
    return  [conn[:, [1,  4,  3,  2]];
             conn[:, [1,  2,  6,  5]];
             conn[:, [2,  3,  7,  6]];
             conn[:, [3,  4,  8,  7]];
             conn[:, [4,  1,  5,  8]];
             conn[:, [6,  7,  8,  5]]];
end

function privboundaryfe(self::FESetH8)
    # Get  the constructor of the class of the  boundary finite element.
    return FESetQ4;
end

function privinparametric(self::FESetH8, param_coords::FFltVec, tolerance::FFlt)
    s = findall(v -> (v >= (-1.0 - tolerance)) && (v <= (+1.0 + tolerance)), param_coords);
    return  (length(s) == length(param_coords));
end

function privcentroidparametric(self::FESetH8)
    return vec([0.0 0.0 0.0])
end

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

"""
    FESetH20

Type for sets of volume-like hexahedral finite elements with 20 nodes.
"""
@define_FESet FESetH20 AbstractFESet3Manifold 20

function privbfun(self::FESetH20,  param_coords::FFltVec)
    # Evaluate the basis function matrix for an 20-node hexahedron.
    xim    = (-1 + param_coords[1]);
    etam   = (-1 + param_coords[2]);
    zetam  = (-1 + param_coords[3]);
    xip    = (1 + param_coords[1]);
    etap   = (1 + param_coords[2]);
    zetap  = (1 + param_coords[3]);
    val = [ 1.0/8*xim*etam*zetam*(2+param_coords[1]+param_coords[2]+param_coords[3]);
           -1.0/8*xip*etam*zetam*(2-param_coords[1]+param_coords[2]+param_coords[3]);
           1.0/8*xip*etap*zetam*(2-param_coords[1]-param_coords[2]+param_coords[3]);
           -1.0/8*xim*etap*zetam*(2+param_coords[1]-param_coords[2]+param_coords[3]);
           1.0/8*xim*etam*zetap*(-2-param_coords[1]-param_coords[2]+param_coords[3]);
           -1.0/8*xip*etam*zetap*(-2+param_coords[1]-param_coords[2]+param_coords[3]);
           1.0/8*xip*etap*zetap*(-2+param_coords[1]+param_coords[2]+param_coords[3]);
           -1.0/8*xim*etap*zetap*(-2-param_coords[1]+param_coords[2]+param_coords[3]);
           -1.0/4*xim*xip*etam*zetam;
           1.0/4*etam*etap*xip*zetam;
           1.0/4*xim*xip*etap*zetam;
           -1.0/4*etam*etap*xim*zetam;
           1.0/4*xim*xip*etam*zetap;
           -1.0/4*etam*etap*xip*zetap;
           -1.0/4*xim*xip*etap*zetap;
           1.0/4*etam*etap*xim*zetap;
           -1.0/4*zetam*zetap*xim*etam;
           1.0/4*zetam*zetap*xip*etam;
           -1.0/4*zetam*zetap*xip*etap;
           1.0/4*zetam*zetap*xim*etap];
    return reshape(val, 20, 1); # Make sure this is a matrix
end

function privbfundpar(self::FESetH20,  param_coords::FFltVec)
    # Evaluate the derivatives of the basis function matrix.
    xim    = -(-1 + param_coords[1]);
    etam   = -(-1 + param_coords[2]);
    zetam  = -(-1 + param_coords[3]);
    xip    = (1 + param_coords[1]);
    etap   = (1 + param_coords[2]);
    zetap  = (1 + param_coords[3]);
    twoppp =(2+param_coords[1]+param_coords[2]+param_coords[3]);
    twompp =(2-param_coords[1]+param_coords[2]+param_coords[3]);
    twopmp =(2+param_coords[1]-param_coords[2]+param_coords[3]);
    twoppm =(2+param_coords[1]+param_coords[2]-param_coords[3]);
    twommp =(2-param_coords[1]-param_coords[2]+param_coords[3]);
    twopmm =(2+param_coords[1]-param_coords[2]-param_coords[3]);
    twompm =(2-param_coords[1]+param_coords[2]-param_coords[3]);
    twommm =(2-param_coords[1]-param_coords[2]-param_coords[3]);
    val =zeros(FFlt, 20, 3)
   val[:, 1]= [1.0/8*etam*zetam*twoppp-1.0/8*xim*etam*zetam;
           -1.0/8*etam*zetam*twompp+1.0/8*xip*etam*zetam;
           -1.0/8*etap*zetam*twommp+1.0/8*xip*etap*zetam;
           1.0/8*etap*zetam*twopmp-1.0/8*xim*etap*zetam;
           1.0/8*etam*zetap*twoppm-1.0/8*xim*etam*zetap;
           -1.0/8*etam*zetap*twompm+1.0/8*xip*etam*zetap;
           -1.0/8*etap*zetap*twommm+1.0/8*xip*etap*zetap;
           1.0/8*etap*zetap*twopmm-1.0/8*xim*etap*zetap;
           -1.0/4*xip*etam*zetam+1.0/4*xim*etam*zetam;
           1.0/4*etam*etap*zetam;
           -1.0/4*xip*etap*zetam+1.0/4*xim*etap*zetam;
           -1.0/4*etam*etap*zetam;
           -1.0/4*xip*etam*zetap+1.0/4*xim*etam*zetap;
           1.0/4*etam*etap*zetap;
           -1.0/4*xip*etap*zetap+1.0/4*xim*etap*zetap;
           -1.0/4*etam*etap*zetap;
           -1.0/4*zetam*zetap*etam;
           1.0/4*zetam*zetap*etam;
           1.0/4*zetam*zetap*etap;
           -1.0/4*zetam*zetap*etap];
    val[:, 2]= [1.0/8*xim*zetam*twoppp-1.0/8*xim*etam*zetam;
             1.0/8*xip*zetam*twompp-1.0/8*xip*etam*zetam;
             -1.0/8*xip*zetam*twommp+1.0/8*xip*etap*zetam;
             -1.0/8*xim*zetam*twopmp+1.0/8*xim*etap*zetam;
             1.0/8*xim*zetap*twoppm-1.0/8*xim*etam*zetap;
             1.0/8*xip*zetap*twompm-1.0/8*xip*etam*zetap;
             -1.0/8*xip*zetap*twommm+1.0/8*xip*etap*zetap;
             -1.0/8*xim*zetap*twopmm+1.0/8*xim*etap*zetap;
             -1.0/4*xim*xip*zetam;
             -1.0/4*xip*etap*zetam+1.0/4*xip*etam*zetam;
             1.0/4*xim*xip*zetam;
             -1.0/4*xim*etap*zetam+1.0/4*xim*etam*zetam;
             -1.0/4*xim*xip*zetap;
             -1.0/4*xip*etap*zetap+1.0/4*xip*etam*zetap;
             1.0/4*xim*xip*zetap;
             -1.0/4*xim*etap*zetap+1.0/4*xim*etam*zetap;
             -1.0/4*zetam*zetap*xim;
             -1.0/4*zetam*zetap*xip;
             1.0/4*zetam*zetap*xip;
             1.0/4*zetam*zetap*xim];
    val[:, 3]= [1.0/8*xim*etam*twoppp-1.0/8*xim*etam*zetam;
             1.0/8*xip*etam*twompp-1.0/8*xip*etam*zetam;
             1.0/8*xip*etap*twommp-1.0/8*xip*etap*zetam;
             1.0/8*xim*etap*twopmp-1.0/8*xim*etap*zetam;
             -1.0/8*xim*etam*twoppm+1.0/8*xim*etam*zetap;
             -1.0/8*xip*etam*twompm+1.0/8*xip*etam*zetap;
             -1.0/8*xip*etap*twommm+1.0/8*xip*etap*zetap;
             -1.0/8*xim*etap*twopmm+1.0/8*xim*etap*zetap;
             -1.0/4*xim*xip*etam;
             -1.0/4*etam*etap*xip;
             -1.0/4*xim*xip*etap;
             -1.0/4*etam*etap*xim;
             1.0/4*xim*xip*etam;
             1.0/4*etam*etap*xip;
             1.0/4*xim*xip*etap;
             1.0/4*etam*etap*xim;
             -1.0/4*xim*etam*zetap+1.0/4*xim*etam*zetam;
             -1.0/4*xip*etam*zetap+1.0/4*xip*etam*zetam;
             -1.0/4*xip*etap*zetap+1.0/4*xip*etap*zetam;
             -1.0/4*xim*etap*zetap+1.0/4*xim*etap*zetam];
    return reshape(val, 20, 3);
end

function privboundaryconn(self::FESetH20)
    conn = connasarray(self)
    return  [conn[:, [1,  4,  3,  2,  12,  11,  10,  9]];
             conn[:, [1,  2,  6,  5,  9,  18,  13,  17]];
             conn[:, [2,  3,  7,  6,  10,  19,  14,  18]];
             conn[:, [3,  4,  8,  7,   11,  20,  15,  19]];
             conn[:, [4,  1,  5,  8,  12,  17,  16,  20]];
             conn[:, [6,  7,  8,  5,  14,  15,  16,  13]]];
end

function privboundaryfe(self::FESetH20)
    # Get  the constructor of the class of the  boundary finite element.
    return FESetQ8;
end

function privinparametric(self::FESetH20, param_coords::FFltVec, tolerance::FFlt)
    s = findall(v -> (v >= (-1.0 - tolerance)) && (v <= (+1.0 + tolerance)), param_coords);
    return  (length(s) == length(param_coords));
end

function privcentroidparametric(self::FESetH20)
    return vec([0.0 0.0 0.0])
end

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

"""
    FESetH27

Type for sets of volume-like hexahedral finite elements with 27 nodes.
"""
@define_FESet FESetH27 AbstractFESet3Manifold 27

function privbfun(self::FESetH27,  param_coords::FFltVec)
    # Evaluate the basis function matrix for an 8-node hexahedron.
    xi=param_coords[1];
    eta=param_coords[2];
    zet=param_coords[3];
    x1 =(xi-1);
    y1=(eta-1);
    z1 =(zet-1);
    x2 =(xi+1);
    y2=(eta+1);
    z2 =(zet+1);
    val =[1.0/8.0*z1*zet*x1*xi*y1*eta
          1.0/8.0*z1*zet*x2*xi*y1*eta
          1.0/8.0*z1*zet*x2*xi*y2*eta
          1.0/8.0*z1*zet*x1*xi*y2*eta
          1.0/8.0*z2*zet*x1*xi*y1*eta
          1.0/8.0*z2*zet*x2*xi*y1*eta
          1.0/8.0*z2*zet*x2*xi*y2*eta
          1.0/8.0*z2*zet*x1*xi*y2*eta
          1.0/4.0*z1*zet*(-x2)*x1*y1*eta
          1.0/4.0*z1*zet*x2*xi*(-y2)*y1
          1.0/4.0*z1*zet*(-x2)*x1*y2*eta
          1.0/4.0*z1*zet*x1*xi*(-y2)*y1
          1.0/4.0*z2*zet*(-x2)*x1*y1*eta
          1.0/4.0*z2*zet*x2*xi*(-y2)*y1
          1.0/4.0*z2*zet*(-x2)*x1*y2*eta
          1.0/4.0*z2*zet*x1*xi*(-y2)*y1
          1.0/4.0*(-z2)*z1*x1*xi*y1*eta
          1.0/4.0*(-z2)*z1*x2*xi*y1*eta
          1.0/4.0*(-z2)*z1*x2*xi*y2*eta
          1.0/4.0*(-z2)*z1*x1*xi*y2*eta
          1.0/2.0*z1*zet*(-x2)*x1*(-y2)*y1
          1.0/2.0*(-z2)*z1*(-x2)*x1*y1*eta
          1.0/2.0*(-z2)*z1*x2*xi*(-y2)*y1
          1.0/2.0*(-z2)*z1*(-x2)*x1*y2*eta
          1.0/2.0*(-z2)*z1*x1*xi*(-y2)*y1
          1.0/2.0*z2*zet*(-x2)*x1*(-y2)*y1
          (-z2)*z1*(-x2)*x1*(-y2)*y1];
    return reshape(val, 27, 1); # Make sure this is a matrix
end

function privbfundpar(self::FESetH27,  param_coords::FFltVec)
    # Evaluate the derivatives of the basis function matrix.
    xi=param_coords[1];
    eta=param_coords[2];
    zet=param_coords[3];
    x1 =(xi-1.0/2.0);
    x2 =(xi+1.0/2.0);
    x3 =(xi-1.0);
    x4 =(xi+1.0);
    z1 =(zet-1.0);
    z2 =(zet-1.0/2.0);
    z3 =(zet+1.0);
    z4 =(zet+1.0/2.0);
    y1 =(eta-1.0);
    y2 =(eta-1.0/2.0);
    y3 =(eta+1.0);
    y4 = (eta+1.0/2.0);
    val = [      1.0/4.0*z1*zet*x1*y1*eta        1.0/4.0*z1*zet*x3*xi*y2        1.0/4.0*z2*x3*xi*y1*eta;
                  1.0/4.0*z1*zet*x2*y1*eta        1.0/4.0*z1*zet*x4*xi*y2        1.0/4.0*z2*x4*xi*y1*eta;
                  1.0/4.0*z1*zet*x2*y3*eta        1.0/4.0*z1*zet*x4*xi*y4        1.0/4.0*z2*x4*xi*y3*eta;
                  1.0/4.0*z1*zet*x1*y3*eta        1.0/4.0*z1*zet*x3*xi*y4        1.0/4.0*z2*x3*xi*y3*eta;
                  1.0/4.0*z3*zet*x1*y1*eta        1.0/4.0*z3*zet*x3*xi*y2        1.0/4.0*z4*x3*xi*y1*eta;
                  1.0/4.0*z3*zet*x2*y1*eta        1.0/4.0*z3*zet*x4*xi*y2        1.0/4.0*z4*x4*xi*y1*eta;
                  1.0/4.0*z3*zet*x2*y3*eta        1.0/4.0*z3*zet*x4*xi*y4        1.0/4.0*z4*x4*xi*y3*eta;
                  1.0/4.0*z3*zet*x1*y3*eta        1.0/4.0*z3*zet*x3*xi*y4        1.0/4.0*z4*x3*xi*y3*eta;
                       -1.0/2.0*z1*zet*xi*y1*eta   1.0/2.0*z1*zet*(-x4)*x3*y2   1.0/2.0*z2*(-x4)*x3*y1*eta;
             1.0/2.0*z1*zet*x2*(-y3)*y1             -1.0/2.0*z1*zet*x4*xi*eta   1.0/2.0*z2*x4*xi*(-y3)*y1;
                       -1.0/2.0*z1*zet*xi*y3*eta   1.0/2.0*z1*zet*(-x4)*x3*y4   1.0/2.0*z2*(-x4)*x3*y3*eta;
             1.0/2.0*z1*zet*x1*(-y3)*y1             -1.0/2.0*z1*zet*x3*xi*eta   1.0/2.0*z2*x3*xi*(-y3)*y1;
                       -1.0/2.0*z3*zet*xi*y1*eta   1.0/2.0*z3*zet*(-x4)*x3*y2   1.0/2.0*z4*(-x4)*x3*y1*eta;
             1.0/2.0*z3*zet*x2*(-y3)*y1             -1.0/2.0*z3*zet*x4*xi*eta   1.0/2.0*z4*x4*xi*(-y3)*y1;
                       -1.0/2.0*z3*zet*xi*y3*eta   1.0/2.0*z3*zet*(-x4)*x3*y4   1.0/2.0*z4*(-x4)*x3*y3*eta;
             1.0/2.0*z3*zet*x1*(-y3)*y1             -1.0/2.0*z3*zet*x3*xi*eta   1.0/2.0*z4*x3*xi*(-y3)*y1;
             1.0/2.0*(-z3)*z1*x1*y1*eta            1.0/2.0*(-z3)*z1*x3*xi*y2             -1.0/2.0*zet*x3*xi*y1*eta;
             1.0/2.0*(-z3)*z1*x2*y1*eta            1.0/2.0*(-z3)*z1*x4*xi*y2             -1.0/2.0*zet*x4*xi*y1*eta;
             1.0/2.0*(-z3)*z1*x2*y3*eta            1.0/2.0*(-z3)*z1*x4*xi*y4             -1.0/2.0*zet*x4*xi*y3*eta;
             1.0/2.0*(-z3)*z1*x1*y3*eta            1.0/2.0*(-z3)*z1*x3*xi*y4             -1.0/2.0*zet*x3*xi*y3*eta;
                      -z1*zet*xi*(-y3)*y1            -z1*zet*(-x4)*x3*eta  z2*(-x4)*x3*(-y3)*y1;
                      -(-z3)*z1*xi*y1*eta  (-z3)*z1*(-x4)*x3*y2            -zet*(-x4)*x3*y1*eta;
            (-z3)*z1*x2*(-y3)*y1            -(-z3)*z1*x4*xi*eta            -zet*x4*xi*(-y3)*y1;
                      -(-z3)*z1*xi*y3*eta  (-z3)*z1*(-x4)*x3*y4            -zet*(-x4)*x3*y3*eta;
            (-z3)*z1*x1*(-y3)*y1            -(-z3)*z1*x3*xi*eta            -zet*x3*xi*(-y3)*y1;
                      -z3*zet*xi*(-y3)*y1            -z3*zet*(-x4)*x3*eta  z4*(-x4)*x3*(-y3)*y1;
               -2.0*(-z3)*z1*xi*(-y3)*y1     -2.0*(-z3)*z1*(-x4)*x3*eta     -2.0*zet*(-x4)*x3*(-y3)*y1;
           ];
    return reshape(val, 27, 3);
end

function privboundaryconn(self::FESetH27)
    conn = connasarray(self)
    return  [conn[:, [1,  4,  3,  2,  12,  11,  10,  9,  21]];
             conn[:, [1,  2,  6,  5,  9,  18,  13,  17,  22]];
             conn[:, [2,  3,  7,  6,  10,  19,  14,  18,  23]];
             conn[:, [3,  4,  8,  7,   11,  20,  15,  19,  24]];
             conn[:, [4,  1,  5,  8,  12,  17,  16,  20,  25]];
             conn[:, [6,  7,  8,  5,  14,  15,  16,  13,  26]]];
end

function privboundaryfe(self::FESetH27)
    return FESetQ9;
end

function privinparametric(self::FESetH27, param_coords::FFltVec, tolerance::FFlt)
    s = findall(v -> (v >= (-1.0 - tolerance)) && (v <= (+1.0 + tolerance)), param_coords);
    return  (length(s) == length(param_coords));
end

function privcentroidparametric(self::FESetH27)
    return vec([0.0 0.0 0.0])
end

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

"""
    FESetT4

Type for sets of volume-like tetrahedral finite elements with four nodes.
"""
@define_FESet FESetT4 AbstractFESet3Manifold 4

function privbfun(self::FESetT4,  param_coords::FFltVec)
    # Evaluate the basis function matrix for an 3-node triangle.
     val = [(1 - param_coords[1] - param_coords[2] - param_coords[3]);
            param_coords[1];
            param_coords[2];
            param_coords[3]];
    return reshape(val, 4, 1)
end

function privbfundpar(self::FESetT4,  param_coords::FFltVec)
    # Evaluate the derivatives of the basis function matrix.
    val = [-1. -1. -1.;
           +1.  0.  0.;
           0. +1.  0.;
           0.  0. +1.];
    return reshape(val, 4, 3)
end

function privboundaryconn(self::FESetT4)
    conn = connasarray(self)
    return [conn[:, [1,  3,  2]]; conn[:, [1,  2,  4]]; conn[:, [2,  3,  4]]; conn[:, [1,  4,  3]]];
end

function privboundaryfe(self::FESetT4)
    return FESetT3;
end

function privinparametric(self::FESetT4, param_coords::FFltVec, tolerance::FFlt)
    s = findall(v -> (v >= (-tolerance)) && (v <= (+1.0 + tolerance)), param_coords);
    return  (length(s) == length(param_coords)) &&
      (-3.0 * tolerance <= sum(param_coords) <= 1.0 + 3.0 * tolerance);
end

function privcentroidparametric(self::FESetT4)
    return vec([1.0/4 1.0/4 1.0/4])
end

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

"""
    FESetT10

Type for sets of volume-like tetrahedral finite elements with 10 nodes.
"""
@define_FESet FESetT10 AbstractFESet3Manifold 10

function privbfun(self::FESetT10,  param_coords::FFltVec)
    # Evaluate the basis function matrix for an 3-node triangle.
    r = param_coords[1];
    s = param_coords[2];
    t = param_coords[3];
    val = [(1-r-s-t) * (2*(1-r-s-t)-1);
           r * (2*r-1);
           s * (2*s-1);
           t * (2*t-1);
           4*(1-r-s-t)*r;
           4*r*s;
           4*s*(1-r-s-t);
           4*(1-r-s-t)*t;
           4*r*t;
           4*s*t;
           ];
    return reshape(val, 10, 1)
end

function privbfundpar(self::FESetT10,  param_coords::FFltVec)
    # Evaluate the derivatives of the basis function matrix.
    r = param_coords[1];
    s = param_coords[2];
    t = param_coords[3];
            val = [-3+4*r+4*s+4*t   4*r-1 0   0  -8*r+4-4*s-4*t  4*s -4*s  -4*t  4*t  0;
                -3+4*r+4*s+4*t  0  4*s-1  0   -4*r   4*r 4-4*r-8*s-4*t -4*t  0 4*t;
                -3+4*r+4*s+4*t  0  0  4*t-1   -4*r   0  -4*s -8*t+4-4*r-4*s  4*r  4*s
                ]';
    return reshape(val, 10, 3)
end

function privboundaryconn(self::FESetT10)
    conn = connasarray(self)
    return [conn[:, [1,  3,  2,  7,  6,  5]];
            conn[:, [1,  2,  4,  5,  9,  8]];
            conn[:, [2,  3,  4,  6,  10,  9]];
            conn[:, [3,  1,  4,  7,  8,  10]]];
end

function privboundaryfe(self::FESetT10)
    return FESetT6;
end

function privinparametric(self::FESetT10, param_coords::FFltVec, tolerance::FFlt)
    s = findall(v -> (v >= (-tolerance)) && (v <= (+1.0 + tolerance)), param_coords);
    return  (length(s) == length(param_coords)) &&
      (-3.0 * tolerance <= sum(param_coords) <= 1.0 + 3.0 * tolerance);
end

function privcentroidparametric(self::FESetT10)
    return vec([1.0/4 1.0/4 1.0/4])
end

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

end


text:
#
# Function that checks if two atoms belong to the same residue
# without, of course, checking the residue counter
#
function same_residue(atom1::Atom, atom2::Atom)
    atom1.resnum != atom2.resnum && return false
    atom1.model != atom2.model && return false
    atom1.chain != atom2.chain && return false
    atom1.resname != atom2.resname && return false
    atom1.segname != atom2.segname && return false
    return true
end

ext: jl

max_stars_repo_name: mcubeg/PDBTools

max_issues_repo_name: mcubeg/PDBTools

max_forks_repo_name: mcubeg/PDBTools


text:
module WinTypes

# References:
# https://docs.microsoft.com/en-us/windows/win32/winprog/windows-data-types
# https://docs.microsoft.com/en-us/windows/win32/learnwin32/windows-coding-conventions?redirectedfrom=MSDN

const FALSE     = Cint(0)
const TRUE      = Cint(1)

const BYTE      = Cuchar
const WORD      = Cushort
const DWORD     = Culong
const DWORDLONG = UInt64
const DWORD32   = UInt32
const DWORD64   = UInt64

const BOOL      = Cint
const BOOLEAN   = BYTE

const CHAR      = Cchar
const CCHAR     = Cchar

const FLOAT     = Cfloat

const LONG      = Culong
const ULONG     = Culong
const LONG32    = Cint
const LONG64    = Int64
const LONG_PTR  = Int

const SHORT     = Cshort
const USHORT    = Cushort

const ATOM      = WORD
const LANGID    = WORD

const COLORREF  = DWORD
const LGRPID    = DWORD
const LCTYPE    = DWORD

const LCID      = DWORD

const INT       = Int32
const INT8      = Int8
const INT16     = Int16
const INT32     = Int32
const INT64     = Int64
const INT_PTR   = Int
const UINT_PTR  = UInt

const WPARAM    = UINT_PTR
const LPARAM    = LONG_PTR

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

const HANDLE = Ptr{Cvoid}

const HACCEL       = HANDLE
const HBITMAP      = HANDLE
const HBRUSH       = HANDLE
const HCOLORSPACE  = HANDLE
const HCONV        = HANDLE
const HCONVLIST    = HANDLE
const HDC          = HANDLE
const HDDEDATA     = HANDLE
const HDESK        = HANDLE
const HDROP        = HANDLE
const HDWP         = HANDLE
const HENHMETAFILE = HANDLE
const HFILE        = Cint
const HFONT        = HANDLE
const HGDIOBJ      = HANDLE
const HGLOBAL      = HANDLE
const HHOOK        = HANDLE
const HICON        = HANDLE
const HCURSOR      = HICON
const HINSTANCE    = HANDLE
const HKEY         = HANDLE
const HKL          = HANDLE
const HLOCAL       = HANDLE
const HMENU        = HANDLE
const HMETAFILE    = HANDLE
const HMODULE      = HANDLE
const HMONITOR     = HANDLE
const HPALETTE     = HANDLE
const HPEN         = HANDLE
const HRESULT      = Clong
const HRGN         = HANDLE
const HRSRC        = HANDLE
const HSZ          = HANDLE
const HWINSTA      = HANDLE
const HWND         = HANDLE

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

const PHANDLE = Ptr{HANDLE}

const UCHAR   = Cuchar
const WCHAR   = Cwchar_t
const PWCHAR  = Ptr{WCHAR}

const PWORD   = Ptr{WORD}
const LPWORD  = Ptr{WORD}
const PDWORD  = Ptr{DWORD}
const LPDWORD = Ptr{DWORD}

const VOID    = Cvoid
const PVOID   = Ptr{Cvoid}
const LPVOID  = Ptr{Cvoid}
const LPCVOID = Ptr{Cvoid}

const LPINT   = Ptr{Cint}
const PINT    = Ptr{Cint}

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

const PSTR    = Ptr{CHAR}
const PCSTR   = Ptr{CHAR}
const LPSTR   = Ptr{CHAR}
const LPCSTR  = Ptr{CHAR}

const PWSTR   = Ptr{WCHAR}
const PCWSTR  = Ptr{WCHAR}
const LPWSTR  = Ptr{WCHAR}
const LPCWSTR = Ptr{WCHAR}

end


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

using Markdown
using InteractiveUtils

# ╔═╡ 0ffc59b2-f8f7-11eb-1d41-d7be3a660bff
begin
	using Pkg
	Pkg.activate(joinpath(@__DIR__, ".."))
end

# ╔═╡ f601a1f4-0da2-457a-8276-97801af40c82
begin
	using LRHousing
	using CategoricalArrays
end

# ╔═╡ 2df81d5b-72e2-44b7-b84b-b1cb43bb05b5
using StatsBase

# ╔═╡ add411c5-4026-4c9c-92a3-e7181ccd6467
using DataFrames

# ╔═╡ 1bcda19e-3bac-48cb-a550-f23a76e72cf2
using Chain: @chain

# ╔═╡ 4049eb62-6ffa-427b-8a35-7f8b67e51f70
train = LRHousing.read_data("../data/train.csv")

# ╔═╡ 24af9da5-8615-4c36-b5df-e10409ca008a
pl = LRHousing.DataPipeline(train; ignore=["Id"])

# ╔═╡ b06fbcf9-6f2b-412e-aaea-3497fae925c3
trainp = select(pl(train), Not(pl.ignore))

# ╔═╡ f7aa44b0-a3c0-470a-a1d6-1274d80ca67e
describe(train)

# ╔═╡ 0fb1ac77-1d12-4ebf-91e1-983bc925d86d
function number_binning(X, k=10)
	minx = minimum(X)
	maxx = maximum(X)
	bins = range(start=minx, stop=maxx, length=k)
	map(x -> findfirst(x .< bins), X)
end

# ╔═╡ e8dd05ec-b8ce-4c84-b613-6fa38b4f0fb0
number_binning(trainp.GrLivArea)

# ╔═╡ 36dc88cc-dbf2-487b-b24a-494c4722ad32
begin
	function crammerv(x::CategoricalArray, y::CategoricalArray)
		@assert length(x) == length(y)
		lx = levelcode.(x)
		ly = levelcode.(y)
		cx = countmap(lx)
		cy = countmap(ly)
		ux = unique(lx)
		uy = unique(ly)
		k = length(unique(x))
		r = length(unique(y))
		cxy = countmap(zip(lx, ly))
		n = length(x)
		chisq = sum(begin
				nume = (get(cxy, (i, j), 0) - cx[i] * cy[j] / n)^2
				deno = (cx[i] * cy[j] / n)
				nume / deno
			end for (i, j) in Iterators.product(ux, uy))
		sqrt(chisq / n / min(k - 1, r - 1))
	end
	crammerv(x::CategoricalArray, y) = crammerv(x, number_binning(y, 10))
	crammerv(x, y::CategoricalArray) = crammerv(number_binning(x, 10), y)
end

# ╔═╡ 37c0a640-bda7-4e25-916c-ef182de83fc9
crammerv(categorical(train.MSZoning), train.LotShape |> categorical)

# ╔═╡ a0ba3080-b1ce-4c63-8480-28594c4a69de
crammerv(trainp.BldgType, trainp.BldgType)

# ╔═╡ f025061b-99ac-4663-a104-30bb9def46ba
function cor_matrix(df)
	result = DataFrame(col=names(df))
	function compute_cor(name_1, name_2)
		col_1 = df[!, name_1]
		col_2 = df[!, name_2]
		if eltype(col_1) <: Number && eltype(col_2) <: Number
			StatsBase.cor(col_1, col_2)
		else
			crammerv(categorical(col_1), categorical(col_2))
		end
	end
	transformations = map(names(df)) do name_1
		:col => ByRow(name_2 -> compute_cor(name_1, name_2)) => Symbol(name_1)
	end
	transformations
	transform(result, transformations...)
end

# ╔═╡ f7643ebf-da43-42fb-bd75-d204b54d0b08
eltype(trainp.MSSubClass) <: Number

# ╔═╡ 7c52649d-77a0-4fda-991f-89ea0ade4211
cordf = cor_matrix(trainp)

# ╔═╡ 4dbe4a01-bbc3-4274-9ca0-03e0a903d1e1
function find_multicolinearity(cordf, ycol; threshold=0.8)
	to_be_removed = String[]
	sizehint!(to_be_removed, size(cordf, 1)÷3)
	considers = setdiff(cordf.col, ["SalePrice"])
	name_index = String.(cordf.col)
	df = select(cordf, Not("col"))
	map(Iterators.product(considers, considers)) do (n1, n2)
		row_ind = findfirst(n1 .== name_index)
		col_ind = findfirst(n2 .== name_index)
		df[row_ind, col_ind]
	end
	for n1 in considers, n2 in considers
		if n1 in to_be_removed || n2 in to_be_removed || n1 == n2
			continue
		end
		row_ind = findfirst(n1 .== name_index)
		col_ind = findfirst(n2 .== name_index)
		if abs(df[row_ind, col_ind]) ≥ threshold
			append!(to_be_removed, [df[row_ind, ycol] > df[col_ind, ycol] ? n2 : n1])
		end
	end
	to_be_removed
end

# ╔═╡ 97d8c5d2-58ac-4bb1-a5e5-4eefcfe1edfd
find_multicolinearity(cordf, "SalePrice"; threshold=0.8)

# ╔═╡ Cell order:
# ╠═0ffc59b2-f8f7-11eb-1d41-d7be3a660bff
# ╠═f601a1f4-0da2-457a-8276-97801af40c82
# ╠═4049eb62-6ffa-427b-8a35-7f8b67e51f70
# ╠═24af9da5-8615-4c36-b5df-e10409ca008a
# ╠═b06fbcf9-6f2b-412e-aaea-3497fae925c3
# ╠═2df81d5b-72e2-44b7-b84b-b1cb43bb05b5
# ╠═f7aa44b0-a3c0-470a-a1d6-1274d80ca67e
# ╠═37c0a640-bda7-4e25-916c-ef182de83fc9
# ╠═e8dd05ec-b8ce-4c84-b613-6fa38b4f0fb0
# ╠═0fb1ac77-1d12-4ebf-91e1-983bc925d86d
# ╠═36dc88cc-dbf2-487b-b24a-494c4722ad32
# ╠═a0ba3080-b1ce-4c63-8480-28594c4a69de
# ╠═f025061b-99ac-4663-a104-30bb9def46ba
# ╠═add411c5-4026-4c9c-92a3-e7181ccd6467
# ╠═f7643ebf-da43-42fb-bd75-d204b54d0b08
# ╠═7c52649d-77a0-4fda-991f-89ea0ade4211
# ╠═97d8c5d2-58ac-4bb1-a5e5-4eefcfe1edfd
# ╠═1bcda19e-3bac-48cb-a550-f23a76e72cf2
# ╠═4dbe4a01-bbc3-4274-9ca0-03e0a903d1e1


text:
reduce1arg = [
(:sum,     :(ones(x))),
(:sumabs,  :(sign(x))),
(:sumabs2, :(2x)),
(:prod,    :(y./x)),
(:maximum, :(y.==x)),
(:minimum, :(y.==x)),
(:maxabs,  :(y.==abs(x))),
(:minabs,  :(y.==abs(x))),
]

for (f,g) in reduce1arg
    @eval @primitive $f(x,i...),dy,y   (dy.*($g))
    addtest(f, rand(2))
    addtest(f, rand(2,2), 1)
    addtest(f, rand(2,2), 2)
    # @eval @primitive $f(x::Tuple),dy,y (x=[x...];tuple((dy.*($g))...))
    # addtest(f, (rand(2)...))
end

# TODO: more general tuple reduction
@primitive  sum(x::Tuple),dy  ntuple(i->dy,length(x))
addtest(sum, (rand(2)...))

# TODO: implement more general sum ops
# TODO: other functions in reduce.jl:
# r_promote: Not exported
# mapfoldl_impl: Not exported
# mapfoldl
# foldl
# mapfoldr_impl: Not exported
# mapfoldr
# foldr
# mapreduce_seq_impl: Not exported
# mapreduce_pairwise_impl: Not exported
# mapreduce
# mapreduce_impl: Not exported
# mr_empty: Not exported
# _mapreduce: Not exported
# reduce
# shortcircuits: Not exported
# shorted: Not exported
# sc_finish: Not exported
# mapreduce_sc_impl: Not exported
# mapreduce_no_sc: Not exported
# mapreduce_sc: Not exported
# sum_pairwise_blocksize: Not exported
# sum
# sumabs
# sumabs2
# sum_kbn
# prod
# maximum
# minimum
# maxabs
# minabs
# extrema
# any
# all
# in
# ∉
# ∋
# ∌
# contains
# count
# call
# countnz


text:
include("dyna_agent.jl")
include("env_models/env_models.jl")


text:
# Intervals testing
I1 = Intervals([(0.0,1.0)])
I2 = Intervals([(0.0,1.0)])
@test I1 == I2
@test (I1 ∪ I2).list == [(0.0,1.0)]
@test (I1 ∩ I2).list == [(0.0,1.0)]
@test I1 ∪ I2 == I1 ∩ I2
@test I1 ∪ I2 == I2 ∪ I1
@test I1 ∩ I2 == I2 ∩ I1

I1 = Intervals([(0.0,1.0),(1.0,2.0)])
I2 = Intervals([(0.0,1.0),(1.0,2.0)])
@test I1 == I2
@test I1 ∪ I2 == I1 ∩ I2
@test (I1 ∩ I2).list == [(0.0,2.0)]
@test (I1 ∪ I2).list == [(0.0,2.0)]

I1 = Intervals([(0.0,1.0)])
I2 = Intervals([(0.0,1.0),(2.0,3.0)])
@test length(I1) == 1.0
@test count(I1) == 1.0
@test length(I2) == 2.0
@test count(I2) == 2.0
@test (I1 ∩ I2).list == [(0.0,1.0)]
@test length(I1 ∩ I2) == 1.0
@test count(I1 ∩ I2) == 1.0
@test (I1 ∪ I2).list == [(0.0,1.0),(2.0,3.0)]
@test length(I1 ∪ I2) == 2.0
@test I1 ⊆ I2
@test !(I2 ⊆ I1)
@test 0.5 ∈ I1
@test 1.0 ∈ I1
@test 2.5 ∈ (I1 ∪ I2)

I1 = Intervals([(0.0,1.0),(2.0,3.0)])
I2 = Intervals([(0.5,4.0)])
@test (I1 ∩ I2).list == [(0.5,1.0),(2.0,3.0)]
@test (I1 ∪ I2).list == [(0.0,4.0)]

I1 = Intervals([(0.0,1.0),(2.0,3.0)])
I2 = Intervals([(1.1,2.1),(2.9,4.0)])
@test (I1 ∩ I2).list == [(2.0,2.1),(2.9,3.0)]
@test (I1 ∪ I2).list == [(0.0,1.0),(1.1,4.0)]

I1 = Intervals([(0.0,1.0),(2.0,3.0)])
I2 = Intervals([(1.1,1.9)])
@test (I1 ∩ I2).list == []
@test (I1 ∪ I2).list == [(0.0,1.0),(1.1,1.9),(2.0,3.0)]

I1 = Intervals([(0.0,1.0),(2.0,3.0)])
I2 = Intervals([(1.0,2.0)])
@test (I1 ∩ I2).list ==[(1.0,1.0),(2.0,2.0)]
@test (I1 ∪ I2).list == [(0.0,3.0)]

I1 = Intervals()
I2 = Intervals([(0.0,1.0),(2.5,3.5)])
@test length(I1) == 0
@test count(I1) == 0
@test I1 ∩ I2 == Intervals()
@test I1 ∪ I2 == I2

I1 = Intervals([(0.2,0.6),(2.0,4.2)])
I2 = Intervals([(0.0,1.0),(2.5,3.5)])
I3 = Intervals([(0.1,3.0),(3.2,3.5)])
@test I1 ∩ I2 == Intervals([(0.2,0.6),(2.5,3.5)])
I4 = Intervals([(0.2,0.6),(2.5,3.5)])
@test I4 ∩ I3 == Intervals([(0.2,0.6),(2.5,3.0),(3.2,3.5)])
@test I1 ∩ I2 ∩ I3 == Intervals([(0.2,0.6),(2.5,3.0),(3.2,3.5)])
@test I1 ∪ I2 ∪ I3 == Intervals([(0.0,4.2)])

I1 = Intervals([(0.0,1.0),(1.0,2.0)])
I2 = Intervals([(0.3,0.6)])
@test I1 ∩ I2 == I2
@test I1 ∪ I2 == I1

I1 = Intervals([(0.0,1.0),(1.0,2.0)])
I2 = Intervals([(0.000001,0.99)])
@test I1 ∩ I2 == I2
@test I1 ∪ I2 == I1


text:
module NES

# Importing our NES Interface
include("NESInterface/NESInterface.jl")

include("SMBEnv/SMBEnv.jl")
export SMBEnv, reset!, step!, render
end # end module


text:
# An archive saves information about interesting points during an optimization
# run. It also saves (at least) the two last best fitness values since they
# can be used to estimate a confidence interval for how much improvement
# potential there is.
abstract Archive

# A top list archive saves a top list of the best performing (best fitness)
# candidates/individuals seen.
type TopListArchive <: Archive
  start_time::Float64   # Time when archive created, we use this to approximate the starting time for the opt...
  num_fitnesses::Int64  # Number of calls to add_candidate

  size::Int64           # Max size of top lists
  count::Int64          # Current number of values in the top lists

  fitnesses::Array{Float64,1}  # Top fitness values
  candidates::Array{Any, 1}    # Top candidates corresponding to top fitness values

  # We save a fitness history that we can later print to a csv file.
  # For each magnitude class (as defined by magnitude_class function below) we
  # we save the first entry of that class. The tuple saved for each magnitude
  # class is: (magnitude_class, time, num_fevals, fitness, width_of_confidence_interval)
  fitness_history::Array{(Float64, Int64, Float64, Float64),1}

  numdims::Int64        # Number of dimensions in opt problem. Needed for confidence interval estimation.

  TopListArchive(numdims, size::Int64 = 10) = begin
    new(time(), 0, size, 0, Float64[], Any[],
      (Float64, Int64, Float64, Float64)[], int(numdims))
  end
end

best_candidate(a::Archive) = a.candidates[1]
best_fitness(a::Archive) = a.fitnesses[1]
last_top_fitness(a::Archive) = a.fitnesses[a.count]
should_enter_toplist(fitness::Float64, a::Archive) = fitness < last_top_fitness(a)

# Delta fitness is the difference between the top two candidates found so far.
#
function delta_fitness(a::Archive)
  if length(a.fitness_history) < 2
    Inf
  else
    abs(a.fitness_history[2][3] - a.fitness_history[1][3])
  end
end

# Add a candidate with a fitness to the archive (if it is good enough).
function add_candidate!(a::TopListArchive, fitness::Number,
  candidate, num_fevals::Int64 = -1)
  a.num_fitnesses += 1

  if a.count < a.size

    if a.count == 0 || fitness < best_fitness(a)
      push_to_fitness_history!(a, fitness, num_fevals)
    end

    push_then_sort_by_fitness!(fitness, candidate, a)
    a.count += 1

  elseif should_enter_toplist(fitness, a)

    if fitness < best_fitness(a)
      push_to_fitness_history!(a, fitness, num_fevals)
    end

    push_then_sort_by_fitness!(fitness, candidate, a)

    # Cut so only top list values are saved
    a.fitnesses = a.fitnesses[1:a.size]
    a.candidates = a.candidates[1:a.size]

  end
end

# The magnitude class of a value is a tuple indicating its sign and scale in a
# tuple. This is used for filtering so that we only need to save one history
# value per magnitude class instead of saving them all.
function magnitude_class(f::Number)
  f = float(f)
  if f == 0.0
    (-1.0, 1e100)
  else
    (sign(f), floor(10.0*log10(abs(f)))/10.0)
  end
end

# Save fitness history so we can reconstruct the most important events later.
function push_to_fitness_history!(a::Archive, fitness::Number, num_fevals::Int64 = -1)
  push!(a.fitness_history, make_fitness_event(a, fitness, num_fevals))
end

function fitness_improvement_ratio(a::Archive, newFitness)
  try
    bestfitness = best_fitness(a)
    return abs( (bestfitness - newFitness) / bestfitness )
  catch
    return 0.0
  end
end

# Calculate the distance from a fitness value to the optimum/best known fitness value.
function distance_to_optimum(fitness, bestfitness)
  abs(fitness - bestfitness)
end

function make_fitness_event(a::Archive, fitness::Number, num_fevals::Int64 = -1)
  nf = num_fevals == -1 ? a.num_fitnesses : num_fevals
  (time(), nf, fitness, fitness_improvement_ratio(a, fitness))
end

function fitness_history_csv_header(a::Archive)
  "Date,Time,ElapsedTime,Magnitude,NumFuncEvals,FitnessImprovementRatio,Fitness"
end

function save_fitness_history_to_csv_file(a::Archive, filename = "fitness_history.csv";
  header_prefix = "", line_prefix = "",
  include_header = true, bestfitness = nothing)

  fh = open(filename, "a+")

  if include_header

    header = [header_prefix, fitness_history_csv_header(a)]

    if bestfitness != nothing
      push!(header, "DistanceToBestFitness")
    end

    println(fh, join(header, ","))

  end

  for(event in a.fitness_history)

    t, nf, f, fir = event

    mc = magnitude_class(f)

    line = [line_prefix, strftime("%Y-%m-%d,%T", t), t-a.start_time,
      mc[1]*mc[2], nf, fir, f]

    if bestfitness != nothing
      push!(line, distance_to_optimum(f, bestfitness))
    end

    println(fh, join(line, ","))

  end

  close(fh)

end

function push_then_sort_by_fitness!(fitness::Number, candidate, a::Archive)
  push!(a.fitnesses, fitness)
  push!(a.candidates, candidate)
  order = sortperm(a.fitnesses)
  a.fitnesses = a.fitnesses[order]
  a.candidates = a.candidates[order]
end

# Merge multiple fitness histories and calculate the min, max, avg and median
# values for fitness and fir at all fitness eval change points.
function merge_fitness_histories(histories)
  numhist = length(histories)
  counts = ones(numhist)
  fitnesses = zeros(numhist)
  firs = zeros(numhist)
  current_feval = 1

  # Find max history length
  histlengths = [length(h) for h in histories]
  maxlen = maximum(histlengths)

  while maximum(counts) < maxlen

    # Find min feval for current events, this is the next current feval
    for i in 1:numhist
      t, nf, f, fir = histories[i][counts[i]]
    end

  end

end

# Merge the fitness histories and save the average values of the fitness,
# and distance to best fitness for each change in any of the histories.
#function merge_fitness_histories_to_csv(archives::Archive[], filename = "fitness_history.csv";
#  header_prefix = "", line_prefix = "",
#  include_header = true, bestfitness = nothing)
#
#end

# Calculate the width of the confidence interval at a certain p-value.
# This is based on the paper:
#  Carvalho (2011), "Confidence intervals for the minimum of a
#    function using extreme value statistics"
#
# This means that the current estimate of the confidence interval for the minimum
# of the optimized function lies within the interval
#
#     ] l1 - w, l1 [
#
# with probability (1-p) as the number of sampled function points goes to infinity,
# where
#   w = width_of_confidence_interval(a, p)
#   l1 = best_fitness(a)
#
function width_of_confidence_interval(a::Archive, p = 0.01)
  if length(a.fitness_history) < 2
    return Inf
  else
    l1 = a.fitnesses[1]
    l2 = a.fitnesses[2]
    # We use abs below so it works also for maximization.
    abs(l2 - l1) / ( (1-p)^(-2/a.numdims) - 1)
  end
end

# We define the improvement potential as the percentage improvement
# that can be expected from the current fitness value at a given p value.
# In theory, an optimization run should be terminated when this value is very
# small, i.e. there is little improvement potential left in the run.
function fitness_improvement_potential(a::Archive, p = 0.01)
  width_of_confidence_interval(a, p) / abs(best_fitness(a))
end