heuristic_functions.py 45.6 KB
Newer Older
cyts's avatar
cyts committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
from operator import itemgetter
import random
import numpy as np
import networkx as nx
import time
from itertools import product
from utils import (
    graph_cost,
    graph_with_calc_edge_capacity,
    create_compl_graph_from_other,
    
)

"""  
This script generates spanning trees and iterates over the heuristic network optimization algorithms
"""
     
def optimal_from_prufer_sequence(input_graph, cost_function): 
    
    """  
    Finds global minimum-cost graph from another graph object,
    by iterating over all possible trees.
    
    Based on the method described in the appendices of Heijenen et al 2020:
    "A method for designing minimum‐cost multisource multisink network layouts"
    https://doi.org/10.1002/sys.21492
    
    Args:
        input_graph: NetworkX Graph 
        cost_function: a cost function of the capacity
        
    Returns:
        optimal networkx Graph
        cost of optimal graph        
    """

    graph = input_graph.copy()
    
    #Number of nodes
    n = len(list(graph))

    nodes_list = graph.nodes()
    #Rename nodes as integers
    dict_lab = dict(zip(range(n), nodes_list))
    graph = nx.convert_node_labels_to_integers(graph)
    
    
    cost_initial = np.inf
    #Generate prufer sequences of all possible trees
    prufer_sequences = list(product(range(n), repeat = n-2))

    #Check every possible tree for a lower cost solution
    for tree in prufer_sequences:
        candidate_graph = graph.copy()
        candidate_graph.remove_edges_from(list(candidate_graph.edges()))
        candidate_graph.add_edges_from(list(nx.from_prufer_sequence(tree).edges()))
        
        #Relabel with original node labels
        candidate_graph = nx.relabel_nodes(candidate_graph, dict_lab)
        candidate_graph = graph_with_calc_edge_capacity(candidate_graph)

        new_cost = graph_cost(cost_function, candidate_graph)

        if new_cost < cost_initial:
            incumbent_graph = candidate_graph.copy()
            cost_initial = new_cost

    return incumbent_graph, cost_initial


def local_search(input_graph, cost_function):
    
    """  
    Heuristic to provide minimum-cost tree from another tree,
    by creating and breaking cycles in the trees before calculating cost.
    Steepest-descent algorithm. Not guaranteed to provide global minima solutions.
    
    Based on the method described in Brimberg 2003
    "An oil pipeline design problem."
    http://dx.doi.org/10.1287/opre.51.2.228.12786
    
    Args:
        input_graph: NetworkX Graph 
        cost_function: a cost function of the capacity
        
    Returns:
        optimized Networkx Graph
        cost of optimized graph        
    """

    complete_graph = create_compl_graph_from_other(input_graph)
    initial_cost = graph_cost(cost_function, input_graph)

    nodes_i = list(input_graph.nodes())

    #This variable serves to stop the while loop if no better solution is found
    better_found = True

    candidate_graph = input_graph.copy()

    while better_found == True:

        initial_cost = graph_cost(cost_function, candidate_graph)
        better_found = False
        
        #Potential list of cheaper trees:
        results = []
        
        #Add current best solution
        results.append([initial_cost, candidate_graph.copy()])

        for node_iter_i in nodes_i:
            
            #Find all potential neighbours to node_i
            nodes_ij = list(complete_graph[node_iter_i])
            
            #Find all existing neighbours to node_i
            existing_neighbours = list(candidate_graph.neighbors(node_iter_i))

            #Get candiate nodes that aren't connected to node_i
            for node_neighbour in existing_neighbours:
                nodes_ij.remove(node_neighbour)
                
            #Create cycle by connecting candidate nodes
            for node_ij_iter in nodes_ij:

                candidate_graph2 = candidate_graph.copy()
                candidate_graph2.add_edge(
                    node_ij_iter,
                    node_iter_i,
                    weight=complete_graph[node_ij_iter][node_iter_i]["weight"],
                )
                edges_cycle = nx.find_cycle(candidate_graph2)

                try:
                    edges_cycle.remove((node_ij_iter, node_iter_i))
                except ValueError:
                    edges_cycle.remove((node_iter_i, node_ij_iter))
                    
                #Test removing each edge in the new cycle and calculate graph cost 

                for candidate_edge_remove in edges_cycle:

                    candidate_graph3 = candidate_graph2.copy()
                    candidate_graph3.remove_edge(*candidate_edge_remove)
                    candidate_graph3 = graph_with_calc_edge_capacity(candidate_graph3)

                    new_cost = graph_cost(cost_function, candidate_graph3)

                    if new_cost < initial_cost:
                        results.append([new_cost, candidate_graph3])
                        
        #If better solution found:
        if len(results) > 1:  
            better_found = True
            results.sort(key=lambda x: x[0])
            candidate_graph = results[0][1]
            
        #Otherwise exit whle loop (better_found is stil True) and return the incumbent solution
        else:  
            candidate_graph = results[0][1]

    return candidate_graph, graph_cost(cost_function, candidate_graph) 


