text:
@test isempty(detect_ambiguities(EllipsisNotation))

A = Array{Int}(undef,2,4,2)

A[..,1] = [2 1 4 5
           2 2 3 6]

A[..,2] = [3 2 6 5
           3 2 6 6]

@test A[:,:,1] == [2 1 4 5
                   2 2 3 6]

@test A[:,:,2] == [3 2 6 5
                   3 2 6 6]

@test A[:,..,1] == [2 1 4 5
                    2 2 3 6]

@test A[:,..,1] == [2 1 4 5
                    2 2 3 6]

A[1,..] = reshape([3 4
                   5 6
                   4 5
                   6 7],1,4,2)

B = [3 4
     5 6
     4 5
     6 7]

@test B == reshape(A[1,..],4,2) == reshape(view(A, 1,..), 4, 2)

@test A[:,1,2] == A[..,1,2] == @view A[..,1,2]

# [..]
C = zero(B)
C[:] = B[..]
@test B == C
C[1,1] += 1
@test B != C

C = zero(B)
C[..] = B[..]
@test B == C
C[1,1] += 1
@test B != C

C = zero(B)
C[..] .= B[:]
@test B == C
C[1,1] += 1
@test B != C

C = zero(B)
C[..] .= B[..]
@test B == C
C[1,1] += 1
@test B != C

@test B[:] == B[..]
@test B[:] !== B[..]

@testset "ArrayInterface" begin
    A = Array{Int}(undef,2,4,2)
    ArrayInterface.setindex!(A, [2 1 4 5; 2 2 3 6], .., 1)
    ArrayInterface.setindex!(A, [3 2 6 5; 3 2 6 6], .., 2)

    @test ArrayInterface.getindex(A, :, :, 1) == [2 1 4 5; 2 2 3 6]
    @test ArrayInterface.getindex(A, :, :, 2) == [3 2 6 5; 3 2 6 6]

    @test ArrayInterface.getindex(A, :, .., 1) == [2 1 4 5; 2 2 3 6]
    @test ArrayInterface.getindex(A, :, .., 2) == [3 2 6 5; 3 2 6 6]

    ArrayInterface.setindex!(A, reshape([3 4; 5 6; 4 5; 6 7],1,4,2), 1, ..)

    B = [3 4
         5 6
         4 5
         6 7]

    @test B == reshape(ArrayInterface.getindex(A, 1, ..),4,2) == reshape(view(A, 1,..), 4, 2)

    @test A[:,1,2] == ArrayInterface.getindex(A,..,1,2)

    # [..]
    C = zero(B)

    C[:] = ArrayInterface.getindex(B, ..)
    @test B == C
    C[1,1] += 1
    @test B != C

    ArrayInterface.setindex!(C, ArrayInterface.getindex(B, ..), ..)
    @test B == C
    C[1,1] += 1
    @test B != C

end


text:
@testset "instantiation" begin
    sandbox = mktempdir(); bk = pwd(); cd(sandbox)

    @testset "newpage" begin
        newpage(path="foo", overwrite=true)
        @test isdir("foo")
        @test isdir(joinpath("foo", "_layout"))
        @test isfile(joinpath("foo", "config.md"))
        @test isfile(joinpath("foo", "index.md"))
        @test isfile(joinpath(".github", "workflows", "DeployPage.yml"))
    end

    @testset "cd" begin
        bk, outside = PkgPage.clever_cd("foo")
        @test endswith(bk, sandbox)
        @test outside
        @test !isdir("foo")
        bk, outside = PkgPage.clever_cd("foo")
        @test !outside
        cd("..")
    end

    @testset "serve" begin
        @test_throws Exception serve(; single=true)
        serve("foo"; single=true)
        @test isdir(joinpath("foo", "__site"))
        @test isfile(joinpath("foo", "__site", "index.html"))
        @test endswith(pwd(), sandbox)
    end

    @testset "optimize" begin
        PkgPage.optimize(; input="foo", output="bar", purge=false,
                           prerender=false)
        @test readdir(joinpath("foo", "__site")) == ["bar"]
        @test isfile(joinpath("foo", "__site", "bar", "index.html"))
    end
    cd(bk)
    Pkg.activate()
end

ext: jl


text:
module Agent
    using Agents: AbstractAgent

    # Strategy
    const C = true   # Cooperate
    const D = false  # Defect

    mutable struct Player <: AbstractAgent
        id::Int
        pos::Int
        strategy::Bool # C or D
        next_strategy::Bool # C or D
        payoff::Float32 # 6 digits

        Player(id::Int, C_rate::Float64) = new(id, id, rand() < C_rate, D, Float32(0))
        Player(id::Int, strategy::Bool) = new(id, id, strategy, D, Float32(0))
        Player(strategy::Bool) = new(-1, strategy)
    end
end


text:
 using Test
 using LogicCircuits
 using LogicCircuits: SddMgrInnerNode, SddMgrLeafNode

for (V, VI, VL) in [(PlainVtree, PlainVtreeInnerNode, PlainVtreeLeafNode),
                    (SddMgr, SddMgrInnerNode, SddMgrLeafNode)]

    @testset "$V IO test" begin

        vtree = zoo_vtree("little_4var.vtree", V)

        function test_vtree(vtree)
            @test num_variables(vtree) == 4
            @test sort(Int.(variables(vtree))) == [1,2,3,4]
            @test vtree isa VI
            @test left_most_descendent(vtree) isa VL
            @test right_most_descendent(vtree) isa VL
        end

        test_vtree(vtree)

        mktempdir() do tmp
            # Save the vtree
            temp_path = "$tmp/little_4var_temp.vtree"
            write(temp_path, vtree)

            # load from file, and then run the same tests
            vtree2 = read(temp_path, V)
            test_vtree(vtree2)
            @test PlainVtree(vtree) == PlainVtree(vtree2) # we can test equality of plain vtrees!

            mapping = read(temp_path, Dict{String,V})
            for i = 0:num_nodes(vtree)-1 # ids of vtree nodes start at 0
                @test haskey(mapping, "$i")
            end
            @test !haskey(mapping, "$(num_nodes(vtree))")

            # Save dot file
            dot_path = "$tmp/little_4var_temp.dot"
            write(dot_path, vtree)

            # Save unsupported format
            dot_path = "$tmp/little_4var_temp.bad_extension"
            @test_throws String write(dot_path, vtree)
        end

        vtree = zoo_vtree("easy/C17_mince.min.vtree", V)
        @test num_nodes(vtree) == 33
        @test num_variables(vtree) == 17
        @test vtree isa V

    end
end

max_issues_repo_name: Juice-jl/LogicCircuits

max_forks_repo_name: Juice-jl/LogicCircuits


text:
function stations(transit_map::TransitMap)
    transit_map.stations
end

function edges(transit_map::TransitMap)
    transit_map.edges
end

function nstations(transit_map::TransitMap)
    length(stations(transit_map))
end

function nedges(transit_map::TransitMap)
    length(edges(transit_map))
end

function degree(transit_map::TransitMap, station::Station)
    is_connected = x -> x.from == station || x.to == station
    length(filter(is_connected, edges(transit_map)))
end

function inbound_edges(transit_map::TransitMap, station::Station)
    is_connected = x -> x.to == station && !x.is_single_label_edge
    filter(is_connected, edges(transit_map))
end

function outbound_edges(transit_map::TransitMap, station::Station)
    is_connected = x -> x.from == station && !x.is_single_label_edge
    filter(is_connected, edges(transit_map))
end

function out_degree(transit_map::TransitMap, station::Station)
    length(outbound_edges(transit_map, station))
end

# returns a list of incident_edges edges u -> v -> w with the same line
function incident_edges(transit_map::TransitMap)
    edge_list = edges(transit_map)
    incident_edges = Vector{Tuple{GenericEdge, GenericEdge}}()
    for i in 1:nedges(transit_map)
        u = edge_list[i]
        for j in 1:nedges(transit_map)
            if i < j
                v = edge_list[j]
                if u.to == v.from && u.line == v.line
                    push!(incident_edges, (u, v))
                end
            end
        end
    end
    incident_edges
end

function edges_share_face(faces::Set{Set{GenericEdge}}, e1, e2)
    !isempty(filter(x -> e1 in x && e2 in x, faces))
end

# returns a list of faces each edge belongs to
function edge_faces(graph::TransitMap)
    faces = Set{Set{GenericEdge}}()
    for edge in graph.edges
        face_counter_clockwise = find_closest_face(graph, edge, false)
        face_clockwise = find_closest_face(graph, edge, true)
        if length(face_counter_clockwise) <= 2 && length(face_clockwise) <= 2
            continue
        end
        if face_counter_clockwise == face_clockwise
            push!(faces, face_counter_clockwise)
        else
            if length(face_counter_clockwise) > 2
                push!(faces, face_counter_clockwise)
            end
            if length(face_clockwise) > 2
                push!(faces, face_clockwise)
            end
        end
    end
    faces
end

function angle_deg(edge)
    # as on https://stackoverflow.com/a/9970297/2798441
    # by John Ericksen https://stackoverflow.com/users/654187/john-ericksen
    res = rad2deg(atan2(edge.to.coordinate.y - edge.from.coordinate.y,
            edge.to.coordinate.x - edge.from.coordinate.x))
    if res < 0
        res + 360
    else
        res
    end
end

# sorts all outbound nodes counter-clockwise
function sorted_outbound_edges(graph::TransitMap, node::EuclideanStation)
    sort(outbound_edges(graph, node), by = angle_deg)
end

struct FaceNode
    element::GenericEdge
    prev::Union{FaceNode,Nothing}
end

# starts traversing the graph clockwise or counterclockwise until
# it finds the edge again
function find_closest_face(graph::TransitMap, edge, counter_clockwise::Bool)
    queue = Queue(FaceNode)
    enqueue!(queue, FaceNode(edge, Union{FaceNode,Nothing}()))
    visited_nodes = Set{GenericEdge}()
    face = Set{GenericEdge}()
    face_found = true
    while !isempty(queue) || !face_found
        current = dequeue!(queue)
        push!(visited_nodes, current.element)
        if current.element.to.id == edge.from.id
            while current != nothing
                push!(face, current.element)
                current = get(current.prev, nothing)
            end
            face_found = true
        else
            neighbors = sorted_outbound_edges(graph, current.element.to)
            if !isempty(neighbors)
                if counter_clockwise
                    el = first(neighbors)
                else
                    el = last(neighbors)
                end
                if !(el in visited_nodes)
                    node = FaceNode(el, Nullable(current))
                    enqueue!(queue, node)
                end
            end
        end
    end
    face
end

function non_incident_edges(transit_map::TransitMap, faces::Set{Set{GenericEdge}})
    non_incident_edges = Vector{Tuple{GenericEdge, GenericEdge}}()
    edge_list = edges(transit_map)
    for i in 1:nedges(transit_map)
        u = edge_list[i]
        for j in 1:nedges(transit_map)
            if i < j
                v = edge_list[j]
                if u.from != v.from && u.from != v.to && u.to != v.from && u.to != v.to &&
                    edges_share_face(faces, u, v)
                    push!(non_incident_edges, (u, v))
                end
            end
        end
    end
    non_incident_edges
end

