text:
################################################################################
#
#   METHODS FOR LUTTINGER TISZA CALCULATION
#
#   STRUCTURE OF THE FILE
#
#   1) TYPE LTBANDSTRUCTURE
#       - type definition
#       - printInfo function
#
#   2) CALCULATION OF LT BAND STRUCTURES OF UNTICELL OBJECTS
#      (- LT constraint and deviation functions (NOT EXPORTED) )
#       - spin interaction matrices
#       - calculation of band structures
#
#   3) PLOTTING OF LT BAND STRUCTURES
#       - plotting of Bandstructure objects
#       - plotting of bandstructures of unitcells along paths
#
#   4) CALCULATION OF LT GROUND STATES (k space manifold)
#
#   5) PLOTTING OF LT GROUND STATES (k space manifold)
#       - plotting from points
#       - plotting from unitcell
#
#   BUG: So far, the constraint in all calculation ONLY works for k=0
#        For all other k values, the constraint is checked also for k=0
#
#   NOTE: Constraint has to be check with exp(ikr) factors on other sites as well
#
#   NOTE: Factor of 1/2 is missing for correct minimal energy
#
################################################################################

################################################################################
#
#   TYPE LTBANDSTRUCTURE
#       - type definition
#       - printInfo function
#
################################################################################
"""
    struct LTBandstructure

The type that contains information on a Luttinger Tisza band structure (energy values in momentum space with spin constraints).
Fields are

    path              :: Path
    bands             :: Array{Array{Array{Float64, 1}, 1}, 1}
    constraint_values :: Array{Array{Array{Float64, 1}, 1}, 1}

Note that the notation of `bands` (and `constraint_values`) is the following:
- `bands[i]`       gives all bands of path segment `i`
- `bands[i][j]`    gives all energy values for band `j` of path segment `i`
- `bands[i][j][k]` gives the energy value at kpoint index `k` of band `j` in path segment `i`

The `constraint_values` carry the minimal deviation from unit spin length within the unitcell, i.e.
value is the minimum of all sum(|s_i - 1.0|^2) for constructed s_i.

New `LTBandstructure` objects can be created only by the default constructor or one of
the several functions to calculate LT band structures.

# Examples

```julia-repl
julia> bandstructure = LTBandstructure(path, bands, constraint_values)
LatticePhysics.LTBandstructure(...)
```
"""
struct LTBandstructure

    # the path along which the band structure is calcualted
    path::Path

    # bands for each segment
    # bands[i] gives all bands of segment i
    # bands[i][j] gives all energy values for band j of segment i
    # bands[i][j][k] gives the energy value at kpoint index k of band j in segment i
    bands::Array{Array{Array{Float64, 1}, 1}, 1}

    # constraint value for all bands
    # value is the minimum of all sum(|s_i - 1.0|^2) for constructed s_i
    constraint_values::Array{Array{Array{Float64, 1}, 1}, 1}

    # ONLY DEFAULT CONSTRUCTOR

end

# export the type
export LTBandstructure

# INFORMATION FUNCTION
"""
    printInfo(bandstructure::LTBandstructure [; constraint::Float64=1e-6])

Prints information about a `LTBandstructure` in terms of
how many eigenvalues of the contained bands satisfy the LT constraint in which region of the path.
The `constraint` parameter allows to sharpen or soften the distinction between constraint
fullfilling and constraint breaking eigenvalues.

# Examples

```julia-repl
julia> printInfo(bandstructure)
...

julia> printInfo(unitcell, constraint=1e-8)
...

```
"""
function printInfo(bandstructure::LTBandstructure; constraint::Float64=1e-6)
    # print the header
    println("Bandstructure (LT), constraint satisfied for var(length) < $(constraint)")
    print("\t$(bandstructure.path.point_names[1])")
    for p in bandstructure.path.point_names[2:end]
        print("\t->-\t$(p)")
    end
    println("")
    # print each band
    for b in 1:length(bandstructure.bands[1])
        # first the point
        print("$(b))\t|")
        # then all segments with closing points
        for s in 1:length(bandstructure.bands)
            # segment s starts
            print("\t")
            # calculate the number of constraint fullfilling values
            print("$(round(100.0*sum([c>constraint ? 0 : 1 for c in bandstructure.constraint_values[s][b]])/length(bandstructure.constraint_values[s][b]),2))%")
            # print the next point
            print("\t|")
        end
        # print a new line
        println("")
    end
end

# INTERNAL FUNCTIONS CONCERNING THE CONSTRAINT IN THE LT CALCULATION
# (not exported)

# function to calculate the deviation from best result
function deviation(spin_eigenvectors::Array{Array{Complex{Float64},1},1}, spin_dimension::Int64, alpha::Array{Float64,1})
    # build up the global spin vector
    spin_vector = zeros(length(spin_eigenvectors[1]))
    for s in 1:length(spin_eigenvectors)
        spin_vector = spin_vector .+ alpha[s].*spin_eigenvectors[s]
    end
    # find out the individual lengths of spins
    spin_lengths = spin_vector .* conj.(spin_vector)
    spin_lengths = [sum(spin_lengths[s:s+spin_dimension-1]) for s in 1:spin_dimension:length(spin_lengths)]
    # find out the deviation from unity
    dl = sum(abs.((spin_lengths .- 1)))
    return dl
end

# Definition of constraint
function getLTConstraint(spin_eigenvectors::Array{Array{Complex{Float64},1},1}, spin_dimension::Int64)
    # detrmine what to do based on the number of eigenvectors
    if length(spin_eigenvectors) == 1
        # find out the individual lengths of spins
        spin_lengths = spin_eigenvectors[1] .* conj.(spin_eigenvectors[1])
        spin_lengths = [sum(spin_lengths[s:s+spin_dimension-1]) for s in 1:spin_dimension:length(spin_lengths)]
        # find out the deviation from unity
        dl = sum(abs.((spin_lengths .- mean(spin_lengths))))
    else
        # optimize the function
        return Optim.minimum(Optim.optimize(x -> deviation(spin_eigenvectors, spin_dimension, x), ones(length(spin_eigenvectors))))
    end
end

################################################################################
#
#   INTERACTION MATRICES FOR SPIN SYSTEMS
#   - bond interaction matrices
#   - global interaction matrix (for the unitcell)
#
################################################################################

# Function to create a bond strength matrix
"""
    getBondInteractionMatrixHeisenbergKitaev(connection::Array{Any,1})

Constructs the spin interaction matrix for a single bond given by `connection`.
The matrix is of size 3x3 and is of type `Array{Float64,2}`, i.e. it has real entries.
Depending on the parameter of the bond strength (i.e. `c[3]`), it chooses to return a different matrix.
The replaced strings contain
- `"J1"`,`"J2"` --> Heisenberg type, strength 1.0
- `"Jx"`,`"Jy"`,`"Jz"`,`"tx"`,`"tx"`,`"tx"` --> Kitaev type, strength 1.0

Otherwise, the matrix has entries `matrix[i,i] = c[3]`

# Examples

```julia-repl
julia> getBondInteractionMatrixHeisenbergKitaev(Any[1,2, "tx", (0,0,0)])
3×3 Array{Float64,2}:
 1.0  0.0  0.0
 0.0  0.0  0.0
 0.0  0.0  0.0

julia> getBondInteractionMatrixHeisenbergKitaev(Any[1,2, "J1", (0,0,0)])
3×3 Array{Float64,2}:
 1.0  0.0  0.0
 0.0  1.0  0.0
 0.0  0.0  1.0

```
"""
function getBondInteractionMatrixHeisenbergKitaev(connection::Array{Any,1})
    # new 3x3 matrix
    bond_matrix = zeros(3,3)
    # get the bond strength
    strength = connection[3]
    # check what type the bond is
    if typeof(strength) == String
        if strength == "J1"
            bond_matrix[1,1] = 1.0
            bond_matrix[2,2] = 1.0
            bond_matrix[3,3] = 1.0
        elseif strength == "J2"
            bond_matrix[1,1] = 1.0
            bond_matrix[2,2] = 1.0
            bond_matrix[3,3] = 1.0
        elseif strength == "Jx" || strength == "tx"
            bond_matrix[1,1] = 1.0
        elseif strength == "Jy" || strength == "ty"
            bond_matrix[2,2] = 1.0
        elseif strength == "Jz" || strength == "tz"
            bond_matrix[3,3] = 1.0
        end
    else
        bond_matrix[1,1] = strength
        bond_matrix[2,2] = strength
        bond_matrix[3,3] = strength
    end
    # return the matrix
    return bond_matrix
end
export getBondInteractionMatrixHeisenbergKitaev

"""
    getBondInteractionMatrixHeisenberg(connection::Array{Any,1})

Constructs the spin interaction matrix for a single bond given by `connection`.
The matrix is of size 1x1 and is of type `Array{Float64,2}`, i.e. it has only one (real) entry.
Depending on the parameter of the bond strength (i.e. `c[3]`), it chooses to return a different matrix.
The replaced strings contain
- `"J1"`,`"J2"` --> Heisenberg type, strength 1.0

Otherwise, the matrix has entry `matrix[i,i] = c[3]`

# Examples

```julia-repl

julia> getBondInteractionMatrixHeisenberg(Any[1,2, "J1", (0,0,0)])
1×1 Array{Float64,2}:
 1.0

julia> getBondInteractionMatrixHeisenberg(Any[1,2, 3.0, (0,0,0)])
1×1 Array{Float64,2}:
 3.0
```
"""
function getBondInteractionMatrixHeisenberg(connection::Array{Any,1})
    # new 3x3 matrix
    bond_matrix = zeros(1,1)
    # get the bond strength
    strength = connection[3]
    # check what type the bond is
    if typeof(strength) == String
        if strength == "J1"
            bond_matrix[1,1] = 1.0
        elseif strength == "J2"
            bond_matrix[1,1] = 1.0
        end
    else
        bond_matrix[1,1] = strength
    end
    # return the matrix
    return bond_matrix
end
export getBondInteractionMatrixHeisenberg

# Function to produce interaction matrices for entire unitcell
"""
    getSpinInteractionMatrixKSpace(
        unitcell::Unitcell,
        k_vector::Array{Float64,1},
     [  bondInteractionMatrix::Function
      ; enforce_hermitian::Bool=false ]
    )

Constructs the spin interaction matrix (in *momentum* space at point `k_vector`)
for use in a Luttinger Tisza calculation of a given `Unitcell` object.
The matrix is a sNxsN matrix where N is the number of sites in the given object and s is the spin dimension.
Entries (i+s,j+s') contain the interactions between sites i and j and spin components s and s' as well as the phase factor exp(i k*delta).

The matrix is of type `Array{Complex,2}`, i.e. it has complex entries to account for the phase factor.

The precise form of the spin interaction along a bond is given by the function `bondInteractionMatrix`
(which has the default of giving 1x1 matrices with strength = c[3] for connections c).
This function can be customised along the format f(connection) = matrix.
The default of this function is `getBondInteractionMatrixHeisenberg`.

Note that it is a possibility to add custom parameters by distinguishing different bond types based on
their `String` valued strength but then returning a matrix of `Float64` entries.

# Examples

```julia-repl
julia> getSpinInteractionMatrixKSpace(unitcell, k)
2×2 Array{Complex,2}:
...

julia> getSpinInteractionMatrixKSpace(unitcell, [pi/2.0, 0.0])
2×2 Array{Complex,2}:
...

julia> getSpinInteractionMatrixKSpace(unitcell, [pi/2.0, 0.0], c->diagm([ c[3], c[3] ]))
4×4 Array{Complex,2}:
...
```
"""
function getSpinInteractionMatrixKSpace(unitcell::Unitcell, k_vector::Array{Float64,1}, bondInteractionMatrix::Function)
    # get the spin dimension
    spin_dimension = size(bondInteractionMatrix(unitcell.connections[1]), 1)
    # create a new matrix
    matrix = zeros(Complex, spin_dimension*length(unitcell.basis), spin_dimension*length(unitcell.basis))
    # iterate over all connections
    for c in unitcell.connections
        # get the indices
        index_from  = Int(c[1])
        index_to    = Int(c[2])
        strength    = c[3]
        wrap        = c[4]
        # get the difference vector
        pos_delta   = unitcell.basis[index_to] .- unitcell.basis[index_from]
        if size(unitcell.lattice_vectors,1) > 0
            for pair in zip(wrap, unitcell.lattice_vectors)
                pos_delta .+= pair[1].*pair[2]
            end
        end
        # get the bond interaction matrix
        bond_interaction_matrix = bondInteractionMatrix(c)
        # add to the complete matrix twice
        for i in 1:spin_dimension
        for j in 1:spin_dimension
            matrix[(index_from-1)*spin_dimension + i, (index_to-1)*spin_dimension + j] += 0.5 * bond_interaction_matrix[i,j] * exp(-im * sum(pos_delta .* k_vector))
            matrix[(index_to-1)*spin_dimension + j, (index_from-1)*spin_dimension + i] += 0.5 * bond_interaction_matrix[i,j] * exp( im * sum(pos_delta .* k_vector))
        end
        end
    end
    # return the matrix
    return matrix
