Graphs, k-centers and recursion


#1

Hello!
I have been recently trying to improve the performance of the classic brute force algorithm for k-centers in a simple graph. My reasons for this are not really that serious, I was only trying to reach the speed of C++ using the same algorithm or see how close I could get.

The definition of problem goes like this: for a given simple graph G, with weighted edges find k vertices K such that the maximum length of the shortest path from a vertex in G to any vertex in K is as short as possible.
For example imagine trying to build k shops on a grid of n lots so that the customers living on those lots will be able to get to a shop as quickly as possible.

The algorithm that solves that and is implemented below does this by recursively going through all the possible combinations of vertices in set K of size k, rating them and then selecting the best one (the one with smallest distance to centers, that appears first).

using StaticArrays

#Input graph as a list - to conserve space
L = [[(2, 30), (29, 6), (95, 31), (100, 88)],[(1, 30), (3, 46)],[(2, 46), (4, 1), (79,57), (99, 49)],[(3, 1), (5, 28), (27, 27), (33, 45), (68, 28), (73, 41)],[(4, 28), (6, 31), (7, 8), (80, 90), (96, 87)],[(5, 31), (7, 69), (63, 83)],[(5, 8), (6, 69), (8, 39), (19, 9), (35, 7), (62, 55)],[(7, 39), (9, 14), (13, 77), (49, 96), (52, 31), (53, 66)],[(8, 14), (10, 84), (48, 92), (83, 75), (99, 28)],[(9, 84), (11, 59)],[(10, 59), (12, 10), (37, 4)],[(11, 10), (13, 28), (44, 56)],[(8, 77), (12, 28), (14, 63), (42, 3), (85, 22)],[(13, 63), (15, 9), (21, 82), (27, 48), (85, 90)],[(14, 9), (16, 100), (50,31), (69, 46)],[(15, 100), (17, 98), (38, 70), (79, 96)],[(16, 98), (18, 70), (54, 45), (58, 60)],[(17, 70),(19, 94), (38, 27), (75, 100)],[(7, 9), (18, 94), (20, 30), (60, 53), (92, 70)],[(19, 30), (21, 14), (93, 35)],[(14, 82), (20, 14), (22, 87), (39, 59), (81, 34)],[(21, 87),(23, 82), (52, 11), (75, 48)],[(22, 82), (24, 55)],[(23, 55), (25, 2), (59, 63)],[(24, 2), (26, 32), (99, 19)],[(25, 32), (27, 77), (29, 20), (32, 26)],[(4, 27), (14, 48), (26, 77), (28, 95), (50, 16), (60, 86)],[(27,95), (29, 29)],[(1, 6), (26, 20), (28, 29), (30, 59)],[(29, 59), (31, 91), (56, 76), (57, 37), (70, 74), (81, 99), (95, 96)],[(30, 91), (32, 89), (63, 82), (69, 53)],[(26, 26), (31, 89), (33, 50), (99, 35)],[(4, 45), (32, 50), (34, 40),(49, 77), (73, 88), (88, 8)],[(33, 40), (35, 88), (94, 52)],[(7, 7), (34, 88), (36, 94), (60, 5), (80, 46)],[(35, 94), (37, 60)],[(11, 4), (36, 60), (38, 21), (81, 14)],[(16, 70), (18, 27), (37, 21), (39, 89), (62, 56)],[(21, 59), (38, 89), (40, 47), (76, 89)],[(39, 47), (41, 63),(54, 52)],[(40, 63), (42, 45), (44, 80), (76, 42)],[(13, 3), (41, 45), (43, 46), (45, 68), (57, 5), (85, 34)],[(42, 46), (44, 24)],[(12, 56), (41, 80), (43, 24), (45, 77)],Tuple{Int64,Int64}[(42, 68), (44, 77), (46, 60), (59, 33), (68, 19)],[(45, 60), (47, 45)],[(46, 45), (48, 50), (58, 68)],[(9, 92), (47, 50), (49, 93), (98, 53)],[(8, 96), (33, 77), (48, 93), (50, 22)],[(15, 31), (27, 16), (49, 22), (51,84)],[(50, 84), (52, 16), (64, 89)],[(8, 31), (22, 11), (51, 16), (53, 85), (60, 34), (95, 96)],[(8, 66), (52, 85), (54, 68), (74, 50)],[(17, 45), (40, 52), (53, 68), (55, 93), (59, 75), (86, 35)],[(54, 93), (56, 37), (95, 47)],[(30, 76), (55, 37), (57, 26)],[(30, 37), (42, 5), (56, 26), (58, 29)],[(17, 60), (47, 68), (57, 29), (59, 38), (73, 37)],[(24, 63), (45, 33), (54, 75), (58, 38), (60, 10)],[(19, 53), (27, 86), (35, 5), (52, 34), (59, 10), (61, 32)],[(60, 32), (62, 67), (65, 73)],[(7, 55), (38, 56), (61, 67), (63, 66), (94, 96)],[(6, 83), (31, 82), (62, 66), (64, 52)],Tuple{Int64,Int64}[(51, 89), (63, 52), (65, 19)],[(61, 73), (64, 19), (66, 39), (94, 61)],[(65, 39), (67, 12), (97, 81)],[(66, 12), (68, 86)],[(4, 28), (45, 19), (67, 86), (69, 72)],[(15, 46), (31, 53), (68, 72), (70, 73), (85, 51)],[(30, 74), (69, 73), (71, 65)],[(70, 65), (72, 2)],[(71, 2), (73, 8)],[(4, 41), (33, 88), (58, 37), (72, 8), (74, 96)],[(53, 50), (73, 96), (75, 43)],[(18, 100), (22, 48), (74, 43), (76, 39), (91, 33)],[(39, 89), (41, 42), (75, 39), (77, 61), (91, 51)],[(76, 61), (78, 90)],[(77, 90), (79, 8), (91, 32)],[(3, 57), (16, 96), (78, 8), (80, 58)],[(5, 90), (35, 46), (79, 58), (81, 91), (91, 84)],[(21, 34), (30, 99), (37, 14), (80, 91), (82, 58), (95, 54)],[(81, 58), (83, 13)],[(9, 75),(82, 13), (84, 79)],[(83, 79), (85, 59)],[(13, 22), (14, 90), (42, 34), (69, 51), (84, 59), (86, 28)],[(54, 35), (85, 28), (87, 46)],[(86, 46), (88, 24), (91, 17)],[(33, 8), (87, 24), (89, 63)],[(88, 63), (90, 81)],[(89, 81), (91, 14), (94, 63), (95, 96)],[(75, 33), (76, 51), (78, 32), (80, 84), (87, 17), (90, 14), (92, 52)],[(19, 70), (91, 52), (93, 64)],[(20, 35), (92, 64), (94, 75)],[(34, 52), (62, 96), (65, 61), (90, 63), (93, 75), (95, 71)],[(1, 31), (30, 96), (52, 96), (55, 47), (81, 54), (90, 96), (94, 71), (96, 51)],[(5, 87), (95, 51), (97, 75)],[(66, 81), (96, 75), (98, 57)],[(48, 53), (97, 57), (99, 31)],[(3, 49), (9, 28), (25, 19), (32, 35), (98, 31), (100, 49)],[(1, 88), (99, 49)]]