function find_deg2_sequences(transit_map::TransitMap, start_edge::ProcessedEdge, visited_edges::Set{ProcessedEdge})
    paths = Set()
    line = start_edge.line ## not support const declaration on local variables
    current_edge = start_edge
    current_path = Vector{ProcessedEdge}([current_edge])
    while true
        out_edges = filter(x -> x.line == line && !(x in visited_edges),
                                outbound_edges(transit_map, current_edge.to))
        is_deg2edge = degree(transit_map, current_edge.to) == 2
        if is_deg2edge
            is_deg2edge = length(out_edges) == 1
        end
        if is_deg2edge
            current_edge = first(out_edges)
            push!(current_path, current_edge)
        else
            # path has ended
            if length(current_path) > 2
                push!(paths, current_path)
            end
            visited_edges = union(visited_edges, Set(current_path))
            new_paths = map(x -> find_deg2_sequences(transit_map, x, visited_edges), out_edges)
            for np in new_paths
                union!(paths, np)
            end
            break
        end
    end
    paths
end

function classify_direction_sector(direction_angle::Real)
    if direction_angle > 90 - 22.5 && direction_angle <= 90 + 22.5
        0
    elseif direction_angle > 135 - 22.5 && direction_angle <= 135 + 22.5
        7
    elseif direction_angle > 180 - 22.5 && direction_angle <= 180 + 22.5
        6
    elseif direction_angle > 225 - 22.5 && direction_angle <= 225 + 22.5
        5
    elseif direction_angle > 270 - 22.5 && direction_angle <= 270 + 22.5
        4
    elseif direction_angle > 315 - 22.5 && direction_angle <= 315 + 22.5
        3
    elseif direction_angle > 360 - 22.5 || direction_angle <= 22.5
        2
    elseif direction_angle > 45 - 22.5 && direction_angle <= 45 + 22.5
        1
    else
        error("Should not happen")
    end
end

# reduce the network so it is easier to solve
# reduces all sequences of deg. 2 nodes to a single, longer edge
function reduce_transitmap(transit_map::InputGraph)
    new_transit_map = transit_map
    # for each line we remove degree 2 edges
    for line in transit_map.lines
        filter_fun = x -> length(inbound_edges(new_transit_map, x.from)) == 0 && x.line == line
        start_edges = filter(filter_fun, edges(new_transit_map))
        if length(start_edges) == 0
            # we have a cycle
            start_edges = filter(x -> x.line == line, edges(new_transit_map))
            if length(start_edges) == 0
                println(STDERR, "Line ", line.id, " does not have any edges.")
                exit()
            end
        end
        start_edge = first(start_edges)
        deg2seqs = find_deg2_sequences(new_transit_map, start_edge, Set{ProcessedEdge}([]))
        for seq in deg2seqs
            first_edge = first(seq)
            last_edge = last(seq)
            from = first_edge.from
            to = last_edge.to
            dir = classify_direction_sector(angle_deg(Edge(from, to, line)))
            min_length = sum(map(x -> x.min_length, seq))
            new_edge = ProcessedEdge(from, last_edge.to, line, dir, min_length, false)
            new_edges = filter(x -> !(x in seq), edges(new_transit_map))
            push!(new_edges, new_edge)
            new_nodes = unique(map(x -> x.from, new_edges) ∪ map(x -> x.to, new_edges))
            new_transit_map = InputGraph(new_nodes, new_edges, transit_map.lines)
        end
    end
    new_transit_map
end

# restores reduced nodes
function restore_transitmap(layout::TransitMapLayout, original_map::InputGraph)
    # find edges that are not part of the original layout
    # sequentially add those edges
    # It was late ... TODO: refactor
    nodes = layout.stations
    layout_node_ids = Set(map(x -> x.id, nodes))
    node_ids = layout_node_ids # remember what nodes we already added
    edges = layout.edges
    missing_edges = filter(x -> length(filter(y -> y.from.id == x.from.id &&
                                                y.to.id == x.to.id && y.line.id == x.line.id, edges)) == 0, original_map.edges)
    while length(missing_edges) > 0
        # find an edge that starts at a layout node
        edge = first(filter(x -> x.from.id in layout_node_ids, missing_edges))

        # now find the path to a known edge along the same line
        path = [edge.to]
        current_node = edge.to
        node_found = false
        rm_edges = Set([edge])
        while !node_found
            next_edges = filter(x -> x.line.id == edge.line.id &&
                                        x.from.id == current_node.id, missing_edges)
            @assert length(next_edges) <= 1
            if length(next_edges) > 0
                current_node = first(next_edges).to
                push!(rm_edges, first(next_edges))
                if !(current_node.id in layout_node_ids)
                    push!(path, current_node)
                else
                    node_found = true
                end
            else
                node_found = true
            end
        end
        filter!(x -> !(x in rm_edges), missing_edges)
        missing_nodes = filter(x -> !(x in node_ids), path)

        # now add the missing nodes and edges in equi distance
        reduced_edge = first(filter(x -> x.from.id == edge.from.id && x.line.id == edge.line.id, layout.edges))
        start_coord = reduced_edge.from.coordinate
        end_coord = reduced_edge.to.coordinate
        edge_dist = [end_coord.x - start_coord.x, end_coord.y - start_coord.y] ./ (length(path) + 1)
        i = 1
        last_station = reduced_edge.from
        for new_node in path
            if new_node.id in node_ids
                new_node = first(filter(x -> x.id == new_node.id, nodes))
            else
                new_node.coordinate = EuclideanCoordinate(start_coord.x + edge_dist[1] * i,
                                            start_coord.y + edge_dist[2] * i)
                push!(nodes, new_node)
                push!(node_ids, new_node.id)
            end
            new_edge = Edge(last_station, new_node, edge.line)
            last_station = new_node
            push!(edges, new_edge)
            i = i + 1
        end
        filter!(x -> x != reduced_edge, edges)
        # now add the final edge
        push!(edges, Edge(last_station, reduced_edge.to, edge.line))
    end
    TransitMapLayout(nodes, edges, layout.lines, layout.faces)
end

max_stars_repo_name: brainfo/MetroMap

max_issues_repo_name: brainfo/MetroMap

max_forks_repo_name: brainfo/MetroMap


text:
using JuliaPetra
using Test
using TypeStability

include("TestUtil.jl")

#function based tests
include("DenseMultiVectorTests.jl")
include("BasicDirectoryTests.jl")
include("LocalCommTests.jl")

@testset "Serial tests" begin

    #a generic serial comm for the tests that need to be called with a comm object
    comm = SerialComm{UInt64, UInt16, UInt32}()

    @testset "Util Tests" begin
        include("UtilsTests.jl")
    end

    @testset "Comm Tests" begin
        include("SerialCommTests.jl")
        include("Import-Export Tests.jl")
        include("BlockMapTests.jl")

        runLocalCommTests(serialComm)

        basicDirectoryTests(serialComm)
    end

    @testset "Data Structure Tests" begin
        denseMultiVectorTests(serialComm)

        include("SparseRowViewTests.jl")
        include("LocalCSRGraphTests.jl")
        include("LocalCSRMatrixTests.jl")

        include("CSRGraphTests.jl")
        include("CSRMatrixTests.jl")
    end
end

code = """
    $(Base.load_path_setup_code(false))
    cd("$(pwd())")
    include("mpi-runtests.jl")
"""

run(```
    mpiexec -n 4
        $(Base.julia_cmd())
        --code-coverage=$(Bool(Base.JLOptions().code_coverage) ? "user" : "none")
        --color=$(Base.have_color ? "yes" : "no")
        --compiled-modules=$(Bool(Base.JLOptions().use_compiled_modules) ? "yes" : "no")
        --check-bounds=yes
        --startup-file=$(Base.JLOptions().startupfile == 1 ? "yes" : "no")
        --track-allocation=$(("none", "user", "all")[Base.JLOptions().malloc_log + 1])
        --eval $code
```)


text:
using CUDA
using GeometricFlux
using Profile
using ProfileView

ENV["JULIA_NUM_THREADS"] = 1

d = 50
nbins = 20
l = 2^20

hist = zeros(Float32, d, nbins)
δ = rand(Float32, d, l)
idx = rand(1:nbins, l)
scatter_add!(hist, δ, idx)

@profile scatter_add!(hist, δ, idx)
Profile.print()

@profview scatter_add!(hist, δ, idx)

# sudo nvprof --profile-from-start off julia benchmark/scatter.jl
# sudo nvprof --profile-from-start off --print-gpu-trace julia --proj benchmark/scatter.jl
# sudo chown $USER -R $HOME/.julia/

# @profview scatter_add!(hist, δ, idx)
# CUDAdrv.@profile scatter_add!(hist_gpu, δ_gpu, idx_gpu)


text:
module SkeletonTriangulations

using Gridap

export SkeletonTriangulation

import Gridap: CellPoints
import Gridap: CellBasis
import Gridap: CellGeomap
import Gridap: CellRefFEs

struct SkeletonTriangulation{Z,D} <: Triangulation{Z,D}
  trian::Triangulation{Z,D}
  descriptor1::BoundaryDescriptor
  descriptor2::BoundaryDescriptor

  function SkeletonTriangulation(
    trian::Triangulation{Z,D},
    descriptor1::BoundaryDescriptor,
    descriptor2::BoundaryDescriptor) where {Z,D}
    @assert D == Z + 1
    new{Z,D}(trian,descriptor1,descriptor2)
  end

end

for op in (:CellPoints,:CellBasis,:CellGeomap,:CellRefFEs)
  @eval begin
    $op(t::SkeletonTriangulation) = $op(t.trian)
  end
end

end # module


text:
export PDataIter
export PDataLDLt, PDataKTR, PDataKARC, PDataST

abstract type TPData{T} end  # Ancestor of all PreProcess data

abstract type PDataFact{T} <: TPData{T} end # Variants using matricial factorization

abstract type PDataIter{T} <: TPData{T} end # Variants using iterative (Krylov) solvers

mutable struct PDataKARC{T} <: PDataIter{T}
    d::Array{T,1}             # (H+λI)\g ; on first call = g
    λ::T                      # "active" value of λ; on first call = 0
    ζ::T                      # Inexact Newton order parameter: stop when ||∇q|| < ξ * ||g||^(1+ζ)
    ξ::T                      # Inexact Newton order parameter: stop when ||∇q|| < ξ * ||g||^(1+ζ)
    maxtol::T                 # Largest tolerance for Inexact Newton
    mintol::T                 # Smallest tolerance for Inexact Newton

    indmin::Int               # index of best shift value  within "positive". On first call = 0

    positives::Array{Bool,1}   # indices of the shift values yielding (H+λI)⪰0
    xShift::Array{Array{T,1},1}        # solutions for each shifted system
    shifts::Array{T,1}        # values of the shifts
    nshifts::Int              # number of shifts
    norm_dirs::Array{T,1}     # norms of xShifts
    OK::Bool                  # preprocess success

    solver::CgLanczosShiftSolver
end

function PDataKARC(
    ::Type{S},
    ::Type{T},
    n;
    ζ = T(0.5),
    ξ = T(0.01),
    maxtol = T(0.01),
    mintol = (1.0e-8),
    shifts = T.(10.0 .^ collect(-20.0:1.0:20.0)),
    kwargs...,
) where {S,T}
    d = S(undef, n)
    λ = zero(T)
    indmin = 1
    nshifts = length(shifts)
    positives = Array{Bool,1}(undef, nshifts)
    xShift = Array{S,1}(undef, nshifts)
    for i = 1:nshifts
        xShift[i] = S(undef, n)
    end
    norm_dirs = S(undef, nshifts)
    OK = true
    solver = CgLanczosShiftSolver(n, n, nshifts, S)
    return PDataKARC(
        d,
        λ,
        ζ,
        ξ,
        maxtol,
        mintol,
        indmin,
        positives,
        xShift,
        shifts,
        nshifts,
        norm_dirs,
        OK,
        solver,
    )