end
export getSpinInteractionMatrixKSpace

################################################################################
#
#   BAND STRUCTURE CALCULATION
#
################################################################################
"""
    getLTBandstructure(
                unitcell::Unitcell,
                path::Path,
                bondInteractionMatrix::Function
             [; resolution::Int64=-1,
                enforce_hermitian::Bool=false ]
            )

Calculates the Luttinger Tisza band struture of a `Unitcell` object
along some path given by a `Path` object and returns a `LTBandstructure` object.

Note 1: Optionally a bond interaction matrix can be given as a function
that constructs a matrix for a passed bond argument. If this function is not passed,
the default function will be used.

Note 2: The bond interaction matrix specifies the dimension of interacting spins.

# Examples

```julia-repl
julia> bandstructure = getLTBandstructure(unitcell, path)
LatticePhysics.LTBandstructure(...)

julia> bandstructure = getLTBandstructure(unitcell, path, resolution=1000)
LatticePhysics.LTBandstructure(...)
```
"""
function getLTBandstructure(
                unitcell::Unitcell,
                path::Path,
                bondInteractionMatrix::Function = getBondInteractionMatrixHeisenberg;
                resolution::Int64=-1,
                enforce_hermitian::Bool=false,
                epsilon_degenerate::Float64=1e-6
            )

    # maybe modify the path resolution
    if resolution > 0
        setTotalResolution!(path, resolution)
    end

    # build up the list of all bands of all segments (and all constraints)
    segments_total    = Array{Array{Float64,1},1}[]
    constraints_total = Array{Array{Float64,1},1}[]

    # get the spin dimension
    spin_dimension = size(bondInteractionMatrix(unitcell.connections[1]), 1)
    # iterate over all path segments and push empty lists into the segments list
    for i in 1:length(path.segment_resolution)
        # build an empty band structure for this segment
        segment     = Array{Float64, 1}[]
        constraints = Array{Float64, 1}[]
        for b in 1:length(unitcell.basis)*spin_dimension
            push!(segment,     zeros(Float64, path.segment_resolution[i]))
            push!(constraints,  ones(Float64, path.segment_resolution[i]))
        end
        # push the segment band structure into the complete segment list
        push!(segments_total,    segment)
        push!(constraints_total, constraints)
    end

    # iterate over all segments
    for s in 1:length(path.segment_resolution)
        # get the grid in between two points
        segment_resolution = path.segment_resolution[s]
        # get all multipliers of k vectors (i.e. all alpha in (1-alpha)*k_1 + alpha*k_2)
        multipliers = linspace(0, 1, segment_resolution)
        # get the local start and end point of the segment
        k1 = path.points[s]
        k2 = path.points[s+1]
        # calculate all energies
        for i in 1:segment_resolution
            # get the current k
            k = (k2 .* multipliers[i]) .+ (k1 .* (1-multipliers[i]))
            # get the interaction matrix for this k
            matrix = getSpinInteractionMatrixKSpace(unitcell, k, bondInteractionMatrix)
            # diagonalize the matrix
            eigenfactorization = eigfact(matrix)
            eigenvalues  = eigenfactorization[:values]
            eigenvectors = eigenfactorization[:vectors]
            # save all the eigenvalues to their lists
            for b in 1:length(eigenvalues)
                segments_total[s][b][i] = eigenvalues[b]
            end
            # compute the constraint, first find out what bands are degenerate
            # list of bands they are degenerate with
            degenerate = zeros(Int64,   length(eigenvalues)) .- 1
            treat      =  ones(Int64,   length(eigenvalues))
            for b in 2:length(eigenvalues)
                # check if degenerate
                if eigenvalues[b] - epsilon_degenerate <= eigenvalues[b-1]
                    # band b is degenerate with band b-1
                    degenerate[b-1] = b
                    # dont treat the band as it is treated before
                    treat[b] = 0
                end
            end
            # iterate over all bands
            for b in 1:length(eigenvalues)
                # if not treated, continue
                if treat[b] == 0
                    continue
                end
                # if treated, compile list of all bands
                degenerate_bands = Int64[b]
                while degenerate[degenerate_bands[end]] != -1
                    push!(degenerate_bands, degenerate[degenerate_bands[end]])
                end
                # get the LT constraint
                LT_constraint = getLTConstraint([eigenvectors[:,j] for j in degenerate_bands], spin_dimension)
                # save the constraint
                for d in degenerate_bands
                    constraints_total[s][d][i] = LT_constraint
                end
            end
        end
    end

    # generate a new LT band structure object
    bandstructure = LTBandstructure(path, segments_total, constraints_total)

    # return the LT band structure
    return bandstructure
end
export getLTBandstructure

################################################################################
#
#   BAND STRUCTURE PLOTTING
#
################################################################################

"""
    plotLTBandstructure(
            bandstructure::LTBandstructure
         [; constraint::Float64=1e-6
            limits_energy="AUTO",
            plot_title::String="",
            plot_color_valid="b",
            plot_color_invalid="r",
            figsize::Tuple=(6,4),
            showPlot::Bool=true,
            save_filename::String="NONE" ]
        )

    plotLTBandstructure(
            unitcell::Unitcell,
            path::Path,
         [  bondInteractionMatrix::Function
          ; resolution::Int64=-1,
            enforce_hermitian::Bool=false,
            ... ]
        )

Plots the Luttinger Tisza band struture of a passed `LTBandstructure` object along some its path
and returns the plot as a `PyPlot.Figure` object.
Alternatively, one can pass a `Unitcell` and `Path` (and optionally a bond matrix function)
to calculate the Luttinger Tisza band structure which is plotted.

Additional options include setting the LT constraint,
plotting related options of `PyPlot` as well as determining if the plot is saved or shown.

# Examples

```julia-repl
julia> plotLTBandstructure(unitcell, path)
PyPlot.Figure(...)

julia> plotLTBandstructure(unitcell, path, c->diagm([ c[3] ]))
PyPlot.Figure(...)

julia> plotLTBandstructure(unitcell, path, showPlot=false)
PyPlot.Figure(...)

julia> plotLTBandstructure(unitcell, save_filename="myplot.pdf")
PyPlot.Figure(...)

julia> plotLTBandstructure(bandstructure)
PyPlot.Figure(...)
```
"""
function plotLTBandstructure(
            bandstructure::LTBandstructure;
            constraint::Float64=1e-6,
            limits_energy="AUTO",
            plot_title::String="",
            plot_color_valid="b",
            plot_color_invalid="r",
            figsize::Tuple=(6,4),
            showPlot::Bool=true,
            save_filename::String="NONE"
        )

    ###########################
    #   INITIAL SETTINGS
    ###########################

    # get the path from the bandstructure
    path = bandstructure.path

    # configure plot environment
    rc("font", family="serif")

    # create a new figure
    fig = figure(figsize=figsize)

    ###########################
    #   PLOT BANDS
    ###########################

    # plot the band structure
    for s in 1:length(bandstructure.bands)
        # plot the segment (only invalid stuff)
        for b in 1:length(bandstructure.bands[s])
            # xvalues
            xvals = collect(1:path.segment_resolution[s]) .+ sum(path.segment_resolution[1:s-1])
            yvals = bandstructure.bands[s][b]
            # check which values satisfy the constraint
            invalid_indices = collect(1:path.segment_resolution[s])[bandstructure.constraint_values[s][b] .>= constraint]
            # if no invalid indices are found, just skip
            if length(invalid_indices) == 0
                continue
            end
            # plot everything
            plot(
                [xvals[i] for i in invalid_indices],
                [yvals[i] for i in invalid_indices],
                ".$(plot_color_invalid)"
            )
        end
        # plot the segment (only valid stuff)
        for b in 1:length(bandstructure.bands[s])
            # xvalues
            xvals = collect(1:path.segment_resolution[s]) .+ sum(path.segment_resolution[1:s-1])
            yvals = bandstructure.bands[s][b]
            # check which values satisfy the constraint
            valid_indices   = collect(1:path.segment_resolution[s])[bandstructure.constraint_values[s][b] .< constraint]
            # if no valid indices are found, just skip
            if length(valid_indices) == 0
                continue
            end
            # plot everything
            plot(
                [xvals[i] for i in valid_indices],
                [yvals[i] for i in valid_indices],
                ".$(plot_color_valid)"
            )
        end
    end

    ###########################
    #   SET ALL TICKS (POINTS)
    ###########################

    # get the current axis
    ax = gca()
    axx = ax[:get_xaxis]()
    # compile tick positions and labels
    point_pos = Int64[]
    push!(point_pos, 1)
    for l in 1:length(path.segment_resolution)
        push!(point_pos, sum(path.segment_resolution[1:l]))
    end
    point_labels = String[path.point_names[i] for i in 1:length(path.points)]
    # configure tick labels
    xticks(point_pos, point_labels)
    # configure ticks
    axx[:set_tick_params](which="both", direction="out")
    axx[:set_tick_params](which="top", color="none")
    axy = ax[:get_yaxis]()
    axy[:set_tick_params](which="both", direction="out")

    # plot vertical lines for each point
    for p in point_pos
        axvline(p,color=[0.6, 0.6, 0.6], linestyle="--")
    end

    ###########################
    #   CONFIGURE AXIS & TITLE
    ###########################

    # label the axis
    xlabel("momentum")
    ylabel("energy")

    # energy limits
    # check if specific boundaries are desired
    if !(limits_energy == "AUTO")
        ylim(limits_energy[1], limits_energy[2])
    end

    # momentum limits (x axis)
    xlim(0, maximum(point_pos)+1)

    # set the title
    if plot_title == "AUTO"
        # set the title to an automatically generated title
        title("Luttinger Tisza spectrum along $(getPathString(path)), constraint $(constraint)")
    elseif plot_title == ""
        # do nothing title related
    else
        # set the title to the given title
        title(plot_title)
    end

    ###########################
    #   FINISH THE PLOT
    ###########################

    # tighten the layout
    tight_layout()

    # save the plot
    if save_filename != "NONE"
        # make sure the directory exists
        if contains(save_filename, "/")
    		# get the containing folder
    		folder = save_filename[1:findlast(save_filename, '/')]
    		# build the path to that folder
    		mkpath(folder)
    	end
        # save the plot
        savefig(save_filename)
    end

    # maybe show the plot
    if showPlot
        show()
    end

    # return the figure object
    return fig
end
function plotLTBandstructure(
            unitcell::Unitcell,
            path::Path,
            bondInteractionMatrix::Function = getBondInteractionMatrixHeisenberg;
            constraint::Float64=1e-6,
            resolution::Int64=-1,
            enforce_hermitian::Bool=false,
            limits_energy="AUTO",
            plot_title::String="",
            plot_color_valid="b",
            plot_color_invalid="r",
            figsize::Tuple=(6,4),
            showPlot::Bool=true,
            save_filename::String="NONE"
        )
    # calculate the bandstructure
    bandstructure = getLTBandstructure(unitcell, path, bondInteractionMatrix, resolution=resolution, enforce_hermitian=enforce_hermitian)
    # call the respective function
    return plotLTBandstructure(
                bandstructure;
                constraint=constraint,
                limits_energy=limits_energy,
                plot_title=plot_title,
                plot_color_valid=plot_color_valid,
                plot_color_invalid=plot_color_invalid,
                figsize=figsize,
                showPlot=showPlot,
                save_filename=save_filename
            )