#A Function to read this list.
function graphListToMatrix(L)
    n = length(L)
    G = zeros(Int, n, n)
    for el in enumerate(L)
        j, links = el
        for link in links
            G[link[1], j] = link[2]
        end
    end
    return G
end

#A classic Floyd Warshall algorithm.
@inbounds function floydWarshall(G)
    const n = size(G, 1)
    D = copy(G)
    for k = 1:n
        for j = 1:(n - 1)
            for i = (j + 1):n
                if (D[i, k] != 0) & (D[j, k] != 0)
                    newDist = D[i, k] + D[j, k]

                    if D[i, j] == 0
                        D[i, j] = newDist
                        D[j, i] = newDist
                    else
                        if D[i, j] > newDist
                            D[i, j] = newDist
                            D[j, i] = newDist
                        end
                    end
                end
            end
        end
    end
    return D
end

#The main function
@inbounds function findBestKCenters(G, k)
    n = size(G, 1)
    bestScore = typemax(Int)
    bestChoice = MVector{k, Int}(zeros(k))
    D = floydWarshall(G)
    bestScore = useBruteForceRec(MVector{k, Int}(zeros( k)), D, 1, 1, k, n, bestScore, bestChoice)
    return bestScore
end

#The recursive function
 @inbounds function useBruteForceRec(vec, D, i, d, k, n, bestScore, bestChoice)
     if d == (k + 1)
         score = rateAChoice(D, vec)
          if bestScore < score
              bestScore = score
              bestChoice = vec
          end
     else
         for j = i:n
             vec[d] = j
             bestScore = useBruteForceRec(vec, D, i + 1, d + 1, k, n, bestScore, bestChoice)
         end
     end
     return bestScore
 end

#Function to evaluate a choice of center
@views @inbounds function rateAChoice(D, vec)
    const n = size(D, 1)
    score = 0
    for j = 1:n
        distanceToCenter = typemax(Int)
        for i in vec
            distanceToCenter = min(distanceToCenter, D[i, j])
        end
        score = max(score, distanceToCenter)
    end
    return score
end

#Benchmark
function repeatedTest(k, r)
    t = 0
    G =   graphListToMatrix(L)

    for j = 1:r
        t += @elapsed findBestKCenters(G, k)
    end
    println("Average time: $(t/r)")
    return t/r
end

repeatedTest(3, 5)

On my computer I get (2nd run and after that) arround 0.27 s.
Similar C++ code however runs 4x as fast.

If you happen to have any ideas on how to improve this code, please share them!
Those without special non-standard packages will be especially welcom!

The C++ code I’m comparing Julia with can be found here.
The same algorithm can be executed by