def delta_change(input_graph, cost_function):
    
    """  
    Heuristic to provide minimum-cost tree from another tree,
    by creating and breaking cycles in the trees before calculating cost.
    First-descent algorithm. Not guaranteed to provide global minima solutions.
    
    Based on the method described by Rothfarb 1970
    "Optimal Design of Offshore Natural-Gas Pipeline Systems"
    https://doi.org/10.1287/opre.18.6.992
    
    with additions (highlighted here) from Andre 2013
    "Design and dimensioning of hydrogen transmission pipeline networks"
    https://doi.org/10.1016/j.ejor.2013.02.036
    
    Args:
        input_graph: NetworkX Graph 
        cost_function: a cost function of the capacity
        
    Returns:
        optimized Networkx Graph
        cost of optimized graph        
    """

    complete_graph = create_compl_graph_from_other(input_graph)


    nodes_i = list(input_graph.nodes())
    
    #Randomisation of node list as indicated in Andre 2013:
    random.shuffle(nodes_i)
    
    #This variable serves to stop the while loop if no better solution is found
    better_found = True

    candidate_graph = input_graph.copy()

    while better_found == True:
        better_found = False
205
        initial_cost = graph_cost(cost_function, candidate_graph)
cyts's avatar
cyts committed
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254

        for node_iter_i in nodes_i:
            
            #The following paragraph of code would correspond to an addition of Andre 2013 
            #In which only the closest nodes to reconnect are considered. 
            #The value of 5 closest nodes is set here. 
            #We choose not to implement this addition
            #Here the weight property of the graph is simply distance
        
            #sorted_dist = [(x[0],x[1]['weight']) for x in dict(complete_graph[node_iter_i]).items()]
            #sorted_dist.sort(key=itemgetter(1))
            #nodes_ij = [x[0]  for x in  sorted_dist][:5]
             
            #Find all potential neighbours to node_i
            nodes_ij = list(complete_graph[node_iter_i])            
            #Find all existing neighbours to node_i
            existing_neighbours = list(candidate_graph.neighbors(node_iter_i))
            
            #Get candiate nodes that aren't connected to node_i
            for node_neighbour in existing_neighbours:
                nodes_ij.remove(node_neighbour)
                
            #Create cycle by connecting candidate nodes
            for node_ij_iter in nodes_ij:

                candidate_graph2 = candidate_graph.copy()
                candidate_graph2.add_edge(
                    node_ij_iter,
                    node_iter_i,
                    weight=complete_graph[node_ij_iter][node_iter_i]["weight"],
                )
                edges_cycle = nx.find_cycle(candidate_graph2)

                try:
                    edges_cycle.remove((node_ij_iter, node_iter_i))
                except ValueError:
                    edges_cycle.remove((node_iter_i, node_ij_iter))
                    
                #Test removing each edge in the new cycle and calculate graph cost
                for candidate_edge_remove in edges_cycle:

                    candidate_graph3 = candidate_graph2.copy()
                    candidate_graph3.remove_edge(*candidate_edge_remove)
                    candidate_graph3 = graph_with_calc_edge_capacity(candidate_graph3)

                    new_cost = graph_cost(cost_function, candidate_graph3)
                    
                    #As soon as better solution is found break all loops and restart:
                    if new_cost < initial_cost:
255
                        
cyts's avatar
cyts committed
256
257
258
259
                        candidate_graph = candidate_graph3.copy()
                        initial_cost = new_cost
                        better_found = True
                        break
260
                    
cyts's avatar
cyts committed
261
262
263
264
265
266
267
                if better_found == True:
                    break
            if better_found == True:
                break

    return candidate_graph, graph_cost(cost_function, candidate_graph)

268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345


def delta_change_recursive_under(input_graph, cost_function,recursion_level):
    
    complete_graph = create_compl_graph_from_other(input_graph)
    nodes_i = list(input_graph.nodes())

    candidate_graph = input_graph.copy()

    better_found = False
    
    initial_cost = graph_cost(cost_function, candidate_graph)

    for node_iter_i in nodes_i:
        
        #Find all potential neighbours to node_i
        nodes_ij = list(complete_graph[node_iter_i])
        #Get candiate nodes that aren't connected to node_i
        
        existing_neighbours = list(candidate_graph.neighbors(node_iter_i))
        
        for node_neighbour in existing_neighbours:
            nodes_ij.remove(node_neighbour)
            
        #Create cycle by connecting candidate nodes
        for node_ij_iter in nodes_ij:
            
            candidate_graph2 = candidate_graph.copy()
            candidate_graph2.add_edge(
                node_ij_iter,
                node_iter_i,
                weight=complete_graph[node_ij_iter][node_iter_i]["weight"],
            )
            edges_cycle = nx.find_cycle(candidate_graph2)

            try:
                edges_cycle.remove((node_ij_iter, node_iter_i))
            except ValueError:
                edges_cycle.remove((node_iter_i, node_ij_iter))
            
            #Test removing each edge in the new cycle and calculate graph cost
            for candidate_edge_remove in edges_cycle:

                candidate_graph3 = candidate_graph2.copy()
                candidate_graph3.remove_edge(*candidate_edge_remove)
                candidate_graph3 = graph_with_calc_edge_capacity(candidate_graph3)

                new_cost = graph_cost(cost_function, candidate_graph3)
                #As soon as better solution is found break all loops and restart:
               
                if new_cost < initial_cost:

                    candidate_graph = candidate_graph3.copy()
                    initial_cost = new_cost
                    better_found = True

                    break
                
                if recursion_level >1:
                    
                    nested_candidate,nested_cost = delta_change_recursive(candidate_graph3, cost_function,recursion_level-1)
                    
                    if nested_cost < initial_cost:

                        candidate_graph = nested_candidate.copy()
                        initial_cost = nested_cost
                        better_found = True
                        
                        break

            if better_found == True:
                break
        if better_found == True:
            break

    return candidate_graph, graph_cost(cost_function, candidate_graph)
    