end
export plotLTBandstructure

################################################################################
#
#   LT GROUNDSTATE MANIFOLDS (K SPACE)
#
#   Algorithm:
#   1) (maybe) find ground state energy by minimizing eigenvalues
#   2) use Newton method to minimize random points
#   3) refold in k space to 1.BZ
#   4) calculate the constraint
#
################################################################################

# Function for calculating LT Groundstate manifold in 2D (not exported)
function getLTGroundstateKSpace2D(
            unitcell::Unitcell,
            N_points::Int64,
            bondInteractionMatrix::Function=getBondInteractionMatrixHeisenberg;
            groundstate_energy::Float64=Inf,
            epsilon::Float64=1e-10,
            epsilon_k::Float64=1e-10,
            epsilon_degenerate::Float64=1e-6,
            bounds_lower::Array{Float64,1}=-2*pi.*ones(2),
            bounds_upper::Array{Float64,1}=2*pi.*ones(2),
            refold_to_first_BZ::Bool=true
        )

    # list of k values that contribute to the BZ
    k_values          = zeros(Float64, N_points, 2)
    constraint_values =  ones(Float64, N_points)
    # get the spin dimension
    spin_dimension = size(bondInteractionMatrix(unitcell.connections[1]), 1)

    # check if the GS energy has to be calculated
    if groundstate_energy == Inf
        # search for 10 random points and compare the found minimum
        for tries in 1:10
            # find a suitable starting point for the Newton algorithm
            k_start = Float64[rand(), rand()]
            for j in 1:length(k_start)
                k_start[j] = k_start[j] * bounds_lower[j] + (1-k_start[j]) * bounds_upper[j]
            end
            # optimize the energy
            groundstate_energy = min(Optim.minimum(Optim.optimize(k->minimum(eigvals(getSpinInteractionMatrixKSpace(unitcell, k, bondInteractionMatrix))), k_start)), groundstate_energy)
        end
        # print the groundstate_energy
        println("Groundstate energy is E_min = $(groundstate_energy)")
    end

    # funtion of energy
    function energy(kvector)
        eigenvalues = eigvals(getSpinInteractionMatrixKSpace(unitcell, kvector, bondInteractionMatrix)) .- groundstate_energy
        eigenvalues = eigenvalues .* eigenvalues
        return minimum(eigenvalues)
    end

    # the current search index
    index = 1
    # search until there are enough points
    while index <= N_points

        # find a suitable starting point for the Newton algorithm
        k = Float64[rand(), rand()]
        for j in 1:length(k)
            k[j] = k[j] * bounds_lower[j] + (1-k[j]) * bounds_upper[j]
        end

        # start with the initial energy
        e0 = energy(k)
        # iterate i over 100 newton steps (maximum)
        for i in 1:100
            # check if the energy is already converged
            if e0 < epsilon
                # save the k vector
                k_values[index,:] = k
                # increment the index
                index = index+1
                # break the newton loop
                break
            end
            # the current energy
            H_0 = e0
            # the gradient of the energy
            H_eps =
            [
                energy(k .+ [epsilon_k, 0]),
                energy(k .+ [0, epsilon_k])
            ]
            dH = (H_eps - H_0) ./ epsilon_k
            # absolute value of the gradient
            dHdH = dot(dH, dH)
            # break the Newton loop if the gradient diverges or is flattened too much
            if abs(dHdH) < 1e-20 || abs(dHdH) > 1e20
                break
            end
            # increment k
            dk = dH * (H_0 / dHdH)
            k -= dk
            # calculate a new energy
            e0 = energy(k)
        end
    end

    # if refolding is enabled, refold all points
    if refold_to_first_BZ
        # get the lattice vectors
        a1 = unitcell.lattice_vectors[1]
        a2 = unitcell.lattice_vectors[2]
        # get the reciprocal lattice vectors
        b1 = [a2[2], -a2[1]]
        b2 = [a1[2], -a1[1]]
        # normalize the vectors
        b1 .*= 2*pi/sum(b1.*a1)
        b2 .*= 2*pi/sum(b2.*a2)

        # iterate over all k points
        for index in 1:N_points
            # assumend to be refolded to begin with
            refolded = true
            # iterate while still refolding possible
            while refolded
                # new iteration, not refolded
                refolded = false
                # try all possible refoldings
                for r in [b1, b2, -b1, -b2, b1+b2, -b2-b1, b1-b2, -b2+b1]
                    if sum((k_values[index,:].-r).*(k_values[index,:].-r)) < sum((k_values[index,:]).*(k_values[index,:]))
                        k_values[index,:] .-= r
                        refolded=true
                    end
                end
            end
        end
    end

    # calculate the constraint
    for index in 1:N_points
        # the current k value
        k = k_values[index,:]
        # get the interaction matrix for this k
        matrix = getSpinInteractionMatrixKSpace(unitcell, k, bondInteractionMatrix)
        # get the matrix dimension and maybe (for d=1) simply set the constraint to be satisfied
        if size(matrix) == (1,1)
            constraint_values[index] = 0
            continue
        end
        # diagonalize the matrix
        eigenfactorization = eigfact(matrix)
        eigenvalues  = eigenfactorization[:values]
        eigenvectors = eigenfactorization[:vectors]
        # compute the constraint, first find out what bands are degenerate
        # list of bands they are degenerate with
        degenerate = zeros(Int64,   length(eigenvalues)) .- 1
        treat      =  ones(Int64,   length(eigenvalues))
        for b in 2:length(eigenvalues)
            # check if degenerate
            if eigenvalues[b] - epsilon_degenerate <= eigenvalues[b-1]
                # band b is degenerate with band b-1
                degenerate[b-1] = b
                # dont treat the band as it is treated before
                treat[b] = 0
            end
        end
        # iterate over all bands but only treat the first one
        treated_first = false
        for b in 1:length(eigenvalues)
            # if not treated, continue
            if treat[b] == 0
                continue
            elseif treated_first
                break
            end
            # if treated, compile list of all bands
            degenerate_bands = Int64[b]
            while degenerate[degenerate_bands[end]] != -1
                push!(degenerate_bands, degenerate[degenerate_bands[end]])
            end
            # get the LT constraint
            LT_constraint = getLTConstraint([eigenvectors[:,j] for j in degenerate_bands], spin_dimension)
            # save the constraint
            constraint_values[index] = LT_constraint
            # treated at least a first band
            treated_first = true
        end
    end

    # return the list
    return k_values, constraint_values

end

# Function for calculating LT Groundstate manifold in 2D (not exported)
function getLTGroundstateKSpace3D(
            unitcell::Unitcell,
            N_points::Int64,
            bondInteractionMatrix::Function=getBondInteractionMatrixHeisenberg;
            groundstate_energy::Float64=Inf,
            epsilon::Float64=1e-10,
            epsilon_k::Float64=1e-10,
            epsilon_degenerate::Float64=1e-6,
            bounds_lower::Array{Float64,1}=-2*pi.*ones(3),
            bounds_upper::Array{Float64,1}=2*pi.*ones(3),
            refold_to_first_BZ::Bool=true
        )

    # list of k values that contribute to the BZ
    k_values          = zeros(Float64, N_points, 3)
    constraint_values =  ones(Float64, N_points)
    # get the spin dimension
    spin_dimension = size(bondInteractionMatrix(unitcell.connections[1]), 1)

    # check if the GS energy has to be calculated
    if groundstate_energy == Inf
        # search for 10 random points and compare the found minimum
        for tries in 1:10
            # find a suitable starting point for the Newton algorithm
            k_start = Float64[rand(), rand(), rand()]
            for j in 1:length(k_start)
                k_start[j] = k_start[j] * bounds_lower[j] + (1-k_start[j]) * bounds_upper[j]
            end
            # optimize the energy
            groundstate_energy = min(Optim.minimum(Optim.optimize(k->minimum(eigvals(getSpinInteractionMatrixKSpace(unitcell, k, bondInteractionMatrix))), k_start)), groundstate_energy)
        end
        # print the groundstate_energy
        println("Groundstate energy is E_min = $(groundstate_energy)")
    end

    # funtion of energy
    function energy(kvector)
        eigenvalues = eigvals(getSpinInteractionMatrixKSpace(unitcell, kvector, bondInteractionMatrix)) .- groundstate_energy
        eigenvalues = eigenvalues .* eigenvalues
        return minimum(eigenvalues)
    end

    # the current search index
    index = 1
    # search until there are enough points
    while index <= N_points

        # find a suitable starting point for the Newton algorithm
        k = Float64[rand(), rand(), rand()]
        for j in 1:length(k)
            k[j] = k[j] * bounds_lower[j] + (1-k[j]) * bounds_upper[j]
        end

        # start with the initial energy
        e0 = energy(k)
        # iterate i over 100 newton steps (maximum)
        for i in 1:100
            # check if the energy is already converged
            if e0 < epsilon
                # save the k vector
                k_values[index,:] = k
                # increment the index
                index = index+1
                # break the newton loop
                break
            end
            # the current energy
            H_0 = e0
            # the gradient of the energy
            H_eps =
            [
                energy(k .+ [epsilon_k, 0, 0]),
                energy(k .+ [0, epsilon_k, 0]),
                energy(k .+ [0, 0, epsilon_k])
            ]
            dH = (H_eps - H_0) ./ epsilon_k
            # absolute value of the gradient
            dHdH = dot(dH, dH)
            # break the Newton loop if the gradient diverges or is flattened too much
            if abs(dHdH) < 1e-20 || abs(dHdH) > 1e20
                break
            end
            # increment k
            dk = dH * (H_0 / dHdH)
            k -= dk
            # calculate a new energy
            e0 = energy(k)
        end
    end

    # if refolding is enabled, refold all points
    if refold_to_first_BZ
        # get the lattice vectors
        a1 = unitcell.lattice_vectors[1]
        a2 = unitcell.lattice_vectors[2]
        a3 = unitcell.lattice_vectors[3]
        # get the reciprocal lattice vectors
        b1 = cross(a2, a3)
        b2 = cross(a3, a1)
        b3 = cross(a1, a2)
        # normalize the vectors
        b1 .*= 2*pi/sum(b1.*a1)
        b2 .*= 2*pi/sum(b2.*a2)
        b3 .*= 2*pi/sum(b3.*a3)
        # put together refolding
        refoldings = Array{Float64,1}[]
        for i in [-1,0,1]
        for j in [-1,0,1]
        for l in [-1,0,1]
            refolding = i*b1 + j*b2 + l*b3
            if sum(abs.(refolding))>1e-8
                push!(refoldings, refolding)
            end
        end
        end
        end
        # iterate over all k points
        for index in 1:N_points
            # assumend to be refolded to begin with
            refolded = true
            # iterate while still refolding possible
            while refolded
                # new iteration, not refolded
                refolded = false
                # try all possible refoldings
                for r in refoldings
                    if sum((k_values[index,:].-r).*(k_values[index,:].-r)) < sum((k_values[index,:]).*(k_values[index,:]))
                        k_values[index,:] .-= r
                        refolded=true
                    end
                end
            end
        end
    end

    # calculate the constraint
    for index in 1:N_points
        # the current k value
        k = k_values[index,:]
        # get the interaction matrix for this k
        matrix = getSpinInteractionMatrixKSpace(unitcell, k, bondInteractionMatrix)
        # get the matrix dimension and maybe (for d=1) simply set the constraint to be satisfied
        if size(matrix) == (1,1)
            constraint_values[index] = 0
            continue
        end
        # diagonalize the matrix
        eigenfactorization = eigfact(matrix)
        eigenvalues  = eigenfactorization[:values]
        eigenvectors = eigenfactorization[:vectors]
        # compute the constraint, first find out what bands are degenerate
        # list of bands they are degenerate with
        degenerate = zeros(Int64,   length(eigenvalues)) .- 1
        treat      =  ones(Int64,   length(eigenvalues))
        for b in 2:length(eigenvalues)
            # check if degenerate
            if eigenvalues[b] - epsilon_degenerate <= eigenvalues[b-1]
                # band b is degenerate with band b-1
                degenerate[b-1] = b
                # dont treat the band as it is treated before
                treat[b] = 0
            end
        end
        # iterate over all bands but only treat the first one
        treated_first = false
        for b in 1:length(eigenvalues)
            # if not treated, continue
            if treat[b] == 0
                continue
            elseif treated_first
                break
            end
            # if treated, compile list of all bands
            degenerate_bands = Int64[b]
            while degenerate[degenerate_bands[end]] != -1
                push!(degenerate_bands, degenerate[degenerate_bands[end]])
            end
            # get the LT constraint
            LT_constraint = getLTConstraint([eigenvectors[:,j] for j in degenerate_bands], spin_dimension)
            # save the constraint
            constraint_values[index] = LT_constraint
            # treated at least a first band
            treated_first = true
        end
    end

    # return the list
    return k_values, constraint_values

