Ticket #13808: trac_13808-rev.patch

File trac_13808-rev.patch, 11.8 KB (added by ncohen, 8 years ago)
  • sage/graphs/hyperbolicity.pyx

    # 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  
    304304
    305305    An elimination ordering of simplicial vertices is an elimination ordering of
    306306    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
    308308    that the induced subgraph of its neighbors is a clique, we remove ``u`` from
    309309    the graph, add it to the elimination ordering (list of vertices), and
    310310    repeat. This method is inspired from the decomposition of a graph by
     
    343343        ...
    344344        ValueError: The parameter max_degree must be > 0.
    345345
    346     Giving a graph build from a bipartite graph plus an edge::
     346    Giving a graph built from a bipartite graph plus an edge::
    347347
    348348        sage: G = graphs.CompleteBipartiteGraph(2,10)
    349349        sage: G.add_edge(0,1)
     
    433433                if j<jmax:
    434434                    _invert_cells(pairs, j, jmax)
    435435
    436                 if jmax>0:
    437                     jmax -= 1
     436                jmax -= 1
    438437
    439438            else: # This pair is at a correct position.
    440439                j += 1
     
    442441        # We record the position of the first pair of vertices in elim
    443442        last_pair[0] = jmax+1
    444443
    445 
    446444cdef tuple __hyperbolicity__(int N,
    447445                             unsigned short **  distances,
    448446                             int D,
     
    504502    cdef dict distr = {}
    505503    cdef list certificate = []
    506504
    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
    508508    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
    509513    for 0 <= i < N:
    510514        for i+1 <= j < N:
    511515            nb_pairs_of_length[ distances[i][j] ] += 1
    512516
    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
    514520    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)
    515528    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
    516535    for 1 <= i <= D:
    517536        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
    519548    for 0 <= i < N:
    520549        for i+1 <= j < N:
    521550            l = distances[i][j]
    522551            pairs_of_length[ l ][ cpt_pairs[ l ] ].s = i
    523552            pairs_of_length[ l ][ cpt_pairs[ l ] ].t = j
    524553            cpt_pairs[ l ] += 1
     554
    525555    sage_free(cpt_pairs)
    526556
    527557    if verbose:
    528558        print "Current 2 connected component has %d vertices and diameter %d" %(N,D)
    529559        print "Paths length distribution:", [ (l, nb_pairs_of_length[l]) for l in range(1, D+1) ]
    530560
     561    # ==> Allocates last_pair
     562    #
    531563    # We improve the ordering of the pairs according the elimination
    532564    # ordering. We store in last_pair[l] the index of the last pair of length l
    533565    # to consider when the lower bound on the hyperbolicity >=1.
    534566    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
    535575    for 1 <= l <= D:
    536576        _order_pairs_according_elimination_ordering(elim, N, pairs_of_length[l], nb_pairs_of_length[l], last_pair+l)
    537577
     
    556596    for S1, l1, l2 in triples:
    557597
    558598        # 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.
    559602        if l2 <= h:
    560603            h_UB = h
    561604            break
     
    564607            h_UB = l2
    565608
    566609            if verbose:
    567                 print "New upper bound:",ZZ(l2)/2
     610                print "New upper bound:",ZZ(h_UB)/2
    568611
    569612            # Termination if required approximation is found
    570613            if certificate and ( (h_UB <= h*approximation_factor) or (h_UB-h <= additive_gap) ):
    571614                break
    572615
    573616        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]
    575621        x = 0
    576622        while x < nb_pairs_of_length_l1:
    577623            a = pairs_of_length_l1[x].s
    578624            b = pairs_of_length_l1[x].t
    579625
     626            # If we cannot improve further, we stop
     627            #
     628            # See the module's documentation for an proof that this cut is
     629            # valid.
    580630            if l2 <= h:
    581                 # We cut current exploration if lower bound cannot be improved
    582631                break
     632
     633            # We do not want to test pairs of pairs twice if l1 == l2
    583634            elif l1 == l2:
    584635                y = x+1
    585636            else:
    586637                y = 0
    587638
    588639            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]
    590644            while y < nb_pairs_of_length_l2:
    591645                c = pairs_of_length_l2[y].s
    592646                d = pairs_of_length_l2[y].t
     
    614668                    # search space
    615669                    h = hh
    616670                    certificate = [a, b, c, d]
    617                     if h>=1:
     671                    if h>=2:
    618672                        nb_pairs_of_length_l2 = last_pair[l2]
    619673                        nb_pairs_of_length_l1 = last_pair[l1]
    620674                        if x >= nb_pairs_of_length_l1:
     
    626680                # We go for the next pair c-d
    627681                y += 1
    628682                # 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.
    629686                if l2 <= h:
    630687                    h_UB = h
    631688                    break
     689
    632690            # We go for the next pair a-b
    633691            x += 1
    634692
    635693
    636     # We now release the memory
     694    # We now free the memory
    637695    sage_free(nb_pairs_of_length)
    638696    for 1 <= i <= D:
    639697        sage_free(pairs_of_length[i])
     
    646704    else:
    647705        return (h, certificate, h_UB)
    648706
    649 
    650 
    651707def hyperbolicity(G, algorithm='cuts', approximation_factor=1.0, additive_gap=0, verbose = False):
    652708    r"""
    653709    Return the hyperbolicity of the graph or an approximation of this value.
     
    660716    hyperbolicity of the graph is the maximum over all possible 4-tuples `(a, b,
    661717    c, d)` divided by 2. The worst case time complexity is in `O( n^4 )`.
    662718
     719    See the documentation of :mod:`sage.graphs.hyperbolicity` for more
     720    information.
     721
    663722    INPUT:
    664723
    665724    - ``G`` -- a Graph
     
    685744      soon as the ratio between the upper bound and the best found solution is
    686745      less than the approximation factor. When the approximation factor is 1.0,
    687746      the problem is solved optimaly. This parameter is used only when the
    688       chosen algorithm is ``'cuts'``.
     747      chosen algorithm is ``'cuts'`` or ``'cuts+'``.
    689748
    690749    - ``additive_gap`` -- (default: 0.0) When sets to a positive number, the
    691750      function stop computations as soon as the difference between the upper
    692751      bound and the best found solution is less than additive gap. When the gap
    693752      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+'``.
    695754
    696755    - ``verbose`` -- (default: ``False``) is a boolean set to True to display
    697756      some information during execution: new upper and lower bounds, etc.
     
    822881    elif G.num_verts() == G.num_edges() + 1:
    823882        # G is a tree
    824883        # 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
    826885
    827886    elif G.is_clique():
    828887        # 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
    831889
    832890    cdef unsigned short * _distances_
    833891    cdef unsigned short ** distances
     
    874932            # in the range [0..N-1].
    875933            mymap = H.relabel( return_map=True )
    876934
    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
    881936            _distances_ = c_distances_all_pairs(H)
    882937            distances = <unsigned short **>sage_malloc(sizeof(unsigned short *)*N)
     938            if distances == NULL:
     939                sage_free(_distances_)
     940                raise MemoryError
     941
    883942            D = 0
    884943            for 0 <= i < N:
    885944                distances[i] = _distances_+i*N
     
    9621021    cdef int i
    9631022
    9641023    cdef uint64_t * hdistr = <uint64_t *>sage_calloc(N+1,sizeof(uint64_t))
     1024    if hdistr == NULL:
     1025        raise MemoryError
    9651026
    9661027    # We now compute the hyperbolicity of each 4-tuple
    9671028    cdef int a, b, c, d
     
    9731034
    9741035    # We prepare the dictionnary of hyperbolicity distribution to return
    9751036    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}
    9771038
    9781039    sage_free(hdistr)
    9791040
     
    10171078    cdef int i, a, b, c, d
    10181079    cdef uint64_t j
    10191080
     1081    if N < 4:
     1082        raise ValueError("N must be at least 4")
     1083
    10201084    # We initialize the table of hyperbolicity. We use an array of unsigned long
    10211085    # int instead of a dictionnary since it is much faster.
    10221086    cdef uint64_t * hdistr = <uint64_t *>sage_calloc(N+1,sizeof(uint64_t))
     1087    if hdistr == NULL:
     1088        raise MemoryError
    10231089
    10241090    # We now compute the hyperbolicity of each quadruple
    10251091    for 0 <= j < sampling_size:
     
    11351201    H = G.relabel( inplace = False )
    11361202    _distances_ = c_distances_all_pairs(H)
    11371203    distances = <unsigned short **>sage_malloc(sizeof(unsigned short *)*N)
     1204    if distances == NULL:
     1205        sage_free(_distances_)
     1206        raise MemoryError
     1207
    11381208    for 0 <= i < N:
    11391209        distances[i] = _distances_+i*N
    11401210