./main -i "samples/pmed/pmed1.txt" -f "pmed" -m "bf" -c 3

#2

Did you try profiling it?


#3

For one thing, @inbounds needs to be in front of the forloop, not in front of the function.

const also doesn’t do anything inside a function.

I would recomend BenchmarkTools for benchmarking. They run repeatedly (changing the number of runs based on how fast the function is, so you get plenty of samples for something fast without waiting an eternity for something slow), track allocations, and give you min/median/mean/max times.

Anyway, I made a few changes, reducing runtime by about a third:

using Compat, StaticArrays

#Input graph as a list - to conserve space
const L = [[(2, 30), (29, 6), (95, 31), (100, 88)],[(1, 30), (3, 46)],[(2, 46), (4, 1), (79,57), (99, 49)],[(3, 1), (5, 28), (27, 27), (33, 45), (68, 28), (73, 41)],[(4, 28), (6, 31), (7, 8), (80, 90), (96, 87)],[(5, 31), (7, 69), (63, 83)],[(5, 8), (6, 69), (8, 39), (19, 9), (35, 7), (62, 55)],[(7, 39), (9, 14), (13, 77), (49, 96), (52, 31), (53, 66)],[(8, 14), (10, 84), (48, 92), (83, 75), (99, 28)],[(9, 84), (11, 59)],[(10, 59), (12, 10), (37, 4)],[(11, 10), (13, 28), (44, 56)],[(8, 77), (12, 28), (14, 63), (42, 3), (85, 22)],[(13, 63), (15, 9), (21, 82), (27, 48), (85, 90)],[(14, 9), (16, 100), (50,31), (69, 46)],[(15, 100), (17, 98), (38, 70), (79, 96)],[(16, 98), (18, 70), (54, 45), (58, 60)],[(17, 70),(19, 94), (38, 27), (75, 100)],[(7, 9), (18, 94), (20, 30), (60, 53), (92, 70)],[(19, 30), (21, 14), (93, 35)],[(14, 82), (20, 14), (22, 87), (39, 59), (81, 34)],[(21, 87),(23, 82), (52, 11), (75, 48)],[(22, 82), (24, 55)],[(23, 55), (25, 2), (59, 63)],[(24, 2), (26, 32), (99, 19)],[(25, 32), (27, 77), (29, 20), (32, 26)],[(4, 27), (14, 48), (26, 77), (28, 95), (50, 16), (60, 86)],[(27,95), (29, 29)],[(1, 6), (26, 20), (28, 29), (30, 59)],[(29, 59), (31, 91), (56, 76), (57, 37), (70, 74), (81, 99), (95, 96)],[(30, 91), (32, 89), (63, 82), (69, 53)],[(26, 26), (31, 89), (33, 50), (99, 35)],[(4, 45), (32, 50), (34, 40),(49, 77), (73, 88), (88, 8)],[(33, 40), (35, 88), (94, 52)],[(7, 7), (34, 88), (36, 94), (60, 5), (80, 46)],[(35, 94), (37, 60)],[(11, 4), (36, 60), (38, 21), (81, 14)],[(16, 70), (18, 27), (37, 21), (39, 89), (62, 56)],[(21, 59), (38, 89), (40, 47), (76, 89)],[(39, 47), (41, 63),(54, 52)],[(40, 63), (42, 45), (44, 80), (76, 42)],[(13, 3), (41, 45), (43, 46), (45, 68), (57, 5), (85, 34)],[(42, 46), (44, 24)],[(12, 56), (41, 80), (43, 24), (45, 77)],Tuple{Int64,Int64}[(42, 68), (44, 77), (46, 60), (59, 33), (68, 19)],[(45, 60), (47, 45)],[(46, 45), (48, 50), (58, 68)],[(9, 92), (47, 50), (49, 93), (98, 53)],[(8, 96), (33, 77), (48, 93), (50, 22)],[(15, 31), (27, 16), (49, 22), (51,84)],[(50, 84), (52, 16), (64, 89)],[(8, 31), (22, 11), (51, 16), (53, 85), (60, 34), (95, 96)],[(8, 66), (52, 85), (54, 68), (74, 50)],[(17, 45), (40, 52), (53, 68), (55, 93), (59, 75), (86, 35)],[(54, 93), (56, 37), (95, 47)],[(30, 76), (55, 37), (57, 26)],[(30, 37), (42, 5), (56, 26), (58, 29)],[(17, 60), (47, 68), (57, 29), (59, 38), (73, 37)],[(24, 63), (45, 33), (54, 75), (58, 38), (60, 10)],[(19, 53), (27, 86), (35, 5), (52, 34), (59, 10), (61, 32)],[(60, 32), (62, 67), (65, 73)],[(7, 55), (38, 56), (61, 67), (63, 66), (94, 96)],[(6, 83), (31, 82), (62, 66), (64, 52)],Tuple{Int64,Int64}[(51, 89), (63, 52), (65, 19)],[(61, 73), (64, 19), (66, 39), (94, 61)],[(65, 39), (67, 12), (97, 81)],[(66, 12), (68, 86)],[(4, 28), (45, 19), (67, 86), (69, 72)],[(15, 46), (31, 53), (68, 72), (70, 73), (85, 51)],[(30, 74), (69, 73), (71, 65)],[(70, 65), (72, 2)],[(71, 2), (73, 8)],[(4, 41), (33, 88), (58, 37), (72, 8), (74, 96)],[(53, 50), (73, 96), (75, 43)],[(18, 100), (22, 48), (74, 43), (76, 39), (91, 33)],[(39, 89), (41, 42), (75, 39), (77, 61), (91, 51)],[(76, 61), (78, 90)],[(77, 90), (79, 8), (91, 32)],[(3, 57), (16, 96), (78, 8), (80, 58)],[(5, 90), (35, 46), (79, 58), (81, 91), (91, 84)],[(21, 34), (30, 99), (37, 14), (80, 91), (82, 58), (95, 54)],[(81, 58), (83, 13)],[(9, 75),(82, 13), (84, 79)],[(83, 79), (85, 59)],[(13, 22), (14, 90), (42, 34), (69, 51), (84, 59), (86, 28)],[(54, 35), (85, 28), (87, 46)],[(86, 46), (88, 24), (91, 17)],[(33, 8), (87, 24), (89, 63)],[(88, 63), (90, 81)],[(89, 81), (91, 14), (94, 63), (95, 96)],[(75, 33), (76, 51), (78, 32), (80, 84), (87, 17), (90, 14), (92, 52)],[(19, 70), (91, 52), (93, 64)],[(20, 35), (92, 64), (94, 75)],[(34, 52), (62, 96), (65, 61), (90, 63), (93, 75), (95, 71)],[(1, 31), (30, 96), (52, 96), (55, 47), (81, 54), (90, 96), (94, 71), (96, 51)],[(5, 87), (95, 51), (97, 75)],[(66, 81), (96, 75), (98, 57)],[(48, 53), (97, 57), (99, 31)],[(3, 49), (9, 28), (25, 19), (32, 35), (98, 31), (100, 49)],[(1, 88), (99, 49)]];