end

# obtain Ground state manifold
"""
    getLTGroundstateKSpace(
                unitcell::Unitcell,
                N_points::Int64,
             [  bondInteractionMatrix::Function=getBondInteractionMatrixHeisenberg;
              ; groundstate_energy::Float64=Inf,
                epsilon::Float64=1e-10,
                epsilon_k::Float64=1e-10,
                epsilon_degenerate::Float64=1e-6,
                bounds_lower::Array{Float64,1}=-2*pi.*ones(4),
                bounds_upper::Array{Float64,1}=2*pi.*ones(4),
                refold_to_first_BZ::Bool=true ]
            )

calculates `N_points` points which belong to the Luttinger Tisza ground state manifold of the spin model given by
the `Unitcell` object `unitcell` and the bond interaction matrix `bondInteractionMatrix` (function).
The energy of the ground state can be either manually given by using `groundstate_energy=...`
or is searched using optimization (by default).

The points are located using Newton's method with random starting positions. The procedure can be
modified by using some of the optional keywords:
- `epsilon` energy threshold for values to be considered close to the Fermi energy
- `epsilon_k` small k increment that is used for calculation of derivatives
- `bounds_upper` / `bounds_lower` the bounds between which the random starting location for the Newton method is searched

After the search has been completed, the values will optionally be refolded into the first brillouin zone by choosing `refold_to_first_BZ=true`.

Also the LT constraint is calculated and returned for all points.

# Examples

```julia-repl
julia> (groundstate, constraint) = getLTGroundstateKSpace(unitcell, 100)
...

julia> (groundstate, constraint) = getLTGroundstateKSpace(unitcell, 100, groundstate_energy=-4.0)
...

julia> (groundstate, constraint) = getLTGroundstateKSpace(unitcell, 100, refold_to_first_BZ=false)
...

```
"""
function getLTGroundstateKSpace(
            unitcell::Unitcell,
            N_points::Int64,
            bondInteractionMatrix::Function=getBondInteractionMatrixHeisenberg;
            groundstate_energy::Float64=Inf,
            epsilon::Float64=1e-10,
            epsilon_k::Float64=1e-10,
            epsilon_degenerate::Float64=1e-6,
            bounds_lower::Array{Float64,1}=-2*pi.*ones(4),
            bounds_upper::Array{Float64,1}=2*pi.*ones(4),
            refold_to_first_BZ::Bool=true
        )

    # check if the unitcell can be put in either method
    if length(unitcell.basis[1]) != length(unitcell.lattice_vectors)
        println("Unitcell has not the same number of lattice vectors as dimensions")
        return zeros(1,1)
    end
    # check which function to pass to
    if length(unitcell.lattice_vectors) == 2
        # 2D case
        return getLTGroundstateKSpace2D(
                unitcell,
                N_points,
                bondInteractionMatrix,
                groundstate_energy=groundstate_energy,
                epsilon=epsilon,
                epsilon_k=epsilon_k,
                epsilon_degenerate=epsilon_degenerate,
                bounds_lower=bounds_lower,
                bounds_upper=bounds_upper,
                refold_to_first_BZ=refold_to_first_BZ
            )
    else length(unitcell.lattice_vectors) == 3
        # 3D case
        return getLTGroundstateKSpace3D(
                unitcell,
                N_points,
                bondInteractionMatrix,
                groundstate_energy=groundstate_energy,
                epsilon=epsilon,
                epsilon_k=epsilon_k,
                epsilon_degenerate=epsilon_degenerate,
                bounds_lower=bounds_lower,
                bounds_upper=bounds_upper,
                refold_to_first_BZ=refold_to_first_BZ
            )
    end

end
export getLTGroundstateKSpace

# plot the Luttinger Tisza ground state manifold in 2D (not exported)
function plotLTGroundstateKSpace2D(
            k_values::Array{Float64, 2},
            constraint_values::Array{Float64, 1};
            brillouin_zone::BrillouinZone=BrillouinZone(),
            plot_title::String="",
            plot_color_valid="b",
            plot_color_invalid="r",
            figsize::Tuple=(6,6),
            showPlot::Bool=true,
            save_filename::String="NONE",
            constraint::Float64=1e-6
        )

    ###########################
    #   INITIAL SETTINGS
    ###########################

    # configure plot environment
    rc("font", family="serif")

    # create a new figure
    fig = figure(figsize=figsize)

    ###########################
    #   PLOT FERMI SURFACE
    ###########################

    # if brillouin zone not empty, plot it as well
    if length(brillouin_zone.points) > 0
        plotBrillouinZone(brillouin_zone, new_figure=false)
    end

    # scatter the points
    invalid_indices = collect(1:length(constraint_values))[constraint_values .>= constraint]
    valid_indices   = collect(1:length(constraint_values))[constraint_values  .< constraint]

    if length(invalid_indices) > 0
        scatter(
                [k_values[i,1] for i in invalid_indices],
                [k_values[i,2] for i in invalid_indices],
                color=plot_color_invalid,
                label="constrained NOT fullfilled"
            )
    end
    if length(valid_indices) > 0
        scatter(
                [k_values[i,1] for i in valid_indices],
                [k_values[i,2] for i in valid_indices],
                color=plot_color_valid,
                label="constrained fullfilled"
            )
    end

    # add a legend
    legend()

    ###########################
    #   CONFIGURE AXIS & TITLE
    ###########################

    # set the title
    if plot_title == "AUTO"
        # set the title to an automatically generated title
        title("Ground state manifold")
    elseif plot_title == ""
        # do nothing title related
    else
        # set the title to the given title
        title(plot_title)
    end

    # equal axis
    gca()[:set_aspect]("equal")

    ###########################
    #   FINISH THE PLOT
    ###########################

    # tighten the layout
    tight_layout()

    # save the plot
    if save_filename != "NONE"
        # make sure the directory exists
        if contains(save_filename, "/")
    		# get the containing folder
    		folder = save_filename[1:findlast(save_filename, '/')]
    		# build the path to that folder
    		mkpath(folder)
    	end
        # save the plot
        savefig(save_filename)
    end

    # maybe show the plot
    if showPlot
        show()
    end

    # return the figure object
    return fig

end

function plotLTGroundstateKSpace2D(
            unitcell::Unitcell,
            N_points::Int64,
            bondInteractionMatrix::Function=getBondInteractionMatrixHeisenberg;
            groundstate_energy::Float64=Inf,
            epsilon::Float64=1e-10,
            epsilon_k::Float64=1e-10,
            epsilon_degenerate::Float64=1e-6,
            bounds_lower::Array{Float64,1}=-2*pi.*ones(2),
            bounds_upper::Array{Float64,1}=2*pi.*ones(2),
            refold_to_first_BZ::Bool=true,
            brillouin_zone::BrillouinZone=BrillouinZone(),
            plot_title::String="",
            plot_color_valid="b",
            plot_color_invalid="r",
            figsize::Tuple=(6,6),
            showPlot::Bool=true,
            save_filename::String="NONE",
            constraint::Float64=1e-6
        )

    # check if BZ should be calculated
    if length(brillouin_zone.points) == 0
        brillouin_zone = createBrillouinZone(unitcell)
    end

    # calculate the Ground state manifold first
    (groundstate_points, groundstate_constraints) = getLTGroundstateKSpace2D(
            unitcell,
            N_points,
            bondInteractionMatrix,
            groundstate_energy=groundstate_energy,
            epsilon=epsilon,
            epsilon_k=epsilon_k,
            epsilon_degenerate=epsilon_degenerate,
            bounds_lower=bounds_lower,
            bounds_upper=bounds_upper,
            refold_to_first_BZ=refold_to_first_BZ
        )

    # plot the Ground state manifold
    plotLTGroundstateKSpace2D(
            groundstate_points,
            groundstate_constraints,
            brillouin_zone=brillouin_zone,
            plot_title=plot_title,
            plot_color_valid=plot_color_valid,
            plot_color_invalid=plot_color_invalid,
            figsize=figsize,
            showPlot=showPlot,
            save_filename=save_filename,
            constraint=constraint
        )

end

# plot the Luttinger Tisza ground state manifold in 3D (not exported)
function plotLTGroundstateKSpace3D(
            k_values::Array{Float64, 2},
            constraint_values::Array{Float64, 1};
            brillouin_zone::BrillouinZone=BrillouinZone(),
            plot_title::String="",
            plot_color_valid="b",
            plot_color_invalid="r",
            figsize::Tuple=(6,6),
            showPlot::Bool=true,
            save_filename::String="NONE",
            constraint::Float64=1e-6
        )

    ###########################
    #   INITIAL SETTINGS
    ###########################

    # configure plot environment
    rc("font", family="serif")

    # create a new figure
    fig = figure(figsize=figsize)

    ###########################
    #   PLOT FERMI SURFACE
    ###########################

    # if brillouin zone not empty, plot it as well
    if length(brillouin_zone.points) > 0
        plotBrillouinZone(brillouin_zone, new_figure=false)
    end

    # scatter the points
    invalid_indices = collect(1:length(constraint_values))[constraint_values .>= constraint]
    valid_indices   = collect(1:length(constraint_values))[constraint_values  .< constraint]

    if length(invalid_indices) > 0
        scatter3D(
                [k_values[i,1] for i in invalid_indices],
                [k_values[i,2] for i in invalid_indices],
                [k_values[i,3] for i in invalid_indices],
                color=plot_color_invalid,
                label="constrained NOT fullfilled"
            )
    end
    if length(valid_indices) > 0
        scatter3D(
                [k_values[i,1] for i in valid_indices],
                [k_values[i,2] for i in valid_indices],
                [k_values[i,3] for i in valid_indices],
                color=plot_color_valid,
                label="constrained fullfilled"
            )
    end

    # add a legend
    legend()

    ###########################
    #   CONFIGURE AXIS & TITLE
    ###########################

    # set the title
    if plot_title == "AUTO"
        # set the title to an automatically generated title
        title("Ground state manifold")
    elseif plot_title == ""
        # do nothing title related
    else
        # set the title to the given title
        title(plot_title)
    end

    # equal axis
    gca()[:set_aspect]("equal")

    ###########################
    #   FINISH THE PLOT
    ###########################

    # tighten the layout
    tight_layout()

    # save the plot
    if save_filename != "NONE"
        # make sure the directory exists
        if contains(save_filename, "/")
    		# get the containing folder
    		folder = save_filename[1:findlast(save_filename, '/')]
    		# build the path to that folder
    		mkpath(folder)
    	end
        # save the plot
        savefig(save_filename)
    end

    # maybe show the plot
    if showPlot
        show()
    end

    # return the figure object
    return fig

end