end

mutable struct PDataKTR{T} <: PDataIter{T}
    d::Array{T,1}             # (H+λI)\g ; on first call = g
    λ::T                      # "active" value of λ; on first call = 0
    ζ::T                      # Inexact Newton order parameter: stop when ||∇q|| < ξ * ||g||^(1+ζ)
    ξ::T                      # Inexact Newton order parameter: stop when ||∇q|| < ξ * ||g||^(1+ζ)
    maxtol::T                 # Largest tolerance for Inexact Newton
    mintol::T                 # Smallest tolerance for Inexact Newton

    indmin::Int               # index of best shift value  within "positive". On first call = 0

    positives::Array{Bool,1}   # indices of the shift values yielding (H+λI)⪰0
    xShift::Array{Array{T,1},1}        # solutions for each shifted system
    shifts::Array{T,1}        # values of the shifts
    nshifts::Int              # number of shifts
    norm_dirs::Array{T,1}     # norms of xShifts
    OK::Bool                  # preprocess success

    solver::CgLanczosShiftSolver
end

function PDataKTR(
    ::Type{S},
    ::Type{T},
    n;
    ζ = T(0.5),
    ξ = T(0.01),
    maxtol = T(0.01),
    mintol = T(1.0e-8),
    shifts = T[0.0; 10.0 .^ (collect(-20.0:1.0:20.0))],
    kwargs...,
) where {S,T}
    d = S(undef, n)
    λ = zero(T)
    indmin = 1
    nshifts = length(shifts)
    positives = Array{Bool,1}(undef, nshifts)
    xShift = Array{S,1}(undef, nshifts)
    for i = 1:nshifts
        xShift[i] = S(undef, n)
    end
    norm_dirs = S(undef, nshifts)
    OK = true
    solver = CgLanczosShiftSolver(n, n, nshifts, S)
    return PDataKTR(
        d,
        λ,
        ζ,
        ξ,
        maxtol,
        mintol,
        indmin,
        positives,
        xShift,
        shifts,
        nshifts,
        norm_dirs,
        OK,
        solver,
    )
end

mutable struct PDataST{S,T} <: PDataIter{T}
    d::S
    λ::T
    ζ::T                      # Inexact Newton order parameter: stop when ||∇q|| < ξ * ||g||^(1+ζ)
    ξ::T                      # Inexact Newton order parameter: stop when ||∇q|| < ξ * ||g||^(1+ζ)
    maxtol::T                 # Largest tolerance for Inexact Newton
    mintol::T                 # Smallest tolerance for Inexact Newton

    OK::Bool    # preprocess success
    solver::CgSolver
end

function PDataST(
    ::Type{S},
    ::Type{T},
    n;
    ζ = T(0.5),
    ξ = T(0.01),
    maxtol = T(0.01),
    mintol = T(1.0e-8),
    kwargs...,
) where {S,T}
    d = S(undef, n)
    λ = zero(T)
    OK = true
    solver = CgSolver(n, n, S)
    return PDataST(d, λ, ζ, ξ, maxtol, mintol, OK, solver)
end

mutable struct PDataLDLt{T} <: PDataFact{T}
    L::Array{T,2}          # could be sparse
    D::Array{T,2}          # block diagonal 1X1 and 2X2
    pp::Array{Int,1}        # permutation vector LDL' = H[pp,pp]
    # P*v = v[pp]
    # P'*v = v[invperm(pp)]
    #s::Array{T,1}                  # Scaling matrix S=diagm(s)
    Δ::Array{T,1}          # diagonal, eigenvalues of D
    Q::Array{T,2}          # orthogonal matrix, eigenvectors of D:  should be sparse
    # QΔQ'  =  D
    g̃::Array{T,1}           # transformed gradient
    λ::T
    success::Bool                 # previous iteration was successfull
    OK::Bool                 # preprocess success

    PDataLDLt() = new{Nothing}()
    PDataLDLt(L, D, pp, Δ, Q, g, l, success, OK) =
        new{eltype(L)}(L, D, pp, Δ, Q, g, l, success, OK)
end

function PDataLDLt(::Type{S}, ::Type{T}, n; kwargs...) where {S,T}
    L = Array{T,2}(undef, n, n) # ::Array{T,2}          # could be sparse
    D = Array{T,2}(undef, n, n) #::Array{T,2}          # block diagonal 1X1 and 2X2
    pp = collect(1:n) #::Array{Int,1}        # permutation vector LDL' = H[pp,pp]
    Δ = Array{T,1}(undef, n) #::Array{T,1}          # diagonal, eigenvalues of D
    Q = Array{T,2}(undef, n, n) #::Array{T,2}          # orthogonal matrix, eigenvectors of D:  should be sparse
    g̃ = Array{T,1}(undef, n) #::Array{T,1}           # transformed gradient
    λ = zero(T) # ::T
    success = true::Bool                 # previous iteration was successfull
    OK = true
    return PDataLDLt(L, D, pp, Δ, Q, g̃, λ, success, OK)
end

mutable struct PDataSpectral{T} <: PDataFact{T}
    V::Array{T,2}           # orthogonal matrix, eigenvectors of Hessian matrix
    Δ::Array{T,1}           # eigenvalues of Hessian matrix
    g̃::Array{T,1}           # transformed gradient
    λ::T
    success::Bool           # previous iteration was successfull
    OK::Bool                # preprocess success

    PDataSpectral() = new{T}()
    PDataSpectral(V, Λ, g, l, success, OK) = new{eltype(V)}(V, Λ, g, l, success, OK)
end

function PDataSpectral(::Type{S}, ::Type{T}, n; kwargs...) where {S,T}
    V = Array{T,2}(undef, n, n) # ::Array{T,2}          # could be sparse
    Δ = Array{T,1}(undef, n) #::Array{T,1}          # diagonal, eigenvalues of D
    g̃ = Array{T,1}(undef, n) #::Array{T,1}           # transformed gradient
    λ = zero(T) # ::T
    success = true::Bool                 # previous iteration was successfull
    OK = true
    return PDataSpectral(V, Δ, g̃, λ, success, OK)
end


text:
#GracePlot demonstration 1: Plot basics
#(Uses template to avoid specifying too many parameters)
#-------------------------------------------------------------------------------

using GracePlot

#==Input data
===============================================================================#
x=collect(-10:0.1:10)
y2 = x.^2
y3 = x.^3

#==Generate plot
===============================================================================#
plot = GracePlot.new(fixedcanvas=false, emptyplot=false)
g = graph(plot, 0)
	set(g, title = "Parabolic Trajectory")
	set(g, subtitle = "(\\f{Times-Italic}y=+/- x\\f{}\\S2\\N)")
	set(g, xlabel = "Time (s)", ylabel = "Normalized height")
	#Add datasets:
		ds = add(g, x, y2)
		ds = add(g, x, -y2)
	autofit(g)

#Finalize:
redraw(plot)

#Last line


text:
using Optim, ForwardDiff, DiffResults, NLopt, NLsolve, LineSearches, PositiveFactorizations
using LinearAlgebra, BenchmarkTools, DiffResults, Roots
abstract type AbstractEquilibriumSpec end

struct AbstractHelmholtzPhases{N,V,T,W}
    phases::Int64
    molar_fractions::Vector{Vector{N}}
    volumes::Vector{V}
    temperature::T
    partition::Vector{W}
    #function AbstractHelmholtzPhases(v,T,x,phi)
     #   phas  = length(v)
     #   return new(phas,v,T,x,phi)
    #nd
end

#nn is the number of phases minus one
function _build_correction_matrix(nn)

    AA =[Diagonal(ones(nn)) Diagonal(zeros(nn))]
    BB =[Transpose(fill(10.0,nn)) Transpose(zeros(nn))]
    CC = [Diagonal(zeros(nn)) Diagonal(ones(nn))]
    DD = hcat(Transpose(zeros(nn)),Transpose(-ones(nn)))
    return dropzeros!(vcat(AA,BB,CC,DD))
end

#basis for stability testing
#TPD2(fn,x,y) = sum(ForwardDiff.gradient(fn,x),(x - y)) - fn(x) - fn(y)
TPD_vt(fn,x,y) = sum(ForwardDiff.gradient(fn,x) .* (x - y)) - fn(x) + fn(y)

dTPDdy_vt(fn,x,y) =   ForwardDiff.gradient(fn,y) .- ForwardDiff.gradient(fn,x)

function dTPDdy_vt!(fn,F,dx,y)
    F .=  ForwardDiff.gradient(fn,y) - dx
end
d2TPDdy2_vt(fn,x,y) = ForwardDiff.hessian(fn,y)

function d2TPDdy2_vt!(fn,J,x,y)
    J .= ForwardDiff.hessian(fn,y)
    #cholesky!(Positive,J)
end

function TPD_allderivatives_vt!(fn,diffresult,F,J,x,y,dx)
    ForwardDiff.hessian!(diffresult,fn,y)
    F .= DiffResults.gradient(diffresult) - dx
    J .= DiffResults.hessian(diffresult)
end

fn_vt(C,model::AbstractHelmholtzModel,T) = begin
    C2= max.(0.0,C)
    _helmholtzd(model,T,C2)
end

#stability in two phases
function st_vt(model::AbstractHelmholtzModel,V0,T0,n0)
    b0 = covolumes(model)
    n = length(n0) #length of compounds
    q = 0 #counter of random aproximations
    s = 10*n #maximum random aproximations
    k_max = 10 #max inner iterations of cholesky minimizations
    ε = 1e-06
    diff_c = 0
    yy = zeros(n) #holder of solutions
    delta = zeros(n) #holder of delta
    xx0 = zeros(n) #holder of x0
    randx0 = zeros(n) #holder of random vector
    F0 = zeros(n) #holder of function gradient
    J0 = zeros(n,n) #holder of function hessian
    yyplus1 = zeros(n) #another holder
    diffresult= DiffResults.HessianResult(yyplus1)

    for i in 1:n
        xx0[i] = n0[i]/V0
    end

    #xx0[end] = 1.0
    fn = z-> fn_vt(z,model,T0)
    dx = ForwardDiff.gradient(fn,xx0)
    tpd = z ->TPD_vt(fn,xx0,z)
    tpd0 = sum(ForwardDiff.gradient(fn,xx0) .* xx0) - fn(xx0)
    f! = (F,z) -> dTPDdy_vt!(fn,F,dx,z)
    j! = (J,z) -> d2TPDdy2_vt!(fn,J,xx0,z)
    fj! = (F,J,z) -> TPD_allderivatives_vt!(fn,diffresult,F,J,xx0,z,dx)
    norm_sol = z -> norm_vt(xx0,z)
    y_c = 0
    while true
        #building q-random aproximation to n0
        q==s && break
        q+=1
        k = 0

        randx0 .=4*rand(n).^2
        yy .= randx0 .* xx0

        normq = 1.0
        iterq = 0

        while true
            k == k_max && break
            k+=1

            fj!(F0,J0,yy)
            diff_c +=1
            cholesky!(Positive,J0)
            delta .= J0\(-F0)
            #println(yy)
            #println(delta)
            #println(yy+delta)

            #initial λ0, so maximum(yy+λ0*delta)>0
            λ0 = 1.0
            for i = 1:n
                if yy[i]+λ0*delta[i] <0
                    λ0 = (-yy[i]/delta[i])
                end
            end
            #println(λ0)
            #println("")
            λ=λ0
            #println("lambda0: ",λ0)
            iter_λ = 0
            tpdk0 = tpd(yy)
            tpdk1 = tpd(yyplus1)
            yyplus1 .= yy+λ*delta

            while iter_λ < 20
                iter_λ +=1
                λ*=0.5
                tpdk1 < tpdk0 && break
                if tpdk1 < 0
                    if (all(yyplus1.>sqrt(eps(eltype(n0)))/V0) && sum(b0.*yyplus1) <1)

                     return (yyplus1,2,q)
                    end
                end
                tpdk1 = tpd(yyplus1)
                yyplus1 .= yy+λ*delta
            end

            tol = sqrt(sum(λ*delta./yy).^2)
            tol< 1e-06 && break
            yy .=yyplus1 #new value of yy

        end #end of qth aproximation

        if tpd(yy) < 0.
            if (all(yy.>sqrt(eps(eltype(n0)))/V0) && sum(b0.*yy) <1)
                return (yy,2,q)
            end

        end
    end
    return(xx0,1,s,-1)
