# Graphs, k-centers and recursion

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)
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
``````
1 Like

Did you try profiling it?

1 Like

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)
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`.

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.

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)
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
``````

``````function graphListToMatrix(L)
n = length(L)
G = zeros(Int, n, n)
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 @anon94023334 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).

1 Like

Thanks Liso for improving my code!

Thanks foobar_lv2 for greatly improving â€śmyâ€ť algorithm!