def delta_change_recursive(input_graph, cost_function,recursion_level = 2):
cyts's avatar
cyts committed
346
347
348
349
350
351
    
    """  
    Heuristic to provide minimum-cost tree from another tree,
    by creating and breaking cycles in the trees before calculating cost.
    First-descent algorithm. Not guaranteed to provide global minima solutions.
    
352
    This algorithm is a recursive version of a simpler algorithm "delta_change".
cyts's avatar
cyts committed
353
    In this algorithm the delta_change is performed,
354
355
    and on each intermediate candidate tree, a nested delta_change is performed, and so on...
    allowing jumps into the neigbourhood solution space of the recursion level in one step.
cyts's avatar
cyts committed
356
357
    In this implementation we remove the Andre's (2013) additions the delta_change algorithm (and references to them),
    although including them (notably only taking closest candidate nodes) could reduce calculation time.
358
359
360
    
    Note that calcaulation time becomes usually becomes impratical for input recusion levels above 3 and node numbers above 10.
    
cyts's avatar
cyts committed
361
362
363
364
365
366
367
368
369
370
371
372

    Args:
        input_graph: NetworkX Graph 
        cost_function: a cost function of the capacity
        
    Returns:
        optimized Networkx Graph
        cost of optimized graph        
    """    
    
    complete_graph = create_compl_graph_from_other(input_graph)
    nodes_i = list(input_graph.nodes())
373
374
375
376
    
    #Randomisation of node list as indicated in Andre 2013:
    random.shuffle(nodes_i)
    
cyts's avatar
cyts committed
377
378
379
380
381
382
383
384
385
386
    #This variable serves to stop the while loop if no better solution is found
    better_found = True

    candidate_graph = input_graph.copy()

    while better_found == True:
        better_found = False
        initial_cost = graph_cost(cost_function, candidate_graph)

        for node_iter_i in nodes_i:
387
            #if node_iter
cyts's avatar
cyts committed
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
            #Find all potential neighbours to node_i
            nodes_ij = list(complete_graph[node_iter_i])
            #Get candiate nodes that aren't connected to node_i
            
            existing_neighbours = list(candidate_graph.neighbors(node_iter_i))
            
            for node_neighbour in existing_neighbours:
                nodes_ij.remove(node_neighbour)
                
            #Create cycle by connecting candidate nodes
            for node_ij_iter in nodes_ij:

                candidate_graph2 = candidate_graph.copy()
                candidate_graph2.add_edge(
                    node_ij_iter,
                    node_iter_i,
                    weight=complete_graph[node_ij_iter][node_iter_i]["weight"],
                )
                edges_cycle = nx.find_cycle(candidate_graph2)

                try:
                    edges_cycle.remove((node_ij_iter, node_iter_i))
                except ValueError:
                    edges_cycle.remove((node_iter_i, node_ij_iter))
                
                #Test removing each edge in the new cycle and calculate graph cost
                for candidate_edge_remove in edges_cycle:

                    candidate_graph3 = candidate_graph2.copy()
                    candidate_graph3.remove_edge(*candidate_edge_remove)
                    candidate_graph3 = graph_with_calc_edge_capacity(candidate_graph3)

                    new_cost = graph_cost(cost_function, candidate_graph3)
                    #As soon as better solution is found break all loops and restart:
                   
                    if new_cost < initial_cost:

                        candidate_graph = candidate_graph3.copy()
                        initial_cost = new_cost
                        better_found = True

                        break
                    
431
                    if recursion_level >1:
cyts's avatar
cyts committed
432
                        
433
434
435
436
437
438
439
440
                        nested_candidate,nested_cost = delta_change_recursive_under(candidate_graph3, cost_function,recursion_level-1)
                        
                        if nested_cost < initial_cost:
    
                            candidate_graph = nested_candidate.copy()
                            initial_cost = nested_cost
                            better_found = True
                            
cyts's avatar
cyts committed
441
                            break
442

cyts's avatar
cyts committed
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
                if better_found == True:
                    break
            if better_found == True:
                break

    return candidate_graph, graph_cost(cost_function, candidate_graph)