end #end of all aproximations

function test_st(model::AbstractHelmholtzModel,v,T,n)
    x = n/v
    res = st_vt(model,v,T,n)
    res[2] == 1 && return 0
    y  = res[1]
    βmax = min(1.0,minimum(x./y))

    f0β = β -> β*_helmholtzd(model,T,max.(0.0,y)) +
    (1-β)*_helmholtzd(model,T,max.(0.0,((x .- β*y)/((1-β))))) - _helmholtzd(model,T,x)
    f0 = z -> f0β(z[1])
    β =  Optim.optimize(f0,0,βmax).minimizer
    x1 = y
    x2 = (x - β*y)/(1-β)
    X0 = vcat(β,y...,(1-β),(x - β*y)/(1-β)...)
    #println((f0(β),βmax))

    return (β,x1,x2)
end

#evaluate a function that accepts a vector of length n, N times. this can be ForwardDiffed and accepted by an
#optimization routine
function multi_fval(fn,X,n,N)
res = 0.0
for i = 1:N
    x0 = (i-1)*n+1
    xn = i*n
    res+=fn(X[x0:xn])
end
return res
end

#evaluates helmholtz equations of the form (X,α)
#and adds the term f(X0-X)*(1-sum(α))
#this enforces the restriction of phases and matter conservation
function _multi_fval_red(fn,X,n,N,bufferₓ,X0)
    res = 0.0
    bufferₓ .= 0.0
    α₀=0.0

    for i = 1:N
        x0 = (i-1)*n+1
        xn = i*n
        bufferₓ .= bufferₓ .+ X[x0:xn]
        α₀+=X[end-n+i]
        res+=fn(X[x0:xn])*X[end-n+i]
    end
    return res + fn(X0 .- bufferₓ)*(1-α₀)
    end


text:
"""
    BlockMethod(blocks, encodings) <: LearningMethod

Learning method based on the `Block` and `Encoding` interfaces.
"""
struct BlockMethod{B, E, O} <: LearningMethod
    blocks::B
    encodings::E
    outputblock::O
end

function BlockMethod(blocks, encodings; outputblock = encodedblock(encodings, blocks[2]))
    return BlockMethod(blocks, encodings, outputblock)
end

# Core interface

function encode(method::BlockMethod, context, sample)
    encode(method.encodings, context, method.blocks, sample)
end

function encodeinput(method::BlockMethod, context, input)
    encode(method.encodings, context, method.blocks[1], input)
end

function encodetarget(method::BlockMethod, context, target)
    encode(method.encodings, context, method.blocks[2], target)
end

function decode(method::BlockMethod, context, xy)
    xyblock = encodedblock(method.encodings, method.blocks)
    decode(method.encodings, context, xyblock, xy)
end

function decodeŷ(method::BlockMethod, context, ŷ)
    decode(method.encodings, context, method.outputblock, ŷ)
end

function decodey(method::BlockMethod, context, y)
    yblock = encodedblock(method.encodings, method.blocks[2])
    decode(method.encodings, context, yblock, y)
end

# Training interface

function methodmodel(method::BlockMethod, backbone)
    xblock = encodedblock(method.encodings, method.blocks[1])
    return blockmodel(xblock, method.outputblock, backbone)
end

function methodlossfn(method::BlockMethod)
    yblock = encodedblock(method.encodings, method.blocks[2])
    return blocklossfn(method.outputblock, yblock)
end

# Testing interface

mocksample(method::BlockMethod) = mockblock(method.blocks)

mockmodel(method::BlockMethod) = mockmodel(
    encodedblock(method.encodings, method.blocks[1]),
    method.outputblock
)

function mockmodel(inblock::Block, outblock::Block)
    return function mockmodel_block(xs)
        out = mockblock(outblock)
        DataLoaders.collate([out])
    end
end

# Pretty-printing

function Base.show(io::IO, method::BlockMethod)
    print(io, "BlockMethod(", typeof(method.blocks[1]), " -> ", typeof(method.blocks[2]), ")")
end


text:
@testset "sparse_to_csr" begin
    ei = [1;2;3]
    ej = [3;4;4]
    ev = [8;9;10]
    rp,ci,ai,m = sparse_to_csr(ei,ej,ev)

    # more tests added here
    for t = 1:100
        A = sprand(100,80,0.01)
        (rp,ci,ai,m)=sparse_to_csr(A)
        ii = zeros(Int64,length(ai))
        j = zeros(Int64,length(ai))
        a = zeros(Float64,length(ai))
        n = length(rp)-1
        nz = 0
        for cr = 1:n
            for ri = rp[cr]:rp[cr+1]-1
                nz=nz+1
                ii[nz]=cr
                j[nz]=ci[ri]
                a[nz]=ai[ri]
            end
        end
        A2 = sparse(ii,j,a,n,m)
        @test A == A2
    end
end


text:
module GoogleCharts

using Mustache
using JSON
using DataArrays
using DataFrames
using Dates

import Base.show, Base.display

include("utils.jl")
include("charts.jl")
include("render.jl")
include("conveniences.jl")

tpl_name = joinpath(dirname(@__FILE__),"../tpl/chart.html")
global chart_tpl = Mustache.template_from_file(tpl_name)

export GoogleChart
export render, plot, scatter, help_on_chart
export area_chart, bar_chart, bubble_chart, candlestick_chart, column_chart, combo_chart,
       gauge_chart, geo_chart, line_chart, pie_chart, scatter_chart, stepped_area_chart,
       table_chart, tree_chart, annotated_time_line, intensity_map, motion_chart, org_chart,
       image_spark_line
export boxplot, histogram
export surfaceplot

end ## module


text:
a  = read_SSR(joinpath(dirname(@__FILE__), "..", "data", "test_Hz19.5-testing.bdf"))

@fact a.processing["Side"] --> "Bilateral"
@fact a.processing["Name"] --> "P2"
@fact a.processing["Amplitude"] --> 70.0
@fact a.processing["Carrier_Frequency"] --> 1000.0
@fact modulationrate(a) --> 40.0390625


text:
################################################################################
#
#   IMPLEMENTATIONS OF DIFFERENT LATTICE MODIFICATION FUNCTIONS
#   (MOST OF THEM ALSO WORK FOR UNITCELLS)
#
#   STRUCTURE OF THE FILE
#
#   1) CONNECTIONS ADDING REMOVING
#       - Add a new connection
#       - remove connections (based on indices)
#       - remove connections (based on strength)
#
#   2) SITES ADDING REMOVING
#       - add a new site
#       - remove a site (based on index)
#
#   3) CONNECTION STRENGTH MODIFICATION
#       - all connections
#       - mapping of strengths
#       - evaluate strengths
#       - optimize connections
#
#   4) REAL SPACE ROTATION / SCALING / SHIFTING
#       - Rotation around some axis
#       - Scaling along some axis
#       - Global (isotropic) scaling
#       - Global (isotropic) scaling to mean / minimum bond length
#       - Shifting along some axis
#       - Shifting along some lattice vector
#       - Shifting center of lattice to origin
#
#   5) LATTICE VECTOR OPTIMIZATION
#
################################################################################

################################################################################
#
#   1) CONNECTIONS ADDING REMOVING
#       - Add a new connection
#       - remove connections (based on indices)
#       - remove connections (based on strength)
#
################################################################################

#-------------------------------------------------------------------------------
#
#  Add connection to given unitcell / lattice objects
#
#-------------------------------------------------------------------------------
"""
    addConnection!(unitcell::Unitcell, index_from::Int64, index_to::Int64, strength [, wrap=-1; overwrite::Bool=false])
    addConnection!(lattice::Lattice,   index_from::Int64, index_to::Int64, strength [, wrap=-1; overwrite::Bool=false])

Function to add a connection to either a `Unitcell` object or a `Lattice` object.
Necessary parameters include the indices of sites between which the connection takes place as well as the strength of the connection.
Additionally (in case of a unitcell or periodic boundaries), a wrap can be specified in terms of lattice vectors.
If no wrap is specified, the wrap will be chosen automatically as without wrapping around periodic boundaries.

The additional option `overwrite` specifies if one wants to simply add the connection
or if one wants to first search if there is a connection that already satisfies the new connection (in which case nothing would be added).

NOTE 1: This function always tries to add two connections, one going forward, one going backward.

NOTE 2: This function modifies the given object and does not create a copy.

# Examples

```julia-repl
julia> addConnection!(unitcell, 1, 3, "tx", (0,1))

julia> addConnection!(lattice, 31, 27, "tx")

```
"""
function addConnection!(unitcell::Unitcell, index_from::Int64, index_to::Int64, strength, wrap=-1; overwrite::Bool=false)

    # maybe default wrap
    if wrap == -1
        if length(unitcell.lattice_vectors) == 3
            wrap = (0,0,0)
        elseif length(unitcell.lattice_vectors) == 2
            wrap = (0,0)
        else
            wrap = (0)
        end
    end

    # construct the returning wrap
    if length(wrap) == 3
        wrap_return = (-wrap[1],-wrap[2],-wrap[3])
    elseif length(wrap) == 2
        wrap_return = (-wrap[1],-wrap[2])
    else
        wrap_return = (-wrap)
    end
    # construct the returning strength
    if typeof(strength) == Complex
        strength_return = deepcopy(conj(strength))
    else
        strength_return = deepcopy(strength)
    end

    # construct the connection as an array
    connection_1 = Any[index_from; index_to; strength; wrap]
    # construct the returning connection as an array
    connection_2 = Any[index_to; index_from; strength_return; wrap_return]

    # if not overwrite, check if the connection exists
    if !overwrite
        # report if found
        found_c1 = false
        found_c2 = false
        # check in all connections
        for c in unitcell.connections
            if  Int(c[1]) == Int(connection_1[1]) && Int(c[2]) == Int(connection_1[2]) && c[3] == connection_1[3] && c[4] == connection_1[4]
                found_c1 = true
            end
            if  Int(c[1]) == Int(connection_2[1]) && Int(c[2]) == Int(connection_2[2]) && c[3] == connection_2[3] && c[4] == connection_2[4]
                found_c2 = true
            end
        end
        # only add if not added already
        if !found_c1
            push!(unitcell.connections, connection_1);
        end
        if !found_c2
            push!(unitcell.connections, connection_2);
        end
    # otherwise, just add the connections
    else
        push!(unitcell.connections, connection_1);
        push!(unitcell.connections, connection_2);
    end

    # return nothing
    return nothing
