text:
#-----------------------------------------------------------------------------------------------------------------------------
#-----------------------------------------------------------------------------------------------------------------------------
#
#   LATTICE OPERATIONS (MODIFY LATTICE CONNECTIONS AND POSITIONS)
#
#-----------------------------------------------------------------------------------------------------------------------------
#-----------------------------------------------------------------------------------------------------------------------------

#-----------------------------------------------------------------------------------------------------------------------------
#
#   Build a lattice-X / unitcell-X version of the lattice / unitcell
#   Procedure replaces all bonds by sites
#
#-----------------------------------------------------------------------------------------------------------------------------
function getTransformedLatticeBondToSite(lattice::Lattice; connection_strength="AUTO")
    # new positions and connections
    positions = copy(lattice.positions)
    connections = Array{Any,1}[]
    connectionsTreated = Array[]
    # generate a new filename
    filename_new = replace(lattice.filename, ".jld", "_X.jld")
    # neutral wrap
    neutral_wrap = (0,0,0)
    if size(lattice.lattice_vectors,1) == 2
   	    neutral_wrap = (0,0)
    elseif size(lattice.lattice_vectors,1) == 1
   	    neutral_wrap = (0)
    end
    # iterate over all connections
    for c in lattice.connections
        # only treat one of the two
        treated = false
        for ct in connectionsTreated
            if ct[1] == c[2] && ct[2] == c[1] && ct[3] == c[3]
                treated = true
                break
            end
        end
        if treated
            continue
        end
        # add a new position
        pnew = positions[Int(c[1])].+positions[Int(c[2])]
        for i in 1:length(c[4])
            pnew .+= c[4][i].*lattice.lattice_vectors[i]
        end
        pnew = pnew .* 0.5
        push!(positions, pnew)
        # add new connections
        if typeof(connection_strength) == String && connection_strength == "AUTO"
            con_strength = "con$(size(connectionsTreated,1))"
        elseif typeof(connection_strength) == String && connection_strength == "SAME"
            con_strength = c[3]
        elseif typeof(connection_strength) == String && connection_strength == "SQRT"
            if typeof(c[3]) == String
                con_strength = "sqrt($(c[3]))"
            else
                con_strength = sqrt(c[3])
            end
        else
            con_strength = connection_strength
        end
        push!(connections, [c[1], size(positions,1), con_strength, neutral_wrap])
        push!(connections, [size(positions,1), c[1], con_strength, neutral_wrap])
        push!(connections, [size(positions,1), c[2], con_strength, c[4]])
        if length(c[4]) == 3
            push!(connections, [c[2], size(positions,1), con_strength, (-c[4][1],-c[4][2],-c[4][3])])
        elseif length(c[4]) == 2
            push!(connections, [c[2], size(positions,1), con_strength, (-c[4][1],-c[4][2])])
        else
            push!(connections, [c[2], size(positions,1), con_strength, (-c[4])])
        end
        # add to treated
        push!(connectionsTreated, c)
    end
    # build a new lattice
    lattice_X = Lattice(lattice.lattice_vectors, positions, connections, filename_new)
    # return the new lattice
    return lattice_X
end
export getTransformedLatticeBondToSite

function getTransformedUnitcellBondToSite(unitcell::Unitcell; connection_strength="AUTO")
    # new positions and connections
    positions = copy(unitcell.basis)
    connections = Array{Any,1}[]
    connectionsTreated = Array[]
    # generate a new filename
    filename_new = replace(unitcell.filename, ".jld", "_X.jld")
    # neutral wrap
    neutral_wrap = (0,0,0)
    if size(unitcell.lattice_vectors,1) == 2
   	    neutral_wrap = (0,0)
    elseif size(unitcell.lattice_vectors,1) == 1
   	    neutral_wrap = (0)
    end
    # iterate over all connections
    for c in unitcell.connections
        # only treat one of the two
        treated = false
        for ct in connectionsTreated
            if ct[1] == c[2] && ct[2] == c[1] && ct[3] == c[3]
                treated = true
                break
            end
        end
        if treated
            continue
        end
        # add a new position
        pnew = positions[Int(c[1])].+positions[Int(c[2])]
        for i in 1:length(c[4])
            pnew .+= c[4][i].*unitcell.lattice_vectors[i]
        end
        pnew = pnew .* 0.5
        push!(positions, pnew)
        # add new connections
        if typeof(connection_strength) == String && connection_strength == "AUTO"
            con_strength = "con$(size(connectionsTreated,1))"
        elseif typeof(connection_strength) == String && connection_strength == "SAME"
            con_strength = c[3]
        elseif typeof(connection_strength) == String && connection_strength == "SQRT"
            if typeof(c[3]) == String
                con_strength = "sqrt($(c[3]))"
            else
                con_strength = sqrt(c[3])
            end
        else
            con_strength = connection_strength
        end
        push!(connections, [c[1], size(positions,1), con_strength, neutral_wrap])
        push!(connections, [size(positions,1), c[1], con_strength, neutral_wrap])
        push!(connections, [size(positions,1), c[2], con_strength, c[4]])
        if length(c[4]) == 3
            push!(connections, [c[2], size(positions,1), con_strength, (-c[4][1],-c[4][2],-c[4][3])])
        elseif length(c[4]) == 2
            push!(connections, [c[2], size(positions,1), con_strength, (-c[4][1],-c[4][2])])
        else
            push!(connections, [c[2], size(positions,1), con_strength, (-c[4])])
        end
        # add to treated
        push!(connectionsTreated, c)
    end
    # build a new unitcell
    uc_X = Unitcell(unitcell.lattice_vectors, positions, connections, filename_new)
    # return the new unitcell
    return uc_X
end
export getTransformedUnitcellBondToSite

#-----------------------------------------------------------------------------------------------------------------------------
#
#   Build a the squareroot version of the unitcell
#   Procedure replaces all fully connected plaquettes by sites and then sets the interaction strengths (TODO)
#
#-----------------------------------------------------------------------------------------------------------------------------
function getTransformedUnitcellFCPToSite(unitcell::Unitcell; printFCP=false, connection_strength="AUTO")
    # look for all fully connected plaquettes
    fcp_list = []
    connectionlist = getConnectionList(unitcell)
    neutral_wrap = (0,0,0)
    if length(connectionlist[1][1][4]) == 2
        neutral_wrap = (0,0)
    elseif length(connectionlist[1][1][4]) == 1
        neutral_wrap = (0)
    end
    # define a recursive looking method
    function getFCPs(fcp_current)
        # a site with index current_index is asked if it is in a fully connected plaquette with more sits than given by fcp
        # returns the fcps that the site is in (minimum: the fcp given)
        fcp_list_new = Array[]
        # iterate over all neighbors of the current site
        for nc in connectionlist[Int(fcp_current[end][1])]
            # current test site and test wrap
            test_site = Int(nc[2])
            test_wrap = collect(nc[4]) .+ fcp_current[end][2]
            # check, if "step back"
            if test_site == fcp_current[end-1][1] && test_wrap == fcp_current[end-1][2]
                continue
            end
            # no step back, check if this site is connected to all fcp sites that are already in the fcp
            completely_connected = true
            for site_in_fcp in fcp_current
                # check if new site is connected to this site
                connected = false
                for c in connectionlist[test_site]
                    # check if the destination and wrap agree
                    if Int(c[2]) == Int(site_in_fcp[1]) && collect(c[4]) .+ test_wrap == site_in_fcp[2]
                        # found a match
                        connected = true
                        break
                    end
                end
                # if not connected: complete connection broken
                if !connected
                    completely_connected = false
                    break
                end
            end
            # if it is connected to everything: make a new fcp out of it and go deper into the recursion
            if completely_connected
                fcp_new = []
                for fcp_current_element in fcp_current
                    push!(fcp_new, fcp_current_element)
                end
                push!(fcp_new, [test_site, test_wrap])
                # recursion
                fcp_list_recursion = getFCPs(fcp_new)
                # add those to the list
                for fcp in fcp_list_recursion
                    push!(fcp_list_new, fcp)
                end
            end
        end
        # if empty, insert the current fcp and return
        if length(fcp_list_new) == 0
            push!(fcp_list_new, fcp_current)
        end
        # return fcp list
        return fcp_list_new
    end
    # define a function to give connections for a certain fcp
    function getFCPConnections(fcp)
        # list
        fcp_connections = []
        # find all connections
        for site1 in fcp
        for site2 in fcp
            # add connection from 1 to 2 if found
            for c in connectionlist[Int(site1[1])]
                if c[1] == Int(site1[1]) && c[2] == Int(site2[1]) && collect(c[4]) == site2[2] .- site1[2]
                    push!(fcp_connections, c)
                end
            end
        end
        end
        # return the list of connections
        return fcp_connections
    end
    # define a function to compare to lists
    function sameConnections(list1, list2)
        # check if every element in list 1 is in list 2
        for element1 in list1
            # not found yet
            found = false
            # check all elements of list2
            for element2 in list2
                # check if equal
                if Int(element1[1]) == Int(element2[1]) && Int(element1[2]) == Int(element2[2]) && element1[4] == element2[4]
                    found = true
                    break
                end
            end
            # if not found, return false
            if !found
                return false
            end
        end
        # check if every element in list 2 is in list 1
        for element2 in list2
            # not found yet
            found = false
            # check all elements of list1
            for element1 in list1
                # check if equal
                if Int(element1[1]) == Int(element2[1]) && Int(element1[2]) == Int(element2[2]) && element1[4] == element2[4]
                    found = true
                    break
                end
            end
            # if not found, return false
            if !found
                return false
            end
        end
        # every test was successfull: 2 lists identical
        return true
    end
    # define a function to get the center of the fcp
    function getFCPCenter(fcp)
        # add all sites
        center = unitcell.basis[Int(fcp[1][1])]
        for i in 2:size(fcp,1)
            center = center .+ unitcell.basis[Int(fcp[i][1])]
        end
        # add all wraps
        for i in 1:size(fcp, 1)
            for l in 1:size(unitcell.lattice_vectors,1)
                center = center .+ (fcp[i][2][l] .* unitcell.lattice_vectors[l])
            end
        end
        # devide by number of sites
        center = center ./ size(fcp,1)
    end
    # iterate for all sites
    for s in 1:size(unitcell.basis,1)
        # looking for fully connected plaquettes of site s
        # check if neighbors are part of bigger plaquettes
        for neighbor_connection in connectionlist[s]
            # ask if the neighbor is part of an FCP
            fcp_new = [
                [s, collect(neutral_wrap)],
                [Int(neighbor_connection[2]), collect(neighbor_connection[4])]
            ]
            fcp_list_of_neighbor_and_site = getFCPs(fcp_new)
            # go through all fcps that have been identified and add new ones if necessary
            for fcp in fcp_list_of_neighbor_and_site
                # just add to global list, filtering later on
                push!(fcp_list, fcp)
            end
        end
    end
    # up to here: found all fcps
    fcp_list_raw = fcp_list
    fcp_list = []
    fcp_list_connections = []
    # now: filter the identical fcps out
    for fcp_test in fcp_list_raw
        # boolean if found in the ready to use list
        found = false
        # get the connections
        fcp_test_connections = getFCPConnections(fcp_test)
        # check all elements of the ready to use list
        for (i,fcp) in enumerate(fcp_list)
            # first test: dimensions identical?
            if size(fcp, 1) != size(fcp_test, 1)
                continue
            end
            # second test: get the connections associated with this certain fcp
            fcp_connections = fcp_list_connections[i]
            # test if the same lists, i.e. if both are in the other
            if sameConnections(fcp_connections, fcp_test_connections)
                found = true
                break
            end
        end
        # if not found, add to the list
        if !found
            push!(fcp_list, fcp_test)
            push!(fcp_list_connections, fcp_test_connections)
        end
    end
    # maybe print the FCP list
    if printFCP
        println("In total: $(size(fcp_list,1)) FCPs found")
        for i in 1:size(fcp_list,1)
            println("FCP $(i) with $(size(fcp_list[i], 1)) sites:")
            println(" - sites: $(fcp_list[i])")
            println(" - connections: $(fcp_list_connections[i])")
        end
    end
    # replace all fcps by inserting new sites
    positions_new = []
    for p in unitcell.basis
        push!(positions_new, p)
    end
    # push new sites
    for fcp in fcp_list
        # get the center of the fcp and add it as a new position
        push!(positions_new, getFCPCenter(fcp))
    end
    # define new connections
    connections_new = []
    # iterate over all sites / all fcps
    for (i, fcp) in enumerate(fcp_list)
        # new position index
        pos_index = size(unitcell.basis, 1) + i
        # define new connections for every site in the FCP (new position wrap is neutral wrap)
        for (s,site) in enumerate(fcp)
            # build the fields of the new 2 connections to insert
            con_field_1 = Int(site[1])
            con_field_2 = pos_index
            if typeof(connection_strength) == String && connection_strength == "AUTO"
                con_field_3 = "fcp$(i)con$(s)"
            else
                con_field_3 = connection_strength
            end
            con_field_4_array = collect(neutral_wrap) .- site[2]
            if length(con_field_4_array) == 3
                con_field_4_1 = (con_field_4_array[1], con_field_4_array[2], con_field_4_array[3])
                con_field_4_2 = (-con_field_4_array[1], -con_field_4_array[2], -con_field_4_array[3])
            elseif length(con_field_4_array) == 2
                con_field_4_1 = (con_field_4_array[1], con_field_4_array[2])
                con_field_4_2 = (-con_field_4_array[1], -con_field_4_array[2])
            elseif length(con_field_4_array) == 1
                con_field_4_1 = (con_field_4_array[1])
                con_field_4_2 = (-con_field_4_array[1])
            else
                println("wrong length for field 4: $(con_field_4_array)")
            end
            # build new connections
            con_1 = [con_field_1; con_field_2; con_field_3; con_field_4_1]
            con_2 = [con_field_2; con_field_1; con_field_3; con_field_4_2]
            push!(connections_new, con_1)
            push!(connections_new, con_2)
        end
    end
    # build new unitcell
    unitcell_fcp_replaced = Unitcell(
        unitcell.lattice_vectors,
        positions_new,
        connections_new,
        replace(unitcell.filename, ".jld", "_fcp_to_site.jld"))
    # return the new unitcell
    return unitcell_fcp_replaced