#A Function to read this list.
function graphListToMatrix(L)
    n = length(L)
    G = zeros(Int, n, n)
    for (j, links) in enumerate(L)
        @inbounds for link in links
            G[link[1], j] = link[2]
        end
    end
    return G
end

#A classic Floyd Warshall algorithm.function floydWarshall(G)
    n = size(G, 1)
    D = copy(G)
    @inbounds for k = 1:n, j = 1:(n - 1), i = (j + 1):n
        if (D[i, k] != 0) & (D[j, k] != 0)
            newDist = D[i, k] + D[j, k]

            if D[i, j] == 0
                D[i, j] = newDist
                D[j, i] = newDist
            else
                if D[i, j] > newDist
                    D[i, j] = newDist
                    D[j, i] = newDist
                end
            end
        end
    end
    return D
end

#The main function
function findBestKCenters(G, ::Val{k}) where k
    n = size(G, 1)
    bestScore = typemax(Int)
    bestChoice = @MVector zeros(Int, k)
    D = floydWarshall(G)
    bestScore = useBruteForceRec( (@MVector zeros(Int, k)) , D, 1, 1, k, n, bestScore, bestChoice)
    return bestScore
end

#The recursive function
function useBruteForceRec(vec, D, i, d, k, n, bestScore, bestChoice)
     if d == (k + 1)
         score = rateAChoice(D, vec)
          if bestScore > score
              bestScore = score
              bestChoice .= vec
          end
     else
         for j = i:n
            @inbounds vec[d] = j
            bestScore = useBruteForceRec(vec, D, i + 1, d + 1, k, n, bestScore, bestChoice)
         end
     end
     return bestScore
 end

#Function to evaluate a choice of center
function rateAChoice(D, vec)
    n = size(D, 1)
    score = 0
    @inbounds for j = 1:n
        distanceToCenter = typemax(Int)
        for i in vec
            distanceToCenter = min(distanceToCenter, D[i, j])
        end
        score = max(score, distanceToCenter)
    end
    return score
end

#Benchmark
function repeatedTest(::Val{k}, r) where k
    t = 0.0
    G =   graphListToMatrix(L)

    for j = 1:r
        t += @elapsed findBestKCenters(G, Val(k))
    end
    println("Average time: $(t/r)")
    return t/r
end

repeatedTest(Val(3), 5)

But, I’m pretty sure there’s something wrong:

G =   graphListToMatrix(L);
findBestKCenters(G, 3) == typemax(Int) #true 
findBestKCenters(G, Val(3)) == typemax(Int) #true

So I switched the direction of the inequality in useBruteForceRec.


#4

One thing that immediately jumped out at me was this:

Use the short-circuiting && for this.

Also, we have a pretty good Floyd Warshall implementation in LightGraphs - you might want to try implementing with SimpleWeightedGraphs.jl to see if you can take advantage of the optimization work we’ve been doing.