end
function addConnection!(lattice::Lattice, index_from::Int64, index_to::Int64, strength, wrap=-1; overwrite::Bool=false)

    # maybe default wrap
    if wrap == -1
        if length(lattice.lattice_vectors) == 3
            wrap = (0,0,0)
        elseif length(lattice.lattice_vectors) == 2
            wrap = (0,0)
        else
            wrap = (0)
        end
    end

    # construct the returning wrap
    if length(wrap) == 3
        wrap_return = (-wrap[1],-wrap[2],-wrap[3])
    elseif length(wrap) == 2
        wrap_return = (-wrap[1],-wrap[2])
    else
        wrap_return = (-wrap)
    end
    # construct the returning strength
    if typeof(strength) == Complex
        strength_return = deepcopy(conj(strength))
    else
        strength_return = deepcopy(strength)
    end

    # construct the connection as an array
    connection_1 = Any[index_from; index_to; strength; wrap]
    # construct the returning connection as an array
    connection_2 = Any[index_to; index_from; strength_return; wrap_return]

    # if not overwrite, check if the connection exists
    if !overwrite
        # report if found
        found_c1 = false
        found_c2 = false
        # check in all connections
        for c in lattice.connections
            if  Int(c[1]) == Int(connection_1[1]) && Int(c[2]) == Int(connection_1[2]) && c[3] == connection_1[3] && c[4] == connection_1[4]
                found_c1 = true
            end
            if  Int(c[1]) == Int(connection_2[1]) && Int(c[2]) == Int(connection_2[2]) && c[3] == connection_2[3] && c[4] == connection_2[4]
                found_c2 = true
            end
        end
        # only add if not added already
        if !found_c1
            push!(lattice.connections, connection_1);
        end
        if !found_c2
            push!(lattice.connections, connection_2);
        end
    # otherwise, just add the connections
    else
        push!(lattice.connections, connection_1);
        push!(lattice.connections, connection_2);
    end

    # return nothing
    return nothing
end
export addConnection!

#-------------------------------------------------------------------------------
#
#   Remove interactions with given strengths of the unitcell / lattice
#
#-------------------------------------------------------------------------------
"""
    removeConnections!(unitcell::Unitcell, index_from::Int64, index_to::Int64 [, strength ])
    removeConnections!(lattice::Lattice,   index_from::Int64, index_to::Int64 [, strength ])

Function to remove connections from either a `Unitcell` object or a `Lattice` object based on the sites that they connect.
Optionally, one can pass a connection strength to also filter based on this strength. If no strength is passed, all connnections
between the two sites will be removed.

    removeConnections!(unitcell::Unitcell  [, strength])
    removeConnections!(lattice::Lattice    [, strength])

Function to remove connections from either a `Unitcell` object or a `Lattice` object only based on the strength of connections.
If no strength is passed, all connnections with Float64 strength `0.0` will be removed.

NOTE 1: This function always tries to delte both connections, one going forward, one going backward.

NOTE 2: This function modifies the given object and does not create a copy.

# Examples

```julia-repl
julia> removeConnections!(unitcell, 1, 3, "tx")

julia> removeConnections!(unitcell, 1, 3)

julia> removeConnections!(lattice, 23, 73)

julia> removeConnections!(lattice, "tx")

```
"""
function removeConnections!(lattice::Lattice,   index_from::Int64, index_to::Int64, strength)
    # use the julia filter function to filter out all elements which have the strength passed and the fitting indices
    filter!(x->!(x[3]==strength && Int(x[1])==index_from && Int(x[2])==index_to) && !(x[3]==strength && Int(x[2])==index_from && Int(x[1])==index_to), lattice.connections)
    # return nothing
    return nothing
end
function removeConnections!(unitcell::Unitcell, index_from::Int64, index_to::Int64, strength)
    # use the julia filter function to filter out all elements which have the strength passed and the fitting indices
    filter!(x->!(x[3]==strength && Int(x[1])==index_from && Int(x[2])==index_to) && !(x[3]==strength && Int(x[2])==index_from && Int(x[1])==index_to), unitcell.connections)
    # return nothing
    return nothing
end
function removeConnections!(lattice::Lattice,   index_from::Int64, index_to::Int64)
    # use the julia filter function to filter out all elements which have the fitting indices
    filter!(x->!(Int(x[1])==index_from && Int(x[2])==index_to) && !(Int(x[2])==index_from && Int(x[1])==index_to), lattice.connections)
    # return nothing
    return nothing
end
function removeConnections!(unitcell::Unitcell, index_from::Int64, index_to::Int64)
    # use the julia filter function to filter out all elements which have the fitting indices
    filter!(x->!(Int(x[1])==index_from && Int(x[2])==index_to) && !(Int(x[2])==index_from && Int(x[1])==index_to), unitcell.connections)
    # return nothing
    return nothing
end
function removeConnections!(lattice::Lattice,   strength=0.0)
    # use the julia filter function to filter out all elements which have the strength passed
    filter!(x->x[3]!=strength, lattice.connections)
    # return nothing
    return nothing
end
function removeConnections!(unitcell::Unitcell, strength=0.0)
    # use the julia filter function to filter out all elements which have the strength passed
    filter!(x->x[3]!=strength, unitcell.connections)
    # return nothing
    return nothing
end

export removeConnections!

################################################################################
#
#   2) SITES ADDING REMOVING
#       - add a new site
#       - remove a site (based on index)
#
################################################################################

"""
    addSite!(unitcell::Unitcell, position::Array{Float64,1})
    addSite!(lattice::Lattice,   position::Array{Float64,1})

Function to add a site to either a `Unitcell` object or a `Lattice` object which is located at
the real space position `position`. The function returns the index of the newly added site for further use.

NOTE 1: The newly added site is not connected to any other sites yet.

NOTE 2: This function modifies the given object and does not create a copy.

# Examples

```julia-repl
julia> addSite!(unitcell, Float64[0.0, 0.5])
5

julia> addSite!(lattice, Float64[0.0, 0.5, 1.0])
187
```
"""
function addSite!(lattice::Lattice, position::Array{Float64,1})
    # push the site into the positions list
    push!(lattice.positions, position)
    # return the site index
    return length(lattice.positions)
end
function addSite!(unitcell::Unitcell, position::Array{Float64,1})
    # push the site into the positions list
    push!(unitcell.basis, position)
    # return the site index
    return length(unitcell.basis)
end
export addSite!

"""
    removeSite!(unitcell::Unitcell, index::Int64)
    removeSite!(lattice::Lattice,   index::Int64)

Function to remove a site from either a `Unitcell` object or a `Lattice` object with index `index`.
The function furhtermore removes all connections that are connected to the site and relabels all
other connections as well to account for the index shift introduced by the removal.

NOTE: This function modifies the given object and does not create a copy.

# Examples

```julia-repl
julia> removeSite!(unitcell, 2)

julia> removeSite!(lattice, 128)

```
"""
function removeSite!(lattice::Lattice, index::Int64)
    # delete the element from the positions list
    deleteat!(lattice.positions, index)
    # delete all connections that touch the removed site
    deleteat!(lattice.connections, [c[1]==index || c[2]==index for c in lattice.connections])
    # now shift all indices in the connections array by 1
    for (i,c) in enumerate(lattice.connections)
        if Int(c[1]) > index
            c[1] = Int(c[1]) - 1
        end
        if Int(c[2]) > index
            c[2] = Int(c[2]) - 1
        end
    end
    # return nothing
    return nothing
end
function removeSite!(unitcell::Unitcell, index::Int64)
    # delete the element from the positions list
    deleteat!(unitcell.basis, index)
    # delete all connections that touch the removed site
    deleteat!(unitcell.connections, [c[1]==index || c[2]==index for c in unitcell.connections])
    # now shift all indices in the connections array by 1
    for (i,c) in enumerate(unitcell.connections)
        if Int(c[1]) > index
            c[1] = Int(c[1]) - 1
        end
        if Int(c[2]) > index
            c[2] = Int(c[2]) - 1
        end
    end
    # return nothing
    return nothing
end
export removeSite!

################################################################################
#
#   3) CONNECTION STRENGTH MODIFICATION
#       - all connections
#       - mapping of strengths
#       - evaluate strengths
#       - optimize connections
#
################################################################################

#-----------------------------------------------------------------------------------------------------------------------------
#
#   Set / Map interaction strengths of the unitcell / lattice
#
#-----------------------------------------------------------------------------------------------------------------------------
"""
    setAllInteractionStrengths!(unitcell::Unitcell, strength)
    setAllInteractionStrengths!(lattice::Lattice,   strength)

Function to overwrite ALL connection strength values of either a `Unitcell` object or a `Lattice` object
with the value `strength`.

NOTE: This function modifies the given object and does not create a copy.

# Examples

```julia-repl
julia> setAllInteractionStrengths!(unitcell, "tx")

julia> setAllInteractionStrengths!(lattice, 1.0)

```
"""
function setAllInteractionStrengths!(unitcell::Unitcell, strength)
    # go through all connections and modify their strength
    for c in unitcell.connections
        c[3] = strength
    end
    # return nothing
    return nothing
end
function setAllInteractionStrengths!(lattice::Lattice,   strength)
    # go through all connections and modify their strength
    for c in lattice.connections
        c[3] = strength
    end
    # return nothing
    return nothing
end
export setAllInteractionStrengths!

"""
    mapInteractionStrengths!(unitcell::Unitcell, mapping::Dict [; replace_in_strings::Bool=false, evaluate::Bool=false])
    mapInteractionStrengths!(lattice::Lattice,   mapping::Dict [; replace_in_strings::Bool=false, evaluate::Bool=false])

Function to map connection strength values of either a `Unitcell` object or a `Lattice` object
by replacing them along keys of a dictonary `mapping`.

Furthermore, to be able to resolve strengths like `"tx*ty"` (`String` valued) to a `Float64` value one can use the following two steps:
1) chosing `replace_in_strings=true` to also search for mapping keys in `String` connection strengths.
2) chosing `evaluate=true` to evaluate every `String` connection strength by `eval(parse(.))` after all replacement has taken place.

NOTE: This function modifies the given object and does not create a copy.

# Examples

Example mapping could read

    mapping = Dict("tx"=>1.0, "ty"=>1.0, "tz"=>2.0)

which modifies all `String` valued connection strengths `"tx"`, `"ty"` and `"tz"` to be `Float64` with values `1.0` and `2.0`

```julia-repl
julia> mapInteractionStrengths!(unitcell, mapping)

julia> mapInteractionStrengths!(lattice, mapping, replace_in_strings=true)

```
"""
function mapInteractionStrengths!(unitcell::Unitcell, mapping::Dict; replace_in_strings::Bool=false, evaluate::Bool=false)
    # get the old strengths
    old_strengths = keys(mapping)
    # iterate for replacement
    for old_strength in old_strengths
        # new strength is given by mapping
        new_strength = mapping[old_strength]
        # check all connections that are identical
        for c in unitcell.connections
            if c[3] == old_strength
                c[3] = new_strength
            end
        end
        # check for string replacement
        if replace_in_strings
            for c in unitcell.connections
                if typeof(c[3]) == String && contains(c[3], old_strength)
                    c[3] = replace(c[3], old_strength, new_strength)
                end
            end
        end
    end
    # maybe even evaluate
    if evaluate
        for c in unitcell.connections
            if typeof(c[3]) == String
                c[3] = eval(parse(c[3]))
            end
        end
    end
    # return nothing
    return nothing