end
export getTransformedUnitcellFCPToSite

#function getTransformedLatticeFCPToSite(lattice::Lattice)
#
#end
#export getTransformedLatticeFCPToSite

#-----------------------------------------------------------------------------------------------------------------------------
#
#   Square the lattice / unitcell
#   This maps all connections to NNN connections by multiplication
#
#-----------------------------------------------------------------------------------------------------------------------------
function getSquaredLattice(lattice::Lattice)
    # get the connectivity matrix of the lattice
    connectivity = getConnectionList(lattice)
    # define a list of new connections
    connections_new = Array[]
    # iterate over all sites and check if they host a NNN connection
    for i in 1:size(lattice.positions,1)
        # check the neighbors of site i
        for (i1,c1) in enumerate(connectivity[i])
        for (i2,c2) in enumerate(connectivity[i])
            # propose a connection c1+c2 if they are different
            if i1 < i2
                # build the offsets
                if length(c1[4]) == 1
                    off_1 = (c2[4] - c1[4])
                    off_2 = (c1[4] - c2[4])
                elseif length(c1[4]) == 2
                    off_1 = (c2[4][1] - c1[4][1], c2[4][2] - c1[4][2])
                    off_2 = (c1[4][1] - c2[4][1], c1[4][2] - c2[4][2])
                elseif length(c1[4]) == 3
                    off_1 = (c2[4][1] - c1[4][1], c2[4][2] - c1[4][2], c2[4][3] - c1[4][3])
                    off_2 = (c1[4][1] - c2[4][1], c1[4][2] - c2[4][2], c1[4][3] - c2[4][3])
                end
                # build two new connections
                if typeof(c1[3]) == String || typeof(c2[3]) == String
                    connection_new_1 = [Int(c1[2]); Int(c2[2]); "$(c1[3])*$(c2[3])"; off_1]
                    connection_new_2 = [Int(c2[2]); Int(c1[2]); "$(c1[3])*$(c2[3])"; off_2]
                else
                    connection_new_1 = [Int(c1[2]); Int(c2[2]); (c1[3])*(c2[3]); off_1]
                    connection_new_2 = [Int(c2[2]); Int(c1[2]); (c1[3])*(c2[3]); off_2]
                end
                # push them to the add list
                push!(connections_new, connection_new_1)
                push!(connections_new, connection_new_2)
            end
        end
        end
    end
    # create new Lattice with only new connections (and leave out old ones)
    return getLatticeWithOptimizedConnections(
        Lattice(
            lattice.unitcell,
            lattice.unitcellRepetitions,
            lattice.lattice_vectors,
            lattice.positions,
            lattice.positions_indices,
            connections_new,
            replace(lattice.filename, ".jld", "_squared.jld"))
    )
end
export getSquaredLattice

function getSquaredUnitcell(unitcell::Unitcell)
    # get the connectivity matrix of the unitcell
    connectivity = getConnectionList(unitcell)
    # define a list of new connections
    connections_new = Array[]
    # iterate over all sites and check if they host a NNN connection
    for i in 1:size(unitcell.basis,1)
        # check the neighbors of site i
        for (i1,c1) in enumerate(connectivity[i])
        for (i2,c2) in enumerate(connectivity[i])
            # propose a connection c1+c2 if the connections c1 and c2 are different
            if i1 < i2
                # build the offsets
                if length(c1[4]) == 1
                    off_1 = (c2[4] - c1[4])
                    off_2 = (c1[4] - c2[4])
                elseif length(c1[4]) == 2
                    off_1 = (c2[4][1] - c1[4][1], c2[4][2] - c1[4][2])
                    off_2 = (c1[4][1] - c2[4][1], c1[4][2] - c2[4][2])
                elseif length(c1[4]) == 3
                    off_1 = (c2[4][1] - c1[4][1], c2[4][2] - c1[4][2], c2[4][3] - c1[4][3])
                    off_2 = (c1[4][1] - c2[4][1], c1[4][2] - c2[4][2], c1[4][3] - c2[4][3])
                end
                # build two new connections
                if typeof(c1[3]) == String || typeof(c2[3]) == String
                    connection_new_1 = [Int(c1[2]); Int(c2[2]); "$(c1[3])*$(c2[3])"; off_1]
                    connection_new_2 = [Int(c2[2]); Int(c1[2]); "$(c1[3])*$(c2[3])"; off_2]
                else
                    connection_new_1 = [Int(c1[2]); Int(c2[2]); (c1[3])*(c2[3]); off_1]
                    connection_new_2 = [Int(c2[2]); Int(c1[2]); (c1[3])*(c2[3]); off_2]
                end
                # push them to the add list
                push!(connections_new, connection_new_1)
                push!(connections_new, connection_new_2)
            end
        end
        end
    end
    # create new Lattice with only new connections (and leave out old ones)
    return getUnitcellWithOptimizedConnections(
        Unitcell(
            unitcell.lattice_vectors,
            unitcell.basis,
            connections_new,
            replace(unitcell.filename, ".jld", "_squared.jld"))
    )
end
export getSquaredUnitcell

#-----------------------------------------------------------------------------------------------------------------------------
#
#   Get a list of all independent sublattices / sublattices of a lattice / unitcell (i.e. not connected parts of a lattice)
#   Not connected parts can occur when "squaring" a bipartite lattice
#
#-----------------------------------------------------------------------------------------------------------------------------
function getIndependentSublattices(lattice::Lattice)
    # save relevant lattice data
    connections = lattice.connections
    positions = lattice.positions
    positions_indices = lattice.positions_indices
    # generate a list for the sublattices
    sublattices = []
    # define an array of labels
    labels = zeros(Int64, size(lattice.positions,1))
    currentlabel = 1
    # get the connectivity list of the lattice
    connectivity = getConnectionList(lattice)
    # iterate over all positions
    for i in 1:size(lattice.positions, 1)
        # get the label of site i
        if labels[i] == 0
            labels[i] = currentlabel
            currentlabel = currentlabel + 1
        end
        # iterate over all connections
        for c in connectivity[i]
            # get the connected site
            j = Int(c[2])
            # check the label of j
            if labels[j] == 0
                # set as connected to i
                labels[j] = labels[i]
            elseif labels[j] == labels[i]
                # just ignore
            else
                # merging of two clusters here!
                newlabel = min(labels[i], labels[j])
                # scan through the whole lattice
                for n in 1:length(labels)
                    if labels[n] == labels[i] || labels[n] == labels[j]
                        labels[n] = newlabel
                    end
                end
                # get a new currentlabel
                currentlabel = Int(maximum(labels)) + 1
            end
        end
    end
    # get a list of labels
    labellist = []
    for l in labels
        if !(l in labellist)
            push!(labellist, l)
        end
    end
    # print how many sublattices found
    println("$(length(labellist)) independent sublattice(s) found")
    # check how many sublattices
    if length(labellist) == 1
        # just push the original lattice into the sublattices list
        sublattice = Lattice(
            lattice.unitcell,
            lattice.unitcellRepetitions,
            lattice.lattice_vectors,
            lattice.positions,
            lattice.positions_indices,
            lattice.connections,
            lattice.filename_output
        )
        push!(sublattices, sublattice)
        # return
        return sublattices
    end
    # iterate over labels
    for l in labellist
        # find the new positions
        positions_new = Array[]
        positions_indices_new = []
        connections_new = Array[]
        mapping_indices_old = []
        # iterate over all connections to identify relevant positions
        for c in connections
            # check if belonging to the sublattice
            if labels[Int(c[1])] != l
                # not of this sublattice
                continue
            end
            # get the starting and finishing position
            index_from = Int(c[1])
            index_to = Int(c[2])
            # check where they map to
            if index_from in mapping_indices_old
                index_from_new = findfirst(mapping_indices_old, index_from)
            else
                # add the position to the list of positions
                push!(positions_new, positions[index_from])
                push!(positions_indices_new, positions_indices[index_from])
                # add the index mapping to the list
                push!(mapping_indices_old, index_from)
                # get the new index
                index_from_new = size(positions_new,1)
            end
            if index_to in mapping_indices_old
                index_to_new = findfirst(mapping_indices_old, index_to)
            else
                # add the position to the list of positions
                push!(positions_new, positions[index_to])
                push!(positions_indices_new, positions_indices[index_to])
                # add the index mapping to the list
                push!(mapping_indices_old, index_to)
                # get the new index
                index_to_new = size(positions_new,1)
            end
            # push the new connection
            # format [<from index>; <to index>; <strength>; (<lattice displaced by lattice vector j>)]
            connection_new = [index_from_new; index_to_new; c[3]; c[4]]
			# register as connection
			push!(connections_new, connection_new)
        end
        # generate a new sublattice
        sublattice = Lattice(
            lattice.unitcell,
            lattice.unitcellRepetitions,
            lattice.lattice_vectors,
            positions_new,
            positions_indices_new,
            connections_new,
            "$(lattice.filename[1:end-4])_sublattice_$(l).jld"
        )
        # push it the list
        push!(sublattices, sublattice)
    end
    # return the list
    return sublattices
end
export getIndependentSublattices

