# HG changeset patch
# User Nathann Cohen <nathann.cohen@gmail.com>
# Date 1356105749 3600
# Node ID fc33c8b2560206cfe0e2b795fb52b2cd43c0868e
# Parent 9df346e99f5e36540342bae6545ad2aa6856e593
Gromov hyperbolicity of graphs  reviewer's patch
diff git a/sage/graphs/hyperbolicity.pyx b/sage/graphs/hyperbolicity.pyx
a

b


304  304  
305  305  An elimination ordering of simplicial vertices is an elimination ordering of 
306  306  the vertices of the graphs such that the induced subgraph of their neighbors 
307   is a clique. More precisely, as long as the graph as a vertex ``u`` such 
 307  is a clique. More precisely, as long as the graph has a vertex ``u`` such 
308  308  that the induced subgraph of its neighbors is a clique, we remove ``u`` from 
309  309  the graph, add it to the elimination ordering (list of vertices), and 
310  310  repeat. This method is inspired from the decomposition of a graph by 
… 
… 

343  343  ... 
344  344  ValueError: The parameter max_degree must be > 0. 
345  345  
346   Giving a graph build from a bipartite graph plus an edge:: 
 346  Giving a graph built from a bipartite graph plus an edge:: 
347  347  
348  348  sage: G = graphs.CompleteBipartiteGraph(2,10) 
349  349  sage: G.add_edge(0,1) 
… 
… 

433  433  if j<jmax: 
434  434  _invert_cells(pairs, j, jmax) 
435  435  
436   if jmax>0: 
437   jmax = 1 
 436  jmax = 1 
438  437  
439  438  else: # This pair is at a correct position. 
440  439  j += 1 
… 
… 