end
function mapInteractionStrengths!(lattice::Lattice,   mapping::Dict; replace_in_strings::Bool=false, evaluate::Bool=false)
    # get the old strengths
    old_strengths = keys(mapping)
    # iterate for replacement
    for old_strength in old_strengths
        # new strength is given by mapping
        new_strength = mapping[old_strength]
        # check all connections that are identical
        for c in lattice.connections
            if c[3] == old_strength
                c[3] = new_strength
            end
        end
        # check for string replacement
        if replace_in_strings
            for c in lattice.connections
                if typeof(c[3]) == String && contains(c[3], old_strength)
                    c[3] = replace(c[3], old_strength, new_strength)
                end
            end
        end
    end
    # maybe even evaluate
    if evaluate
        for c in lattice.connections
            if typeof(c[3]) == String
                c[3] = eval(parse(c[3]))
            end
        end
    end
    # return nothing
    return nothing
end
export mapInteractionStrengths!

"""
    evaluateInteractionStrengths!(unitcell::Unitcell [, parameters::Dict])
    evaluateInteractionStrengths!(lattice::Lattice   [, parameters::Dict])

Function to evaluate connection strength values of either a `Unitcell` object or a `Lattice` object
by first optionally replacing them along keys of a dictonary `parameters` and then calling `eval(parse(.))`.

NOTE: This function modifies the given object and does not create a copy.

# Examples

Example parameters could read

    parameters = Dict("tx"=>1.0, "ty"=>1.0, "tz"=>2.0)

which modifies all `String` valued connection strengths `"tx"`, `"ty"` and `"tz"` to be `Float64` with values `1.0` and `2.0`

```julia-repl
julia> evaluateInteractionStrengths!(unitcell)

julia> evaluateInteractionStrengths!(lattice, parameters)

```
"""
function evaluateInteractionStrengths!(unitcell::Unitcell, parameters::Dict=Dict())
    # just call the mapping functions
    mapInteractionStrengths!(unitcell, parameters, replace_in_strings=true, evaluate=true)
end
function evaluateInteractionStrengths!(lattice::Lattice,   parameters::Dict=Dict())
    # just call the mapping functions
    mapInteractionStrengths!(lattice, parameters, replace_in_strings=true, evaluate=true)
end
export evaluateInteractionStrengths!

#-----------------------------------------------------------------------------------------------------------------------------
#
#   Optimize the lattice / unitcell by collapsing redundant connections into fewer connections and modify given object
#
#-----------------------------------------------------------------------------------------------------------------------------

"""
    optimizeConnections!(unitcell::Unitcell)
    optimizeConnections!(lattice::Lattice)

Function to optimize connections of either a `Unitcell` object or a `Lattice` object
by checking which ones can be compressed and put together.

NOTE: This function modifies the given object and does not create a copy.

# Examples

```julia-repl
julia> optimizeConnections!(unitcell)

```
"""
function optimizeConnections!(lattice::Lattice)
    # build up new connections
    connections_new = Array{Any,1}[]
    # go through all old connections to check for redundant entries
    for c in lattice.connections
        # check if the connection already is present in the new list
        found = false
        for c_new in connections_new
            if c[1] == c_new[1] && c[2] == c_new[2] && c[4] == c_new[4]
                # found a redundent connection
                found = true
                # add connection strength
                if typeof(c[3]) == String || typeof(c_new[3]) == String
                    c_new[3] = "$(c_new[3])+$(c[3])"
                else
                    c_new[3] = c_new[3] + c[3]
                end
                # break the inner loop
                break
            end
        end
        # if not, add the current connection
        if !found
            push!(connections_new, c)
        end
    end
    # check if the connection strength is 0
    connections_new_2 = Array{Any,1}[]
    for c in connections_new
        # check if close to 0 connection strength
        if typeof(c[3]) == Float64 || typeof(c[3]) == Int64
            if abs(c[3]) > 1e-18
                push!(connections_new_2, c)
            else
                # do not push! c[3] too small!
            end
        elseif isnull(tryparse(Float64,c[3]))
            # not a float and cannot be parsed to float, just push in list as it is
            push!(connections_new_2, c)
        else
            # check if the strength is non-zero
            if abs(parse(Float64,c[3])) > 1e-18
                push!(connections_new_2, c)
            else
                # do not push! c[3] too small!
            end
        end
    end
    # overwrite the connections in the given lattice
    lattice.connections = connections_new_2
    # return nothing
    return nothing
end
function optimizeConnections!(unitcell::Unitcell)
    # build up new connections
    connections_new = Array{Any,1}[]
    # go through all old connections to check for redundant entries
    for c in unitcell.connections
        # check if the connection already is present in the new list
        found = false
        for c_new in connections_new
            if c[1] == c_new[1] && c[2] == c_new[2] && c[4] == c_new[4]
                # found a redundent connection
                found = true
                # add connection strength
                if typeof(c[3]) == String || typeof(c_new[3]) == String
                    c_new[3] = "$(c_new[3])+$(c[3])"
                else
                    c_new[3] = c_new[3] + c[3]
                end
                # break the inner loop
                break
            end
        end
        # if not, add the current connection
        if !found
            push!(connections_new, c)
        end
    end
    # check if the connection strength is 0
    connections_new_2 = Array{Any,1}[]
    for c in connections_new
        # check if close to 0 connection strength
        if typeof(c[3]) == Float64 || typeof(c[3]) == Int64
            if abs(c[3]) > 1e-18
                push!(connections_new_2, c)
            else
                # do not push! c[3] too small!
            end
        elseif isnull(tryparse(Float64,c[3]))
            # not a float and cannot be parsed to float, just push in list as it is
            push!(connections_new_2, c)
        else
            # check if the strength is non-zero
            if abs(parse(Float64,c[3])) > 1e-18
                push!(connections_new_2, c)
            else
                # do not push! c[3] too small!
            end
        end
    end
    # overwrite the connections in the given lattice
    unitcell.connections = connections_new_2
    # return nothing
    return nothing
end
export optimizeConnections!

################################################################################
#
#   4) REAL SPACE ROTATION / SCALING / SHIFTING
#       - Rotation around some axis
#       - Scaling along some axis
#       - Global (isotropic) scaling
#       - Global (isotropic) scaling to mean / minimum bond length
#       - Shifting along some axis
#       - Shifting along some lattice vector
#       - Shifting center of lattice to origin
#
################################################################################

###################
#   ROTATION
###################

# ROTATION AROUND X AXIS
function rotateAroundXAxis!(unitcell::Unitcell, angle::Float64)
    # rotate the lattice vectors
    for lv in unitcell.lattice_vectors
        # rotate and overwrite simulateously all components
        lv[2], lv[3] = cos(angle)*lv[2] + sin(angle)*lv[3]  ,   -sin(angle)*lv[2] + cos(angle)*lv[3]
    end
    # rotate the basis sites
    for bs in unitcell.basis
        # rotate and overwrite simulateously all components
        bs[2], bs[3] = cos(angle)*bs[2] + sin(angle)*bs[3]  ,   -sin(angle)*bs[2] + cos(angle)*bs[3]
    end
end
function rotateAroundXAxis!(lattice::Lattice, angle::Float64)
    # rotate the lattice vectors
    for lv in lattice.lattice_vectors
        # rotate and overwrite simulateously all components
        lv[2], lv[3] = cos(angle)*lv[2] + sin(angle)*lv[3]  ,   -sin(angle)*lv[2] + cos(angle)*lv[3]
    end
    # rotate the basis sites
    for p in lattice.positions
        # rotate and overwrite simulateously all components
        p[2], p[3] = cos(angle)*p[2] + sin(angle)*p[3]  ,   -sin(angle)*p[2] + cos(angle)*p[3]
    end
end
export rotateAroundXAxis!

function rotateAroundXAxisDeg!(unitcell::Unitcell, angle::Float64)
    rotateAroundXAxis!(unitcell, angle*pi/180.0)
end
function rotateAroundXAxisDeg!(lattice::Lattice, angle::Float64)
    rotateAroundXAxis!(lattice, angle*pi/180.0)
end
export rotateAroundXAxisDeg!

# ROTATION AROUND Y AXIS
function rotateAroundYAxis!(unitcell::Unitcell, angle::Float64)
    # rotate the lattice vectors
    for lv in unitcell.lattice_vectors
        # rotate and overwrite simulateously all components
        lv[3], lv[1] = cos(angle)*lv[3] + sin(angle)*lv[1]  ,   -sin(angle)*lv[3] + cos(angle)*lv[1]
    end
    # rotate the basis sites
    for bs in unitcell.basis
        # rotate and overwrite simulateously all components
        bs[3], bs[1] = cos(angle)*bs[3] + sin(angle)*bs[1]  ,   -sin(angle)*bs[3] + cos(angle)*bs[1]
    end
end
function rotateAroundYAxis!(lattice::Lattice, angle::Float64)
    # rotate the lattice vectors
    for lv in lattice.lattice_vectors
        # rotate and overwrite simulateously all components
        lv[3], lv[1] = cos(angle)*lv[3] + sin(angle)*lv[1]  ,   -sin(angle)*lv[3] + cos(angle)*lv[1]
    end
    # rotate the basis sites
    for p in lattice.positions
        # rotate and overwrite simulateously all components
        p[3], p[1] = cos(angle)*p[3] + sin(angle)*p[1]  ,   -sin(angle)*p[3] + cos(angle)*p[1]
    end
end
export rotateAroundYAxis!

function rotateAroundYAxisDeg!(unitcell::Unitcell, angle::Float64)
    rotateAroundYAxis!(unitcell, angle*pi/180.0)
end
function rotateAroundYAxisDeg!(lattice::Lattice, angle::Float64)
    rotateAroundYAxis!(lattice, angle*pi/180.0)
end
export rotateAroundYAxisDeg!

# ROTATION AROUND Z AXIS
function rotateAroundZAxis!(unitcell::Unitcell, angle::Float64)
    # rotate the lattice vectors
    for lv in unitcell.lattice_vectors
        # rotate and overwrite simulateously all components
        lv[1], lv[2] = cos(angle)*lv[1] + sin(angle)*lv[2]  ,   -sin(angle)*lv[1] + cos(angle)*lv[2]
    end
    # rotate the basis sites
    for bs in unitcell.basis
        # rotate and overwrite simulateously all components
        bs[1], bs[2] = cos(angle)*bs[1] + sin(angle)*bs[2]  ,   -sin(angle)*bs[1] + cos(angle)*bs[2]
    end
end
function rotateAroundZAxis!(lattice::Lattice, anglTODOe::Float64)
    # rotate the lattice vectors
    for lv in lattice.lattice_vectors
        # rotate and overwrite simulateously all components
        lv[1], lv[2] = cos(angle)*lv[1] + sin(angle)*lv[2]  ,   -sin(angle)*lv[1] + cos(angle)*lv[2]
    end
    # rotate the basis sites
    for p in lattice.positions
        # rotate and overwrite simulateously all components
        p[1], p[2] = cos(angle)*p[1] + sin(angle)*p[2]  ,   -sin(angle)*p[1] + cos(angle)*p[2]
    end
end
export rotateAroundZAxis!

function rotateAroundZAxisDeg!(unitcell::Unitcell, angle::Float64)
    rotateAroundZAxis!(unitcell, angle*pi/180.0)
end
function rotateAroundZAxisDeg!(lattice::Lattice, angle::Float64)
    rotateAroundZAxis!(lattice, angle*pi/180.0)
end
export rotateAroundZAxisDeg!

###################
#   SCALING
###################