function getIndependentSubunitcells(unitcell::Unitcell)
    # save relevant lattice data
    connections = unitcell.connections
    positions = unitcell.basis
    # generate a list for the sublattices
    subunitcells = []
    # define an array of labels
    labels = zeros(Int64, size(positions,1))
    currentlabel = 1
    # get the connectivity list of the unitcell
    connectivity = getConnectionList(unitcell)
    # iterate over all positions
    for i in 1:size(positions, 1)
        # get the label of site i
        if labels[i] == 0
            labels[i] = currentlabel
            currentlabel = currentlabel + 1
        end
        # iterate over all connections
        for c in connectivity[i]
            # get the connected site
            j = Int(c[2])
            # check the label of j
            if labels[j] == 0
                # set as connected to i
                labels[j] = labels[i]
            elseif labels[j] == labels[i]
                # just ignore
            else
                # merging of two clusters here!
                newlabel = min(labels[i], labels[j])
                # scan through the whole lattice
                for n in 1:length(labels)
                    if labels[n] == labels[i] || labels[n] == labels[j]
                        labels[n] = newlabel
                    end
                end
                # get a new currentlabel
                currentlabel = Int(maximum(labels)) + 1
            end
        end
    end
    # get a list of labels
    labellist = []
    for l in labels
        if !(l in labellist)
            push!(labellist, l)
        end
    end
    # print how many sublattices found
    println("$(length(labellist)) independent sublattice(s) found")
    # check how many sublattices
    if length(labellist) == 1
        # just push the original lattice into the sublattices list
        subunitcell = Unitcell(
            unitcell.lattice_vectors,
            unitcell.basis,
            unitcell.connections,
            unitcell.filename
        )
        push!(subunitcells, subunitcell)
        # return
        return subunitcells
    end
    # iterate over labels
    for l in labellist
        # find the new positions
        positions_new = Array[]
        connections_new = Array[]
        mapping_indices_old = []
        # iterate over all connections to identify relevant positions
        for c in connections
            # check if belonging to the sublattice
            if labels[Int(c[1])] != l
                # not of this sublattice
                continue
            end
            # get the starting and finishing position
            index_from = Int(c[1])
            index_to = Int(c[2])
            # check where they map to
            if index_from in mapping_indices_old
                index_from_new = findfirst(mapping_indices_old, index_from)
            else
                # add the position to the list of positions
                push!(positions_new, positions[index_from])
                # add the index mapping to the list
                push!(mapping_indices_old, index_from)
                # get the new index
                index_from_new = size(positions_new,1)
            end
            if index_to in mapping_indices_old
                index_to_new = findfirst(mapping_indices_old, index_to)
            else
                # add the position to the list of positions
                push!(positions_new, positions[index_to])
                # add the index mapping to the list
                push!(mapping_indices_old, index_to)
                # get the new index
                index_to_new = size(positions_new,1)
            end
            # push the new connection
            # format [<from index>; <to index>; <strength>; (<lattice displaced by lattice vector j>)]
            connection_new = [index_from_new; index_to_new; c[3]; c[4]]
			# register as connection
			push!(connections_new, connection_new)
        end
        # generate a new sublattice
        subunitcell = Unitcell(
            unitcell.lattice_vectors,
            positions_new,
            connections_new,
            replace(unitcell.filename, ".jld", "_subunitcell_$(l).jld")
        )
        # push it the list
        push!(subunitcells, subunitcell)
    end
    # return the list
    return subunitcells
end
export getIndependentSubunitcells

#-----------------------------------------------------------------------------------------------------------------------------
#
#   Optimize the lattice / unitcell by collapsing redundant connections into fewer connections and return a new object
#
#-----------------------------------------------------------------------------------------------------------------------------
function getLatticeWithOptimizedConnections(lattice::Lattice)
    # build up new connections
    connections_new = Array[]
    # go through all old connections
    for c in lattice.connections
        # check if the connection already is present in the 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[]
    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
    # return a new lattice
    return Lattice(
        lattice.unitcell,
        lattice.unitcellRepetitions,
        lattice.lattice_vectors,
        lattice.positions,
        lattice.positions_indices,
        connections_new_2,
        lattice.filename)
end
export getLatticeWithOptimizedConnections

function getUnitcellWithOptimizedConnections(unitcell::Unitcell)
    # build up new connections
    connections_new = Array[]
    # go through all old connections
    for c in unitcell.connections
        # check if the connection already is present in the 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[]
    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
    # return a new unitcell
    return Unitcell(
        unitcell.lattice_vectors,
        unitcell.basis,
        connections_new_2,
        unitcell.filename)
end
export getUnitcellWithOptimizedConnections

#-----------------------------------------------------------------------------------------------------------------------------
#
#   Add interactions along next-nearest neighbors
#
#-----------------------------------------------------------------------------------------------------------------------------
function addNextNearestNeighborsToConnections!(lattice::Lattice; strengthNNN="AUTO", connectionsNN=["ALL"])
    # get the connectivity matrix of the lattice
    connectivity = getConnectionList(lattice)
    # define a list of new connections
    connections_new = Array[]
    # iterate over all sites and check if they host a NNN connection
    for i in 1:size(lattice.positions,1)
        # check the neighbors of site i
        for (i1,c1) in enumerate(connectivity[i])
        for (i2,c2) in enumerate(connectivity[i])
            # check if the connections are valid
            if !(in("ALL",connectionsNN) || (in(c1[3],connectionsNN) && in(c2[3],connectionsNN)))
                continue
            end
            # propose a connection c1+c2 if they are different
            if i1 < i2
                # build the offsets
                if length(c1[4]) == 1
                    off_1 = (c2[4] - c1[4])
                    off_2 = (c1[4] - c2[4])
                elseif length(c1[4]) == 2
                    off_1 = (c2[4][1] - c1[4][1], c2[4][2] - c1[4][2])
                    off_2 = (c1[4][1] - c2[4][1], c1[4][2] - c2[4][2])
                elseif length(c1[4]) == 3
                    off_1 = (c2[4][1] - c1[4][1], c2[4][2] - c1[4][2], c2[4][3] - c1[4][3])
                    off_2 = (c1[4][1] - c2[4][1], c1[4][2] - c2[4][2], c1[4][3] - c2[4][3])
                end
                # build two new connections
                if strengthNNN == "AUTO"
                    connection_new_1 = [Int(c1[2]); Int(c2[2]); "NNN$(c1[2])to$(c2[2])"; off_1]
                    connection_new_2 = [Int(c2[2]); Int(c1[2]); "NNN$(c1[2])to$(c2[2])"; off_2]
                elseif strengthNNN == "MULTIPLY"
                    connection_new_1 = [Int(c1[2]); Int(c2[2]); "$(c1[3])*$(c2[3])"; off_1]
                    connection_new_2 = [Int(c2[2]); Int(c1[2]); "$(c1[3])*$(c2[3])"; off_2]
                else
                    connection_new_1 = [Int(c1[2]); Int(c2[2]); strengthNNN; off_1]
                    connection_new_2 = [Int(c2[2]); Int(c1[2]); strengthNNN; off_2]
                end
                # push them to the add list
                push!(connections_new, connection_new_1)
                push!(connections_new, connection_new_2)
            end
        end
        end
    end
    # push all connections to the lattice
    for c in connections_new
        push!(lattice.connections, c)
    end
end
function addNextNearestNeighborsToConnections!(unitcell::Unitcell; strengthNNN="AUTO", connectionsNN=["ALL"])
    # get the connectivity matrix of the unitcell
    connectivity = getConnectionList(unitcell)
    # define a list of new connections
    connections_new = Array[]
    # iterate over all sites and check if they host a NNN connection
    for i in 1:size(unitcell.basis,1)
        # check the neighbors of site i
        for (i1,c1) in enumerate(connectivity[i])
        for (i2,c2) in enumerate(connectivity[i])
            # check if the connections are valid
            if !(in("ALL",connectionsNN) || (in(c1[3],connectionsNN) && in(c2[3],connectionsNN)))
                continue
            end
            # propose a connection c1+c2 if they are different
            if i1 < i2
                # build the offsets
                if length(c1[4]) == 1
                    off_1 = (c2[4] - c1[4])
                    off_2 = (c1[4] - c2[4])
                elseif length(c1[4]) == 2
                    off_1 = (c2[4][1] - c1[4][1], c2[4][2] - c1[4][2])
                    off_2 = (c1[4][1] - c2[4][1], c1[4][2] - c2[4][2])
                elseif length(c1[4]) == 3
                    off_1 = (c2[4][1] - c1[4][1], c2[4][2] - c1[4][2], c2[4][3] - c1[4][3])
                    off_2 = (c1[4][1] - c2[4][1], c1[4][2] - c2[4][2], c1[4][3] - c2[4][3])
                end
                # build two new connections
                if strengthNNN == "AUTO"
                    connection_new_1 = [Int(c1[2]); Int(c2[2]); "NNN$(c1[2])to$(c2[2])"; off_1]
                    connection_new_2 = [Int(c2[2]); Int(c1[2]); "NNN$(c1[2])to$(c2[2])"; off_2]
                elseif strengthNNN == "MULTIPLY"
                    connection_new_1 = [Int(c1[2]); Int(c2[2]); "$(c1[3])*$(c2[3])"; off_1]
                    connection_new_2 = [Int(c2[2]); Int(c1[2]); "$(c1[3])*$(c2[3])"; off_2]
                else
                    connection_new_1 = [Int(c1[2]); Int(c2[2]); strengthNNN; off_1]
                    connection_new_2 = [Int(c2[2]); Int(c1[2]); strengthNNN; off_2]
                end
                # push them to the add list
                push!(connections_new, connection_new_1)
                push!(connections_new, connection_new_2)
            end
        end
        end
    end
    # push all connections to the unitcell
    for c in connections_new
        push!(unitcell.connections, c)
    end
end

# export addNextNearestNeighborsToConnections!