442  441  # We record the position of the first pair of vertices in elim 
443  442  last_pair[0] = jmax+1 
444  443  
445   
446  444  cdef tuple __hyperbolicity__(int N, 
447  445  unsigned short ** distances, 
448  446  int D, 
… 
… 

504  502  cdef dict distr = {} 
505  503  cdef list certificate = [] 
506  504  
507   # We count the number of pairs of vertices at distance l for every l 
 505  # ==> Allocates and fills nb_pairs_of_length 
 506  # 
 507  # nb_pairs_of_length[d] is the number of pairs of vertices at distance d 
508  508  cdef uint32_t * nb_pairs_of_length = <uint32_t *>sage_calloc(D+1,sizeof(uint32_t)) 
 509  if nb_pairs_of_length == NULL: 
 510  sage_free(nb_pairs_of_length) 
 511  raise MemoryError 
 512  
509  513  for 0 <= i < N: 
510  514  for i+1 <= j < N: 
511  515  nb_pairs_of_length[ distances[i][j] ] += 1 
512  516  
513   # We organize the pairs by length in an array of pairs 
 517  # ==> Allocates pairs_of_length 
 518  # 
 519  # pairs_of_length[d] is the list of pairs of vertices at distance d 
514  520  cdef pair ** pairs_of_length = <pair **>sage_malloc(sizeof(pair *)*(D+1)) 
 521  if pairs_of_length == NULL: 
 522  sage_free(nb_pairs_of_length) 
 523  raise MemoryError 
 524  
 525  # ==> Allocates cpt_pairs 
 526  # 
 527  # (temporary variable used to fill pairs_of_length) 
515  528  cdef uint32_t * cpt_pairs = <uint32_t *>sage_calloc(D+1,sizeof(uint32_t)) 
 529  if cpt_pairs == NULL: 
 530  sage_free(nb_pairs_of_length) 
 531  sage_free(pairs_of_length) 
 532  raise MemoryError 
 533  
 534  # ==> Allocates pairs_of_length[d] for all d 
516  535  for 1 <= i <= D: 
517  536  pairs_of_length[i] = <pair *>sage_malloc(sizeof(pair)*nb_pairs_of_length[i]) 
518   # cpt_pairs[i] = 0 
 537  
 538  if nb_pairs_of_length[i] > 0 and pairs_of_length[i] == NULL: 
 539  while i>0: 
 540  i = 1 
 541  sage_free(pairs_of_length[i]) 
 542  sage_free(nb_pairs_of_length) 
 543  sage_free(pairs_of_length) 
 544  sage_free(cpt_pairs) 
 545  raise MemoryError 
 546  
 547  # ==> Fills pairs_of_length[d] for all d 
519  548  for 0 <= i < N: 
520  549  for i+1 <= j < N: 
521  550  l = distances[i][j] 
522  551  pairs_of_length[ l ][ cpt_pairs[ l ] ].s = i 
523  552  pairs_of_length[ l ][ cpt_pairs[ l ] ].t = j 
524  553  cpt_pairs[ l ] += 1 
 554  
525  555  sage_free(cpt_pairs) 
526  556  
527  557  if verbose: 
528  558  print "Current 2 connected component has %d vertices and diameter %d" %(N,D) 
529  559  print "Paths length distribution:", [ (l, nb_pairs_of_length[l]) for l in range(1, D+1) ] 
530  560  
 561  # ==> Allocates last_pair 
 562  # 
531  563  # We improve the ordering of the pairs according the elimination 
532  564  # ordering. We store in last_pair[l] the index of the last pair of length l 
533  565  # to consider when the lower bound on the hyperbolicity >=1. 
534  566  cdef uint32_t * last_pair = <uint32_t *>sage_malloc(sizeof(uint32_t)*(D+1)) 
 567  if last_pair == NULL: 
 568  for 1 <= i <= D: 
 569  sage_free(pairs_of_length[i]) 
 570  sage_free(nb_pairs_of_length) 
 571  sage_free(pairs_of_length) 
 572  sage_free(cpt_pairs) 
 573  raise MemoryError 
 574  
535  575  for 1 <= l <= D: 
536  576  _order_pairs_according_elimination_ordering(elim, N, pairs_of_length[l], nb_pairs_of_length[l], last_pair+l) 
537  577  
… 
… 

556  596  for S1, l1, l2 in triples: 
557  597  
558  598  # If we cannot improve further, we stop 
 599  # 
 600  # See the module's documentation for an proof that this cut is 
 601  # valid. Remember that the triples are sorted in a specific order. 
559  602  if l2 <= h: 
560  603  h_UB = h 
561  604  break 
… 
… 

564  607  h_UB = l2 
565  608  
566  609  if verbose: 
567   print "New upper bound:",ZZ(l2)/2 
 610  print "New upper bound:",ZZ(h_UB)/2 
568  611  
569  612  # Termination if required approximation is found 
570  613  if certificate and ( (h_UB <= h*approximation_factor) or (h_UBh <= additive_gap) ): 
571  614  break 
572  615  
573  616  pairs_of_length_l1 = pairs_of_length[l1] 
574   nb_pairs_of_length_l1 = last_pair[l1] if h>=1 else nb_pairs_of_length[l1] 
 617  
 618  # We only need test the pairs after last_pair[l1] if we do not have a 
 619  # proof that h>=2. 
 620  nb_pairs_of_length_l1 = last_pair[l1] if h>=2 else nb_pairs_of_length[l1] 
575  621  x = 0 
576  622  while x < nb_pairs_of_length_l1: 
577  623  a = pairs_of_length_l1[x].s 
578  624  b = pairs_of_length_l1[x].t 
579  625  
 626  # If we cannot improve further, we stop 
 627  # 
 628  # See the module's documentation for an proof that this cut is 
 629  # valid. 
580  630  if l2 <= h: 
581   # We cut current exploration if lower bound cannot be improved 
582  631  break 
 632  
 633  # We do not want to test pairs of pairs twice if l1 == l2 
583  634  elif l1 == l2: 
584  635  y = x+1 
585  636  else: 
586  637  y = 0 
587  638  
588  639  pairs_of_length_l2 = pairs_of_length[l2] 
589   nb_pairs_of_length_l2 = last_pair[l2] if h>=1 else nb_pairs_of_length[l2] 
 640  
 641  # We only need test the pairs after last_pair[l2] if we do not have 
 642  # a proof that h>=2. 
 643  nb_pairs_of_length_l2 = last_pair[l2] if h>=2 else nb_pairs_of_length[l2] 
590  644  while y < nb_pairs_of_length_l2: 
591  645  c = pairs_of_length_l2[y].s 
592  646  d = pairs_of_length_l2[y].t 
… 
… 

614  668  # search space 
615  669  h = hh 
616  670  certificate = [a, b, c, d] 
617   if h>=1: 
 671  if h>=2: 
618  672  nb_pairs_of_length_l2 = last_pair[l2] 
619  673  nb_pairs_of_length_l1 = last_pair[l1] 
620  674  if x >= nb_pairs_of_length_l1: 
… 
… 

626  680  # We go for the next pair cd 
627  681  y += 1 
628  682  # We cut current exploration if we know we can not improve lower bound 
 683  # 
 684  # See the module's documentation for an proof that this cut is 
 685  # valid. 
629  686  if l2 <= h: 
630  687  h_UB = h 
631  688  break 
 689  
632  690  # We go for the next pair ab 
633  691  x += 1 
634  692  
635  693  
636   # We now release the memory 
 694  # We now free the memory 
637  695  sage_free(nb_pairs_of_length) 
638  696  for 1 <= i <= D: 
639  697  sage_free(pairs_of_length[i]) 
… 
… 

646  704  else: 
647  705  return (h, certificate, h_UB) 
648  706  
649   
650   
651  707  def hyperbolicity(G, algorithm='cuts', approximation_factor=1.0, additive_gap=0, verbose = False): 
652  708  r""" 
653  709  Return the hyperbolicity of the graph or an approximation of this value. 
… 
… 

660  716  hyperbolicity of the graph is the maximum over all possible 4tuples `(a, b, 
661  717  c, d)` divided by 2. The worst case time complexity is in `O( n^4 )`. 
662  718  
 719  See the documentation of :mod:`sage.graphs.hyperbolicity` for more 
 720  information. 
 721  
663  722  INPUT: 
664  723  
665  724   ``G``  a Graph 
… 
… 

685  744  soon as the ratio between the upper bound and the best found solution is 
686  745  less than the approximation factor. When the approximation factor is 1.0, 
687  746  the problem is solved optimaly. This parameter is used only when the 
688   chosen algorithm is ``'cuts'``. 
 747  chosen algorithm is ``'cuts'`` or ``'cuts+'``. 
689  748  
690  749   ``additive_gap``  (default: 0.0) When sets to a positive number, the 
691  750  function stop computations as soon as the difference between the upper 
692  751  bound and the best found solution is less than additive gap. When the gap 
693  752  is 0.0, the problem is solved optimaly. This parameter is used only when 
694   the chosen algorithm is ``'cuts'``. 
 753  the chosen algorithm is ``'cuts'`` or ``'cuts+'``. 
695  754  
696  755   ``verbose``  (default: ``False``) is a boolean set to True to display 
697  756  some information during execution: new upper and lower bounds, etc. 
… 
… 

822  881  elif G.num_verts() == G.num_edges() + 1: 
823  882  # G is a tree 
824  883  # Any set of 4 vertices is a valid certificate 
825   return 0, [G.vertices()[x] for x in range(4)], 0 
 884  return 0, G.vertices()[:4], 0 
826  885  
827  886  elif G.is_clique(): 
828  887  # Any set of 4 vertices is a valid certificate 
829   return 0, [G.vertices()[z] for z in xrange(4)], 0 
830   
 888  return 0, G.vertices()[:4], 0 
831  889  
832  890  cdef unsigned short * _distances_ 
833  891  cdef unsigned short ** distances 
… 
… 

874  932  # in the range [0..N1]. 
875  933  mymap = H.relabel( return_map=True ) 
876  934  
877   # We compute the distances and store the results in a 2D array 
878   # Although it seems irrelevant to use a 2D array instead of a 1D 
879   # array, it seems to be faster in practice. We also compute the 
880   # diameter. 
 935  # We compute the distances and store the results in a 2D array, and the diameter 
881  936  _distances_ = c_distances_all_pairs(H) 
882  937  distances = <unsigned short **>sage_malloc(sizeof(unsigned short *)*N) 
 938  if distances == NULL: 
 939  sage_free(_distances_) 
 940  raise MemoryError 
 941  
883  942  D = 0 
884  943  for 0 <= i < N: 
885  944  distances[i] = _distances_+i*N 
… 
… 

962  1021  cdef int i 
963  1022  
964  1023  cdef uint64_t * hdistr = <uint64_t *>sage_calloc(N+1,sizeof(uint64_t)) 
 1024  if hdistr == NULL: 
 1025  raise MemoryError 
965  1026  
966  1027  # We now compute the hyperbolicity of each 4tuple 
967  1028  cdef int a, b, c, d 
… 
… 

973  1034  
974  1035  # We prepare the dictionnary of hyperbolicity distribution to return 
975  1036  Nchoose4 = binomial(N,4) 
976   cdef dict hdict = dict( [ (ZZ(i)/2, ZZ(hdistr[i])/Nchoose4) for 0 <= i <= N if hdistr[i] > 0 ] ) 
 1037  cdef dict hdict = {ZZ(i)/2: (ZZ(hdistr[i])/Nchoose4) for 0 <= i <= N if hdistr[i] > 0} 
977  1038  
978  1039  sage_free(hdistr) 
979  1040  
… 
… 

1017  1078  cdef int i, a, b, c, d 
1018  1079  cdef uint64_t j 
1019  1080  
 1081  if N < 4: 
 1082  raise ValueError("N must be at least 4") 
 1083  
1020  1084  # We initialize the table of hyperbolicity. We use an array of unsigned long 
1021  1085  # int instead of a dictionnary since it is much faster. 
1022  1086  cdef uint64_t * hdistr = <uint64_t *>sage_calloc(N+1,sizeof(uint64_t)) 
 1087  if hdistr == NULL: 
 1088  raise MemoryError 
1023  1089  
1024  1090  # We now compute the hyperbolicity of each quadruple 
1025  1091  for 0 <= j < sampling_size: 
… 
… 

1135  1201  H = G.relabel( inplace = False ) 
1136  1202  _distances_ = c_distances_all_pairs(H) 
1137  1203  distances = <unsigned short **>sage_malloc(sizeof(unsigned short *)*N) 
 1204  if distances == NULL: 
 1205  sage_free(_distances_) 
 1206  raise MemoryError 
 1207  
1138  1208  for 0 <= i < N: 
1139  1209  distances[i] = _distances_+i*N 
1140  1210  