# SCALING ALONG X AXIS
function scaleAlongXAxis!(unitcell::Unitcell, factor::Float64)
    # scale all lattice vectors
    for l in unitcell.lattice_vectors
        l[1] *= factor
    end
    # scale all positions
    for p in unitcell.basis
        p[1] *= factor
    end
end
function scaleAlongXAxis!(lattice::Lattice, factor::Float64)
    # scale all lattice vectors
    for l in lattice.lattice_vectors
        l[1] *= factor
    end
    # scale all positions
    for p in lattice.positions
        p[1] *= factor
    end
end
export scaleAlongXAxis!

# SCALING ALONG Y AXIS
function scaleAlongYAxis!(unitcell::Unitcell, factor::Float64)
    # scale all lattice vectors
    for l in unitcell.lattice_vectors
        l[2] *= factor
    end
    # scale all positions
    for p in unitcell.basis
        p[2] *= factor
    end
end
function scaleAlongYAxis!(lattice::Lattice, factor::Float64)
    # scale all lattice vectors
    for l in lattice.lattice_vectors
        l[2] *= factor
    end
    # scale all positions
    for p in lattice.positions
        p[2] *= factor
    end
end
export scaleAlongYAxis!

# SCALING ALONG Z AXIS
function scaleAlongZAxis!(unitcell::Unitcell, factor::Float64)
    # scale all lattice vectors
    for l in unitcell.lattice_vectors
        l[3] *= factor
    end
    # scale all positions
    for p in unitcell.basis
        p[3] *= factor
    end
end
function scaleAlongZAxis!(lattice::Lattice, factor::Float64)
    # scale all lattice vectors
    for l in lattice.lattice_vectors
        l[3] *= factor
    end
    # scale all positions
    for p in lattice.positions
        p[3] *= factor
    end
end
export scaleAlongZAxis!

# SCALING ISOTROPIC
function scaleIsotropic!(unitcell::Unitcell, factor::Float64)
    # scale all lattice vectors
    for l in unitcell.lattice_vectors
        l .*= factor
    end
    # scale all positions
    for p in unitcell.basis
        p .*= factor
    end
end
function scaleIsotropic!(lattice::Lattice, factor::Float64)
    # scale all lattice vectors
    for l in lattice.lattice_vectors
        l .*= factor
    end
    # scale all positions
    for p in lattice.positions
        p .*= factor
    end
end
export scaleIsotropic!

# SCALING TO MEAN BOND LENGTH
function scaleToMeanBondLength!(unitcell::Unitcell, bond_length::Float64=1.0)
    # get all bond lengths
    bond_lengths = zeros(length(unitcell.connections))
    for i in 1:length(bond_lengths)
        # get the connection
        c = unitcell.connections[i]
        # get the connecting vector
        c_vector = unitcell.basis[Int(c[2])] .- unitcell.basis[Int(c[1])]
        for (j,l) in enumerate(unitcell.lattice_vectors)
            c_vector .+= l.*c[4][j]
        end
        # get how long it is
        bond_lengths[i] = sqrt(sum(c_vector.*c_vector))
    end
    # calculate the mean bond length
    mean_bond_length = mean(bond_lengths)
    # calculate how the scale factor looks like
    factor = bond_length / mean_bond_length
    # scale the object
    scaleIsotropic!(unitcell, factor)
end
function scaleToMeanBondLength!(lattice::Lattice, bond_length::Float64=1.0)
    # get all bond lengths
    bond_lengths = zeros(length(lattice.connections))
    for i in 1:length(bond_lengths)
        # get the connection
        c = lattice.connections[i]
        # get the connecting vector
        c_vector = lattice.positions[Int(c[2])] .- lattice.positions[Int(c[1])]
        for (j,l) in enumerate(lattice.lattice_vectors)
            c_vector .+= l.*c[4][j]
        end
        # get how long it is
        bond_lengths[i] = sqrt(sum(c_vector.*c_vector))
    end
    # calculate the mean bond length
    mean_bond_length = mean(bond_lengths)
    # calculate how the scale factor looks like
    factor = bond_length / mean_bond_length
    # scale the object
    scaleIsotropic!(lattice, factor)
end
export scaleToMeanBondLength!

# SCALING TO MINIMUM BOND LENGTH
function scaleToMinimumBondLength!(unitcell::Unitcell, bond_length::Float64=1.0)
    # get all bond lengths
    bond_lengths = zeros(length(unitcell.connections))
    for i in 1:length(bond_lengths)
        # get the connection
        c = unitcell.connections[i]
        # get the connecting vector
        c_vector = unitcell.basis[Int(c[2])] .- unitcell.basis[Int(c[1])]
        for (j,l) in enumerate(unitcell.lattice_vectors)
            c_vector .+= l.*c[4][j]
        end
        # get how long it is
        bond_lengths[i] = sqrt(sum(c_vector.*c_vector))
    end
    # calculate the minimum bond length
    minimum_bond_length = minimum(bond_lengths)
    # calculate how the scale factor looks like
    factor = bond_length / minimum_bond_length
    # scale the object
    scaleIsotropic!(unitcell, factor)
end
function scaleToMinimumBondLength!(lattice::Lattice, bond_length::Float64=1.0)
    # get all bond lengths
    bond_lengths = zeros(length(lattice.connections))
    for i in 1:length(bond_lengths)
        # get the connection
        c = lattice.connections[i]
        # get the connecting vector
        c_vector = lattice.positions[Int(c[2])] .- lattice.positions[Int(c[1])]
        for (j,l) in enumerate(lattice.lattice_vectors)
            c_vector .+= l.*c[4][j]
        end
        # get how long it is
        bond_lengths[i] = sqrt(sum(c_vector.*c_vector))
    end
    # calculate the minimum bond length
    minimum_bond_length = minimum(bond_lengths)
    # calculate how the scale factor looks like
    factor = bond_length / minimum_bond_length
    # scale the object
    scaleIsotropic!(lattice, factor)
end
export scaleToMinimumBondLength!

###################
#   SHIFTING
###################

# SHIFT ALONG X AXIS
function shiftAlongXAxis!(unitcell::Unitcell, offset::Float64)
    # shift all positions
    for p in unitcell.basis
        p[1] += offset
    end
end
function shiftAlongXAxis!(lattice::Lattice, offset::Float64)
    # shift all positions
    for p in lattice.positions
        p[1] += offset
    end
end
export shiftAlongXAxis!

# SHIFT ALONG Y AXIS
function shiftAlongYAxis!(unitcell::Unitcell, offset::Float64)
    # shift all positions
    for p in unitcell.basis
        p[2] += offset
    end
end
function shiftAlongYAxis!(lattice::Lattice, offset::Float64)
    # shift all positions
    for p in lattice.positions
        p[2] += offset
    end
end
export shiftAlongYAxis!

# SHIFT ALONG Z AXIS
function shiftAlongZAxis!(unitcell::Unitcell, offset::Float64)
    # shift all positions
    for p in unitcell.basis
        p[3] += offset
    end
end
function shiftAlongZAxis!(lattice::Lattice, offset::Float64)
    # shift all positions
    for p in lattice.positions
        p[3] += offset
    end
end
export shiftAlongZAxis!

# SHIFT BY VECTOR
function shiftByVector!(unitcell::Unitcell, vector::Array{Float64,1})
    # shift all positions
    for p in unitcell.basis
        p .+= vector
    end
end
function shiftByVector!(lattice::Lattice, vector::Array{Float64,1})
    # shift all positions
    for p in lattice.positions
        p .+= vector
    end
end
export shiftByVector!

# SHIFT ALONG ANY LATTICE VECTOR
function shiftAlongLatticeVector!(unitcell::Unitcell, lattice_vector::Int64, offset::Float64)
    # define a vector for shifting
    vector = unitcell.lattice_vectors[lattice_vector] .* (offset / sqrt(sum(unitcell.lattice_vectors[lattice_vector].*unitcell.lattice_vectors[lattice_vector])))
    # shift by that vector
    shiftByVector!(unitcell, vector)
end
function shiftAlongLatticeVector!(lattice::Lattice, lattice_vector::Int64, offset::Float64)
    # define a vector for shifting
    vector = lattice.lattice_vectors[lattice_vector] .* (offset / sqrt(sum(lattice.lattice_vectors[lattice_vector].*lattice.lattice_vectors[lattice_vector])))
    # shift by that vector
    shiftByVector!(lattice, vector)
end
export shiftAlongLatticeVector!

# SHIFT ALONG A1
function shiftAlongA1!(unitcell::Unitcell, offset::Float64)
    # shift along the first lattice vector
    shiftAlongLatticeVector!(unitcell, 1, offset)
end
function shiftAlongA1!(lattice::Lattice, offset::Float64)
    # shift along the first lattice vector
    shiftAlongLatticeVector!(lattice, 1, offset)
end
export shiftAlongA1!

# SHIFT ALONG A2
function shiftAlongA2!(unitcell::Unitcell, offset::Float64)
    # shift along the second lattice vector
    shiftAlongLatticeVector!(unitcell, 2, offset)
end
function shiftAlongA2!(lattice::Lattice, offset::Float64)
    # shift along the second lattice vector
    shiftAlongLatticeVector!(lattice, 2, offset)
end
export shiftAlongA2!

# SHIFT ALONG A3
function shiftAlongA3!(unitcell::Unitcell, offset::Float64)
    # shift along the third lattice vector
    shiftAlongLatticeVector!(unitcell, 3, offset)
end
function shiftAlongA3!(lattice::Lattice, offset::Float64)
    # shift along the third lattice vector
    shiftAlongLatticeVector!(lattice, 3, offset)
end
export shiftAlongA3!

# SHIFT ALONG ANY LATTICE VECTOR (Relative)
function shiftAlongLatticeVectorRelative!(unitcell::Unitcell, lattice_vector::Int64, offset::Float64)
    # define a vector for shifting
    vector = unitcell.lattice_vectors[lattice_vector] .* offset
    # shift by that vector
    shiftByVector!(unitcell, vector)
end
function shiftAlongLatticeVectorRelative!(lattice::Lattice, lattice_vector::Int64, offset::Float64)
    # define a vector for shifting
    vector = lattice.lattice_vectors[lattice_vector] .* offset
    # shift by that vector
    shiftByVector!(lattice, vector)
end
export shiftAlongLatticeVectorRelative!

# SHIFT ALONG A1 (Relative)
function shiftAlongA1Relative!(unitcell::Unitcell, offset::Float64)
    # shift along the first lattice vector
    shiftAlongLatticeVectorRelative!(unitcell, 1, offset)
end
function shiftAlongA1Relative!(lattice::Lattice, offset::Float64)
    # shift along the first lattice vector
    shiftAlongLatticeVectorRelative!(lattice, 1, offset)
end
export shiftAlongA1Relative!

# SHIFT ALONG A2 (Relative)
function shiftAlongA2Relative!(unitcell::Unitcell, offset::Float64)
    # shift along the second lattice vector
    shiftAlongLatticeVectorRelative!(unitcell, 2, offset)
end
function shiftAlongA2Relative!(lattice::Lattice, offset::Float64)
    # shift along the second lattice vector
    shiftAlongLatticeVectorRelative!(lattice, 2, offset)
end
export shiftAlongA2Relative!

# SHIFT ALONG A3 (Relative)
function shiftAlongA3Relative!(unitcell::Unitcell, offset::Float64)
    # shift along the third lattice vector
    shiftAlongLatticeVectorRelative!(unitcell, 3, offset)