function plotLTGroundstateKSpace3D(
            unitcell::Unitcell,
            N_points::Int64,
            bondInteractionMatrix::Function=getBondInteractionMatrixHeisenberg;
            groundstate_energy::Float64=Inf,
            epsilon::Float64=1e-10,
            epsilon_k::Float64=1e-10,
            epsilon_degenerate::Float64=1e-6,
            bounds_lower::Array{Float64,1}=-2*pi.*ones(4),
            bounds_upper::Array{Float64,1}=2*pi.*ones(4),
            refold_to_first_BZ::Bool=true,
            brillouin_zone::BrillouinZone=BrillouinZone(),
            plot_title::String="",
            plot_color_valid="b",
            plot_color_invalid="r",
            figsize::Tuple=(6,6),
            showPlot::Bool=true,
            save_filename::String="NONE",
            constraint::Float64=1e-6
        )

    # check if BZ should be calculated
    if length(brillouin_zone.points) == 0
        brillouin_zone = createBrillouinZone(unitcell)
    end

    # calculate the Ground state manifold first
    (groundstate_points, groundstate_constraints) = getLTGroundstateKSpace3D(
            unitcell,
            N_points,
            bondInteractionMatrix,
            groundstate_energy=groundstate_energy,
            epsilon=epsilon,
            epsilon_k=epsilon_k,
            epsilon_degenerate=epsilon_degenerate,
            bounds_lower=bounds_lower,
            bounds_upper=bounds_upper,
            refold_to_first_BZ=refold_to_first_BZ
        )

    # plot the Ground state manifold
    plotLTGroundstateKSpace3D(
            groundstate_points,
            groundstate_constraints,
            brillouin_zone=brillouin_zone,
            plot_title=plot_title,
            plot_color_valid=plot_color_valid,
            plot_color_invalid=plot_color_invalid,
            figsize=figsize,
            showPlot=showPlot,
            save_filename=save_filename,
            constraint=constraint
        )

end

"""
    plotLTGroundstateKSpace(
                k_values::Array{Float64, 2},
                constraint_values::Array{Float64, 1}
             [; brillouin_zone::BrillouinZone=BrillouinZone(),
                plot_title::String="",
                plot_color_valid="b",
                plot_color_invalid="r",
                figsize::Tuple=(6,6),
                showPlot::Bool=true,
                save_filename::String="NONE",
                constraint::Float64=1e-6 ]
            )

    plotLTGroundstateKSpace(
                unitcell::Unitcell,
                N_points::Int64,
             [  bondInteractionMatrix::Function=getBondInteractionMatrixHeisenberg
              ; groundstate_energy::Float64=Inf,
                epsilon::Float64=1e-10,
                epsilon_k::Float64=1e-10,
                epsilon_degenerate::Float64=1e-6,
                bounds_lower::Array{Float64,1}=-2*pi.*ones(4),
                bounds_upper::Array{Float64,1}=2*pi.*ones(4),
                refold_to_first_BZ::Bool=true,
                brillouin_zone::BrillouinZone=BrillouinZone(),
                plot_title::String="",
                plot_color_valid="b",
                plot_color_invalid="r",
                figsize::Tuple=(6,6),
                showPlot::Bool=true,
                save_filename::String="NONE",
                constraint::Float64=1e-6
            )

Plots (or before that, calculates) `N_points` points which belong to the Luttinger Tisza ground state manifold
of the spin model given by the `Unitcell` object `unitcell` and the bond interaction matrix `bondInteractionMatrix` (function)
and plots them using `PyPlot`.

When the points are calculated, the ground state energy can be searched using optimization (DEFAULT)
or can be set by using `groundstate_energy=...`.

When calculating, the points are located using Newton's method with random starting positions. The procedure can be
modified by using some of the optional keywords:
- `enforce_hermitian` determines if the matrix used for energy calculation is made hermitian by construction
- `epsilon` energy threshold for values to be considered close to the Fermi energy
- `epsilon_k` small k increment that is used for calculation of derivatives
- `bounds_upper` / `bounds_lower` the bounds between which the random starting location for the Newton method is searched

After the search has been completed, the values will optionally be refolded into the first brillouin zone by choosing `refold_to_first_BZ=true`.

For plotting, several more options can be used.

# Examples

```julia-repl
julia> plotLTGroundstateKSpace(unitcell, 100)
PyPlot.Figure(...)

julia> plotLTGroundstateKSpace(unitcell, 100, groundstate_energy=-2.0, showPlot=false)
PyPlot.Figure(...)

```
"""
function plotLTGroundstateKSpace(
            k_values::Array{Float64, 2},
            constraint_values::Array{Float64, 1};
            brillouin_zone::BrillouinZone=BrillouinZone(Array{Float64,1}[], Array{Int64,1}[], Array{Int64,1}[]),
            plot_title::String="",
            plot_color_valid="b",
            plot_color_invalid="r",
            figsize::Tuple=(6,6),
            showPlot::Bool=true,
            save_filename::String="NONE",
            constraint::Float64=1e-6
        )

    # check which function to pass to
    if size(k_values,2) == 2
        # 2D case
        return plotLTGroundstateKSpace2D(
                k_values,
                constraint_values,
                brillouin_zone=brillouin_zone,
                plot_title=plot_title,
                plot_color_valid=plot_color_valid,
                plot_color_invalid=plot_color_invalid,
                figsize=figsize,
                showPlot=showPlot,
                save_filename=save_filename,
                constraint=constraint
            )
    else size(k_values,2) == 3
        # 3D case
        return plotLTGroundstateKSpace3D(
                k_values,
                constraint_values,
                brillouin_zone=brillouin_zone,
                plot_title=plot_title,
                plot_color_valid=plot_color_valid,
                plot_color_invalid=plot_color_invalid,
                figsize=figsize,
                showPlot=showPlot,
                save_filename=save_filename,
                constraint=constraint
            )
    end

end
function plotLTGroundstateKSpace(
            unitcell::Unitcell,
            N_points::Int64,
            bondInteractionMatrix::Function=getBondInteractionMatrixHeisenberg;
            groundstate_energy::Float64=Inf,
            epsilon::Float64=1e-10,
            epsilon_k::Float64=1e-10,
            epsilon_degenerate::Float64=1e-6,
            bounds_lower::Array{Float64,1}=-2*pi.*ones(4),
            bounds_upper::Array{Float64,1}=2*pi.*ones(4),
            refold_to_first_BZ::Bool=true,
            brillouin_zone::BrillouinZone=BrillouinZone(Array{Float64,1}[], Array{Int64,1}[], Array{Int64,1}[]),
            plot_title::String="",
            plot_color_valid="b",
            plot_color_invalid="r",
            figsize::Tuple=(6,6),
            showPlot::Bool=true,
            save_filename::String="NONE",
            constraint::Float64=1e-6
        )

    # check if BZ should be calculated
    if length(brillouin_zone.points) == 0
        brillouin_zone = createBrillouinZone(unitcell)
    end
    # check if the unitcell can be put in either method
    if length(unitcell.basis[1]) != length(unitcell.lattice_vectors)
        println("Unitcell has not the same number of lattice vectors as dimensions")
        return zeros(1,1)
    end
    # check which function to pass to
    if length(unitcell.lattice_vectors) == 2
        # 2D case
        return plotLTGroundstateKSpace2D(
                unitcell,
                N_points,
                bondInteractionMatrix,
                groundstate_energy=groundstate_energy,
                epsilon=epsilon,
                epsilon_k=epsilon_k,
                epsilon_degenerate=epsilon_degenerate,
                bounds_lower=bounds_lower,
                bounds_upper=bounds_upper,
                refold_to_first_BZ=refold_to_first_BZ,
                brillouin_zone=brillouin_zone,
                plot_title=plot_title,
                plot_color_valid=plot_color_valid,
                plot_color_invalid=plot_color_invalid,
                figsize=figsize,
                showPlot=showPlot,
                save_filename=save_filename,
                constraint=constraint
            )
    else length(unitcell.lattice_vectors) == 3
        # 3D case
        return plotLTGroundstateKSpace3D(
                unitcell,
                N_points,
                bondInteractionMatrix,
                groundstate_energy=groundstate_energy,
                epsilon=epsilon,
                epsilon_k=epsilon_k,
                epsilon_degenerate=epsilon_degenerate,
                bounds_lower=bounds_lower,
                bounds_upper=bounds_upper,
                refold_to_first_BZ=refold_to_first_BZ,
                brillouin_zone=brillouin_zone,
                plot_title=plot_title,
                plot_color_valid=plot_color_valid,
                plot_color_invalid=plot_color_invalid,
                figsize=figsize,
                showPlot=showPlot,
                save_filename=save_filename,
                constraint=constraint
            )
    end

end
export plotLTGroundstateKSpace


text:
#########################################################################
#    testing script for gradients
#########################################################################

#####       testing functions and definitions    ######

## variables of different dimension for testing
v0ref = 2.
v1ref = [2., 3, 0.1, 0, -5]
v2ref = [-1. 3 0 ; 0 5 -2]

## error thresholds
DIFF_DELTA = 1e-9
ERROR_THRESHOLD = 2e-2

good_enough(x,y) = isfinite(x) ? (abs(x-y) / max(ERROR_THRESHOLD, abs(x))) < ERROR_THRESHOLD : isequal(x,y)
good_enough(t::Tuple) = good_enough(t[1], t[2])

##  gradient check by comparing numerical gradient to automated gradient
function deriv1(ex::Expr, x0::Union(Float64, Vector{Float64}, Matrix{Float64})) #  ex= :(log(x)) ; x0 = [1.0, 2,3]
	println("testing gradient of $ex at x = $x0")

	nx = length(x0)
	model = Expr(:block, :(sum($ex)) )  # sum is there to ensure that model produces a scalar
	myf, dummy = MCMC.generateModelFunction(model, gradient=true, x=x0)

	__beta = vec([x0])
	l0, grad0 = myf(__beta)
	gradn = zeros(nx)
	for i in 1:nx
		__beta = vec([x0])
		__beta[i] += DIFF_DELTA
		l, grad = myf(__beta)
		gradn[i] = (l-l0)/DIFF_DELTA
	end

	@assert all(good_enough, zip([grad0], [gradn]))
		"Gradient false for $ex at x=$x0, expected $(round(gradn,5)), got $(round(grad0,5))"
end

## argument pattern generation for testing
# all args can be scalar, vector or matrices, but with compatible dimensions (i.e same size for arrays)
function testpattern1(parnames, rules)
	arity = length(parnames)
	ps = [ ifloor((i-1) / 2^(j-1)) % 2 for i=1:2^arity, j=1:arity]
	vcat(ps, 2*ps[2:end,:])
end

# all args are scalar
function testpattern2(parnames, rules)
	arity = length(parnames)
	zeros(Int32, 1, arity)
end

# all args are vectors
function testpattern3(parnames, rules)
	arity = length(parnames)
	ones(Int32, 1, arity)
end

# all args are matrices
function testpattern4(parnames, rules)
	arity = length(parnames)
	ones(Int32, 1,arity)*2
end

# all args can be scalar, vector or matrices, but with no more than one array
function testpattern5(parnames, rules)
	arity = length(parnames)
	ps = [ ifloor((i-1) / 2^(j-1)) % 2 for i=1:2^arity, j=1:arity]
	ps = vcat(ps, 2*ps[2:end,:])
	ps[ps * ones(arity) .<= 1, :]
end

## runs arg dimensions combinations
function runpattern(fsym, parnames, rules, combin)
	arity = length(parnames)

	for ic in 1:size(combin,1)  # try each arg dim in combin
		c = combin[ic,:]
		par = [ symbol("arg$i") for i in 1:arity]

		# create variables
		for i in 1:arity  # generate arg1, arg2, .. variables
			vn = par[i]
			vref = [:v0ref, :v1ref, :v2ref][c[i]+1]
			eval(:( $vn = copy($vref)))
		end

		# apply transformations on args
		for r in rules # r = rules[1]
			if isa(r, Expr) && r.head == :(->)
				pos = find(parnames .== r.args[1]) # find the arg the rules applies to
				@assert length(pos)>0 "arg of rule ($r.args[1]) not found among $parnames"
				vn = symbol("arg$(pos[1])")
				eval(:( $vn = map($r, $vn)))
			end
		end

		# now run tests
		prange = [1:arity]
		prange = any(rules .== :(:exceptLast)) ? prange[1:end-1] : prange
		prange = any(rules .== :(:exceptFirstAndLast)) ? prange[2:end-1] : prange
		# println("$prange - $(length(rules)) - $rules")
		for p in prange  # try each argument as parameter
			tpar = copy(par)
			tpar[p] = :x  # replace tested args with parameter symbol :x for deriv1 testing func
			f = Expr(:call, [fsym, tpar...]...)
			vn = symbol("arg$p")
			x0 = eval(vn)  # set x0 for deriv 1
			deriv1(f, x0+0.001)  # shift sightly (to avoid numerical derivation pb on max() and min())
		end
	end