text:
# ----------------------------------------------------------------------------------- #
# Copyright (c) 2019 Varnerlab
# Robert Frederick Smith School of Chemical and Biomolecular Engineering
# Cornell University, Ithaca NY 14850
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.
# ----------------------------------------------------------------------------------- #
#
# ----------------------------------------------------------------------------------- #
# Function: calculate_transcription_control_array
# Description: Calculate the transcriptional control array at time t
# Generated on: 2019-10-02T16:15:21.68
#
# Input arguments:
# t::Float64 => Current time value (scalar)
# x::Array{Float64,1} => State array (number_of_species x 1)
# data_dictionary::Dict{String,Any} => Dictionary holding model parameters
#
# Output arguments:
# control_array::Array{Float64,1} => Transcriptional control array (number_of_genes x 1) at time t
# ----------------------------------------------------------------------------------- #
function calculate_transcription_control_array(t::Float64,x::Array{Float64,1},data_dictionary::Dict{String,Any})

	# initialize the control -
	control_array = zeros(4)

	# hack -
	x = abs.(x)

	# Alias the species -
	clssrA = x[1]
	deGFPssrA = x[2]
	sigma_28 = x[3]
	sigma_70 = x[4]
	mRNA_clssrA = x[5]
	mRNA_deGFPssrA = x[6]
	mRNA_sigma_28 = x[7]
	mRNA_sigma_70 = x[8]
	protein_clssrA = x[9]
	protein_deGFPssrA = x[10]
	protein_sigma_28 = x[11]
	protein_sigma_70 = x[12]

	# Alias the binding parameters -
	binding_parameter_dictionary = data_dictionary["binding_parameter_dictionary"]
	n_clssrA_sigma_28 = binding_parameter_dictionary["n_clssrA_sigma_28"]
	K_clssrA_sigma_28 = binding_parameter_dictionary["K_clssrA_sigma_28"]
	n_deGFPssrA_sigma_70 = binding_parameter_dictionary["n_deGFPssrA_sigma_70"]
	K_deGFPssrA_sigma_70 = binding_parameter_dictionary["K_deGFPssrA_sigma_70"]
	n_deGFPssrA_clssrA = binding_parameter_dictionary["n_deGFPssrA_clssrA"]
	K_deGFPssrA_clssrA = binding_parameter_dictionary["K_deGFPssrA_clssrA"]
	n_sigma_28_sigma_70 = binding_parameter_dictionary["n_sigma_28_sigma_70"]
	K_sigma_28_sigma_70 = binding_parameter_dictionary["K_sigma_28_sigma_70"]
	n_sigma_28_clssrA = binding_parameter_dictionary["n_sigma_28_clssrA"]
	K_sigma_28_clssrA = binding_parameter_dictionary["K_sigma_28_clssrA"]

	# Alias the control function parameters -
	control_parameter_dictionary = data_dictionary["control_parameter_dictionary"]
	W_clssrA_RNAP = control_parameter_dictionary["W_clssrA_RNAP"]
	W_clssrA_sigma_28 = control_parameter_dictionary["W_clssrA_sigma_28"]
	W_deGFPssrA_RNAP = control_parameter_dictionary["W_deGFPssrA_RNAP"]
	W_deGFPssrA_sigma_70 = control_parameter_dictionary["W_deGFPssrA_sigma_70"]
	W_deGFPssrA_clssrA = control_parameter_dictionary["W_deGFPssrA_clssrA"]
	W_sigma_28_RNAP = control_parameter_dictionary["W_sigma_28_RNAP"]
	W_sigma_28_sigma_70 = control_parameter_dictionary["W_sigma_28_sigma_70"]
	W_sigma_28_clssrA = control_parameter_dictionary["W_sigma_28_clssrA"]
	W_sigma_70_RNAP = control_parameter_dictionary["W_sigma_70_RNAP"]

	# Transfer function target:clssrA actor:sigma_28
	actor_set_clssrA_sigma_28 = [
		protein_sigma_28
	]
	actor = (prod(actor_set_clssrA_sigma_28))
	b_clssrA_sigma_28 = ((actor)^(n_clssrA_sigma_28))/((K_clssrA_sigma_28)^(n_clssrA_sigma_28)+(actor)^(n_clssrA_sigma_28))

	# Control function for clssrA -
	control_array[1] = (W_clssrA_RNAP+W_clssrA_sigma_28*b_clssrA_sigma_28)/(1+W_clssrA_RNAP+W_clssrA_sigma_28*b_clssrA_sigma_28)

	# Transfer function target:deGFPssrA actor:sigma_70
	actor_set_deGFPssrA_sigma_70 = [
		protein_sigma_70
	]
	actor = prod(actor_set_deGFPssrA_sigma_70)
	b_deGFPssrA_sigma_70 = ((actor)^(n_deGFPssrA_sigma_70))/((K_deGFPssrA_sigma_70)^(n_deGFPssrA_sigma_70)+(actor)^(n_deGFPssrA_sigma_70))

	# Transfer function target:deGFPssrA actor:clssrA
	actor_set_deGFPssrA_clssrA = [
		protein_clssrA
	]
	actor = prod(actor_set_deGFPssrA_clssrA)
	b_deGFPssrA_clssrA = (actor^(n_deGFPssrA_clssrA))/(K_deGFPssrA_clssrA^(n_deGFPssrA_clssrA)+actor^(n_deGFPssrA_clssrA))

	# Control function for deGFPssrA -
	control_array[2] = (W_deGFPssrA_RNAP+W_deGFPssrA_sigma_70*b_deGFPssrA_sigma_70)/(1+W_deGFPssrA_RNAP+W_deGFPssrA_sigma_70*b_deGFPssrA_sigma_70+W_deGFPssrA_clssrA*b_deGFPssrA_clssrA)

	# Transfer function target:sigma_28 actor:sigma_70
	actor_set_sigma_28_sigma_70 = [
		protein_sigma_70
	]
	actor = prod(actor_set_sigma_28_sigma_70)
	b_sigma_28_sigma_70 = (actor^(n_sigma_28_sigma_70))/(K_sigma_28_sigma_70^(n_sigma_28_sigma_70)+actor^(n_sigma_28_sigma_70))

	# Transfer function target:sigma_28 actor:clssrA
	actor_set_sigma_28_clssrA = [
		protein_clssrA
	]
	actor = prod(actor_set_sigma_28_clssrA)
	b_sigma_28_clssrA = (actor^(n_sigma_28_clssrA))/(K_sigma_28_clssrA^(n_sigma_28_clssrA)+actor^(n_sigma_28_clssrA))

	# Control function for sigma_28 -
	control_array[3] = (W_sigma_28_RNAP+W_sigma_28_sigma_70*b_sigma_28_sigma_70)/(1+W_sigma_28_RNAP+W_sigma_28_sigma_70*b_sigma_28_sigma_70+W_sigma_28_clssrA*b_sigma_28_clssrA)

	# Control function for sigma_70 -
	control_array[4] = (W_sigma_70_RNAP)/(1+W_sigma_70_RNAP)

	# return -
	return control_array
end

#
# ----------------------------------------------------------------------------------- #
# Function: calculate_translation_control_array
# Description: Calculate the translation control array at time t
# Generated on: 2019-10-02T16:15:21.684
#
# Input arguments:
# t::Float64 => Current time value (scalar)
# x::Array{Float64,1} => State array (number_of_species x 1)
# data_dictionary::Dict{String,Any} => Dictionary holding model parameters
#
# Output arguments:
# control_array::Array{Float64,1} => Translation control array (number_of_genes x 1) at time t
# ----------------------------------------------------------------------------------- #
function calculate_translation_control_array(t::Float64,x::Array{Float64,1},data_dictionary::Dict{String,Any})

	# initialize the control -
	control_array = ones(4)

	# for S70 - no translation
	control_array[4] = 0.0

	# correct for "resource"?
    correction_term = (x[13]/100.0)
    control_array = control_array*correction_term

	# return -
	return control_array
end

ext: jl

max_stars_repo_name: varnerlab/Biophysical-TXTL-Model-Code

max_issues_repo_name: varnerlab/Biophysical-TXTL-Model-Code

max_forks_repo_name: varnerlab/Biophysical-TXTL-Model-Code


text:
export export_png

function _nop_yes_no(flag::Union{Nothing,Bool},yes_no::Tuple{AbstractString,AbstractString})
    if flag !== nothing
        if flag
            return first(yes_no)
        else
            return last(yes_no)
        end
    end
    ""
end

function _enhanced(flag::Union{Nothing,Bool})
    _nop_yes_no(flag,("enhanced","noenhanced"))
end

# ================================================================

"""
```julia
export_png(gp::GnuplotScript,
           filename::AbstractString;
           enhanced::Bool)
```

Export the current plots as an `PNG` image
"""
function export_png(gp::GnuplotScript,filename::AbstractString;
                    enhanced::Union{Nothing,Bool}=nothing)

    # add .png ext to filename
    #
    filename = splitext(filename)
    if length(filename)>1 && filename[2]!=".png"
        @warn "File extension $(filename[2]) replaced by \".png\""
    end
    filename = first(filename)*".png"

    # command
    #
    command = ""
    command *= "set terminal push\n"
    command *= "set terminal png\n"

    command *= _enhanced(enhanced) * "\n"

    command *= "set output '$filename'\n"
    command *= "replot\n"
    command *= "set terminal pop\n"

    # append gnuplot instructions
    #
    _append_to_script(gp,command)

    filename
end


text:
# Canonical Correlation Analysis

#### CCA Type
"""
Canonical Correlation Analysis Model
"""
struct CCA{T<:Real} <: RegressionModel
    xmean::Vector{T}  # sample mean of X: of length dx (can be empty)
    ymean::Vector{T}  # sample mean of Y: of length dy (can be empty)
    xproj::Matrix{T}  # projection matrix for X, of size (dx, p)
    yproj::Matrix{T}  # projection matrix for Y, of size (dy, p)
    corrs::Vector{T}  # correlations, of length p

    function CCA(xm::Vector{T},
                 ym::Vector{T},
                 xp::Matrix{T},
                 yp::Matrix{T},
                 crs::Vector{T}) where T<:Real

        dx, px = size(xp)
        dy, py = size(yp)

        isempty(xm) || length(xm) == dx ||
            throw(DimensionMismatch("Incorrect length of xmean."))

        isempty(ym) || length(ym) == dy ||
            throw(DimensionMismatch("Incorrect length of ymean."))

        px == py ||
            throw(DimensionMismatch("xproj and yproj should have the same number of columns."))

        length(crs) == px ||
            throw(DimensionMismatch("Incorrect length of corrs."))

        new{T}(xm, ym, xp, yp, crs)
    end
end

## properties

"""
    size(M:CCA)

Return a tuple with the dimension of `X`, `Y`, and the output dimension.
"""
size(M::CCA) = (size(M.xproj, 1), size(M.yproj, 1), size(M.xproj, 2))

"""
    mean(M::CCA, c::Symbol)

Get the mean vector for the component `c` of the model `M`.
The component parameter can be `:x` or `:y`.
"""
function mean(M::CCA, c::Symbol)
    xi, yi, o = size(M)
    if c == :x
        fullmean(xi, M.xmean)
    elseif c == :y
        fullmean(yi, M.ymean)
    else
        throw(ArgumentError("Unknown component $c"))
    end
end

"""
    projection(M::CCA, c::Symbol)

Get the projection matrix for the component `c` of the model `M`.
The component parameter can be `:x` or `:y`.
"""
function projection(M::CCA, c::Symbol)
    if c == :x
        M.xproj
    elseif c == :y
        M.yproj
    else
        throw(ArgumentError("Unknown component $c"))
    end
end

"""
    cor(M::CCA)

The correlations of the projected componnents (a vector of length `p`).
"""
cor(M::CCA) = M.corrs

## use
"""
    predict(M::CCA, Z::AbstractVecOrMat{<:Real}, c::Symbol)

Given a [`CCA`](@ref) model, one can transform observations into both spaces into a common space, as

```math
\\mathbf{z}_x = \\mathbf{P}_x^T (\\mathbf{x} - \\boldsymbol{\\mu}_x) \\\\
\\mathbf{z}_y = \\mathbf{P}_y^T (\\mathbf{y} - \\boldsymbol{\\mu}_y)
```

Here, ``\\mathbf{P}_x`` and ``\\mathbf{P}_y`` are projection matrices for ``X`` and ``Y``;
``\\boldsymbol{\\mu}_x`` and ``\\boldsymbol{\\mu}_y`` are mean vectors.

Parameter `Z` can be either a vector of length `dx`, `dy`, or a matrix where each column is an observation. The component parameter `c` can be `:x` or `:y`.
"""
function predict(M::CCA, Z::AbstractVecOrMat{<:Real}, c::Symbol)
    if c == :x
        transpose(M.xproj) * centralize(Z, M.xmean)
    elseif c == :y
        transpose(M.yproj) * centralize(Z, M.ymean)
    else
        throw(ArgumentError("Unknown component $c"))
    end
end

## show

function show(io::IO, M::CCA)
    xi, yi, o = size(M)
    print(io, "CCA (xindim = $xi, yindim = $yi, outdim = $o)")
end

#### Perform CCA on data

"""
    ccacov(Cxx, Cyy, Cxy, xmean, ymean, p)

Compute CCA based on analysis of the given covariance matrices, using generalized
eigenvalue decomposition, and return [`CCA`](@ref) model.

Parameters:
- `Cxx`: The covariance matrix of `X`.
- `Cyy`: The covariance matrix of `Y`.
- `Cxy`: The covariance matrix between `X` and `Y`.
- `xmean`: The mean vector of the **original** samples of `X`, which can be
a vector of length `dx`, or an empty vector indicating a zero mean.
- `ymean`: The mean vector of the **original** samples of `Y`, which can be
a vector of length `dy`, or an empty vector indicating a zero mean.
- `p`: The output dimension, *i.e* the dimension of the common space.
"""
function ccacov(Cxx::DenseMatrix{T},
                Cyy::DenseMatrix{T},
                Cxy::DenseMatrix{T},
                xmean::Vector{T},
                ymean::Vector{T},
                p::Int) where T<:Real

    # argument checking
    dx, dx2 = size(Cxx)
    dy, dy2 = size(Cyy)
    dx == dx2 || error("Cxx must be a square matrix.")
    dy == dy2 || error("Cyy must be a square matrix.")
    size(Cxy) == (dx, dy) ||
        throw(DimensionMismatch("size(Cxy) should be equal to (dx, dy)"))

    isempty(xmean) || length(xmean) == dx ||
        throw(DimensionMismatch("Incorrect length of xmean."))

    isempty(ymean) || length(ymean) == dy ||
        throw(DimensionMismatch("Incorrect length of ymean."))

    1 <= p <= min(dx, dy) ||
        throw(DimensionMismatch(""))

    _ccacov(Cxx, Cyy, Cxy, xmean, ymean, p)