end
function shiftAlongA3Relative!(lattice::Lattice, offset::Float64)
    # shift along the third lattice vector
    shiftAlongLatticeVectorRelative!(lattice, 3, offset)
end
export shiftAlongA3Relative!

# SHIFT CENTER TO ORIGIN
function shiftCenterToOrigin!(unitcell::Unitcell)
    # look for the central point
    center = unitcell.basis .* 0.0
    for p in unitcell.basis
        center .+= p
    end
    center ./= length(unitcell.basis)
    # shift by negative of that vector
    shiftByVector!(unitcell, center.*(-1))
end
function shiftCenterToOrigin!(lattice::Lattice)
    # look for the central point
    center = lattice.positions .* 0.0
    for p in lattice.positions
        center .+= p
    end
    center ./= length(lattice.positions)
    # shift by negative of that vector
    shiftByVector!(lattice, center.*(-1))
end
export shiftCenterToOrigin!

################################################################################
#
#   5) LATTICE VECTOR OPTIMIZATION
#
################################################################################

# Function to optimize lattice vectors
function optimizeLatticeVectors2D!(unitcell::Unitcell; verbose::Bool=false)
    # init the loop by assuming one refolded a vector
    refolded_one_vector = true
    # as long as one vector was refolded, repeat this loop
    while refolded_one_vector
        # init inside the loop by assuming one did not refold a vector
        refolded_one_vector = false
        # try a1
        if dot(unitcell.lattice_vectors[1],unitcell.lattice_vectors[1]) > dot(unitcell.lattice_vectors[1].-unitcell.lattice_vectors[2], unitcell.lattice_vectors[1].-unitcell.lattice_vectors[2])
            # print
            if verbose
                println("a1 --> a1 - a2")
            end
            # change the lattice vector
            unitcell.lattice_vectors[1] = unitcell.lattice_vectors[1] .- unitcell.lattice_vectors[2]
            # change all connections
            for c in unitcell.connections
                if c[4][1] != 0
                    c[4] = (c[4][1], c[4][2]+c[4][1])
                end
            end
            # did refold
            refolded_one_vector = true
        end
        if dot(unitcell.lattice_vectors[1],unitcell.lattice_vectors[1]) > dot(unitcell.lattice_vectors[1].+unitcell.lattice_vectors[2], unitcell.lattice_vectors[1].+unitcell.lattice_vectors[2])
            # print
            if verbose
                println("a1 --> a1 + a2")
            end
            # change the lattice vector
            unitcell.lattice_vectors[1] = unitcell.lattice_vectors[1] .+ unitcell.lattice_vectors[2]
            # change all connections
            for c in unitcell.connections
                if c[4][1] != 0
                    c[4] = (c[4][1], c[4][2]-c[4][1])
                end
            end
            # did refold
            refolded_one_vector = true
        end
        # try a2
        if dot(unitcell.lattice_vectors[2],unitcell.lattice_vectors[2]) > dot(unitcell.lattice_vectors[2].-unitcell.lattice_vectors[1], unitcell.lattice_vectors[2].-unitcell.lattice_vectors[1])
            # print
            if verbose
                println("a2 --> a2 - a1")
            end
            # change the lattice vector
            unitcell.lattice_vectors[2] = unitcell.lattice_vectors[2] .- unitcell.lattice_vectors[1]
            # change all connections
            for c in unitcell.connections
                if c[4][2] != 0
                    c[4] = (c[4][1]+c[4][2], c[4][2])
                end
            end
            # did refold
            refolded_one_vector = true
        end
        if dot(unitcell.lattice_vectors[2],unitcell.lattice_vectors[2]) > dot(unitcell.lattice_vectors[2].+unitcell.lattice_vectors[1], unitcell.lattice_vectors[2].+unitcell.lattice_vectors[1])
            # print
            if verbose
                println("a2 --> a2 + a1")
            end
            # change the lattice vector
            unitcell.lattice_vectors[2] = unitcell.lattice_vectors[2] .+ unitcell.lattice_vectors[1]
            # change all connections
            for c in unitcell.connections
                if c[4][2] != 0
                    c[4] = (c[4][1]-c[4][2], c[4][2])
                end
            end
            # did refold
            refolded_one_vector = true
        end
    end
end

# Function to optimize lattice vectors
function optimizeLatticeVectors3D!(unitcell::Unitcell; verbose::Bool=false)
    # init the loop by assuming one refolded a vector
    refolded_one_vector = true
    # as long as one vector was refolded, repeat this loop
    while refolded_one_vector
        # init inside the loop by assuming one did not refold a vector
        refolded_one_vector = false
        # try all combinations of the remaining two lattice vectors
        for comb in Array{Int64,1}[[1,0], [-1,0], [0,1], [0,-1], [1,1], [1,-1], [-1,1], [-1,-1]]
            # try a1
            a1p = unitcell.lattice_vectors[1] .+ comb[1].*unitcell.lattice_vectors[2] .+ comb[2].*unitcell.lattice_vectors[3]
            # check if this is closer
            if dot(unitcell.lattice_vectors[1],unitcell.lattice_vectors[1]) > dot(a1p, a1p)
                # print
                if verbose
                    println("a1 --> a1 + $(comb[1])*a2 + $(comb[2])*a3")
                end
                # change the lattice vector
                unitcell.lattice_vectors[1] = a1p
                # change all connections
                for c in unitcell.connections
                    c[4] = (c[4][1], c[4][2]-c[4][1]*comb[1], c[4][3]-c[4][1]*comb[2])
                end
                # did refold
                refolded_one_vector = true
            end
            # try a2
            a2p = unitcell.lattice_vectors[2] .+ comb[1].*unitcell.lattice_vectors[1] .+ comb[2].*unitcell.lattice_vectors[3]
            # check if this is closer
            if dot(unitcell.lattice_vectors[2],unitcell.lattice_vectors[2]) > dot(a2p, a2p)
                # print
                if verbose
                    println("a2 --> a2 + $(comb[1])*a1 + $(comb[2])*a3")
                end
                # change the lattice vector
                unitcell.lattice_vectors[2] = a2p
                # change all connections
                for c in unitcell.connections
                    c[4] = (c[4][1]-c[4][2]*comb[1], c[4][2], c[4][3]-c[4][2]*comb[2])
                end
                # did refold
                refolded_one_vector = true
            end
            # try a3
            a3p = unitcell.lattice_vectors[3] .+ comb[1].*unitcell.lattice_vectors[1] .+ comb[2].*unitcell.lattice_vectors[2]
            # check if this is closer
            if dot(unitcell.lattice_vectors[3],unitcell.lattice_vectors[3]) > dot(a3p, a3p)
                # print
                if verbose
                    println("a3 --> a3 + $(comb[1])*a1 + $(comb[2])*a2")
                end
                # change the lattice vector
                unitcell.lattice_vectors[3] = a3p
                # change all connections
                for c in unitcell.connections
                    c[4] = (c[4][1]-c[4][3]*comb[1], c[4][2]-c[4][3]*comb[2], c[4][3])
                end
                # did refold
                refolded_one_vector = true
            end
        end
    end
end

# general function
function optimizeLatticeVectors!(unitcell::Unitcell; verbose::Bool=false)
    if length(unitcell.lattice_vectors) == 2
        return optimizeLatticeVectors2D!(unitcell, verbose=verbose)
    elseif length(unitcell.lattice_vectors) == 3
        return optimizeLatticeVectors3D!(unitcell, verbose=verbose)
    else
        # DO NOTHING
    end
end
export optimizeLatticeVectors!


text:
# WARNING:
# ========
# This blog post https://julialang.org/blog/2019/01/fluxdiffeq/ gives a high level
# description of this file.

# Many function calls are fully qualified. This is not necessary. The package name is
# included to get a better sense of the role of the different packages.

# The Lotka-Volterra differential equations model populations of predators
# and preys over time.

# This example starts with an ODE setup with a given set of parameters.
# Then, using the same equations, optimises starting with a different set of parameters

## Import the required packages
using DifferentialEquations, DiffEqSensitivity
using Flux, DiffEqFlux
using Plots

##
## Setup ODE to optimize
##

# The populations are representated as a vector of 2 numbers:
#    U = [numbers of rabbits, number of wolves]
#
# ODE problems expect the equations to be presented with standard parameters:
#   dU: the derivatives of U. This vector is modified in place for speed purposes
#    U: the functions
#    p: a vector of paramters
#    t: the functions' variable (here time)
function lotka_volterra!(du, u, p, t)
  x, y = u
  α, β, δ, γ = p
  du[1] = dx = α * x - β * x * y
  du[2] = dy = -δ * y + γ * x * y
end

# Initial conditions
u0 = [1.0, 1.0]

# Time over which to simulate
t_span = (0.0, 10.0)

# Initial parameters
p0 = [1.5, 1.0, 3.0, 1.0]

# Setup the ODE
prob0 = DiffEqBase.ODEProblem(lotka_volterra!, u0, t_span, p0)

# Solve the ODE
sol0 = DiffEqBase.solve(prob0, Tsit5())

# Plot it
plot(sol0)

##
## Generate data that will be used to train the network
## The data is the number of rabbits/wolves at set times
##

# Solve the ODE (again) and collect solutions at fixed intervals
target_data = DiffEqBase.solve(prob0, Tsit5(), saveat = 0.1)
rabbits = target_data[1, :]                     # vector of 101 data points
wolves  = target_data[2, :]                    # vector of 101 data points

# Plot the data on top of the full solution
t_steps = 0:0.1:10.0
scatter!(t_steps, rabbits, color = [1], label = "rabbits")
scatter!(t_steps, wolves,  color = [2], label = "wolves")

##
## Parameter optimisation
##

# Loss function is the total squared error
# (the number of points is constant - mean is not necessary)
loss_function = function()
    prediction = DiffEqBase.concrete_solve(prob0, Tsit5(), u0, p;
    sensealg = TrackerAdjoint(), saveat = 0.0:0.1:10.0)

    # Calculate squared error
  return sum(abs2, prediction - target_data)
end

# Callback function to observe training
iter = 0
callback = function ()
  global iter += 1
  if iter % 10 == 1
    # Plot the training data
    scatter(t_steps,  rabbits, color = [1], label = "rabbit data")
    scatter!(t_steps, wolves,  color = [2], label = "wolves data")

    # Use `remake` to re-create the original `prob0` with different current parameters `p`
    remade_sol = DiffEqBase.solve(remake(prob0, p = p), Tsit5(), saveat = 0.1)
    remade_rabbits = remade_sol[1, :]
    remade_wolves  = remade_sol[2, :]
    plot!(t_steps, remade_rabbits, ylim = (0, 6), labels = "rabbit model", color = 1)
    display(plot!(t_steps, remade_wolves, ylim = (0, 6), labels = "wolf model", color = 2))
  end
end

# Vector of new parameters different from p0
p = [4.0, 1.0, 2.0, 0.4]
ODEparams = Flux.params(p)

# Display the ODE with the initial parameter values.
callback()

# The parameter `loss_data` in `train!` (below) contains a list of data that are provided
# to the loss function. `train!`` works by iterating over that list. In our case, the
# calculation of the loss function does not require any additional data.
# Therefore `loss_data` is filled with a whole bunch of nothing.
loss_data = Iterators.repeated((), 1000)

# ADAM is a great default optimiser
optimiser = Flux.ADAM(0.1)

# train! is a function that hides some of the complexity of the library
Flux.train!(loss_function, ODEparams, loss_data, optimiser; cb = callback)