def edge_turn(input_graph, cost_function):
    
    """  
    Heuristic to provide minimum-cost tree from another tree,
    by breaking the tree in two,
    and reconnecting the two connected components elswehere in the trees before calculating cost.
    
    Steepest-descent algorithm. Not guaranteed to provide global minima solutions.
    
    Based on the edge-turn method described in Heijenen et al 2019:
    "A method for designing minimum‐cost multisource multisink network layouts"
    https://doi.org/10.1002/sys.21492
    Args:
        input_graph: NetworkX Graph 
        cost_function: a cost function of the capacity
        
    Returns:
        optimized Networkx Graph
        cost of optimized graph        
    """
    
    #This variable serves to stop the while loop if no better solution is found
    better_found = True
    candidate_graph = input_graph.copy()

    while better_found == True:
        better_found = False
        initial_cost = graph_cost(cost_function, candidate_graph)
        
        #Potential list of cheaper trees:
        results = []
        
        #Add current best solution        
        results.append([initial_cost, candidate_graph.copy()])

        #Test removing edges and reconnecting elswehere
        for edge in candidate_graph.edges:
      
            candidate_graph2 = candidate_graph.copy()
            candidate_graph2.remove_edge(*edge)
            
            #For each removed edge, all nodes of either connected component are re-attached to the node in the oher connected compoenent from the removed edge
            # the following codeblock in then repeated only twice in the for loop: once for each connected component
            
            #Iterates over both connected components:
            for component in list(nx.connected_components(candidate_graph2)):
                
                connected = list(component)
    
                #Finding the terminal in the removed edge that is not in the connected component considered
                if edge[0] not in connected:
                    term = edge[0]
                else:
                    term = edge[1]
    
                for comp in connected:
                    candidate_graph_test = candidate_graph2.copy()
                    candidate_graph_test.add_edge(term, comp)
                    candidate_graph_test = graph_with_calc_edge_capacity(candidate_graph_test)
                    new_cost = graph_cost(cost_function, candidate_graph_test)
    
                    if new_cost < initial_cost:
                        results.append([new_cost, candidate_graph_test])
                            
        #If better solution found:
        if len(results) > 1:
            better_found = True

            results.sort(key=lambda x: x[0])
            candidate_graph = results[0][1]

        #Otherwise exit whle loop (check is stil True) and return the incumbent solution
        else:
            candidate_graph = results[0][1]
        input_graph = candidate_graph.copy()

    return input_graph, graph_cost(cost_function, input_graph)

528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595

def edge_turn_recursive_under(input_graph, cost_function,recursion_level):
    
        
    #This variable serves to stop the while loop if no better solution is found
    candidate_graph = input_graph.copy()

    better_found = False
    
    initial_cost = graph_cost(cost_function, candidate_graph)
    
    #Potential list of cheaper trees:
    results = []
    
    #Add current best solution        
    results.append([initial_cost, candidate_graph.copy()])
    
    #Iterates over the edges to remove
    for edge in candidate_graph.edges:
        
        candidate_graph2 = candidate_graph.copy()
        candidate_graph2 = graph_with_calc_edge_capacity(candidate_graph2)
        candidate_graph2.remove_edge(*edge)
        
        #For each removed edge, all nodes of either connected component are re-attached 
        #to the node in the oher connected compoenent from the removed edge.
        #The following codeblock in then repeated only twice in the for loop: once for each connected component
        
        #Iterates over both connected components:
        for component in list(nx.connected_components(candidate_graph2)):
            connected = list(component)
            
            #Finding the terminal in the removed edge that is not in the connected component considered
            if edge[0] not in connected:
                term = edge[0]
            else:
                term = edge[1]

            for comp in connected:
                candidate_graph_test = candidate_graph2.copy()
                candidate_graph_test.add_edge(term, comp)
                candidate_graph_test = graph_with_calc_edge_capacity(candidate_graph_test)
                new_cost = graph_cost(cost_function, candidate_graph_test)

                if new_cost < initial_cost:
                    results.append([new_cost, candidate_graph_test])
                    
                if recursion_level >1:                    
                    nested_candidate,nested_cost = edge_turn_recursive_under(candidate_graph_test, cost_function,recursion_level-1)
                    
                    if nested_cost < initial_cost:
                        results.append([nested_cost, nested_candidate])
           
    #If better solution found:
    if len(results) > 1:
        better_found = True
        results.sort(key=lambda x: x[0])
        candidate_graph = results[0][1]
        
    #Otherwise return the incumbent solution
    else:
        candidate_graph = results[0][1]
    
    input_graph = candidate_graph.copy()
        
    return input_graph, graph_cost(cost_function, input_graph)