#5

I am doing some similar stuff. Just for context: best k-center is NP complete, so we are ultimately putting lipstick on a pig. Whatever we do, this will be slow; and at some point you want to use “real” algorithms.

Second. You actually have a giant n choose k decision tree. You want two things: Fast rejection of failed candidates, and early pruning of entire branches. Improvements are supposed to be rare.

So, we don’t compute diameters; we adjust our graph after each improvement by dropping edges, then then check whether we have a k-dominator. For proper speed, we do this via bitwise operations. Not using recursion would probaby be better, but I am too lazy right now to write a direct version.

So, first some convenience function for bit matrices (speed doesn’t matter here, but I had them lying arround):

@inline function getbit(B::NTuple{N,UInt64},i) where N
    i1, i2 = Base.get_chunks_id(i)
    u = 0x8000000000000000 >> i2
    @inbounds r = (B[i1] & u) != 0
    return r
end
@inline function setbit!(Bptr::Ptr{NTuple{N,UInt64}}, i, b) where N 
    i1, i2 = Base.get_chunks_id(i)
   # @boundscheck i1 <= N || error("oob")
    u = 0x8000000000000000 >> i2
    ptr = convert(Ptr{UInt64}, Bptr)
    bt = unsafe_load(ptr,i1)
    oldbit = bt & u 
    val = ifelse(b, bt | u, bt & (~u))
    unsafe_store!(ptr, val, i1)
    return oldbit!=0
end

@inline function getbit(B::Vector{NTuple{N,UInt64}},i,j) where N
    i1, i2 = Base.get_chunks_id(i)
    u = 0x8000000000000000 >> i2
    @inbounds r = (B[j][i1] & u) != 0
    return r
end
@inline function setbit!(B::Vector{NTuple{N,UInt64}}, i, j, b) where N 
    Bptr=pointer(B,j)
    setbit!(Bptr,i,b)
end

Everything becomes nicer when we write our bitmatrices as the complement of the adjaceny matrix. A set vec is a dominator if &(comp_adj[n] for n in vec)==0. For pruning, we also compute a trail_mat that tells us which points would become connected if all remaining vertices are taken. This allows us to quit the recursion early if some early vertices are required. But we should really do proper graph restructuring after every update. So, our bitmatrices look like

function getbitmats!(distmat, ANB_mat, Trail_mat, thresh)
    #basic checks
    n,m = size(distmat)
    assert(n==m)
    k=length(ANB_mat)-n
    assert(k>0)
    assert(length(ANB_mat)==length(Trail_mat))
    
    
    fill!(reinterpret(Int,ANB_mat),0)
    fill!(reinterpret(Int,Trail_mat),0)
    
    for i=1:n
        for j=1:n            
            b = (distmat[i,j] >= thresh)
            setbit!(ANB_mat, i, j, b)
        end
    end
    for i=1:n
        setbit!(ANB_mat, i, n+1, true)
    end
    for j=n+2:length(ANB_mat)
        ANB_mat[j] = ANB_mat[n+1]
    end
    
    cm = ANB_mat[n+1]
    for i=length(Trail_mat):-1:1
        cm = cm .& ANB_mat[i]
        Trail_mat[i]= cm
    end
    nothing
end

And our outer code looks like

function kcenter_entry(distmat, k, verbose=false)
   n,m=size(distmat)
   assert(n==m)
   sz = 1+ ( ((n-1)>>6) )
   ANB_mat = Vector{NTuple{sz,UInt64}}(n+k)
   Trail_mat= Vector{NTuple{sz,UInt64}}(n+k)
   return _kcenter(distmat,ANB_mat,Trail_mat,k, verbose)
end

function _kcenter(distmat, ANB_mat::Vector{NTuple{N,UInt64}}, Trail_mat, k, verbose) where {N}
   n,m=size(distmat)
   assert(n==m)
   best::Int = typemax(Int)
   output_tape = Int[]
   best_tape = Int[]
   input_tape = Bool[]
   no_miss=ntuple(i->reinterpret(UInt64,-1), N)::NTuple{N,UInt64}
   getbitmats!(distmat, ANB_mat, Trail_mat, best)
   while true
        
        res = check_nxt(ANB_mat, Trail_mat, output_tape, input_tape, no_miss, 1,k)::Bool
        if res==false
            return (best,best_tape)
        end
        #found solution, need to reduce it.
        while output_tape[1]>n
            shift!(output_tape)
        end
        best =  rateAChoice(distmat, output_tape)
        best_tape = copy(output_tape)
        verbose && println("Found new max with $(best): $(best_tape)")
        resize!(input_tape,output_tape[1])
        fill!(input_tape, false)
        for i in output_tape
            input_tape[i]=true
        end 
        resize!(output_tape,0)
        getbitmats!(distmat, ANB_mat, Trail_mat, best)
   end
   
end

Now you see that I the bitmatrices are too long; that is because I permit “nonexisting” vertices to fill up a small dominator up to k. This should never happen in reality, so feel free to refactor this out.