end

function _ccacov(Cxx, Cyy, Cxy, xmean, ymean, p::Int)
    dx = size(Cxx, 1)
    dy = size(Cyy, 1)

    # solve Px and Py

    if dx <= dy
        # solve Px: (Cxy * inv(Cyy) * Cyx) Px = λ Cxx * Px
        # compute Py: inv(Cyy) * Cyx * Px

        G = cholesky(Cyy) \ Cxy'
        Ex = eigen(Symmetric(Cxy * G), Symmetric(Cxx))
        ord = sortperm(Ex.values; rev=true)
        vx, Px = extract_kv(Ex, ord, p)
        Py = qnormalize!(G * Px, Cyy)
    else
        # solve Py: (Cyx * inv(Cxx) * Cxy) Py = λ Cyy Py
        # compute Px: inv(Cx) * Cxy * Py

        H = cholesky(Cxx) \ Cxy
        Ey = eigen(Symmetric(Cxy'H), Symmetric(Cyy))
        ord = sortperm(Ey.values; rev=true)
        vy, Py = extract_kv(Ey, ord, p)
        Px = qnormalize!(H * Py, Cxx)
    end

    # compute correlations
    # Note: Px' * Cxx * Px == I
    #       Py' * Cyy * Py == I
    crs = coldot(Px, Cxy * Py)

    # construct CCA model
    CCA(xmean, ymean, Px, Py, crs)
end

"""
    ccasvd(Zx, Zy, xmean, ymean, p)

Compute CCA based on singular value decomposition of centralized sample matrices `Zx` and `Zy`, and return [`CCA`](@ref) model[^1].

Parameters:
- `Zx`: The centralized sample matrix for `X`.
- `Zy`: The centralized sample matrix for `Y`.
- `xmean`: The mean vector of the **original** samples of `X`, which can be
a vector of length `dx`, or an empty vector indicating a zero mean.
- `ymean`: The mean vector of the **original** samples of `Y`, which can be
a vector of length `dy`, or an empty vector indicating a zero mean.
- `p`: The output dimension, *i.e* the dimension of the common space.
"""
function ccasvd(Zx::DenseMatrix{T},
                Zy::DenseMatrix{T},
                xmean::Vector{T},
                ymean::Vector{T},
                p::Int) where T<:Real

    dx, n = size(Zx)
    dy, n2 = size(Zy)
    n == n2 ||
        throw(DimensionMismatch("Zx and Zy must have the same number of columns."))

    isempty(xmean) || length(xmean) == dx ||
        throw(DimensionMismatch("Incorrect length of xmean."))

    isempty(ymean) || length(ymean) == dy ||
        throw(DimensionMismatch("Incorrect length of ymean."))

    1 <= p <= min(dx, dy) ||
        throw(DimensionMismatch(""))

    _ccasvd(Zx, Zy, xmean, ymean, p)
end

# The implementation is partly based on:
#
#   David Weenink.
#   Canonical Correlation Analysis.
#   Institute of Phonetic Sciences, Univ. of Amsterdam,
#   Proceedings 25 (2003), 81-99.
#
#   Note: in this paper, each row is considered as an observation.
#   The algorithm is adpated to the column-major format here.
#
function _ccasvd(Zx::DenseMatrix{T}, Zy::DenseMatrix{T}, xmean::Vector{T}, ymean::Vector{T}, p::Int) where T<:Real
    # svd factorization of Z

    n = size(Zx, 2)

    # svd decomposition
    Sx = svd(Zx)
    Sy = svd(Zy)
    S = svd!(Sx.Vt * transpose(Sy.Vt)) # svd of Vx * Vy'

    # compute Px and Py
    ord = sortperm(S.S; rev=true)
    si = ord[1:p]
    Px = rmul!(Sx.U, Diagonal(1.0 ./ Sx.S)) * S.U[:, si]
    Py = rmul!(Sy.U, Diagonal(1.0 ./ Sy.S)) * S.V[:, si]

    # scale so that Px' * Cxx * Py == I
    #           and Py' * Cyy * Py == I,
    #
    # with Cxx = Zx * Zx' / (n - 1)
    #      Cyy = Zy * Zy' / (n - 1)
    #
    rmul!(Px, sqrt(n-1))
    rmul!(Py, sqrt(n-1))

    # compute correlations
    crs = rmul!(coldot(Zx'Px, Zy'Py), one(T)/(n-1))

    # construct CCA model
    CCA(xmean, ymean, Px, Py, crs)
end

## interface functions

"""
    fit(CCA, X, Y; ...)

Perform CCA over the data given in matrices `X` and `Y`.
Each column of `X` and `Y` is an observation.

`X` and `Y` should have the same number of columns (denoted by `n` below).

This method returns an instance of [`CCA`](@ref).

**Keyword arguments:**
- `method`: The choice of methods:
    - `:cov`: based on covariance matrices
    - `:svd`: based on SVD of the input data (*default*)
- `outdim`: The output dimension, *i.e* dimension of the common space (*default*: `min(dx, dy, n)`)
- `mean`: The mean vector, which can be either of:
    - `0`: the input data has already been centralized
    - `nothing`: this function will compute the mean (*default*)
    - a pre-computed mean vector

**Notes:** This function calls [`ccacov`](@ref) or [`ccasvd`](@ref) internally, depending on the choice of method.
"""
function fit(::Type{CCA}, X::AbstractMatrix{T}, Y::AbstractMatrix{T};
             outdim::Int=min(min(size(X)...), min(size(Y)...)),
             method::Symbol=:svd,
             xmean=nothing,
             ymean=nothing) where T<:Real

    dx, n = size(X)
    dy, n2 = size(Y)

    n2 == n ||
        throw(DimensionMismatch("X and Y should have the same number of columns."))

    (n >= dx && n >= dy) ||
        @warn("CCA would be numerically instable when n < dx or n < dy.")

    xmv = preprocess_mean(X, xmean)
    ymv = preprocess_mean(Y, ymean)

    Zx = centralize(X, xmv)
    Zy = centralize(Y, ymv)

    if method == :cov
        Cxx = rmul!(Zx*transpose(Zx), inv(n - 1))
        Cyy = rmul!(Zy*transpose(Zy), inv(n - 1))
        Cxy = rmul!(Zx*transpose(Zy), inv(n - 1))
        M = ccacov(Cxx, Cyy, Cxy, xmv, ymv, outdim)
    elseif method == :svd
        M = ccasvd(Zx, Zy, xmv, ymv, outdim)
    else
        error("Invalid method name $(method)")
    end

    return M::CCA
end


text:
#    Problem : GROUP A
#    *********
#	 The extended Woods problem.
#	 This problem is a sum of n/4 sets of 6 terms, each of which is
#    assigned its own group.  For a given set i, the groups are
#    A(i), B(i), C(i), D(i), E(i) and F(i). Groups A(i) and C(i) contain 1
#    nonlinear element each, denoted Y(i) and Z(i).
#
#    The problem dimension is defined from the number of these sets.
#    The number of problem variables is then 4 times larger.
#
#	 This version uses a slightly unorthodox expression of Woods
#    function as a sum of squares (see Buckley)
#
#    Origonal SIF Source: problem 14 in
#    J.J. More', B.S. Garbow and K.E. Hillstrom,
#    "Testing Unconstrained Optimization Software",
#    ACM Transactions on Mathematical Software, vol. 7(1), pp. 17-41, 1981.
#
#    See also Toint#27, Buckley#17 (p. 101), Conn, Gould, Toint#7
#
#    WOODS.SIF classification SUR2-AN-V-0
#
#	 This problem is decomposed in n linear groups, the last n-1 of which
#    are 2 x 2 and singular.
#
#    NS is the number of sets (= n/4)
#
# Daniel Henderson, 08/2021

f = (x) -> begin
	fx = 0.0
	for i in 1:Int(lastindex(x)/4)
		fx += (100*(x[4*i-2]-x[4*i-3]^2)^2 + (1-x[4*i-3])^2 + 90*(x[4*i]-x[4*i-1]^2)^2 + (1-x[4*i-1])^2 + 10*(x[4*i-2]+x[4*i]-2)^2 + 0.1*(x[4*i-2]-x[4*i])^2)
	end
    return fx
end

g! = (g, x) -> begin
	for i in 1:Int(lastindex(x)/4)
		g[4i-3] = -2(1-x[4i-3]) - 400*x[4i-3]*(x[4i-2]-x[4i-3]^2)
		g[4i-2] = 200(x[4i-2]-x[4i-3]^2) + 0.2(x[4i-2]-x[4i]) + 20(x[4i-2] + x[4i]-2)
		g[4i-1] = -2(1-x[4i-1]) - 360*x[4i-1]*(x[4i]-x[4i-1]^2)
		g[4i] = -0.2(x[4i-2]-x[4i]) + 20(x[4i-2] + x[4i]-2) + 180(x[4i]-x[4i-1]^2)
	end
    return g
end

fg! = (g, x) -> begin
	fx = 0.0
	for i in 1:Int(lastindex(x)/4)
		fx += (100*(x[4*i-2]-x[4*i-3]^2)^2 + (1-x[4*i-3])^2 + 90*(x[4*i]-x[4*i-1]^2)^2 + (1-x[4*i-1])^2 + 10*(x[4*i-2]+x[4*i]-2)^2 + 0.1*(x[4*i-2]-x[4*i])^2)
		g[4i-3] = -2(1-x[4i-3]) - 400*x[4i-3]*(x[4i-2]-x[4i-3]^2)
		g[4i-2] = 200(x[4i-2]-x[4i-3]^2) + 0.2(x[4i-2]-x[4i]) + 20(x[4i-2] + x[4i]-2)
		g[4i-1] = -2(1-x[4i-1]) - 360*x[4i-1]*(x[4i]-x[4i-1]^2)
		g[4i] = -0.2(x[4i-2]-x[4i]) + 20(x[4i-2] + x[4i]-2) + 180(x[4i]-x[4i-1]^2)
	end
    return fx, g
end

init = (n::Int=4000) -> begin
	mod(n, 4) > 0 && @warn "WOODS: number of variables must be divisible by 4"
	q = max(1, div(n, 4))
	n = 4*q

	x0 = [j % 2 == 1 ? -3.0 : -1.0 for j in 1:n]
    return n, x0
end

TestSet["WOODS"] = UncProgram("WOODS", f, g!, fg!, init)


text:
using Compose: x_measure, y_measure, Rotation, Vec, Vec2, svg_print_float, resolve, Form, FormPrimitive,
    @makeform, Transform, UnitBox, AbsoluteBox, empty_tag, indent, print_vector_properties

export mathjax

# MathJaxPrimitive form
struct MathJaxPrimitive{P<:Vec, S<:Vec, R<:Rotation, O<:Vec} <: FormPrimitive
    position::P
    size::S
    value::AbstractString
    rot::R
    offset::O
end

const MathJax{P<:MathJaxPrimitive} = Form{P}
Compose.form_string(::MathJax) = "MATH"

"""
    mathjax(x, y, width, heigh, value [,rot::Rotation])

Draw the text `value` at the position (`x`,`y`) relative to the current context.
Parameters `width` and `height` specify the bounding box size.
"""
function mathjax(x, y, width, height, value::AbstractString,
              rot=Rotation(), offset::Vec2=(0mm,0mm);
              tag::Symbol=empty_tag)
    moffset = (x_measure(offset[1]), y_measure(offset[2]))
    prim = MathJaxPrimitive((x_measure(x), y_measure(y)), (x_measure(width), y_measure(height)), string(value), rot, moffset)
    MathJax{typeof(prim)}([prim], tag)
end

"""
    mathjax(xs::AbstractArray, ys::AbstractArray, ws::AbstractArray, hs::AbstractArray, values::AbstractArray [,rots::Rotation])

Arguments can be passed in arrays in order to perform multiple drawing operations at once.
"""
mathjax(xs::AbstractArray, ys::AbstractArray, ws::AbstractArray, hs::AbstractArray, values::AbstractArray{<:AbstractString},
              rots::AbstractArray=[Rotation()], offsets::AbstractArray=[(0mm,0mm)];
              tag::Symbol=empty_tag) =
    @makeform (x in xs, y in ys, w in ws, h in hs, value in values, rot in rots, offset in offsets),
    MathJaxPrimitive((x_measure(x), y_measure(y)), (x_measure(w), y_measure(h)), value, rot, (x_measure(offset[1]), y_measure(offset[2]))) tag

function Compose.resolve(box::AbsoluteBox, units::UnitBox, t::Transform, p::MathJaxPrimitive)
    rot = resolve(box, units, t, p.rot)
    return MathJaxPrimitive(
                resolve(box, units, t, p.position),
                resolve(box, units, t, p.size),
                p.value, rot, p.offset)
end

function Compose.boundingbox(form::MathJaxPrimitive, linewidth::Measure,
                     font::AbstractString, fontsize::Measure)

    width, height = form.size
    x0 = form.position.x
    y0 = form.position.y

    return BoundingBox(x0 - linewidth + form.offset[1],
                       y0 - linewidth + form.offset[2],
                       width + linewidth,
                       height + linewidth)
end

function Compose.draw(img::SVG, prim::MathJaxPrimitive, idx::Int)
    indent(img)

    img.indentation += 1
    print(img.out, "<g transform=\"translate(")
    svg_print_float(img.out, prim.position[1].value)
    print(img.out, ",")
    svg_print_float(img.out, prim.position[2].value)
    print(img.out, ")\"")
    print_vector_properties(img, idx, true)
    print(img.out, ">\n")
    indent(img)

    print(img.out, "<g class=\"primitive\">\n")
    img.indentation += 1
    indent(img)
    print(img.out, """<foreignObject x="0" y="0" """)

    if abs(prim.rot.theta) > 1e-4 || sum(abs.(prim.offset)) > 1e-4mm
        print(img.out, " transform=\"")
        if abs(prim.rot.theta) > 1e-4
            print(img.out, "rotate(")
            svg_print_float(img.out, rad2deg(prim.rot.theta))
            print(img.out, ",")
            svg_print_float(img.out, prim.rot.offset[1].value-prim.position[1].value)
            print(img.out, ", ")
            svg_print_float(img.out, prim.rot.offset[2].value-prim.position[2].value)
            print(img.out, ")")
        end
        if sum(abs.(prim.offset)) > 1e-4mm
            print(img.out, "translate(")
            svg_print_float(img.out, prim.offset[1].value)
            print(img.out, ",")
            svg_print_float(img.out, prim.offset[2].value)
            print(img.out, ")")
        end
        print(img.out, "\"")
    end

    print(img.out, """ width="$(prim.size[1])" height="$(prim.size[2])">
          <div xmlns="http://www.w3.org/1999/xhtml">
    \\(\\displaystyle{$(prim.value)}\\)
    </div>
</foreignObject>""")

    img.indentation -= 1
    indent(img)
    print(img.out, "</g>\n")
    img.indentation -= 1
    indent(img)
    print(img.out, "</g>\n")
end


text:
include("benchmarks/sort/compare.jl")
isdir(".mempool") && rm(".mempool"; recursive=true)
include("benchmarks/sort2/compare.jl")
isdir(".mempool") && rm(".mempool"; recursive=true)
include("benchmarks/deepdag/compare.jl")
isdir(".mempool") && rm(".mempool"; recursive=true)


text:
""" Helper function for merging learner arguments. """
function merge_kwargs(valid_keys; kwargs...)
    nkwargs = Dict{Symbol, Any}()
    for item in kwargs
        if item.first in valid_keys
            nkwargs[item.first] = item.second
        end
    end
    return nkwargs
end

""" Preprocesses merging of kwargs for IAI.fit!, to avoid errors!
TODO: Add nkwargs that can be changed! """
function fit_regressor_kwargs(; kwargs...)
    nkwargs = Dict{Symbol, Any}()
    for item in kwargs
        if item.first in keys(nkwargs)
            nkwargs[item.first] = item.second
        end
    end
    return nkwargs
end

""" Preprocesses merging of kwargs for IAI.fit!, to avoid errors!
TODO: Add nkwargs that can be changed! """
function fit_classifier_kwargs(; kwargs...)
    nkwargs = Dict{Symbol, Any}(:sample_weight => :autobalance)
    for item in kwargs
        if item.first in keys(nkwargs)
            nkwargs[item.first] = item.second
        end
    end
    return nkwargs
end

""" Function that preprocesses merging of kwargs for IAI.OptimalTreeClassifier, to avoid errors! """
function classifier_kwargs(; kwargs...)
    valid_keys = [:random_seed, :max_depth, :cp, :minbucket,
                  :fast_num_support_restarts, :localsearch,
                  :ls_num_hyper_restarts, :ls_num_tree_restarts,
                  :hyperplane_config, :criterion]
    merge_kwargs(valid_keys; kwargs...)
end

""" Function that preprocesses merging of kwargs for IAI.OptimalTreeRegressor, to avoid errors! """
function regressor_kwargs(; kwargs...)
    valid_keys = [:random_seed, :max_depth, :cp, :minbucket,
                  :fast_num_support_restarts, :localsearch,
                  :ls_num_hyper_restarts, :ls_num_tree_restarts,
                  :hyperplane_config,
                  :regression_sparsity, :regression_weighted_betas,
                  :regression_lambda, :criterion]
    return merge_kwargs(valid_keys; kwargs...)
end

"""
    learn_from_data!(X::DataFrame, Y::AbstractArray, lnr::IAI.OptimalTreeLearner,
                          idxs::Union{Nothing, Array}=nothing; kwargs...)

Wrapper around IAI.fit! for constraint learning.
Arguments:
-lnr: OptimalTreeClassifier or OptimalTreeRegressor
-X: matrix of feature data
-Y: matrix of constraint data.
Returns fitted Learner corresponding to the data.
NOTE: kwargs from `learn_constraint!` get unpacked here.
"""
function learn_from_data!(X::DataFrame, Y::AbstractArray, lnr::Union{IAI.OptimalTreeLearner,
                                                                     IAI.Heuristics.RandomForestLearner},
                          idxs::Union{Nothing, Array}=nothing; kwargs...)
    n_samples, n_features = size(X);
    @assert n_samples == length(Y);
    # Making sure that we only consider relevant features.
    if !isnothing(idxs)
        IAI.set_params!(lnr, split_features = idxs)
        if typeof(lnr) in [IAI.OptimalTreeRegressor, IAI.Heuristics.RandomForestRegressor]
            IAI.set_params!(lnr, regression_features = idxs)
        end
    else
        IAI.set_params!(lnr, split_features = All())
        if typeof(lnr) in [IAI.OptimalTreeRegressor, IAI.Heuristics.RandomForestRegressor]
            IAI.set_params!(lnr, regression_features = All())
        end
    end
    IAI.fit!(lnr, X, Y; kwargs...)
    return lnr
end

""" Checks that a BlackBoxClassifier has enough feasible/infeasible samples. """
function check_feasibility(bbc::BlackBoxClassifier)
    return bbc.feas_ratio >= get_param(bbc, :threshold_feasibility)
end

check_feasibility(bbr::BlackBoxRegressor) = true

function check_feasibility(gm::GlobalModel)
    return [check_feasibility(bbl) for bbl in gm.bbls]
end

""" Checks that a BlackBoxLearner.learner has adequate accuracy."""
function check_accuracy(bbc::BlackBoxClassifier)
    return bbc.accuracies[end] >= get_param(bbc, :threshold_accuracy)
end

check_accuracy(bbr::BlackBoxRegressor) = 1. # TODO: use MSE?

function check_sampled(bbl::BlackBoxLearner)
    size(bbl.X, 1) == 0 && throw(OCTHaGOnException(string("BlackBoxLearner ", bbl.name, " must be sampled first.")))
    return
end

"""
    get_random_trees(lnr::IAI.Heuristics.RandomForestLearner)
Returns one of the trees of RandomForestLearner.
"""
function get_random_trees(lnr::IAI.Heuristics.RandomForestLearner)
    trees = []
    for i = 1:lnr.num_trees
        tree = IAI.clone(lnr.inner)
        tree.prb_ = lnr.prb_
        tree.tree_ = lnr.forest_.trees[i]
        push!(trees, tree)
    end
    return trees
end

"""
    boundify(lnr::OptimalTreeRegressor, X::DataFrame, Y, solver::MOI.OptimizerWithAttributes)
    boundify(lnr::IAI.Heuristics.RandomForestRegressor, X::DataFrame, Y, solver::MOI.OptimizerWithAttributes)
    boundify(lnr::OptimalTreeClassifier, X::DataFrame, Y, hypertype::String = "lower", solver::MOI.OptimizerWithAttributes)
Returns bounding hyperplanes of an OptimalTreeLearner.
"""
function boundify(lnr::IAI.OptimalTreeRegressor, X::DataFrame, Y, solver::MOI.OptimizerWithAttributes)
    feas_leaves = find_leaves(lnr)
    leaf_idx = IAI.apply(lnr, X)
    ul_data = Dict()
    for leaf in feas_leaves
        idx = findall(x -> x == leaf, leaf_idx)
        ul_data[-leaf] = u_regress(X[idx,:], Y[idx], solver) # negative leaf shows upper bounding
        ul_data[leaf] = l_regress(X[idx,:], Y[idx], solver)
    end
    return ul_data
end

function boundify(lnr::IAI.OptimalTreeClassifier, X::DataFrame, Y, solver::MOI.OptimizerWithAttributes, hypertype::String = "lower")
    all_leaves = find_leaves(lnr)
    leaf_idx = IAI.apply(lnr, X)
    data = Dict()
    feas_leaves = [i for i in all_leaves if Bool(IAI.get_classification_label(lnr, i))]
    for leaf in feas_leaves
        idx = findall(x -> x == leaf, leaf_idx)
        if hypertype == "upper"
            data[-leaf] = u_regress(X[idx,:], Y[idx], solver)
        elseif hypertype == "lower"
            data[leaf] = l_regress(X[idx,:], Y[idx], solver)
        end
    end
    return data
end

function boundify(lnr::IAI.Heuristics.RandomForestRegressor,
                  X::DataFrame, Y, solver::MOI.OptimizerWithAttributes)
    trees = get_random_trees(lnr)
    data = Dict()
    for i=1:length(trees)
        ul_data = Dict()
        tree = trees[i]
        all_leaves = find_leaves(tree)
        leaf_idx = IAI.apply(tree, X)
        for leaf in all_leaves
            idx = findall(x -> x == leaf, leaf_idx)
            ul_data[-leaf] = u_regress(X[idx,:], Y[idx], solver) # negative leaf shows upper bounding
            ul_data[leaf] = l_regress(X[idx,:], Y[idx], solver)
        end
        data[i] = ul_data
    end
    return data
end

function boundify(lnr::IAI.Heuristics.RandomForestClassifier,
                    X::DataFrame, Y, solver::MOI.OptimizerWithAttributes, hypertype::String = "lower")
    trees = get_random_trees(lnr)
    data = Dict()
    for i=1:length(trees)
        ul_data = Dict()
        tree = trees[i]
        all_leaves = find_leaves(tree)
        leaf_idx = IAI.apply(tree, X)
        feas_leaves = [j for j in all_leaves if Bool(IAI.get_classification_label(tree, j))]
        for leaf in feas_leaves
            idx = findall(x -> x == leaf, leaf_idx)
            if hypertype == "upper"
                ul_data[-leaf] = u_regress(X[idx,:], Y[idx], solver)
            elseif hypertype == "lower"
                ul_data[leaf] = l_regress(X[idx,:], Y[idx], solver)
            end
        end
        data[i] = ul_data
    end
    return data
end

"""
    learn_constraint!(bbl::Union{GlobalModel, BlackBoxLearner, Array}; kwargs...)

Constructs a constraint tree from a BlackBoxLearner and dumps in bbo.learners.
Arguments:
    bbl: OCT structs (GM, bbl, Array{bbl})
    kwargs: arguments for learners and fits.
"""
function learn_constraint!(bbc::BlackBoxClassifier; kwargs...)
    check_sampled(bbc)
    set_param(bbc, :reloaded, false) # Makes sure that we know trees are retrained.
    lnr = base_classifier()
    IAI.set_params!(lnr; minbucket =
        maximum([2*length(bbc.vars)/length(bbc.Y), lnr.minbucket]))
    if bbc.equality # Equalities are approximated more aggressively.
        IAI.set_params!(lnr; max_depth = 6, ls_num_tree_restarts = 30, fast_num_support_restarts = 15)
    end
    IAI.set_params!(lnr; classifier_kwargs(; kwargs...)...)
    if check_feasibility(bbc) || get_param(bbc, :ignore_feasibility)
        lnr = learn_from_data!(bbc.X, bbc.Y .>= 0, lnr; fit_classifier_kwargs(; kwargs...)...)
        push!(bbc.learners, lnr)
        push!(bbc.accuracies, IAI.score(lnr, bbc.X, bbc.Y .>= 0)) # TODO: add ability to specify criterion.
        push!(bbc.learner_kwargs, Dict(kwargs))
    else
        throw(OCTHaGOnException("Not enough feasible samples for BlackBoxClassifier " * string(bbc.name) * "."))
    end
    return
end

function learn_constraint!(bbr::BlackBoxRegressor, threshold::Pair = Pair("reg", nothing);  kwargs...)
    classifications = ["upper", "lower", "upperlower"]
    check_sampled(bbr)
    get_param(bbr, :reloaded) && set_param(bbr, :reloaded, false) # Makes sure that we know trees are retrained.
    if bbr.convex && !bbr.equality
        return # If convex, don't train a tree!
    elseif threshold.first in classifications
        lnr = base_classifier()
        IAI.set_params!(lnr; minbucket =
            maximum([2*length(bbr.vars)/length(bbr.Y), lnr.minbucket]),
            classifier_kwargs(; kwargs...)...)
        ul_data = Dict()
        if threshold.first == "upper" # Upper bounding classifier with upper bounds in leaves
            lnr = learn_from_data!(bbr.X, bbr.Y .<= threshold.second, lnr; fit_classifier_kwargs(; kwargs...)...)
            merge!(ul_data, boundify(lnr, bbr.X, bbr.Y, SOLVER_SILENT, "upper"))
        elseif threshold.first == "lower" # Lower bounding classifier with lower bounds in leaves
            lnr = learn_from_data!(bbr.X, bbr.Y .>= threshold.second, lnr; fit_classifier_kwargs(; kwargs...)...)
            merge!(ul_data, boundify(lnr, bbr.X, bbr.Y, SOLVER_SILENT, "lower"))
        elseif threshold.first == "upperlower" # Upper bounding classifier with bounds in each leaf
            lnr = learn_from_data!(bbr.X, bbr.Y .<= threshold.second, lnr; fit_classifier_kwargs(; kwargs...)...)
            merge!(ul_data, boundify(lnr, bbr.X, bbr.Y, SOLVER_SILENT, "lower"))
            merge!(ul_data, boundify(lnr, bbr.X, bbr.Y, SOLVER_SILENT, "upper"))
        end
        push!(bbr.learners, lnr);
        push!(bbr.learner_kwargs, Dict(kwargs))
        push!(bbr.thresholds, threshold)
        push!(bbr.ul_data, ul_data)
        push!(bbr.accuracies, IAI.score(lnr, bbr.X, bbr.Y .>= 0))
    elseif threshold.first == "reg"
        lnr = base_regressor()
        IAI.set_params!(lnr; minbucket =
            maximum([2*length(bbr.vars)/length(bbr.Y), lnr.minbucket]),
            regressor_kwargs(; kwargs...)...)
        if bbr.equality # Equalities cannot leverage convexity unfortunately...
            if threshold.second == nothing
                lnr = learn_from_data!(bbr.X, bbr.Y, lnr; fit_regressor_kwargs(; kwargs...)...)
            else
                idxs = findall(y -> y <= threshold.second, bbr.Y)
                lnr = learn_from_data!(bbr.X[idxs, :], bbr.Y[idxs], lnr; fit_regressor_kwargs(; kwargs...)...)
            end
        elseif threshold.second == nothing && bbr.local_convexity < 0.75
            lnr = learn_from_data!(bbr.X, bbr.Y, lnr; fit_regressor_kwargs(; kwargs...)...)
        elseif threshold.second == nothing && bbr.local_convexity >= 0.75
            lnr = learn_from_data!(bbr.X, bbr.curvatures .> 0, lnr; fit_regressor_kwargs(; kwargs...)...)
        elseif bbr.local_convexity < 0.75
            idxs = findall(y -> y <= threshold.second, bbr.Y)
            lnr = learn_from_data!(bbr.X[idxs, :], bbr.Y[idxs], lnr; fit_regressor_kwargs(; kwargs...)...)
        elseif bbr.local_convexity >= 0.75
            idxs = findall(y -> y <= threshold.second, bbr.Y)
            lnr = learn_from_data!(bbr.X[idxs, :], bbr.curvatures[idxs] .> 0, lnr; fit_regressor_kwargs(; kwargs...)...)
        end
        push!(bbr.learners, lnr);
        push!(bbr.learner_kwargs, Dict(kwargs))
        push!(bbr.thresholds, threshold)
        push!(bbr.ul_data, boundify(lnr, bbr.X, bbr.Y, SOLVER_SILENT))
        push!(bbr.accuracies, IAI.score(lnr, bbr.X, bbr.Y))
    elseif threshold.first == "rfreg"
        lnr = base_rf_regressor()
        bbr.local_convexity < 0.75 || throw(OCTHaGOnException("Cannot use RandomForestRegressor " *
        "on BBR $(bbr.name) since it is almost convex."))
        IAI.set_params!(lnr; minbucket =
            maximum([2*length(bbr.vars)/length(bbr.Y), lnr.minbucket]),
            regressor_kwargs(; kwargs...)...)
        if threshold.second == nothing
            lnr = learn_from_data!(bbr.X, bbr.Y, lnr; fit_regressor_kwargs(; kwargs...)...)
        else
            idxs = findall(y -> y <= threshold.second, bbr.Y)
            lnr = learn_from_data!(bbr.X[idxs, :], bbr.Y[idxs], lnr; fit_regressor_kwargs(; kwargs...)...)
        end
        push!(bbr.learners, lnr);
        push!(bbr.learner_kwargs, Dict(kwargs))
        push!(bbr.thresholds, threshold)
        push!(bbr.ul_data, boundify(lnr, bbr.X, bbr.Y, SOLVER_SILENT))
        push!(bbr.accuracies, IAI.score(lnr, bbr.X, bbr.Y))
    else
        throw(OCTHaGOnException("$(threshold.first) is not a valid learner type for" *
            " thresholded learning of BBR $(bbr.name)."))
    end
    return
end

function learn_constraint!(bbl::Array; kwargs...)
   for fn in bbl
        learn_constraint!(fn; kwargs...)
   end
end

learn_constraint!(gm::GlobalModel; kwargs...) =
    learn_constraint!(gm.bbls; kwargs...)


text:
module TestGlobal

using Recyclers
using Test

Recyclers.@global RE1 = Recyclers.ShardedRecycler(() -> [1, 2, 3])
atexit() do
    Recyclers.unsafe_destroy!(RE1)
end

function test_init_sharded()
    @test length(RE1.cache.pools) == Threads.nthreads()
end

end  # module


text:
using BenchmarkTools

# preallocate arrays
A = rand(100,100)
B = rand(100,100)

function inner_rows(A,B)
    for i in 1:100, j in 1:100
    C[i,j] = A[i,j] + B[i,j]
    end
end
@btime inner_rows(A,B)

function inner_cols(A,B)
    for j in 1:100, i in 1:100
    C[i,j] = A[i,j] + B[i,j]
    end
end
@btime inner_cols(A,B)

function inner_cols(A,B)
    for j in 1:100, i in 1:100
    val = [A[i,j] + B[i,j]]
    C[i,j] = val[1]
    end
end
@btime inner_cols(A,B)

function inner_cols(A,B)
    for j in 1:100, i in 1:100
    val = A[i,j] + B[i,j]
    C[i,j] = val
    end
end
@btime inner_cols(A,B)

using StaticArrays
val = SVector{3,Float64}(1.0,2.0,3.0)
typeof(val)

function static_inner_alloc(A,B)
  for j in 1:100, i in 1:100
    val = @SVector [A[i,j] + B[i,j]]
    C[i,j] = val[1]
  end
end
@btime static_inner_alloc(A,B)

function no_prealloc(A,B)
    C = rand(100,100)
    for j in 1:100, i in 1:100
    val = A[i,j] + B[i,j]
    C[i,j] = val
    end
end
@btime no_prealloc(A,B)

C = rand(100,100)
function w_prealloc!(C,A,B)
    for j in 1:100, i in 1:100
    val = A[i,j] + B[i,j]
    C[i,j] = val
    end
end
@btime w_prealloc!(C,A,B)

function no_prealloc(A,B)
    C = rand(100,100)
    for j in 1:100, i in 1:100
    val = A[i,j] + B[i,j]
    C[i,j] = val
    end
end
@btime no_prealloc(A,B)

function f(A,B)
    # 10*(A+B)
    sum([A + B for k in 1:10]) # calculate 10 arrays and sum them
end
@btime f(A,B)

function h(A,B)
    # 10*(A+B)
    C = similar(A) # create array for output
    for k in 1:10
        C .+= (A + B) # add all of A+B first then add this temp array to C
    end
end
@btime h(A,B)

function g(A,B)
    # 10*(A+B)
    C = similar(A) # create array for output
    for k in 1:10
        C .+= A.+ B # reuse memory in each iteration
    end
end
@btime g(A,B)

function g!(C,A,B)
    # 10*(A+B)
    for k in 1:10
        C .+= A.+ B # reuse memory in each iteration
    end
end
@btime g!(C,A,B)

function g!(C,A,B)
    # 10*(A+B)
    for k in 1:10
        C .+= A.+ B # reuse memory in each iteration
    end
    return nothing
end
C = zeros(100,100)
g!(C,A,B)
C

D = Matrix{Float64}(undef, 100, 100)
function f_vector!(D,C,A,B)
    D .*= A .* B.* C
end
@btime f_vector!(D,C,A,B)

function f_loop!(D,C,A,B)
    for i in 1:length(D)
        @inbounds D[i] = A[i] * B[i] * C[i]
    end
end
@btime f_loop!(D,C,A,B)

function f_alloc(A)
    B = A[1:5,1:5]
end

function f_point(A)
    B = @view A[1:5,1:5]
end

@btime f_alloc(A)
@btime f_point(A)

f_alloc(A)
B = f_alloc(A)

function basic()
    A = zeros(10,10)
    B = A .- 2
    C = A + B
    D = A - B

    return B, D
end

B, D = basic()

ff(x,y) = 2x + 2y
ff(x::Int,y::Int) = 2x + y
ff(x::Float64,y::Float64) = x/y

@show ff(2,5)
@show ff(2.0,5.0)

ff(2.0,5)

ff(x::Float64,y::Number) = 2x

ff(2.0,5)

ff(x::Float64,y::Number) = 5x + 2y
ff(x::Number,y::Int) = x - y
ff(2.0,5)

f(x,y) = x + y
x = Number[1.0,3]
function q(x)
  a = 4
  b = 2
  c = f(x[1],a)
  d = f(b,c)
  f(d,x[2])
  return nothing
end

using BenchmarkTools
@btime q(x)

x = [1.0,3.0]
@btime q(x)

x = Number[1.0,3]
function r(x)
  a = 4
  b = 2
  for i in 1:100
    c = f(x[1],a)
    d = f(b,c)
    a = f(d,x[2])
  end
  a
end
@btime r(x)

s(x) = _s(x[1],x[2])
function _s(x1,x2)
  a = 4
  b = 2
  for i in 1:100
    c = f(x1,a)
    d = f(b,c)
    a = f(d,x2)
  end
  a
end
@btime s(x)

typeof(x)
typeof([x[1],x[2]])

A = rand(100,100)
B = rand(100,100)
C = rand(100,100)
@btime for j in 1:100, i in 1:100
  global A,B,C
  C[i,j] = A[i,j] + B[i,j]
end

function f(A,B,C)
    for j in 1:100, i in 1:100
        C[i,j] = A[i,j] + B[i,j]
    end
end
@btime f(A,B,C)


text:
using JuMP, Combinatorics, Gurobi

include("./algorithms/TEAVAR.jl")
include("./algorithms/FFC.jl")
include("./algorithms/ARROW.jl")
include("./algorithms/OPTIMAL.jl")
include("./algorithms/HYPERTHETICAL.jl")

## locate the cut fiberlinks, get affected IP links from cut fiber
## for example, cutfiber = [6,14], which is (1,6), (6,1), the down ip links are (1,5) (1,6) (5,1) (6,1)
function FailureLocator(OpticalTopo, optical_failure_scenario)
    damage_links_index = []

    for l in 1:length(optical_failure_scenario)
        if optical_failure_scenario[l] == 0
            push!(damage_links_index, l)
        end
    end
    return damage_links_index
end

## Traffic engineering module
function TrafficEngineering(GRB_ENV, IPTopo, OpticalTopo, algorithm, links, capacity, demand, flows, T, Tf, scenarios, fiber_scenarios, scenario_probs, scenario_restored_bw, optical_rerouting_K, tunnelK, beta, purb, verbose, solve_or_not)
    best_options = ones(size(scenarios,1))  # only randomized rouding will return this other just 1
    best_scenario_resored_bw = []

    nflows = size(flows,1)
    nedges = size(links,1)
    ntunnels = tunnelK
    TunnelBw = zeros(nflows, ntunnels)
    FlowBw = zeros(nflows)
    var = 0
    throughput = 0
    solve_runtime = 0
    best_options = []
    best_scenario_resored_bw = 0

    if algorithm == "TEAVAR"
        var, cvar, TunnelBw, scenario_loss, loss, solve_t, opt_t = TEAVAR(GRB_ENV, links, capacity, flows, demand, beta, T, Tf, scenarios, scenario_probs, average=true)
        for f in 1:nflows
            FlowBw[f] = demand[f] * (1-var)
        end
        throughput = sum(FlowBw)
        solve_runtime = solve_t

    elseif algorithm == "FFC2"
        IPScenarios, OpticalScenarios = GetAllScenarios(IPTopo, OpticalTopo, 0, true, 2)  # all double fiber cut scenarios
        scenarios = IPScenarios["code"][2:end]  # FFC only take failure scenarios except healthy state
        TunnelBw, FlowBw, obj, solve_t, opt_t = FFC(GRB_ENV, links, capacity, flows, demand, scenarios, T, Tf) # here the ffc TunnelBw is the reserved bw of tunnels > actual demand
        throughput = sum(FlowBw)
        solve_runtime = solve_t

    elseif algorithm == "FFC1"
        IPScenarios, OpticalScenarios = GetAllScenarios(IPTopo, OpticalTopo, 0, true, 1)  # all single fiber cut scenarios
        scenarios = IPScenarios["code"][2:end]  # FFC only take failure scenarios except healthy state
        TunnelBw, FlowBw, obj, solve_t, opt_t = FFC(GRB_ENV, links, capacity, flows, demand, scenarios, T, Tf) # here the ffc TunnelBw is the reserved bw of tunnels > actual demand
        throughput = sum(FlowBw)
        solve_runtime = solve_t

    elseif algorithm == "ECMP"
        SplitRatio = ones(size(Tf,1), ntunnels) ./ ntunnels   # ECMP is equal SplitRatio
        TunnelBw = SplitRatio .* demand
        FlowBw = demand
        throughput = sum(FlowBw)
        solve_runtime = 0

    elseif algorithm == "ARROW_NAIVE"   # one restoration option made at planning stage
        # here we only put one option inside the TE and do not make selections
        TunnelBw, FlowBw, obj, solve_t, opt_t = ARROW_ONE(GRB_ENV, links, capacity, scenario_restored_bw, flows, demand, scenarios, T, Tf, verbose)
        throughput = sum(FlowBw)
        solve_runtime = solve_t

    elseif algorithm == "ARROW"   # ARROW TE
        TunnelBw, FlowBw, obj, v_restore, solve_t_1, opt_t_1 = ARROW(GRB_ENV, links, capacity, scenario_restored_bw, flows, demand, scenarios, T, Tf, purb, verbose)
        best_options, virtual_restore, best_option_diviation = RestoreMatchMinLoss(v_restore, scenario_restored_bw) # find best option for each failure scenario
        if verbose println("best_options: ", best_options) end
        if verbose println("virtual_restore: ", virtual_restore) end
        if verbose println("purb for sum delta: ", purb) end
        if verbose println("best_option_diviation ", best_option_diviation) end

        best_scenario_resored_bw = []
        nscenarios = size(scenario_restored_bw,1)
        nedges = size(links,1)
        for s in 1:nscenarios
            best_restore_bw = []
            for e in 1:nedges
                push!(best_restore_bw, scenario_restored_bw[s][Int(best_options[s])][e])
            end
            push!(best_scenario_resored_bw, best_restore_bw)
        end
        if verbose
            print("ARROW best_scenario_resored_bw: sum=")
            println(sum(best_scenario_resored_bw))
            for rr in 1:size(best_scenario_resored_bw, 1)
                println(best_scenario_resored_bw[rr])
            end
            println("\n")
        end

        TunnelBw, FlowBw, obj, solve_t_2, opt_t_2 = ARROW_ONE(GRB_ENV, links, capacity, best_scenario_resored_bw, flows, demand, scenarios, T, Tf, verbose)

        throughput = sum(FlowBw)
        solve_runtime = solve_t_1 + solve_t_2

    elseif algorithm == "ARROW_BIN"   # ARROW BINARY TE
        TunnelBw, FlowBw, obj, ticket_selection, solve_t, opt_t = ARROW_BINARY(GRB_ENV, links, capacity, scenario_restored_bw, flows, demand, scenarios, T, Tf)
        best_options = []
        for s in 1:size(ticket_selection,1)
            for z in 1:size(ticket_selection,2)
                if ticket_selection[s,z] == 0
                    push!(best_options, (1-ticket_selection[s,z])*z)
                    break
                end
            end
        end
        best_scenario_resored_bw = []
        nscenarios = size(scenario_restored_bw,1)
        nedges = size(links,1)
        for s in 1:nscenarios
            best_restore_bw = []
            for e in 1:nedges
                push!(best_restore_bw, scenario_restored_bw[s][Int(best_options[s])][e])
            end
            push!(best_scenario_resored_bw, best_restore_bw)
        end
        if verbose
            print("ARROW best_scenario_resored_bw: sum=")
            println(sum(best_scenario_resored_bw))
            for rr in 1:size(best_scenario_resored_bw, 1)
                println(best_scenario_resored_bw[rr])
            end
            println("\n")
        end

        throughput = sum(FlowBw)
        println("best_options ", best_options)
        solve_runtime = solve_t

    elseif algorithm == "OPTIMAL"   # optimal cross-layer TE (super ILP with optical restoration RWA and IP TE)
        ILP_LP = 1
        time_limit = 36000  # 10 hours maximum optimization runtime
        if solve_or_not
            TunnelBw, FlowBw, obj, solve_t, opt_t, var_num, sum_constraint = OPTIMAL(GRB_ENV, links, capacity, flows, demand, scenarios, fiber_scenarios, T, Tf, OpticalTopo, optical_rerouting_K, ILP_LP, time_limit, solve_or_not)
            throughput = sum(FlowBw)
            solve_runtime = solve_t
        else
            var_num, sum_constraint = OPTIMAL(GRB_ENV, links, capacity, flows, demand, scenarios, fiber_scenarios, T, Tf, OpticalTopo, optical_rerouting_K, ILP_LP, time_limit, solve_or_not)
        end

    elseif algorithm == "OPTIMAL_LP"   # optimal cross-layer TE (super ILP with optical restoration RWA and IP TE)
        ILP_LP = 0
        time_limit = 36000  # 10 hours maximum optimization runtime
        if solve_or_not
            TunnelBw, FlowBw, obj, solve_t, opt_t, var_num, sum_constraint = OPTIMAL(GRB_ENV, links, capacity, flows, demand, scenarios, fiber_scenarios, T, Tf, OpticalTopo, optical_rerouting_K, ILP_LP, time_limit, solve_or_not)
            throughput = sum(FlowBw)
            solve_runtime = solve_t
        else
            var_num, sum_constraint = OPTIMAL(GRB_ENV, links, capacity, flows, demand, scenarios, fiber_scenarios, T, Tf, OpticalTopo, optical_rerouting_K, ILP_LP, time_limit, solve_or_not)
        end

    elseif algorithm == "HYPERTHETICAL"  # a hyperthetical TE that has no failure scenarios (or all failures are restored)
        TunnelBw, FlowBw, solve_t, opt_t = HYPERTHETICAL(GRB_ENV, links, capacity, flows, demand, T, Tf) # here the ffc TunnelBw is the reserved bw of tunnels = actual demand (because no failure)
        throughput = sum(FlowBw)
        solve_runtime = solve_t

    else
        printstyled("Error: TE Algorithm not indentified.", color=:red)
    end

    return TunnelBw, FlowBw, var, throughput, solve_runtime, best_options, best_scenario_resored_bw
end

function RestoreMatchMinLoss(v_restore, scenario_restored_bw)
    noptions = size(scenario_restored_bw[1], 1)
    nscenarios = size(scenario_restored_bw,1)
    nedges = size(scenario_restored_bw[1][1],1)
    best_option = ones(nscenarios)
    virtual_restore = zeros(nscenarios)

    for s in 1:nscenarios
        min_loss = 9999999999
        max_sum_restorable = 0
        for z in 1:noptions
            current_loss = 0
            current_sum_restorable = sum(scenario_restored_bw[s][z])
            ## choose the option with the min positive virtual restorable capacity
            for e in 1:nedges
                if v_restore[s,z,e] > 0
                    current_loss += v_restore[s,z,e]
                end
            end

            if current_loss < min_loss  # optimize for minimum virtual restorable capacity
                best_option[s] = z
                virtual_restore[s] = current_loss
                min_loss = current_loss
            elseif current_loss == min_loss  # if virtual restorable capacity to be the same, max restorable capacity
                if current_sum_restorable > max_sum_restorable
                    best_option[s] = z
                    max_sum_restorable = current_sum_restorable
                end
            end
        end
    end
    best_option_diviation = []
    return best_option, virtual_restore, best_option_diviation
end

max_stars_repo_name: hipersys-team/arrow

max_issues_repo_name: arrow-bow/arrow

max_forks_repo_name: arrow-bow/arrow


text:
using TwitterUploadMedia
using Test
using Base64

TwitterUploadMedia.twitterauth(ENV["CONSUMER_KEY"],
                ENV["CONSUMER_SECRET"],
                ENV["ACCESS_TOKEN"],
                ENV["ACCESS_TOKEN_SECRET"])

@testset "POST media/upload" begin
    postmediaupload = post_media_upload(media_data=base64encode(open("lena.png")))
    postmedia_id = postmediaupload["media_id"]
    @test typeof(postmeidaupload) == Media
    @test postmediaupload["image_type"] == "image/jpeg"
    @test typeof(postmedia_id) == Int
end

@testset "GET media/upload (STATUS)" begin
    getmediaupload = get_media_upload(command="STATUS", media_id=postmedia_id)
    @test getmediaupload != Nothing
    # WIP
end