def edge_turn_recursive(input_graph, cost_function,recursion_level = 2):
cyts's avatar
cyts committed
596
597
598
599
600
601
    
    """  
    Heuristic to provide minimum-cost tree from another tree,
    by creating and breaking cycles in the trees before calculating cost.
    Steepest-descent algorithm. Not guaranteed to provide global minima solutions.
    
602
    This algorithm is a recursive version of a simpler algorithm "edge_turn".
cyts's avatar
cyts committed
603
    In this algorithm the edge turn algorithm is performed,
604
605
    and on each intermediate candidate tree, a nested turn algorithm is performed, and so on...
    allowing jumps into the neigbourhood solution space of the recursion level in one step
cyts's avatar
cyts committed
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
    In our implementation, the steepest descent applies to both the initial 
    and to the nested edge turn algorithm, and the chosen best tree is chosen 
    over all the initial and nested edge turns.
    
    Args:
        input_graph: NetworkX Graph 
        cost_function: a cost function of the capacity
        
    Returns:
        optimized Networkx Graph
        cost of optimized graph        
    """
        
    #This variable serves to stop the while loop if no better solution is found
    better_found = True
    candidate_graph = input_graph.copy()

    while better_found == True:
        
        better_found = False
        initial_cost = graph_cost(cost_function, candidate_graph)
        
        #Potential list of cheaper trees:
        results = []
        
        #Add current best solution        
        results.append([initial_cost, candidate_graph.copy()])
        
        #Iterates over the edges to remove
        for edge in candidate_graph.edges:
            
            candidate_graph2 = candidate_graph.copy()
            candidate_graph2 = graph_with_calc_edge_capacity(candidate_graph2)
            candidate_graph2.remove_edge(*edge)
            
            #For each removed edge, all nodes of either connected component are re-attached 
            #to the node in the oher connected compoenent from the removed edge.
            #The following codeblock in then repeated only twice in the for loop: once for each connected component
            
            #Iterates over both connected components:
            for component in list(nx.connected_components(candidate_graph2)):
                connected = list(component)
                
                #Finding the terminal in the removed edge that is not in the connected component considered
                if edge[0] not in connected:
                    term = edge[0]
                else:
                    term = edge[1]
    
                for comp in connected:
                    candidate_graph_test = candidate_graph2.copy()
                    candidate_graph_test.add_edge(term, comp)
                    candidate_graph_test = graph_with_calc_edge_capacity(candidate_graph_test)
                    new_cost = graph_cost(cost_function, candidate_graph_test)
    
                    if new_cost < initial_cost:
                        results.append([new_cost, candidate_graph_test])
                        
664
665
                    if recursion_level >1:                    
                        nested_candidate,nested_cost = edge_turn_recursive_under(candidate_graph_test, cost_function,recursion_level-1)
cyts's avatar
cyts committed
666
                        
667
668
669
                        if nested_cost < initial_cost:
                            results.append([nested_cost, nested_candidate])
               
cyts's avatar
cyts committed
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
        #If better solution found:
        if len(results) > 1:
            better_found = True
            results.sort(key=lambda x: x[0])
            candidate_graph = results[0][1]
            
        #Otherwise exit whle loop (better_found is stil False) and return the incumbent solution
        else:
            candidate_graph = results[0][1]
        
        input_graph = candidate_graph.copy()
        
    return input_graph, graph_cost(cost_function, input_graph)

def vns(input_graph, cost_function):
    
    """  
    Metaheuristic to provide minimum-cost tree from another tree,
    by exploring higher-cost solutions as ways out of local minima.
    Makes jumps to neighbourhoods by creating and breaking 
    cycles a given number of times before performing local_search on the
    cheapest found tree in the new optimization space.
    
    
    This algorithm is based on the description for VNS given by 
    Brimberg in 2003:
    
    "An oil pipeline design problem."
    http://dx.doi.org/10.1287/opre.51.2.228.12786
    
    The two parameters  "k_max" and "time_stopped" are fixed here at values of 5 and 10 times 
    the time required for an initial local_search on the graph, as set by Brimberg.
    These parameters correspond to the kmax and stopping criterion as described by Brimberg.
    
    Args:
        input_graph: NetworkX Graph 
        cost_function: a cost function of the capacity
        
    Returns:
        optimized Networkx Graph
        cost of optimized graph        
    """
    

    #Initiate initial solution from local_search
    start = time.time()
    candidate_graph, cost_orig = local_search(input_graph, cost_function)
    end = time.time()
    local_search_time = end - start
    
    #Set stoppping and kmax parameters
    k_max = 5
    time_stopped = 10*local_search_time

    complete_graph = create_compl_graph_from_other(input_graph)
    nodes_i = list(candidate_graph.nodes())
    
    #Initiate stopping timer and neighbourhood distance 
    continue_k = False
    start_full = time.time()

    #While stopping criterion not reached
    while time.time() - start_full < time_stopped:
        
        #Start k at 1 and increase as long as kmax not reached
        if continue_k == False:
            k = 1
        if continue_k == True:
            k = k + 1

        candidate_graph2 = candidate_graph.copy()

        nodes_done = []
        
        #Modify create and break a cycle in initial graph k times 
        for dummy in range(k):
            works = True
            #This while loop ensures that the node pairs to connect
            #are not the same as old ones

            while works == True:
                nodes_i_it = nodes_i.copy()
                node1_rand = random.choice(nodes_i_it)
                existing_neighbours = list(candidate_graph.neighbors(node1_rand))

                for node_neighbour in existing_neighbours:
                    nodes_i_it.remove(node_neighbour)

                node2_rand = random.choice(nodes_i_it)

                while (
                    node2_rand == node1_rand
                    or set([node1_rand, node2_rand]) in nodes_done
                ):
                    node2_rand = random.choice(nodes_i_it)
                    if node2_rand == node1_rand and [node2_rand] == nodes_i_it:
                        works = False
                        break
                if works == False:
                    works = True
                    continue
                #If new node pair is valid, connect them to create a cycle
                candidate_graph2.add_edge(
                    node1_rand,
                    node2_rand,
                    weight=complete_graph[node1_rand][node2_rand]["weight"],
                )

                edges_cycle = nx.find_cycle(candidate_graph2)
                #Dont allow the algorithm to remove the newly created edge
                try:
                    edges_cycle.remove((node1_rand, node2_rand))
                except ValueError:
                    edges_cycle.remove((node2_rand, node1_rand))

                pot_graphs = [] 
                #Find and keep the cheapest graph from the new cycle
                for candidate_edge_remove in edges_cycle:

                    candidate_graph3 = candidate_graph2.copy()
                    candidate_graph3.remove_edge(*candidate_edge_remove)
                    candidate_graph3 = graph_with_calc_edge_capacity(candidate_graph3)
                    new_cost = graph_cost(cost_function, candidate_graph3)
                    pot_graphs.append((new_cost, candidate_graph3))

                pot_graphs.sort(key=lambda x: x[0])

                candidate_graph2 = pot_graphs[0][1]

                nodes_done.append(set([node1_rand, node2_rand]))
                works = False
                
        #Do local search on the graph moved out of local minimum k times
        candidate_graph_end, new_cost = local_search(candidate_graph2, cost_function)
        
        #If a better solution found set as the incumbent solution and reinitialize k to 1
        if new_cost < cost_orig:
            candidate_graph = candidate_graph_end.copy()
            continue_k = False
            cost_orig = new_cost
            
        #Else explore further away from incumbent solution
        if k + 1 >= k_max:
            continue_k = False

    return candidate_graph, cost_orig