Now the inner loop (as a recursion). We need to restart from an initial state (input_tape), and then be as fast as we can, only writing outputs onto the output_tape if we actually find a dominator (very rare!). So, here we go:

function check_nxt(ANB_mat, Trail_mat, output_tape, cur_miss, pos, rem) 
    #first, attempt to set the next bit
    rem == 0 &&
        return all(iszero.(cur_miss))
    pos > length(ANB_mat) && 
        error("should not happen")

    res = check_nxt(ANB_mat, Trail_mat, output_tape, cur_miss .& ANB_mat[pos], pos+1 ,rem-1)
    if res
        push!(output_tape, pos)
        return true
    elseif all(iszero.(cur_miss .& Trail_mat[pos]))
        return check_nxt(ANB_mat, Trail_mat, output_tape, cur_miss, pos+1, rem)
    else        
        return false
    end    
end

function check_nxt(ANB_mat, Trail_mat, output_tape, input_tape, cur_miss, pos, rem) 
    if length(input_tape)==0
        return check_nxt(ANB_mat, Trail_mat, output_tape, cur_miss, pos, rem)
    end
    tapeval = pop!(input_tape)    
    if tapeval
        res = check_nxt(ANB_mat, Trail_mat, output_tape, input_tape, cur_miss.& ANB_mat[pos] , pos+1, rem-1)
        if res
            push!(output_tape, pos)
            return true
        end
    end
    return check_nxt(ANB_mat, Trail_mat, output_tape, input_tape, cur_miss, pos+1, rem)  
end

Now, using your setup code:

function graphListToMatrix(L)
    n = length(L)
    G = zeros(Int, n, n)
    for (j, links) in enumerate(L)
        @inbounds for link in links
            G[link[1], j] = link[2]
        end
    end
    return G
end

#A classic Floyd Warshall algorithm.
function floydWarshall(G)
    n = size(G, 1)
    D = copy(G)
    @inbounds for k = 1:n, j = 1:(n - 1), i = (j + 1):n
        if (D[i, k] != 0) & (D[j, k] != 0)
            newDist = D[i, k] + D[j, k]

            if D[i, j] == 0
                D[i, j] = newDist
                D[j, i] = newDist
            else
                if D[i, j] > newDist
                    D[i, j] = newDist
                    D[j, i] = newDist
                end
            end
        end
    end
    return D
end

function rateAChoice(D, vec)
    n = size(D, 1)
    score = 0
    @inbounds for j = 1:n
        distanceToCenter = typemax(Int)
        for i in vec
            distanceToCenter = min(distanceToCenter, D[i, j])
        end
        score = max(score, distanceToCenter)
    end
    return score
