Ticket #4859: 11113.patch

File 11113.patch, 12.1 KB (added by dgordon, 11 years ago)
  • sage/combinat/designs/all.py

    # HG changeset patch
    # User Daniel Gordon <gordon@ccrwest.org>
    # Date 1230048660 28800
    # Node ID a65f968c0f2e9eff0022377d8d5cb9ba200f2017
    # Parent  5be1d5ad833948143908b3a989524c960e8a6f8c
    new code to handle covering designs, using incidence structure machinery
    
    diff -r 5be1d5ad8339 -r a65f968c0f2e sage/combinat/designs/all.py
    a b  
    55
    66from incidence_structures import (IncidenceStructure,
    77                          IncidenceStructureFromMatrix)
     8
     9from covering_design import (CoveringDesign,
     10                             schonheim,
     11                             trivial_covering_design,
     12                             best_known_covering_design_www)
  • new file sage/combinat/designs/covering_design.py

    diff -r 5be1d5ad8339 -r a65f968c0f2e sage/combinat/designs/covering_design.py
    - +  
     1r"""
     2Covering designs: coverings of $t$-element subsets of a $v$-set by $k$-sets
     3
     4A $(v,k,t)$ covering design $C$ is an incidence structure consisting of a
     5set of points $P$ of order $v$, and a set of blocks $B$, where each
     6block contains $k$ points of $P$.  Every $t$-element subset of $P$
     7must be contained in at least one block.
     8
     9If every $t$-set is contained in exactly one block of $C$, then we
     10have a block design.  Following the block design implementation, the
     11standard representation of a covering design uses $P = [0,1,..., v-1]$.
     12
     13There is an online database of the best known covering designs, the La
     14Jolla Covering Repository (LJCR), at [1].  This is maintained with an SQL
     15database, also in Sage, but since it changes frequently, and is over
     1660MB, that code is not included here.  A user may get individual
     17coverings from the LJCR using best_known_covering_design_www.
     18
     19In addition to the parameters and incidence structure for a covering
     20design from this database, we include extra information:
     21
     22\begin{itemize}
     23\item Best known lower bound on the size of a (v,k,t)-covering design
     24\item Name of the person(s) who produced the design
     25\item Method of construction used
     26\item Date when the design was added to the database
     27\end{itemize}
     28
     29REFERENCES
     30  [1] La Jolla Covering Repository, http://www.ccrwest.org/cover.html
     31  [2] Coverings, by Daniel Gordon and Douglas Stinson,
     32      http://www.ccrwest.org/gordon/hcd.pdf, in Handbook of Combinatorial Designs
     33
     34
     35AUTHORS:
     36    -- Daniel M. Gordon (2008-12-22): initial version
     37
     38"""       
     39
     40#*****************************************************************************
     41#      Copyright (C) 2008 Daniel M. Gordon <dmgordo@gmail.com>
     42#
     43# Distributed  under  the  terms  of  the  GNU  General  Public  License (GPL)
     44#                         http://www.gnu.org/licenses/
     45#*****************************************************************************
     46
     47import types
     48import urllib
     49from sage.misc.sage_eval import sage_eval
     50from sage.structure.sage_object import SageObject
     51from sage.rings.integer import Integer
     52from sage.rings.integer_ring import ZZ
     53from sage.rings.arith import binomial
     54from sage.combinat.combination import Combinations
     55from sage.rings.real_double import RDF
     56from sage.combinat.designs.incidence_structures import IncidenceStructure
     57
     58###################### covering design functions ##############################
     59
     60
     61def schonheim(v,k,t):
     62    r"""
     63    Schonheim lower bound for size of covering design
     64   
     65
     66    INPUT:
     67        v -- integer, size of point set
     68        k -- integer, cardinality of each block
     69        t -- integer, cardinality of sets being covered
     70
     71    OUTPUT:
     72        The Schonheim lower bound for such a covering design's size:
     73        $C(v,k,t) \leq \lceil(\frac{v}{k}  \lceil \frac{v-1}{k-1} \cdots \lceil \frac{v-t+1}{k-t+1} \rceil \cdots \rceil \rceil$
     74       
     75    EXAMPLES:
     76        sage: from sage.combinat.designs.covering_design import schonheim
     77        sage: schonheim(10,3,2)
     78        17
     79        sage: schonheim(32,16,8)
     80        930
     81    """
     82    bound = 1.0
     83    for i in range(t-1,-1,-1):
     84        bound = (bound*RDF(v-i)/RDF(k-i)).ceiling()
     85
     86    return bound
     87
     88
     89def trivial_covering_design(v,k,t):
     90    r"""
     91    Construct a trivial covering design.
     92
     93    INPUT:
     94        v -- integer, size of point set
     95        k -- integer, cardinality of each block
     96        t -- integer, cardinality of sets being covered
     97
     98    OUTPUT:
     99        (v,k,t) covering design
     100
     101    EXAMPLE:
     102        sage: C = trivial_covering_design(8,3,1)
     103        sage: C.show()
     104        C(8,3,1) = 3
     105        Method: Trivial
     106        0   1   2   
     107        0   6   7   
     108        3   4   5   
     109        sage: C = trivial_covering_design(5,3,2)
     110        sage: C.show()
     111        4 <= C(5,3,2) <= 10
     112        Method: Trivial
     113        0   1   2   
     114        0   1   3   
     115        0   1   4   
     116        0   2   3   
     117        0   2   4   
     118        0   3   4   
     119        1   2   3   
     120        1   2   4   
     121        1   3   4   
     122        2   3   4   
     123
     124    NOTES:
     125        Cases are:
     126        \begin{enumerate}
     127        \item $t=0$: This could be empty, but it's a useful convention to
     128        have one block (which is empty if $k=0$).
     129        \item $t=1$: This contains $\lceil v/k \rceil$ blocks:
     130        $[0,...,k-1]$,[k,...,2k-1],...$.  The last block
     131        wraps around if $k$ does not divide $v$.
     132        \item anything else:  Just use every $k$-subset of $[0,1,...,v-1]$.
     133        \end{enumerate}
     134
     135    """
     136    if t==0:     #single block [0,...,k-1]
     137        blk=[]
     138        for i in range(k):
     139            blk.append(i)
     140        return CoveringDesign(v,k,t,1,range(v),[blk],1,"Trivial")
     141
     142    if t==1:     #blocks [0,...,k-1],[k,...,2k-1],...
     143        size = (RDF(v)/RDF(k)).ceiling()
     144        blocks=[]
     145        for i in range(size-1):
     146            blk=[]
     147            for j in range(i*k,(i+1)*k):
     148                blk.append(j)
     149            blocks.append(blk)
     150
     151        blk=[]   # last block: if k does not divide v, wrap around
     152        for j in range((size-1)*k,v):
     153            blk.append(j)
     154        for j in range(k-len(blk)):
     155            blk.append(j)
     156        blk.sort()
     157        blocks.append(blk)
     158        return CoveringDesign(v,k,t,size,range(v),blocks,size,"Trivial")
     159
     160    # default case, all k-subsets
     161    return CoveringDesign(v,k,t,binomial(v,k),range(v),Combinations(range(v),k),schonheim(v,k,t),"Trivial")
     162
     163
     164class CoveringDesign(SageObject):
     165    """
     166    Covering design.
     167
     168    INPUT:
     169      data -- parameters (v,k,t)
     170              size
     171              incidence structure (points and blocks) -- default points are $[0,...v-1]$
     172              lower bound for such a design
     173              database information: creator, method, timestamp
     174    """
     175
     176    def __init__(self, v=0, k=0, t=0, size=0, points=[], blocks=[], low_bd=0, method='', created_by ='',timestamp=''):
     177        """
     178        EXAMPLES:
     179            sage: C=CoveringDesign(7,3,2,7,range(7),[[0, 1, 2], [0, 3, 4], [0, 5, 6], [1, 3, 5], [1, 4, 6], [2, 3, 6], [2, 4, 5]],0, 'Projective Plane')
     180            sage: C.show()
     181            C(7,3,2) = 7
     182            Method: Projective Plane
     183            0   1   2   
     184            0   3   4   
     185            0   5   6   
     186            1   3   5   
     187            1   4   6   
     188            2   3   6   
     189            2   4   5   
     190        """   
     191        self.v = v
     192        self.k = k
     193        self.t = t
     194        self.size = size
     195        if low_bd > 0:
     196            self.low_bd = low_bd
     197        else:
     198            self.low_bd = schonheim(v,k,t)
     199        self.method = method
     200        self.created_by = created_by
     201        self.timestamp = timestamp
     202        self.incidence_structure = IncidenceStructure(points, blocks)
     203
     204
     205    def show(self, max_size=100):
     206        """
     207        Displays a covering design.
     208
     209        INPUT:
     210            max_size -- maximum number of blocks (to avoid trying to show huge ones)
     211
     212        OUTPUT:
     213            covering design parameters and blocks, in a readable form
     214           
     215        EXAMPLES:
     216            sage: C=CoveringDesign(5,4,3,4,range(5),[[0,1,2,3],[0,1,2,4],[0,1,3,4],[0,2,3,4]],4, 'Lexicographic Covering')
     217            sage: C.show()
     218            C(5,4,3) = 4
     219            Method: Lexicographic Covering
     220            0   1   2   3   
     221            0   1   2   4   
     222            0   1   3   4   
     223            0   2   3   4   
     224        """
     225        if self.size == self.low_bd:    # check if the covering is known to be optimal
     226            print 'C(%d,%d,%d) = %d'%(self.v,self.k,self.t,self.size)
     227        else:
     228            print '%d <= C(%d,%d,%d) <= %d'%(self.low_bd,self.v,self.k,self.t,self.size);
     229        if self.created_by != '':
     230            print 'Created by: %s'%(self.created_by)
     231        if self.method != '':
     232            print 'Method: %s'%(self.method)
     233        if self.timestamp != '':
     234            print 'Submitted on: %s'%(self.timestamp)
     235        for i in range(min(max_size, self.size)):
     236            for j in range(self.k):
     237                print self.incidence_structure.blocks()[i][j]," ",
     238            print "\n",
     239        if(max_size < self.size):   # if not showing all blocks, indicate it with ellipsis
     240            print "..."
     241
     242
     243                   
     244    def is_covering(self):
     245        """
     246        Checks that all t-sets are in fact covered.
     247
     248        INPUT:
     249            putative covering design
     250
     251        OUTPUT:
     252            True if all t-sets are in at least one block
     253
     254
     255        NOTES:
     256            This is very slow and wasteful of memory.  A faster cython
     257            version will be added soon.
     258
     259        EXAMPLES:
     260            sage: C=CoveringDesign(7,3,2,7,range(7),[[0, 1, 2], [0, 3, 4], [0, 5, 6], [1, 3, 5], [1, 4, 6], [2, 3, 6], [2, 4, 5]],0, 'Projective Plane')
     261            sage: C.is_covering()
     262            True
     263            sage: C=CoveringDesign(7,3,2,7,range(7),[[0, 1, 2], [0, 3, 4], [0, 5, 6], [1, 3, 5], [1, 4, 6], [2, 3, 6], [2, 4, 6]],0, 'not a covering')   # last block altered
     264            sage: C.is_covering()
     265            False
     266        """
     267        v = self.v
     268        k = self.k
     269        t = self.t
     270        Svt = Combinations(range(v),t)
     271        Skt = Combinations(range(k),t)
     272        tset = {}       # tables of t-sets: False = uncovered, True = covered
     273        for i in Svt:
     274            tset[tuple(i)] = False
     275
     276        # mark all t-sets covered by each block
     277        for a in self.incidence_structure.blocks():
     278            for z in Skt:
     279                y = [a[x] for x in z]
     280                tset[tuple(y)] = True
     281
     282        for i in Svt:
     283            if tset[tuple(i)] == False:     # uncovered
     284                return False
     285       
     286        return True                  # everything was covered
     287
     288
     289
     290
     291   
     292
     293def best_known_covering_design_www(v, k, t, verbose=False):
     294    r"""
     295    Gives the best known (v,k,t) covering design, using database at
     296         \url{http://www.ccrwest.org/}.
     297
     298    INPUT:
     299        v -- integer, the size of the point set for the design
     300        k -- integer, the number of points per block
     301        t -- integer, the size of sets covered by the blocks
     302        verbose -- bool (default=False), print verbose message
     303
     304    OUTPUT:
     305        CoveringDesign -- (v,k,t) covering design with smallest number of blocks
     306
     307    EXAMPLES:
     308        sage: C = best_known_covering_design_www(7, 3, 2)   # requires internet, optional
     309        sage: C.show()                                      # requires internet, optional
     310        C(7,3,2) = 7
     311        Method: lex covering
     312        Submitted on: 1996-12-01 00:00:00
     313        0   1   2   
     314        0   3   4   
     315        0   5   6   
     316        1   3   5   
     317        1   4   6   
     318        2   3   6   
     319        2   4   5   
     320       
     321    This function raises a ValueError if the (v,k,t) parameters are
     322    not found in the database.
     323    """
     324    v = int(v)
     325    k = int(k)
     326    t = int(t)
     327
     328    param = ("?v=%s&k=%s&t=%s"%(v,k,t))
     329
     330    url = "http://www.ccrwest.org/cover/get_cover.php"+param
     331    if verbose:
     332        print "Looking up the bounds at %s"%url
     333    f = urllib.urlopen(url)
     334    s = f.read()
     335    f.close()
     336
     337    if 'covering not in database' in s:   #not found
     338        str = "no (%d,%d,%d) covering design in database\n"%(v,k,t)
     339        raise ValueError, str
     340
     341    return sage_eval(s)
     342