def tabu(input_graph, cost_function):
    """  
    Metaheuristic to provide minimum-cost tree from another tree,
    by exploring higher-cost solutions as ways out of local minima.
    Makes jumps to different solutions using by creating and breaking 
    cycles, and pursuing potentially higher-cost solutions in hope of finding
    a lower minimum in the wider optimization space.
    Once a transformation is made the reverse transformation is not allowed and
    is added to the Tabu list.
    
    
    This algorithm is based on the description for Tabu search given by 
    Brimberg in 2003:
    
    "An oil pipeline design problem."
    http://dx.doi.org/10.1287/opre.51.2.228.12786
    
    The parameter "stopped" is fixed here at a value of 10.
    This parameter corresponds to the stopping criterion as described by Brimberg.
    
    Args:
        s_ref: NetworkX Graph 
        func: a cost function of the capacity
        
    Returns:
        optimized Networkx Graph
        cost of optimized graph        .
    """
    
    start = time.time()
    
    #Initiate initial solution from local_search
    best_graph, best_cost = local_search(input_graph, cost_function)
    end = time.time()
    local_time = end - start
      
    #Set stopping parameter
    time_stopped = 10*local_time

    length_tabu = 7

    candidate_graph = best_graph.copy()

    #Complete graph for edges to add    
    complete_graph = create_compl_graph_from_other(candidate_graph)

    #Initial Tabu list and stopping timer
    tabu_list = []
    start_full = time.time()
    
    while time.time() - start_full < time_stopped:


        sets_complete = [set(x) for x in complete_graph.edges()]
        sets_candidate = [set(x) for x in candidate_graph.edges()]

        #Potential edges to add
        edges_not_there = [tuple(x) for x in sets_complete if x not in sets_candidate]
        
        #Initiate potential solutions list
        results = []
        
        #Iterate over neighbourhood
        for edge_test in edges_not_there:

            candidate_graph2 = candidate_graph.copy()

            candidate_graph2.add_edge(
                edge_test[0],
                edge_test[1],
                weight=complete_graph[edge_test[0]][edge_test[1]]["weight"],)

            edges_cycle = nx.find_cycle(candidate_graph2)
            
            #Dont allow the algorithm to remove the newly created edge

            try:
                edges_cycle.remove((edge_test[0], edge_test[1]))
            except ValueError:
                edges_cycle.remove((edge_test[1], edge_test[0]))
              
            #Break cycle and add solutions
            for edge_rem in edges_cycle:
                
                #Don't add solutions from Tabu list
                if (set(edge_test),set(edge_rem)) in tabu_list or (set(edge_rem),set(edge_test)) in tabu_list:
                    continue

                candidate_graph3 = candidate_graph2.copy()
                candidate_graph3.remove_edge(*edge_rem)
                candidate_graph3 = graph_with_calc_edge_capacity(candidate_graph3)
                new_cost = graph_cost(cost_function, candidate_graph3)
                results.append([new_cost, candidate_graph3, edge_test, edge_rem]) 
                
        #Choose current "exploratory" solution as best from neighbourhood
        results.sort(key=lambda x: x[0])
        candidate_graph = results[0][1].copy()
        
        #Add the transformation to Tabu list
        if len(tabu_list)<length_tabu:  
            tabu_list.insert(0,(set(results[0][2]), set(results[0][3])))
        
        else:
            tabu_list.insert(0,(set(results[0][2]), set(results[0][3])))
            tabu_list = tabu_list[:length_tabu]

        #If best "exploratory" solution is better than incumbent,
        # it becomes new incumbent solution
        if results[0][0] < best_cost:

            best_graph = results[0][1].copy()
            best_cost = results[0][0]

    return best_graph, best_cost