end
L = [[(2, 30), (29, 6), (95, 31), (100, 88)],[(1, 30), (3, 46)],[(2, 46), (4, 1), (79,57), (99, 49)],[(3, 1), (5, 28), (27, 27), (33, 45), (68, 28), (73, 41)],[(4, 28), (6, 31), (7, 8), (80, 90), (96, 87)],[(5, 31), (7, 69), (63, 83)],[(5, 8), (6, 69), (8, 39), (19, 9), (35, 7), (62, 55)],[(7, 39), (9, 14), (13, 77), (49, 96), (52, 31), (53, 66)],[(8, 14), (10, 84), (48, 92), (83, 75), (99, 28)],[(9, 84), (11, 59)],[(10, 59), (12, 10), (37, 4)],[(11, 10), (13, 28), (44, 56)],[(8, 77), (12, 28), (14, 63), (42, 3), (85, 22)],[(13, 63), (15, 9), (21, 82), (27, 48), (85, 90)],[(14, 9), (16, 100), (50,31), (69, 46)],[(15, 100), (17, 98), (38, 70), (79, 96)],[(16, 98), (18, 70), (54, 45), (58, 60)],[(17, 70),(19, 94), (38, 27), (75, 100)],[(7, 9), (18, 94), (20, 30), (60, 53), (92, 70)],[(19, 30), (21, 14), (93, 35)],[(14, 82), (20, 14), (22, 87), (39, 59), (81, 34)],[(21, 87),(23, 82), (52, 11), (75, 48)],[(22, 82), (24, 55)],[(23, 55), (25, 2), (59, 63)],[(24, 2), (26, 32), (99, 19)],[(25, 32), (27, 77), (29, 20), (32, 26)],[(4, 27), (14, 48), (26, 77), (28, 95), (50, 16), (60, 86)],[(27,95), (29, 29)],[(1, 6), (26, 20), (28, 29), (30, 59)],[(29, 59), (31, 91), (56, 76), (57, 37), (70, 74), (81, 99), (95, 96)],[(30, 91), (32, 89), (63, 82), (69, 53)],[(26, 26), (31, 89), (33, 50), (99, 35)],[(4, 45), (32, 50), (34, 40),(49, 77), (73, 88), (88, 8)],[(33, 40), (35, 88), (94, 52)],[(7, 7), (34, 88), (36, 94), (60, 5), (80, 46)],[(35, 94), (37, 60)],[(11, 4), (36, 60), (38, 21), (81, 14)],[(16, 70), (18, 27), (37, 21), (39, 89), (62, 56)],[(21, 59), (38, 89), (40, 47), (76, 89)],[(39, 47), (41, 63),(54, 52)],[(40, 63), (42, 45), (44, 80), (76, 42)],[(13, 3), (41, 45), (43, 46), (45, 68), (57, 5), (85, 34)],[(42, 46), (44, 24)],[(12, 56), (41, 80), (43, 24), (45, 77)],Tuple{Int64,Int64}[(42, 68), (44, 77), (46, 60), (59, 33), (68, 19)],[(45, 60), (47, 45)],[(46, 45), (48, 50), (58, 68)],[(9, 92), (47, 50), (49, 93), (98, 53)],[(8, 96), (33, 77), (48, 93), (50, 22)],[(15, 31), (27, 16), (49, 22), (51,84)],[(50, 84), (52, 16), (64, 89)],[(8, 31), (22, 11), (51, 16), (53, 85), (60, 34), (95, 96)],[(8, 66), (52, 85), (54, 68), (74, 50)],[(17, 45), (40, 52), (53, 68), (55, 93), (59, 75), (86, 35)],[(54, 93), (56, 37), (95, 47)],[(30, 76), (55, 37), (57, 26)],[(30, 37), (42, 5), (56, 26), (58, 29)],[(17, 60), (47, 68), (57, 29), (59, 38), (73, 37)],[(24, 63), (45, 33), (54, 75), (58, 38), (60, 10)],[(19, 53), (27, 86), (35, 5), (52, 34), (59, 10), (61, 32)],[(60, 32), (62, 67), (65, 73)],[(7, 55), (38, 56), (61, 67), (63, 66), (94, 96)],[(6, 83), (31, 82), (62, 66), (64, 52)],Tuple{Int64,Int64}[(51, 89), (63, 52), (65, 19)],[(61, 73), (64, 19), (66, 39), (94, 61)],[(65, 39), (67, 12), (97, 81)],[(66, 12), (68, 86)],[(4, 28), (45, 19), (67, 86), (69, 72)],[(15, 46), (31, 53), (68, 72), (70, 73), (85, 51)],[(30, 74), (69, 73), (71, 65)],[(70, 65), (72, 2)],[(71, 2), (73, 8)],[(4, 41), (33, 88), (58, 37), (72, 8), (74, 96)],[(53, 50), (73, 96), (75, 43)],[(18, 100), (22, 48), (74, 43), (76, 39), (91, 33)],[(39, 89), (41, 42), (75, 39), (77, 61), (91, 51)],[(76, 61), (78, 90)],[(77, 90), (79, 8), (91, 32)],[(3, 57), (16, 96), (78, 8), (80, 58)],[(5, 90), (35, 46), (79, 58), (81, 91), (91, 84)],[(21, 34), (30, 99), (37, 14), (80, 91), (82, 58), (95, 54)],[(81, 58), (83, 13)],[(9, 75),(82, 13), (84, 79)],[(83, 79), (85, 59)],[(13, 22), (14, 90), (42, 34), (69, 51), (84, 59), (86, 28)],[(54, 35), (85, 28), (87, 46)],[(86, 46), (88, 24), (91, 17)],[(33, 8), (87, 24), (89, 63)],[(88, 63), (90, 81)],[(89, 81), (91, 14), (94, 63), (95, 96)],[(75, 33), (76, 51), (78, 32), (80, 84), (87, 17), (90, 14), (92, 52)],[(19, 70), (91, 52), (93, 64)],[(20, 35), (92, 64), (94, 75)],[(34, 52), (62, 96), (65, 61), (90, 63), (93, 75), (95, 71)],[(1, 31), (30, 96), (52, 96), (55, 47), (81, 54), (90, 96), (94, 71), (96, 51)],[(5, 87), (95, 51), (97, 75)],[(66, 81), (96, 75), (98, 57)],[(48, 53), (97, 57), (99, 31)],[(3, 49), (9, 28), (25, 19), (32, 35), (98, 31), (100, 49)],[(1, 88), (99, 49)]];
G =   graphListToMatrix(L);
distmat = floydWarshall(G);

we can time this beast:

using BenchmarkTools

@btime kcenter_entry(distmat, $1) 
@btime kcenter_entry(distmat, $2) 
@btime kcenter_entry(distmat, $3) 
 @btime kcenter_entry(distmat, $4) 
@time kcenter_entry(distmat,1,true)
@time kcenter_entry(distmat,2,true)
@time kcenter_entry(distmat,3,true)
@time kcenter_entry(distmat,4,true)
@time kcenter_entry(distmat,5,true)
@time kcenter_entry(distmat,6,true)

for an output of:

 100.353 μs (40 allocations: 5.41 KiB)
  335.961 μs (81 allocations: 7.55 KiB)
  2.066 ms (68 allocations: 6.95 KiB)
  38.938 ms (93 allocations: 8.38 KiB)