end

## macro to simplify tests expression
macro mtest(pattern, func::Expr, constraints...)
	tmp = [ isa(e, Symbol) ? Expr(:quote, e) : e for e in constraints]
	quote
		local fsym = $(Expr(:quote, func.args[1]))
		local pars = $(Expr(:quote, [func.args[2:end]...]) )
		local rules = $(Expr(:quote, [tmp...]) )

		combin = ($pattern)(pars, rules)
		runpattern(fsym, pars, rules, combin)
	end
end

#########################################################################
#  Tests on functions
#########################################################################

## regular functions
@mtest testpattern1 x+y
@mtest testpattern1 x+y+z
@mtest testpattern1 sum(x)
@mtest testpattern1 x-y
@mtest testpattern1 x.*y
@mtest testpattern1 x./y  y -> y==0 ? 0.1 : y
@mtest testpattern1 x.^y  x -> x<=0 ? 0.2 : x
@mtest testpattern1 sin(x)
@mtest testpattern1 abs(x)
@mtest testpattern1 cos(x)
@mtest testpattern1 exp(x)
@mtest testpattern1 log(x) x -> x<=0 ? 0.1 : x

@mtest testpattern1 transpose(x)
deriv1(:(x'), [-3., 2, 0])

@mtest testpattern1 max(x,y)
@mtest testpattern1 min(x,y)

@mtest testpattern2 x^y

@mtest testpattern5 x/y   y->y==0 ? 0.1 : y

@mtest testpattern5 x*y
tz = transpose(v1ref)
deriv1(:(x*tz), [-3., 2, 0])
deriv1(:(tz*x), v1ref)
deriv1(:(v2ref*x), [-3., 2, 0])
deriv1(:(v2ref[:,1:2]*x), [-3. 2 0 ; 1 1 -2])

@mtest testpattern2 dot(x,y)
@mtest testpattern3 dot(x,y)

## continuous distributions
@mtest testpattern1 logpdfNormal(mu,sigma,x)  sigma->sigma<=0?0.1:sigma
@mtest testpattern1 logpdfUniform(a,b,x)      a->a-10 b->b+10
@mtest testpattern1 logpdfBeta(a,b,x)         x->clamp(x, 0.01, 0.99) a->a<=0?0.1:a b->b<=0?0.1:b
@mtest testpattern1 logpdfTDist(df,x)         df->df<=0?0.1:df
@mtest testpattern1 logpdfExponential(sc,x)   sc->sc<=0?0.1:sc  x->x<=0?0.1:x
@mtest testpattern1 logpdfGamma(sh,sc,x)      sh->sh<=0?0.1:sh  sc->sc<=0?0.1:sc  x->x<=0?0.1:x
@mtest testpattern1 logpdfCauchy(mu,sc,x)      sc->sc<=0?0.1:sc
@mtest testpattern1 logpdfLogNormal(lmu,lsc,x)  lsc->lsc<=0?0.1:lsc x->x<=0?0.1:x
@mtest testpattern1 logpdfWeibull(sh,sc,x)    sh->sh<=0?0.1:sh  sc->sc<=0?0.1:sc  x->x<=0?0.1:x

## discrete distributions
#@mtest testpattern1 logpdfBernoulli(prob,x)   exceptLast prob->clamp(prob, 0.01, 0.99) x->(x>0)+0.
# note for Bernoulli : having prob=1 or 0 is ok but will make the numeric differentiator fail => not tested

# @mtest testpattern1 logpdfPoisson(l,x)   exceptLast l->l<=0?0.1:l x->iround(abs(x))
# @mtest testpattern1 logpdfBinomial(n, prob,x)   exceptFirstAndLast prob->clamp(prob, 0.01, 0.99) x->iround(abs(x)) n->iround(abs(n)+10)

#########################################################################
#   misc. tests
#########################################################################

# Parsing should throw an error when model parameter is used as an integer variable
# try
# 	deriv1(:(logpdfBernoulli(1, x)), [0.])
# 	deriv1(:(logpdfPoisson(1, x)), [0.])
# 	deriv1(:(logpdfBinomial(3, 0.5, x)), [0.])
# 	deriv1(:(logpdfBinomial(x, 0.5, 2)), [0.])
# 	throw("no error !!")
# catch e
# 	@assert e != "no error !!"
# 		"parser not throwing error when discrete distribution has a parameter dependant integer argument"
# end

##  ref  testing
deriv1(:(x[2]),              v1ref)
deriv1(:(x[2:3]),            v1ref)
deriv1(:(x[2:end]),          v1ref)

deriv1(:(x[2:end]),          v2ref)
deriv1(:(x[2]),              v2ref)
deriv1(:(x[2:4]),            v2ref)
deriv1(:(x[:,2]),            v2ref)
deriv1(:(x[1,:]),            v2ref)
deriv1(:(x[2:end,:]),        v2ref)
deriv1(:(x[:,2:end]),        v2ref)

deriv1(:(x[2]+x[1]),          v2ref)
deriv1(:(log(x[2]^2+x[1]^2)), v2ref)

# fail case when individual elements of an array are set several times
# FIXME : correct var renaming step in unfold...
# model = :(x::real(3); y=x; y[2] = x[1] ; y ~ TestDiff())

ext: jl


text:
function binarytree!(grid::AbstractGrid)
    for cell in grid
        neighbors = Cell[]
        n = north(grid, cell)
        e = east(grid, cell)
        isnothing(n) || push!(neighbors, n)
        isnothing(e) || push!(neighbors, e)
        if !isempty(neighbors)
            target = rand(neighbors)
            link!(cell, target)
        end
    end
    return grid
end

function sidewinder!(grid::AbstractGrid)
    for row in eachrow(grid)
        run = Cell[]
        for cell in row
            push!(run, cell)
            atebound = isnothing(east(grid, cell))
            atnbound = isnothing(north(grid, cell))
            closerun = atebound || (!atnbound && iszero(rand((0, 1))))
            if closerun
                target = rand(run)
                n = north(grid, target)
                !isnothing(n) && link!(target, n)
                empty!(run)
            else
                link!(cell, east(grid, cell))
            end
        end
    end
    return grid
end

max_stars_repo_name: basp/ariadne

max_issues_repo_name: basp/ariadne

max_forks_repo_name: basp/ariadne


text:
module NNlib
using Requires

# Include APIs
include("dim_helpers.jl")

# NNPACK support
include(joinpath(@__DIR__, "..", "deps", "deps.jl"))
if check_deps() == nothing
  include("nnpack/NNPACK.jl")
else
  is_nnpack_available() = false
end

include("activation.jl")
include("softmax.jl")
include("gemm.jl")
include("conv.jl")
include("pooling.jl")

## Include implementations
include("impl/padding_edges.jl")

# Direct implementations of convolutional and depthwise-convolutional algorithms
include("impl/conv_direct.jl")
include("impl/depthwiseconv_direct.jl")
# im2col implementations of convolutional and depthwise-convolutional algorithms
include("impl/conv_im2col.jl")
include("impl/depthwiseconv_im2col.jl")

# Direct implementations of pooling
include("impl/pooling_direct.jl")

end # module NNlib


text:
#!/usr/bin/env julia

using SEM
using LinearAlgebra,Plots,UnPack

#----------------------------------#
# Set up grid
#----------------------------------#
Ex = 8; nr1 = 8;
Ey = 8; ns1 = 8;

function deform(x,y)
    x = @. 0.5*(x+1)
    y = @. 0.5*(y+1)
    return x,y
end

ifperiodic = [false,false] # [xdir, ydir]

m1 = Mesh(nr1,ns1,Ex,Ey,deform,ifperiodic)

#----------------------------------#
# case setup
#----------------------------------#
# solve with Newton Raphson approach
#
# ν∇²u + G(u,x) = f(x), x∈Ω
#     u = ub,         , x∈∂Ω
#
#----------------------------------#
x1 = m1.x
y1 = m1.y

function G(u,msh::Mesh) # G(u,x⃗)
    @unpack x,y = msh
#   Gu = @. u^3
    Gu = @. -exp.(-u)
    return Gu
end
function Gdu(u,msh::Mesh) # d/du G(u,x⃗)
    @unpack x,y = msh
#   Gdu = @. 3u^2
    Gdu = @. exp.(-u)
    return Gdu
end

ν  = @. 1+0*x1
u  = @. 0+0*x1 # initial guess (must agree with boundary data)
δu = @. 0+0*x1 # initialize
f  = @. 1+0*x1

# contrived example
kx = 3
ky = 3
ut = @. 0 + sin(kx*pi*x1)*sin(ky*pi*y1) # true solution
f  = @. ut*((kx^2+ky^2)*pi^2)         # forcing/RHS
f += G(ut,m1)

# boundary condition
# 'N': Neumann (homogeneous)
# 'D': Dirichlet (inhomogeniety applied through boundary data)
M1 = generateMask(['D','D','D','D'],m1) # [xmin,xmax,ymin,ymax]

#----------------------------------#
# solve linearized homogeneous equation for δu
#----------------------------------#
function residual(v,msh::Mesh)
    resi  = mass(f,msh)
    resi -= mass(G(v,msh),msh) .+ ν.*lapl(v,msh)
    resi .= mask(resi,M1)
    resi .= gatherScatter(resi,m1)
    return resi
end

function opHlmz(v,k,msh::Mesh) # LHS op
    Hu = hlmz(v,ν,k,msh)
    Hu = gatherScatter(Hu,msh)
    Hu = mask(Hu,M1)
    return Hu
end

for i=1:100
    local g = G(u,m1)
    local k = Gdu(u,m1)

    # LHS
    opLHS(v) = opHlmz(v,k,m1)

    # RHS
    b  = mass(f,m1)                  # forcing
    b -= ν.*lapl(u,m1) .+ mass(g,m1) # from LHS
    b .= mask(b,M1)
    b .= gatherScatter(b,m1)

    # solve
    pcg!(δu,b,opLHS,mult=m1.mult,ifv=false)
    u .+= δu

    r  = residual(u,m1);rr = norm(r,Inf);println("iter $i, res norm $rr")
    if(rr < 1e-8) break end
end
#----------------------------------#
er=ut-u; er=norm(er,Inf); println("Error ∞ norm: $er")
plt = meshplt(u,m1)
display(plt)
#----------------------------------#
nothing


text:
# plot factor graph using Makie for belief space representation

using Distributed

using Dates
using Caesar, RoME, DistributedFactorGraphs
@everywhere using Dates
@everywhere using Caesar, RoME, DistributedFactorGraphs
@everywhere using KernelDensityEstimate
@everywhere using ProgressMeter

@everywhere import RoME: Point2, Point3
using Makie
using MakieLayout
using DocStringExtensions

"""
    $SIGNATURES

Get the cartesian range over which the factor graph variables span.

Notes:
- Optional `regexFilter` can be used to subselect, according to label, which variable IDs to use.

DevNotes
- TODO, allow `tags` as filter too.
"""
function getRangeCartesian(dfg::AbstractDFG,
                           regexFilter::Union{Nothing, Regex}=nothing;
                           extend::Float64=0.2,
                           digits::Int=6,
                           xmin::Real=99999999,
                           xmax::Real=-99999999,
                           ymin::Real=99999999,
                           ymax::Real=-99999999,
                           force::Bool=false )
  #

  if !force
  # which variables to consider
  vsyms = getVariableIds(dfg, regexFilter)

  # find the cartesian range over all the vsyms variables

  for vsym in vsyms
    lran = getKDERange(getVariable(dfg, vsym) |> getKDE)
    xmin = lran[1,1] < xmin ? lran[1,1] : xmin
    ymin = lran[2,1] < ymin ? lran[2,1] : ymin
    xmax = xmax < lran[1,2] ? lran[1,2] : xmax
    ymax = ymax < lran[2,2] ? lran[2,2] : ymax
  end

  # extend the range for looser bounds on plot
  xra = xmax-xmin; xra *= extend
  yra = ymax-ymin; yra *= extend
  xmin -= extend; xmax += extend
  ymin -= extend; ymax += extend

  # clamp to nearest integers
  xmin = floor(xmin, digits=digits); xmax = ceil(xmax, digits=digits)
  ymin = floor(ymin, digits=digits); ymax = ceil(ymax, digits=digits)
  end

  return [xmin xmax; ymin ymax]
end

"""
    $SIGNATURES

2D plot of variable marginal belief estimates.

Notes:
- Uses `Makie.contour` as backend
- Disable with `fadeFloor=0`, must be ∈ [0,1].
- `autoStride::Int` if at default `varStride==-1`, else ignore `autoStride`

DevNotes
- TODO, allow `tags` as filter too.

Example
-------
```julia
fg = generateCanonicalFG_Hexagonal()
pl, Z = plotVariableBeliefs(fg, r"x\\d") # using optional Regex filter
```

Related

getRangeCartesian,
"""
function plotVariableBeliefs(dfg::AbstractDFG,
                             regexFilter::Union{Nothing, Regex}=nothing;
                             solvable::Int=0,
                             vsyms::Vector{Symbol}=listVariables(dfg, regexFilter, solvable=solvable),
                             extras::Vector{Symbol}=Symbol[],
                             N::Int=100,
                             minColorBase::Float64=-0.3,
                             maxColorBase::Float64=1.0,
                             sortVars::Bool=false,
                             varStride::Int=-1,
                             autoStride::Int=300,
                             fade::Int=0,
                             fadeExtras::Bool=false,
                             fadeFloor::Real=0.3,
                             # fadeClamp::Bool=true,
                             tail::Int=-1,
                             digits::Int=-1,
                             extend::Real=0.2,
                             resolution::Union{Nothing, Tuple{Int, Int}}=nothing,
                             xmin::Real=99999999,
                             xmax::Real=-99999999,
                             ymin::Real=99999999,
                             ymax::Real=-99999999,
                             force::Bool=false,
                             scale::Float64=1.0,
                             origin=(0,0),
                             scene=resolution==nothing ? Scene() : Scene(resolution=resolution),
                             colormap=:viridis  )
  #
  # get range over which to plot
  dfgran = getRangeCartesian(dfg, regexFilter, digits=digits, extend=extend,
                              xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, force=force)

  x = LinRange(dfgran[1,1], dfgran[1,2], N)
  y = LinRange(dfgran[2,1], dfgran[2,2], N)
  Z = zeros(N,N)
  zz = zeros(N,N)
  xy = zeros(Threads.nthreads(), 2,N)
  for i in 1:Threads.nthreads()
    xy[i,2,:] .= y
  end

  ## vsyms are variables for plotting, while applying available filters

  # specialty feature
  sortVars ? (vsyms .= vsyms |> sortDFG) : nothing
  # add additional variables that dont fit regex (AFTER SORTING)
  fadeExtras ? reverse!(vsyms) : nothing
  # add extras to front or back depending on fading request
  union!(vsyms, extras)
  fadeExtras ? reverse!(vsyms) : nothing
  varStride = varStride != -1 ? varStride : Int(floor(length(vsyms)/autoStride))+1
  sortVars && varStride != 1 ? (vsyms = vsyms[1:varStride:end]) : nothing
  sortVars && 0 < tail ? (vsyms = vsyms[end-tail:end]) : nothing
  !sortVars && varStride != 1 ? @warn("set sortVars=true to use varStride") : nothing
  !sortVars && 0 < fade ? @warn("set sortVars=true to use fade > 0") : nothing
  # walk through all variables and plotting accordingly
  len = length(vsyms)
  count = 0

  @showprogress "Evaluating symbols" for vsym in vsyms[1:end]
    if !isInitialized(getVariable(dfg, vsym))
      @warn "skipping belief plot of $vsym since not initialized"
      continue
    end
    count += 1
    XY = marginal(getVariable(dfg, vsym) |> getKDE, [1;2])
    Threads.@threads for i in 1:N
      xy[Threads.threadid(), 1,:] .= x[i]
      zz[i,:] = XY(xy[Threads.threadid(), :,:])
    end
    # normalize all beliefs to same scope
    zz ./= maximum(zz)
    # do the requested fading
    fadeval = maximum( [(fade-(len-count))/fade; 0.0]) * (1-fadeFloor) + fadeFloor
    zz .*= fadeval

    if isnan(maximum(zz)) || isnan(minimum(zz))
      @warn "skipping $vsym"
      continue
    end

    # Accumulate and clamp max value of accumulated beliefs
    Z .+= zz
    if count < len-fade && 0 < fadeFloor
      # i.e. count still in fade tail, must now do fade clamping
      Z[Z .> fadeFloor] .= fadeFloor
    else
      Z[Z .> 1] .= 1
    end
    # fade == -1 ? (Z .*= count/len) : nothing # alternatively
  end

  # set the base "background" color level by dropping one element to the desired minimum
  Z[1,end] += minColorBase
  Z[1,end-1] = maxColorBase

  # finally use Makie to draw the figure
  Makie.contour!(scene, scale.*x.+origin[1], scale.*y.+origin[2], Z, levels = 0, linewidth = 0, fillrange = true, colormap=colormap), Z
end

function addLinesBelief!(fg, pl, TTm;
                         scale::Float64=1.0,
                         origin=(0,0),
                         ppe::Bool=true, drt::Bool=true, ref::Bool=true,
                         maskRef::Tuple{Second, Second}=(Second(0),Second(0))   )
  ## This is a little excessive, but doesnt really matter
  # set all main variables and factors solvable
  setSolvable!(fg, :l1, 1)
  drtf = lsf(fg, r"drt")
  drts = ls(fg, r"drt")
  # set all necessary variables solvable
  map(v->setSolvable!(fg, v, 1), ls(fg, r"x\d"))
  vsyms = ls(fg, r"x\d", solvable=1) |> sortDFG
  # set all necessary factors solvable
  fcSo = union(map(v->setdiff(ls(fg, v), drtf), vsyms)...)
  map(v->setSolvable!(fg, v, 1), fcSo)
  # disallow factors connected to not solvables
  filter!(x->isSolved(fg, x), vsyms)
  vsyNo = setdiff(ls(fg, r"x\d",solvable=0), vsyms)
  fcNo = union(map(v->ls(fg, v), vsyNo)...)
  map(v->setSolvable!(fg, v, 0), fcNo)

  # time mask the syms
  asyms = ls(fg, r"x\d") |> sortDFG
  atims = getTimestamp.(map(x->getVariable(fg, x),asyms))
  tmask = atims .< getTimestamp(getVariable(fg, asyms[end]))

  # now draw vsyms
  # pl = plotVariableBeliefs(fg, r"x\d", vsyms=vsyms, sortVars=true, fade=15, fadeFloor=0.2)

  # overlay lines
  XYTv = map(x->getVariablePPE(getVariable(fg, x)).suggested, vsyms)
  XYT = hcat((v->v[1]).(XYTv), (v->v[2]).(XYTv), (v->v[3]).(XYTv))

  drttm = getTimestamp(getVariable(fg, asyms[1])) .< TTm .< getTimestamp(getVariable(fg, asyms[end]))
  XXmm, YYmm = XXm[drttm], YYm[drttm]

  posesyms = ls(fg, r"x\d") |> sortDFG
  filter!(x->isInitialized(fg, x), posesyms)
  filter!(x->getSolverData(getVariable(fg, x), :lbl) != nothing, posesyms)
  XXlbl = (x->(getSolverData(getVariable(fg, x), :lbl).val[1,1])).(posesyms)
  YYlbl = (x->(getSolverData(getVariable(fg, x), :lbl).val[2,1])).(posesyms)

  ts = (x->getTimestamp(getVariable(fg, x))).(posesyms) .|> datetime2unix
  T0 = ts[1]
  ts .-= ts[1]
  mask = maskRef[1].value .< ts .< maskRef[2].value
  imask = xor.(mask, true)
  @show sum(mask)

  XXlbl = XXlbl[imask]
  YYlbl = YYlbl[imask]

  # XXlbl[mask] .= Inf
  # YYlbl[mask] .= Inf

  # draw slam PPE suggested solution
  try
     !ppe ? nothing : lines!(pl, scale.*XYT[:,1].+origin[1], scale.*XYT[:,2].+origin[2], color=:black)
     !drt ? nothing : lines!(pl, scale.*XXmm.+origin[1], scale.*YYmm.+origin[2], color=:red)
     !ref ? nothing : lines!(pl, scale.*XXlbl.+origin[1], scale.*YYlbl.+origin[2], color=:black)
  catch ex
    @error ex
  end
  pl
end

0

# fg = generateCanonicalFG_Hexagonal()
# pl = plotVariableBeliefs(fg, r"x\d") # using optional Regex filter
# pl = plotVariableBeliefs(fg, r"x\d", sortVars=true, fade=3) # using optional Regex filter
# pl = plotVariableBeliefs(fg, r"x\d", sortVars=true, fade=2, tail=4) # using optional Regex filter

# targetResultsDir = "2020-01-08T17:57:25.612"
# targetResultsDir = "2020-01-15T18:30:40.98"
# fg = LightDFG{SolverParams}(params=SolverParams())
# loadDFG("/tmp/caesar/$targetResultsDir/fg_final.tar.gz", Main, fg)
# getSolverParams(fg).logpath = "/tmp/caesar/$targetResultsDir"
# dontMarginalizeVariablesAll!(fg)

# scene = plotVariableBeliefs(fg, r"x\d", sortVars=true, varStride=2, fade=20, digits=-2, fadeFloor=0.1)

# mask = YYf .< -32
# drt_data = readdlm(joinLogPath(fg, "DRT.csv"), ',')
# ... from GenerateResults.jl
# XXfm = XXf[mask]; YYfm = YYf[mask];
# lines!(scene, XXfm, YYfm, color=:red)

# xyt = getPPESuggestedAll(fg2, r"x\d")
# lines!(scene, xyt[2][:,1], xyt[2][:,2], color=:black)

# N = 20
# x = LinRange(-0.3, 1, N)
# y = LinRange(-1, 0.5, N)
# z = x .* y'
# hbox(
#    vbox(
#        contour(x, y, z, levels = 20, linewidth =3),
#        contour(x, y, z, levels = 0, linewidth = 0, fillrange = true),
#        heatmap(x, y, z),
#    ),
#    vbox(
#        image(x, y, z, colormap = :viridis),
#        surface(x, y, fill(0f0, N, N), color = z, shading = false),
#        image(-0.3..1, -1..0.5, AbstractPlotting.logo())
#    )
# )


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE
g = TypeDBClient

# Scenario: create one database
@when("connection create database: alice") do context
    client = context[:client]
    db = g.create_database(client, "alice")
    @expect db
end

@then("connection has database: alice") do context
    result = g.contains_database(context[:client], "alice")
    @expect result
end

# Scenario: create many databases
@when("connection create databases:") do context
    db_names = context.datatable
    all_db = Bool[]
    for db in db_names
        push!(all_db, g.create_database(context[:client], db[1]))
    end
    @expect all(all_db)
end

@then("connection has databases:") do context
    db_names = [db[1] for db in context.datatable]
    server_dbs = [item.name for item in g.get_all_databases(context[:client])]

    @expect (Set(db_names) == Set(server_dbs))
end

# Scenario: create many databases in parallel
@when("connection create databases in parallel:") do context
    db_names = context.datatable
    @sync @async for db in db_names
        g.create_database(context[:client], db[1])
    end
end

# Scenario: delete one database
@when("connection delete database: alice") do context
    g.delete_database(context[:client],"alice")
end

@then("connection does not have database: alice") do context
    @expect g.contains_database(context[:client], "alice") == false
end

# Scenario: connection can delete many databases
@when("connection delete databases:") do context
    db_names = [db[1] for db in context.datatable]
    for db in db_names
        g.delete_database(context[:client], db)
    end
end

@then("connection does not have databases:") do context
    db_names = [db[1] for db in context.datatable]
    db_there = Bool[]
    for db in db_names
        push!(db_there, g.contains_database(context[:client], db))
    end
    @expect all(db_there) !== true
end

# Scenario: delete many databases in parallel
@when("connection delete databases in parallel:") do context
    db_names = [db[1] for db in context.datatable]
    @sync @async for db in db_names
        g.delete_database(context[:client], db)
    end
    @expect isempty(g.get_all_databases(context[:client]))
end

# Scenario: delete a database causes open sessions to fail
@when("connection create database: typedb") do context
    g.create_database(context[:client], "typedb")
    context[:db_name] = "typedb"
end

@when("connection open session for database: typedb") do context
    context[:session] = g.CoreSession(context[:client], context[:db_name] , g.Proto.Session_Type.DATA, request_timeout=Inf)
end

@when("connection delete database: typedb") do context
    g.delete_database(context[:client], context[:db_name])
end

@then("connection does not have database: typedb") do context
    @expect g.contains_database(context[:client], context[:db_name]) === false
end

@then("session open transaction of type; throws exception: write") do context
    try
        transaction(context[:session], g.Proto.Transaction_Type.WRITE)
    catch ex
        @expect ex !== nothing
    end
end

@then("graql define; throws exception containing \"transaction has been closed\"") do context
    define_string = "define person sub entity;"
    try
        g.define(context[:transaction], define_string)
    catch ex
        @expect ex !== nothing
    end
end

@then("connection delete database; throws exception: typedb") do context
    try
        g.delete_database(context[:client], "typedb")
    catch ex
        @expect (typeof(ex) == g.TypeDBClientException)
        @expect occursin("The database typedb does not exist", string(ex))
    end
end


text:
import Base.Iterators
function load( filename )

    open(filename) do f
        return [parse(Int,x) for x in readlines(f)]
    end

end

function check( numbers, item )

    for (i, n1) in enumerate(numbers[1:end-1])
        for n2 in numbers[i+1:end]
	    n1 + n2 == item && ( return true )
        end
    end
    return false

end

function part1(puzzle, preamble)
    numbers = load(puzzle)
    size = length(numbers)
    index = preamble + 1

    while index < size
        item = numbers[index]
        previous_numbers = numbers[index-preamble:index-1]
        if !check(previous_numbers, item)
            return item
        else
            index += 1
        end
    end
end

function part2(puzzle, preamble)

    numbers = load(puzzle)
    size = length(numbers)
    invalid = part1(puzzle, preamble)
    deque = numbers[1:2]
    s = sum(deque)
    index = 3
    while s !== invalid && index < size
        if s < invalid
           push!(deque, numbers[index])
           index += 1
        elseif s > invalid
           popfirst!(deque)
        end
        s = sum(deque)
    end
    return maximum(deque) + minimum(deque)

end
println(part1("input09.txt", 25))
println(part2("input09.txt", 25))

max_stars_repo_name: pnavaro/AdventOfCode

max_issues_repo_name: pnavaro/AdventOfCode

max_forks_repo_name: pnavaro/AdventOfCode


text:
"""
`module  JAC.CoulombIonization`
    ... a submodel of JAC that contains all methods for computing Coulomb ionization cross sections and alignment parameters
        for the excitation of target or projectile electrons by fast ion impact.
"""
module CoulombIonization

    using ..Basics, ..ManyElectron, ..Radial

    """
    `struct  CoulombIonization.Settings`  ... defines a type for the details and parameters of computing Coulomb-ionization lines.

        + multipoles              ::Array{EmMultipole}           ... Specifies the multipoles of the radiation field that are to be included.
        + gauges                  ::Array{UseGauge}              ... Specifies the gauges to be included into the computations.
        + energies                ::Array{Float64,1}             ... List of ... energies.
        + calcAlignment           ::Bool                         ... True, if alignment parameters to be calculated and false otherwise.
        + printBefore  ::Bool                         ... True, if all energies and lines are printed before their evaluation.
        + selectLines             ::Bool                         ... True, if lines are selected individually for the computations.
        + selectedLines           ::Array{Tuple{Int64,Int64},1}  ... List of lines, given by tupels (inital-level, final-level).
    """
    struct Settings
        multipoles                ::Array{EmMultipole}
        gauges                    ::Array{UseGauge}
        photonEnergies            ::Array{Float64,1}
        calcAlignment             ::Bool
        printBefore    ::Bool
        selectLines               ::Bool
        selectedLines             ::Array{Tuple{Int64,Int64},1}
    end

    """
    `CoulombIonization.Settings()`  ... constructor for the default values of Coulomb-excitation line computations.
    """
    function Settings()
        Settings(EmMultipole[], UseGauge[], Float64[], false, false, false, Tuple{Int64,Int64}[])
    end

    # `Base.show(io::IO, settings::CoulombIonization.Settings)`  ... prepares a proper printout of the variable settings::CoulombIonization.Settings.
    function Base.show(io::IO, settings::CoulombIonization.Settings)
        println(io, "multipoles:               $(settings.multipoles)  ")
        println(io, "gauges:                   $(settings.gauges)  ")
        println(io, "photonEnergies:           $(settings.photonEnergies)  ")
        println(io, "calcAlignment:            $(settings.calcAlignment)  ")
        println(io, "printBefore:   $(settings.printBefore)  ")
        println(io, "selectLines:              $(settings.selectLines)  ")
        println(io, "selectedLines:            $(settings.selectedLines)  ")
    end

    """
    `struct  CoulombIonization.Channel`
        ... defines a type for a Coulomb-excitation channel to help characterize a single magnetic substate. !!! This need to adapted !!!

        + multipole      ::EmMultipole          ... Multipole of the photon absorption.
        + gauge          ::EmGauge              ... Gauge for dealing with the (coupled) radiation field.
        + kappa          ::Int64                ... partial-wave of the free electron
        + symmetry       ::LevelSymmetry        ... total angular momentum and parity of the scattering state
        + phase          ::Float64              ... phase of the partial wave
        + amplitude      ::Complex{Float64}     ... Photoionization amplitude associated with the given channel.
   """
    struct  Channel
        multipole        ::EmMultipole
        gauge            ::EmGauge
        kappa            ::Int64
        symmetry         ::LevelSymmetry
        phase            ::Float64
        amplitude        ::Complex{Float64}
    end

    """
    `struct  Line`  ... defines a type for a Coulomb-excitation line that may include the definition of channels.

        + initialLevel   ::Level                  ... initial-(state) level
        + finalLevel     ::Level                  ... final-(state) level
        + crossSection   ::EmProperty             ... Cross section for this photoionization.
        + alignmentA2    ::EmProperty             ... Alignment A_2 parameter.
        + hasChannels    ::Bool                   ... Determines whether the individual (sub-) channels are defined in terms of their
                                                      magnetic substates, etc., or not.
        + channels       ::Array{CoulombIonization.Channel,1}  ... List of CoulombIonization.Channels of this line.
    """
    struct  Line
        initialLevel     ::Level
        finalLevel       ::Level
        crossSection     ::EmProperty
        alignmentA2      ::EmProperty
        hasChannels      ::Bool
        channels         ::Array{CoulombIonization.Channel,1}
    end

    # `Base.show(io::IO, line::CoulombIonization.Line)`  ... prepares a proper printout of the variable line::CoulombIonization.Line.
    function Base.show(io::IO, line::CoulombIonization.Line)
        println(io, "initialLevel:      $(line.initialLevel)  ")
        println(io, "finalLevel:        $(line.finalLevel)  ")
        println(io, "crossSection:      $(line.crossSection)  ")
    end

end # module


text:
using TimeseriesPrediction
using Test
using Statistics, LinearAlgebra

@testset "SymmetricEmbedding" begin
    @testset "Ordering in β_groups" begin
        γ = 0; τ = 1; r = 4; c = 0; bc = ConstantBoundary(10.);
        dummy_data = [rand(10,10,10) for _ in 1:10]
        em = light_cone_embedding(dummy_data, γ, τ, r, c, bc)

        symmetries = [(Reflection(1),), (Reflection(2),), (Rotation(1,3),),
        (Rotation(1,2,3),), (Rotation(1,3), Reflection(2))]

        # symmetries = [ [[1]], [[2]], [[1,3]], [[1,2,3]], [[3,1], [2]] ]
        @testset "Symmetry $sym" for sym in symmetries
            sem = SymmetricEmbedding(em, sym)
            # Check that first entry is always origin
            @test sem.β_groups[1] == [CartesianIndex(0,0,0)]

            # Check that all points within group have the same distance to origin
            # Check that β groups are sorted
            let dist = -Inf
                for group in sem.β_groups
                    @test dist <= norm(group[1].I)
                    dist = norm(group[1].I)
                    @test 1 == length(unique(norm(getproperty.(group, :I))))
                end
            end
        end

    end
    @testset "Ordering in τ (γ=$γ)" for γ ∈ [5,10]
        #Entries are ordered with increasing τ
        τ = 1; r = 4; c = 0; bc = ConstantBoundary(10.);
        dummy_data = [rand(10,10,10,10) for _ in 1:10]
        em = light_cone_embedding(dummy_data, γ, τ, r, c, bc)

        symmetries = [(Rotation(1,4),), (Rotation(1,2,3,4),),
                    (Rotation(1,3), Reflection(2))]

        # symmetries = [ [[1,4]], [[1,2,3,4]], [[3,1], [2]] ]
        @testset "Symmetry $sym" for sym in symmetries
            sem = SymmetricEmbedding(em, sym)

            @test issorted(sem.τ)
            @test all( sem.τ .∈ Ref(em.τ))

            #Check that points are not accidentally moved to wrong timstep
            for (t,group) in zip(sem.τ, sem.β_groups)
                for g in group
                    @test (t,g) in zip(em.τ, em.β)
                end
            end
        end
    end
    @testset "None are lost" begin
        #Check that no points go missing in the reduction
        γ = 1; τ = 1; r = 3; c = 0; bc = PeriodicBoundary();
        dummy_data = [rand(10,10,10) for _ in 1:10]
        em = light_cone_embedding(dummy_data, γ, τ, r, c, bc)

        symmetries = [(Reflection(1),), (Reflection(2),), (Rotation(1,3),),
                      (Rotation(1,2,3),), (Rotation(1,3), Reflection(2))]
        # symmetries = [ [[1]], [[2]], [[1,3]], [[1,2,3]], [[3,1], [2]] ]
        @testset "Symmetry $sym" for sym in symmetries
            sem = SymmetricEmbedding(em, sym)

            all_points = vcat(sem.β_groups...)
            @test all(all_points .∈ Ref(em.β))
            @test all(em.β .∈ Ref(all_points))
        end
    end

    @testset "Parameter checking" begin
        γ = 0; τ = 1; r = 2; c = 0; bc = ConstantBoundary(10.);
        dummy_data = [rand(10,10) for _ in 1:10]
        em = light_cone_embedding(dummy_data, γ, τ, r, c, bc)

        symmetries = [(Reflection(0),), (Reflection(3),), (Rotation(1,3),),
                      (Rotation(2,0),), (Rotation(1,2), Reflection(2))]
        # symmetries = [ [[0]], [[3]], [[1,3]], [[2,0]], [[2,1], [2]] ]

        @testset "Symmetry $sym" for sym in symmetries
            @test_throws ArgumentError SymmetricEmbedding(em, sym)
        end
    end

    @testset "Reconstruction" begin
        #Quality of results depends on the system..
        @testset "$bc" for bc in [ConstantBoundary(10.), PeriodicBoundary()]
            γ = 2; τ = 5; r = 1; c = 0;
            data = [rand(25,25) for _ in 1:50]
            em = light_cone_embedding(data, γ, τ, r, c, bc)
            sem = SymmetricEmbedding(em, (Rotation(1,2),))
            # sem = SymmetricEmbedding(em, [[1,2]])

            #Check wether this works in principle
            r1 = reconstruct(data, sem);
            @test outdim(sem) == size(r1,2)
            #Check wether PCA works
            pcaem = PCAEmbedding(data, sem)
            r2 = reconstruct(data, sem);
            @test outdim(pcaem) == size(r2,2)

            @test length(r1) == length(r2)
        end
    end
end