def high_valency_shuffle_local_search(input_graph, cost_function):
    
    """  
    Metaheuristic to provide minimum-cost tree from another tree,
    by switching high valency nodes with neighbours.
    
    After a local minimum is achieved with a lower-level heuristic,
    high valency nodes are identified and their edes are transferered to 
    a nearby node. Cycles are broken in case any are created.
    
    The lower-level heuristic used here is local search.
    
    This metaheuristic algorithm is original work,
    currently under review for publication.
    
    Args:
        input_graph: NetworkX Graph 
        cost_function: a cost function of the capacity
        
    Returns:
        optimized Networkx Graph
        cost of optimized graph        .
    """
    
    #Initiate initial solution from local_search 
    candidate_graph, orig_cost = local_search(input_graph, cost_function)
    
    #Complete graph is used to identify closest nodes for node shuffle    
    complete_graph = create_compl_graph_from_other(input_graph)

    better_found = True

    while better_found == True:

        better_found = False

        #Initiate results list
        results = []
        results.append([orig_cost, candidate_graph])

        #Find high valency (here chosen as <2) nodes
        candidates_high_valency = [(x[0], x[1]) for x in list(candidate_graph.degree) if x[1] > 2]

        #Iterating over initial high valency nodes
        for cand in candidates_high_valency:

            sorted_dist = [
                (x[0], x[1]["weight"]) for x in dict(complete_graph[cand[0]]).items()
            ]
            
            #Find closest neighbours to transfer edges to (here closest 4)
            sorted_dist.sort(key=itemgetter(1))
            nodes_ij = [x[0] for x in sorted_dist][:4]
            
            #Iterating over nodes to transfer edges to
            for node_ij in nodes_ij:

                candidate_graph2 = candidate_graph.copy()
                
                #Edges to transfer
                edges = list(candidate_graph2[cand[0]])

                for to_reconnect in edges:

                    if to_reconnect != node_ij:

                        candidate_graph2.add_edge(node_ij, to_reconnect)
                        candidate_graph2.remove_edge(cand[0], to_reconnect)
                        
                #After edges removal, if node is not connected, reconnect to candidate
                if len(list(candidate_graph2[cand[0]])) == 0:
                    candidate_graph2.add_edge(nodes_ij[0], cand[0])
                    
                
                #Break cycle (if found) and add solutions
                try:
                    edges_cycle = nx.find_cycle(candidate_graph2)

                    for candidate_edge_remove in edges_cycle:

                        candidate_graph3 = candidate_graph2.copy()
                        candidate_graph3.remove_edge(*candidate_edge_remove)
                        candidate_graph3 = graph_with_calc_edge_capacity(candidate_graph3)
                        #Local Search is the lower-level heuristic
                        candidate_graph4, new_cost = local_search(candidate_graph3, cost_function)

                        results.append([new_cost, candidate_graph4])

                except nx.exception.NetworkXNoCycle:

                    candidate_graph2 = graph_with_calc_edge_capacity(candidate_graph2)
                    candidate_graph3, new_cost = local_search(candidate_graph2, cost_function)
                    results.append([new_cost, candidate_graph3])

        results.sort(key=lambda x: x[0])

        if results[0][0] < orig_cost:

            orig_cost = results[0][0]
            candidate_graph = results[0][1].copy()
            better_found = True

    return candidate_graph, orig_cost