Found new max with 231: [1]
Found new max with 197: [3]
Found new max with 196: [4]
Found new max with 186: [5]
  0.000536 seconds (241 allocations: 12.047 KiB)
Found new max with 215: [2, 1]
Found new max with 191: [3, 1]
Found new max with 186: [5, 1]
Found new max with 184: [68, 1]
Found new max with 183: [58, 2]
Found new max with 178: [7, 3]
Found new max with 173: [11, 3]
Found new max with 172: [33, 3]
Found new max with 171: [35, 3]
Found new max with 167: [45, 3]
Found new max with 163: [59, 3]
Found new max with 162: [13, 5]
  0.001233 seconds (652 allocations: 27.141 KiB)
Found new max with 191: [3, 2, 1]
Found new max with 186: [5, 2, 1]
Found new max with 178: [7, 4, 1]
Found new max with 172: [11, 4, 1]
Found new max with 165: [19, 13, 1]
Found new max with 157: [42, 7, 1]
Found new max with 155: [11, 5, 3]
Found new max with 151: [13, 5, 3]
Found new max with 148: [42, 5, 3]
  0.004046 seconds (515 allocations: 23.297 KiB)
Found new max with 191: [4, 3, 2, 1]
Found new max with 186: [5, 3, 2, 1]
Found new max with 178: [7, 4, 3, 1]
Found new max with 158: [8, 5, 4, 1]
Found new max with 157: [9, 5, 4, 1]
Found new max with 154: [13, 8, 5, 1]
Found new max with 152: [64, 13, 7, 1]
Found new max with 147: [75, 64, 57, 1]
Found new max with 144: [41, 13, 5, 3]
Found new max with 143: [90, 13, 5, 3]
Found new max with 140: [91, 13, 6, 3]
Found new max with 138: [41, 13, 7, 3]
Found new max with 137: [75, 13, 7, 3]
Found new max with 133: [75, 65, 13, 3]
  0.048291 seconds (816 allocations: 36.000 KiB)
Found new max with 186: [5, 4, 3, 2, 1]
Found new max with 178: [7, 4, 3, 2, 1]
Found new max with 158: [8, 6, 5, 4, 1]
Found new max with 157: [9, 5, 4, 3, 1]
Found new max with 154: [13, 8, 6, 5, 1]
Found new max with 151: [68, 13, 8, 6, 1]
Found new max with 143: [98, 91, 61, 56, 1]
Found new max with 139: [87, 68, 38, 8, 1]
Found new max with 136: [75, 61, 42, 24, 1]
Found new max with 130: [97, 78, 42, 35, 1]
Found new max with 129: [99, 62, 41, 7, 3]
Found new max with 128: [97, 75, 62, 13, 3]
Found new max with 127: [78, 63, 57, 9, 5]
  0.466803 seconds (3.56 k allocations: 128.547 KiB)
Found new max with 186: [6, 5, 4, 3, 2, 1]
Found new max with 178: [7, 5, 4, 3, 2, 1]
Found new max with 158: [8, 6, 5, 4, 3, 1]
Found new max with 157: [9, 6, 5, 4, 3, 1]
Found new max with 156: [9, 7, 6, 5, 4, 1]
Found new max with 155: [11, 6, 5, 4, 3, 1]
Found new max with 153: [27, 13, 8, 7, 6, 1]
Found new max with 150: [87, 59, 24, 20, 15, 1]
Found new max with 148: [98, 81, 75, 64, 29, 1]
Found new max with 141: [94, 87, 35, 24, 18, 1]
Found new max with 132: [94, 85, 78, 60, 8, 1]
Found new max with 130: [97, 78, 42, 35, 10, 1]
Found new max with 127: [99, 78, 63, 58, 20, 1]
Found new max with 126: [91, 63, 48, 42, 37, 1]
Found new max with 124: [97, 91, 85, 60, 38, 1]
Found new max with 120: [78, 63, 13, 9, 5, 3]
Found new max with 117: [78, 63, 13, 9, 7, 3]
Found new max with 116: [91, 63, 42, 11, 9, 3]
Found new max with 115: [99, 91, 81, 63, 57, 3]
Found new max with 114: [99, 91, 81, 63, 57, 4]
Found new max with 113: [99, 91, 81, 63, 57, 45]
  8.273631 seconds (4.23 k allocations: 212.766 KiB)

Almost all of the run-time should be function-call overhead (pushing stuff on the stack and returning). A loop would have been better.

PS @sbromberger Is there a very fast way of detecting isomorphism type of tiny graphs? That would be quite handy (my graphs are minimal under the following operation: If v and w are connected, and all N(v) is a subset of N(w), then I contract v->w, i.e. remove v; I compute this reduction from a bitmatrix graph, which is quite fast for up to 1000 vertex graphs, but scales like a good polynomial anyway).


#6

Thanks Liso for improving my code!


#7

Thanks foobar_lv2 for greatly improving “my” algorithm!