def high_valency_shuffle_edge_turn(input_graph, cost_function):
    
    """  
    Metaheuristic to provide minimum-cost tree from another tree,
    by switching high valency nodes with neighbours.
    
    After a local minimum is achieved with a lower-level heuristic,
    high valency nodes are identified and their edes are transferered to 
    a nearby node. Cycles are broken in case any are created.
    
    The lower-level heuristic used here is the edge-turn algorithm.
    
    This metaheuristic algorithm is original work,
    currently under review for publication.
    
    Args:
        input_graph: NetworkX Graph 
        cost_function: a cost function of the capacity
        
    Returns:
        optimized Networkx Graph
        cost of optimized graph        .
    """
    
    #Initiate initial solution from edge_turn
    candidate_graph, orig_cost = edge_turn(input_graph, cost_function)
    
    #Complete graph is used to identify closest nodes for node shuffle    
    complete_graph = create_compl_graph_from_other(input_graph)

    better_found = True

    while better_found == True:
        better_found = False
        
        #Initiate results list
        results = []
        results.append([orig_cost, candidate_graph])
        #Find high valency (here chosen as <2) nodes
        candidates = [(x[0], x[1]) for x in list(candidate_graph.degree) if x[1] > 2]

        #Iterating over initial high valency nodes
        for cand in candidates:
            sorted_dist = [
                (x[0], x[1]["weight"]) for x in dict(complete_graph[cand[0]]).items()
            ]
            #Find closest neighbours to transfer edges to (here closest 4)
            sorted_dist.sort(key=itemgetter(1))
            nodes_ij = [x[0] for x in sorted_dist][:4]
            
            #Iterating over nodes to transfer edges to
            for node_ij in nodes_ij:

                candidate_graph2 = candidate_graph.copy()
                
                #Edges to transfer
                edges = list(candidate_graph2[cand[0]])

                for to_reconnect in edges:

                    if to_reconnect != node_ij:

                        candidate_graph2.add_edge(node_ij, to_reconnect)
                        candidate_graph2.remove_edge(cand[0], to_reconnect)
                        
                #After edges removal, if node is not connected, reconnect to candidate
                if len(list(candidate_graph2[cand[0]])) == 0:
                    candidate_graph2.add_edge(nodes_ij[0], cand[0])
                    
                
                #Break cycle (if found) and add solutions
                try:
                    edges_cycle = nx.find_cycle(candidate_graph2)

                    for candidate_edge_remove in edges_cycle:

                        candidate_graph3 = candidate_graph2.copy()
                        candidate_graph3.remove_edge(*candidate_edge_remove)
                        candidate_graph3 = graph_with_calc_edge_capacity(candidate_graph3)
                        #The Edge turn algorithm is the lower-level heuristic
                        candidate_graph4, new_cost = edge_turn(candidate_graph3, cost_function)

                        results.append([new_cost, candidate_graph4])

                except nx.exception.NetworkXNoCycle:

                    candidate_graph2 = graph_with_calc_edge_capacity(candidate_graph2)
                    candidate_graph3, new_cost = edge_turn(candidate_graph2, cost_function)
                    results.append([new_cost, candidate_graph3])

        results.sort(key=lambda x: x[0])

        if results[0][0] < orig_cost:

            orig_cost = results[0][0]
            candidate_graph = results[0][1].copy()
            better_found = True

    return candidate_graph, orig_cost

def high_valency_shuffle_delta_change(input_graph, cost_function):
    
    """  
    Metaheuristic to provide minimum-cost tree from another tree,
    by switching high valency nodes with neighbours.
    
    After a local minimum is achieved with a lower-level heuristic,
    high valency nodes are identified and their edes are transferered to 
    a nearby node. Cycles are broken in case any are created.
    
    The lower-level heuristic used here is the delta change algorithm.
    
    This metaheuristic algorithm is original work,
    currently under review for publication.
    
    Args:
        input_graph: NetworkX Graph 
        cost_function: a cost function of the capacity
        
    Returns:
        optimized Networkx Graph
        cost of optimized graph        .
    """
    
    #Initiate initial solution from delta_change
    candidate_graph, orig_cost = delta_change(input_graph, cost_function)
    
    #Complete graph is used to identify closest nodes for node shuffle    
    complete_graph = create_compl_graph_from_other(input_graph)

    better_found = True

    while better_found == True:

        better_found = False

        #Initiate results list
        results = []

        results.append([orig_cost, candidate_graph])

        #Find high valency (here chosen as >2) nodes
        candidates = [(x[0], x[1]) for x in list(candidate_graph.degree) if x[1] > 2]

        #Iterating over initial high valency nodes
        for cand in candidates:

            sorted_dist = [
                (x[0], x[1]["weight"]) for x in dict(complete_graph[cand[0]]).items()
            ]
            
            #Find closest neighbours to transfer edges to (here closest 4)
            sorted_dist.sort(key=itemgetter(1))
            nodes_ij = [x[0] for x in sorted_dist][:4]
            
            #Iterating over nodes to transfer edges to
            for node_ij in nodes_ij:

                candidate_graph2 = candidate_graph.copy()
                
                #Edges to transfer
                edges = list(candidate_graph2[cand[0]])

                for to_reconnect in edges:

                    if to_reconnect != node_ij:

                        candidate_graph2.add_edge(node_ij, to_reconnect)
                        candidate_graph2.remove_edge(cand[0], to_reconnect)
                        
                #After edges removal, if node is not connected, reconnect to candidate
                if len(list(candidate_graph2[cand[0]])) == 0:
                    candidate_graph2.add_edge(nodes_ij[0], cand[0])
                    
                
                #Break cycle (if found) and add solutions
                try:
                    edges_cycle = nx.find_cycle(candidate_graph2)

                    for candidate_edge_remove in edges_cycle:

                        candidate_graph3 = candidate_graph2.copy()
                        candidate_graph3.remove_edge(*candidate_edge_remove)
                        candidate_graph3 = graph_with_calc_edge_capacity(candidate_graph3)
                        #The delta change algorithm is the lower-level heuristic
                        candidate_graph4, new_cost = delta_change(candidate_graph3, cost_function)

                        results.append([new_cost, candidate_graph4])

                except nx.exception.NetworkXNoCycle:

                    candidate_graph2 = graph_with_calc_edge_capacity(candidate_graph2)
                    candidate_graph3, new_cost = delta_change(candidate_graph2, cost_function)
                    results.append([new_cost, candidate_graph3])

        results.sort(key=lambda x: x[0])

        if results[0][0] < orig_cost:

            orig_cost = results[0][0]
            candidate_graph = results[0][1].copy()
            better_found = True

    return candidate_graph, orig_cost