CirclePacking.py

From CSclasswiki
Jump to: navigation, search
# pack.py
# D. Thiebaut
#
# reference: circle intersection
#            http://local.wasp.uwa.edu.au/~pbourke/geometry/2circle/
# reference: 4th inner circle adjacent to 3 touching circles:
#            http://en.wikipedia.org/wiki/Descartes'_theorem
#
from graphics import *
from math import sqrt, pi
import time
import sys
import random
import svg


DEBUG =  True      # if true the circles are labeled as they are added
SAMPLESIZE = 20    # the number of circle-i circle-j circles to pick
SCREENRADIUSRATIO =1.1
WIDTH = 700
HEIGHT = 700

def is_nan(num):
    """a NaN testing function"""
    return str(num) == "nan"


def rgb2hex(rgb_tuple):
    """ convert an (R, G, B) tuple to #RRGGBB """
    hexcolor = '#%02x%02x%02x' % rgb_tuple
    # that's it! '%02x' means zero-padded, 2-digit hex values
    return hexcolor

def hex2rgb(colorstring):
    """ convert #RRGGBB to an (R, G, B) tuple """
    colorstring = colorstring.strip()
    if colorstring[0] == '#': colorstring = colorstring[1:]
    if len(colorstring) != 6:
        raise ValueError, "input #%s is not in #RRGGBB format" % colorstring
    r, g, b = colorstring[:2], colorstring[2:4], colorstring[4:]
    r, g, b = [int(n, 16) for n in (r, g, b)]
    return (r, g, b)

"""
  ____ _          _       ____ _
 / ___(_)_ __ ___| | ___ / ___| | __ _ ___ ___
| |   | | '__/ __| |/ _ \ |   | |/ _` / __/ __|
| |___| | | | (__| |  __/ |___| | (_| \__ \__ \
 \____|_|_|  \___|_|\___|\____|_|\__,_|___/___/
                                              """
class circleClass:
    """ a class to implement a circle that can be moved and drawn
    on the screen """
    def __init__( self, xx, yy, rr, win, color=None ):
           self.x         = xx
           self.y         = yy
           self.r         = rr
           self.c         = None       # will hold the graphics object
           self.color     = color
           self.win       = win
           self.label     = None
           self.R0        = 200
           self.caption   = None       # will hold the words shown inside 
           self.t         = None       # anchor for the text object
           self.pointSize = None
           self.scale     = 1.0        # what we mult real radius by to
                                       # get # of pixels
           self.isDecorated = False    # if true, show color and caption
                                       # of circle
           
    def setCaption( self, C ):
        self.caption = C.replace( ' ', '\n' )
        self.computePointSize()

    def setPointSize( self, p ):
        self.pointSize = p

    def setDecorated( self, b ):
        self.isDecorated = b
        
    def setScale( self, s ):
        """ set the scale, define as the real radius in pixels divided
        by the real radius"""
        self.scale = float( s )
        
    def getCaption( self ):
        if self.caption is None:
            return ""
        return self.caption.replace( '\n', ' ' )
    
    def setLabel( self, label, R0=200 ):
        self.label = label
        self.R0    = R0
        
    def printStatus( self, caption="" ):
        print "circle %s: x=%d y=%d r=%d" % ( caption, self.x,
                                              self.y, self.r )

    def undraw( self ):
        if self.c is not None:
            self.c.undraw()
        if self.t is not None:
            self.t.undraw()
        
    def clear( self ):
        self.undraw()
        self.c = None
        self.t = None

    def display( self ):
        if self.c is None:          
            self.c = Circle( Point( self.x, self.y ), self.r )

        if DEBUG and self.label is not None and self.label!=0:
            t= Text( Point( self.x, self.y ), str( self.label ) )
            t.draw( self.win )
            
        if not self.isDecorated:
            self.c.draw( self.win )
            return
        
        if self.color is not None:
            self.c.setFill( self.color )
            self.c.setOutline( self.color )
            
        self.c.draw( self.win )
        
        if self.t is not None:
            self.t.draw( self.win )
        elif self.caption is not None:
            try:
                words = self.caption.split()
                maxLen= max( [len(k) for k in words] )
                if self.t is None:
                    self.t = Text( Point( self.x, self.y ), self.caption )
                    self.t.setSize( self.pointSize )
                    self.t.setFace( 'arial' )
                self.t.draw( self.win )
            except:
                pass
            
    def computePointSize( self ):
        if self.caption is not None:
            try:
                words = self.caption.split()
                maxLen= max( [len(k) for k in words] )
                self.pointSize = self.r*2 * self.scale / maxLen
                self.pointSize = int( min( 32, max( 5, self.pointSize ) ) )
            except:
                self.pointSize = None
            
    def setColor( self, col ):
         self.color = col
         if self.c is not None:
             self.c.setFill( col )
             self.c.setOutline( col )

    def setx( self, xx ):
         self.x = xx

    def sety( self, yy ):
           self.y = yy

    def getx( self ):
           return self.x

    def gety( self ):
           return self.y

    def getr( self ):
           return self.r


"""
  ____             __ _        ____ _
 / ___|___  _ __  / _(_) __ _ / ___| | __ _ ___ ___
| |   / _ \| '_ \| |_| |/ _` | |   | |/ _` / __/ __|
| |__| (_) | | | |  _| | (_| | |___| | (_| \__ \__ \
 \____\___/|_| |_|_| |_|\__, |\____|_|\__,_|___/___/
                        |___/
"""
class configurationClass:
    """class for a configuration of N disks, where the first disk is
       the containing disk, and the others are tangential """
   
    def __init__( self, win ):
         self.circleList = {}     # dico. key=cicle index, value=circle
         self.xoffset       = 0
         self.yoffset       = 0
         self.N             = 0
         self.win           = win
         self.dist          = []
         self.r0            = 200
         self.innerCircles  = []
         self.unusedCircles = []
         self.cornerList    = []
         self.degreeList    = []
         self.dico1         = {}      # used to cache adjacent circles to (p,q)
         self.dico2         = {}      # used to cache adjacent circles to (p,q)
         self.distDico      = {}      # cache for distance computation
         self.color         = []      # array of colors for the circles (except C0)
         self.labels        = []      # labels for each circle
         self.show          = True
         self.maxRadius     = None
         self.lastGoodConf  = []
         self.scale         = 1.0     # what real radius must be x by to get
                                      # number of pixels
         self.dataFileName  = None
         self.height        = 1000
         self.width         = 1000
         
    def setDataFileName( self, f ):
        self.dataFileName = f
        
    def setCoords( self, x, y ):
        self.win.setCoords( -x, -y, x, y )
        self.xoffset = x
        self.yoffset = y

    def setWidthHeight( self, w, h ):
        self.height = h
        self.width  = w
        
    def setShow( self, b ):
        self.show = b
        
    def setLabels( self, L ):
        self.labels = L
        
    def complete( self ):
        return len( self.innerCircles )==len( self.dist )

    def degreeListEmpty( self ):
         return len( self.degreeList ) == 0

    def setDecorated( self, b ):
        for key in self.circleList:
            self.circleList[key].setDecorated( b )
            
    def clear( self ):
         #self.win.clear()
         self.circleList   = {}
         self.innerCircles = []
         self.unusedCircles= range( 1, len( self.dist ) )
         self.degreeList   = []
         self.dico1        = {}
         self.dico2        = {}
         self.distDico     = {}
         
    def listOfUnusedCircles( self ):
        return self.unusedCircles
        
    def setOffsets( self, W, H ):
         """ makes the middle point of the window win the origin
         of the system of axes"""
         self.xoffset = W
         self.yoffset = H
         
    def getNoCircles( self ):
         """ returns the number of circles in the distribution
         of disks, including C0"""
         return len( self.dist )

    def appendCircle( self, k, xk, yk, rk ):
         c = circleClass( xk, yk, rk, self.win, self.color[ k ] )
         c.setScale( self.scale )
         #print "setting scale of circle %d to %f" % ( k, self.scale )
         
         if len( self.labels ) > k:
             c.setCaption( self.labels[k] )
             #print "setting caption of circle",k,"to", self.labels[k]
             
         self.innerCircles.append( k )
         self.unusedCircles.remove( k )
         #print "removing %d: unusedCircles = %s" % ( k, str( self.unusedCircles ) )
         self.circleList[k] = c
         c.setLabel( k )
         if self.show:
             c.display()
             self.win.update()

    def setDist( self, dist ):
         self.dist = dist[:]
         self.unusedCircles = range( 1, len( dist )+1 )
         N = len( self.dist )
         print "length of unusedCircles = ", len( self.unusedCircles )
         print "max unused index = ", max( self.unusedCircles )
         
         #--- set the largest circle after the containing one ---
         self.maxRadius = max( self.dist )
         
         #--- add containing circle at index 0 ---
         self.dist = [self.r0] + self.dist

         #--- define unknown colors for each ---
         for i in range( len( self.dist ) ):
             self.color.append( None )
             
    def setLabels( self, L ):
        self.labels = [""] + L # add fake label for C0
        
    def colorCircles( self ):
        """ assign a color to each of the circles """
        N = len( self.dist )
         
        #--- set color range ---
        r1, g1, b1 = hex2rgb( "f00080" )
        r2, g2, b2 = hex2rgb( "f0f080" )
        dr = (r2-r1)*1.0/ N
        dg = (g2-g1)*1.0/ N
        db = (b2-b1)*1.0/ N
        self.color[0] = "#f0f0B0"
        for i in range( 1, len( self.dist ) ):
             self.color[i] = rgb2hex( ( int(r1+i*dr), int(g1+i*dg), int(b1+i*db) ) )

    def scaleCircles( self ):
        """assign a scale factor to all the circles so that the point-size used
        to display the caption will be right """
        global SCREENRADIUSRATIO
        pixelWidth = self.win.getWidth()
        realWidth = SCREENRADIUSRATIO * self.dist[ 0 ] * 2
        scale = pixelWidth * 1.0 / realWidth
        for i in range( 1, len( self.dist ) ):
            self.circleList[ i ].setScale( scale )
        
    def setContainingRadius( self, r ):
         self.r0 = r
         if len( self.dist ) > 0:
               self.dist[0] = r

         # compute the scaling factor
         #  |<--------- Win.width() ---------->|
         #  = 2 * containingRadius * SCREENRADIUSRATIO
         #
         self.scale = self.win.getWidth() / 2.0 / r / SCREENRADIUSRATIO

    def printStatus( self, caption="Configuration" ):
         print "\n\n========================================"
         print caption
         print "========================================"
         print "Dist = ", str( self.dist )
         print "# circles = ", len( self.dist )
         print "InnerCircles = ", str( self.innerCircles )
         print "# inner circles = ", len( self.innerCircles )
         print "# circles in circleList = ", len( self.circleList )
         print "circleList = "
         for key in self.circleList:
               c= self.circleList[ key ]
               print "%d:(%d,%d, %d) " % ( key, c.x, c.y, c.r ),
         print "\nR0 = ", self.r0
         print "CornerList = ", str( self.cornerList )
         print "DegreeList = "
         for tuple in self.degreeList:
                 ( degree, k, xk, yk, rk, p, q ) = tuple
                 print "  --- k=%d d=%5.2f (%d, %d) " % ( k, degree, p, q )
         print "========================================"
         print "\n\n"

    def saveCurrentConfiguration( self ):
        """ save the current configuration so that we can redraw it later"""
        self.goodConf = []
        for key in self.circleList:
            c = self.circleList[ key ]
            caption = c.getCaption()
            self.goodConf.append( ( key, c.x, c.y, c.r, c.color, c.pointSize, caption ) )
        self.goodConf.sort()
        self.goodConf = [ len( self.circleList ) ] + self.goodConf
        
        if self.dataFileName is not None:
            FILE = open( self.dataFileName + ".out.xml", "wt" )
            FILE.write( "\n".join( self.parametersInXML() ) )
            FILE.close()
                
    def displayGoodConf( self ):
        """ display the current configuration """
        N = 0
        for i, tuple in enumerate( self.goodConf ):
            if i==0:
                N = tuple
            else:
                ( k, x, y, r, color, pointSize, caption ) = tuple
                c = circleClass( x, y, r, self.win, color )
                c.setDecorated( True )
                c.setCaption( caption )
                c.setPointSize( pointSize )
                c.display()

    def printGoodParameters( self ):
        print '\n'.join( self.parametersInXML() )
        
    def parametersInXML( self ):
        """ display all that the user might be interested in"""
        #print "parametersInXML: ", str( self.goodConf )
        lines = ["<circlepacking>"]
        N = 0
        for i, tuple in enumerate( self.goodConf ):
            #print str( tuple )
            if i==0:
                N = tuple
                lines.append( "  <noOfDisks>%d</noOfDisks>" % N )
            elif i==1:
                ( k, x, y, r, color, pointSize, caption ) = tuple
                lines.append( "  <containingRadius>%d</containingRadius>" % r )
            else:
                ( k, x, y, r, color, pointSize, caption ) = tuple
                lines.append( "  <circle>" )
                lines.append( "    <no>%d</no>" % k )
                lines.append( "    <x>%f</x>" % x )
                lines.append( "    <y>%f</y> " % y )
                lines.append( "    <r>%f</r> " % r )
                lines.append( "    <color>%s</color>" % color )
                lines.append( "    <pointsize>%s</pointsize>" % str( pointSize ) )
                lines.append( "    <caption>%s</caption>" % caption )
                lines.append( "  </circle>" )
                
        lines.append( "</circlepacking>" )
        return lines
    
    def generateSVGFile( self, fileName = None ):
        """ generate an svg file """
        global WIDTH, HEIGHT
        
        #print "generateSVGFile: ", str( self.goodConf )
        
        if fileName==None:
            fileName = self.dataFileName + "_%dx%d.svg" % ( WIDTH, HEIGHT )

        scene = svg.Scene( fileName, self.width, self.height )
        scene.setScale( self.scale )
        scene.setOffset( self.r0 * SCREENRADIUSRATIO )
        scene.createGradients( "f00080", "f0f080", len( self.dist ) )
        
        for i, tuple in enumerate( self.goodConf ):
            #print str( tuple )
            if i>=1:
                ( k, x, y, r, color, pointSize, caption ) = tuple                    
                scene.add( svg.Circle( ( x, y ), r, color, k ) )
                if pointSize is not None:
                    pointSize = int( ( 1.5 * pointSize ) / self.scale )
                    #pointSize = int( pointSize  )
                scene.add( svg.Text( ( x, y ), caption, pointSize ) )
        scene.write_svg()
                           
    def saveParameters( self, dataFileName ):
        """ save all that the user might be interested in"""
        
        file = open( dataFileName+".txt", "wt" )
        file.write( "Circle Packing             : %s\n" % dataFileName )
        file.write( "Containing Radius          : %10.3f\n", self.dist[0] )
        file.write( "Center of containing circle: 0, 0\n" )
        file.write( "Coordinate system          : -%d, -%d, %d, %d \n" 
              % ( self.xoffset, self.yoffset, self.xoffset, self.yoffset ) )
        file.write( "Number of inner circles    :", len( self.dist )-1 )
        file.write( "%10s  %10s  %10s   %8s %5s %s\n" %
                    ("x", "y", "radius", "color", "point", "caption" ) )
        for key in self.circleList:
            c = self.circleList[ key ]
            caption = c.getCaption()
            file.write( '%10.4f  %10.4f   %10.4f  %8s %5d %1s \n' 
                  % ( c.x, c.y, c.r, c.color, c.pointSize, caption ) )
        file.close()
        
    def initialize( self, i, j, containingRadius=None ):
        self.clear()
        if len( self.dist )==0:
             print "configuration.initialize() error. 0-lengt dist"
             sys.exit()

        if containingRadius is not None:
             self.setContainingRadius( containingRadius )

        r0 = self.r0
        c0 = circleClass( 0, 0, r0, self.win, self.color[0] )
        c0.setScale( self.scale )
        c0.setLabel( 0 )
        self.circleList[0] = c0
        try:
               ri = self.dist[ i ]
               rj = self.dist[ j ]
        except:
               print "*** ERROR, initial i or j outside list boundaries ***"
               sys.exit()

        #--- no need to explore this is circles cannot fit
        if ri> r0 or rj>r0 or ri+rj > r0:
            return False

        #--- create first inner circle, Circle i                         
        #--- put it on the positive x-axis, against the containing circle
        ci = circleClass( 0, 0, ri, self.win )
        ci.setScale( self.scale )
        ci.setLabel( i )
        ci.setColor( self.color[i] )
        ci.setx( r0-ri )
        ci.sety( 0 )
        if len( self.labels ) > i:
            ci.setCaption( self.labels[i] )
        self.circleList[i] = ci
        
        # create second circle, Circle j above circle i, tangent to 
        # containing circle
        a = r0 - ri
        b = ri + rj
        c = r0 - rj
        x = ( a*a +c*c - b*b )/ 2 / a
        cj = circleClass( x, sqrt( c**2 - x**2 ), rj, self.win )
        cj.setScale( self.scale )
        cj.setLabel( j )
        cj.setColor( self.color[j] )
        if len( self.labels ) > j:
            cj.setCaption( self.labels[j] )
        self.circleList[j] = cj

        #--- keep track of indexes of already placed circles ---
        self.innerCircles = [ 0, i, j ]
        self.unusedCircles.remove( i )
        self.unusedCircles.remove( j )
        
        #--- create a new corner list ---
        self.createCornerList( i, j )
        return True
    
    def createCornerList( self, i, j ):
        """create a list of pairs of the form [ (0,i), (0,j), (i,j) ]
        where 0, i, and j are the indexes of the container and first 2
        circles in the configuration"""
        self.cornerList = [ ( 0, i ), ( 0, j ), ( min(i,j), max(i,j) ) ]
        return

    def addToCornerList( self, tuple ):
        """ given that k is the index of the newly added circle to the
        configuration, generate the list of all new corners it forms with
        the circles already in the configuration """
        #--- add corner only if two circles closer than radius of largest
        #--- remaining circle
        listRadius = []
        for j in self.unusedCircles:
            listRadius.append( self.dist[j] )
        try:
            maxUnusedRadius = max( listRadius )
        except:
            maxUnusedRadius = 0
            
        ( hd, k, xk, yk, rk, p, q ) = tuple
        for i in self.circleList:
            if i==k:
                continue
            ci = self.circleList[i]
            distance_i_k = self.distance( ci.x, ci.y, xk, yk )
            if distance_i_k < ci.r + rk + 2* maxUnusedRadius:               
                self.cornerList.append( ( min( i, k ), max( i, k ) ) )

    def validateCornerList( self ):
        """prune corners that are not any longer viable.  Two circles form
        and invalid corner if the distance between their centers plus their
        radius is larger than the largest radius remaining"""

        #--- compute largest radius remaining ---
        listRadius = []
        for j in self.unusedCircles:
            listRadius.append( self.dist[j] )
        try:
            maxUnusedRadius = max( listRadius )
        except:
            maxUnusedRadius = 0

        toRemove = []
        for ( p, q ) in self.cornerList:
            cp = self.circleList[ p ]
            cq = self.circleList[ q ]
            distance_p_q = self.distance( cp.x, cp.y, cq.x, cq.y )
            if distance_p_q < cp.r + cq.r + 2 * maxUnusedRadius:
                toRemove.append( ( p, q ) )

        for ( p, q ) in toRemove:
            self.cornerList.remove( (p,q) )
            
    def getCornerList( self ):
        return self.cornerList

    def getLengthCornerList( self ):
        return len( self.cornerList )
    
    def holeDegree( self, xk, yk, rk, i, j ):
       """computes the hole degree, i.e. shortest distance from the k-Circle
       to the closest circle that is not Circle i, or Circle j"""
       degreeMin = 999999999.0
       for n, cn in self.circleList.items():
           if n==i or n==j: continue
           if n==0:
               r0 = self.circleList[0].getr()
               #d = r0 - sqrt(xk**2 + yk**2 )
               d = r0 - self.distance( 0, 0, xk, yk ) 
               if d < degreeMin:
                   degreeMin = d
           else:
               cn = self.circleList[n]
               xn = cn.getx()
               yn = cn.gety()
               rn = cn.getr()
               #d = sqrt( (xn-xk)**2 + (yn-yk)**2 ) - rk - rn
               d = self.distance( xn, yn, xk, yk ) - rk - rn
               if d < degreeMin:
                   degreeMin = d
       return 1 - degreeMin/rk
    
    def isGood( self, j, xj, yj, rj ):
       """returns False if circle (xj, yj, rj) overlaps with some of the circles
       or extends out the containing circle.  Return True otherwise"""
       #print "isGood(%d,%d,%f)" % (xj,yj,rj)
       #print "circleList = ",str( self.circleList )
        
       for i, ci in  self.circleList.items():
           #print " isGood: i=", i
           if i==j:
               continue
            
           xi = ci.getx()
           yi = ci.gety()
           ri = ci.getr()
           epsillon = 0.001
           if i==0:
               d2 = xj**2 + yj**2
               if d2  > (ri-rj)**2 + epsillon:
                   #print "     case i==0 d2=%f (ri-rj)**2=%f" \
                   #      % ( d2, (ri-rj)**2 )
                   return False
           else:
               d2 = (xj-xi)**2 + (yj-yi)**2
               if d2 + epsillon < (ri+rj)**2:
                   #print "     case d2=%f<(ri+rj)**2=%f" \
                   #      % ( d2, (ri+rj)**2 )
                   return False
       return True
    
    def addHoleDegree( self, tuple ):
        """ add the tuple ( degree, k, xk, yk, rk, p, q ) to the list, where
        p and q are the indices of the two circles forming the corner"""        
        self.degreeList.append( tuple )
        self.degreeList.sort()
        self.degreeList.reverse()

    def updateDegreeList( self, plast, qlast ):
        """all the tuples associated to the corner (p,q) where we just put the
        k-th circle should have their degree recomputed and tested.
        All the tuples who have one circle in common with p or q should
        also have their degree recomputed"""

        L = []          # list of new tuples
        toRemove = []   # list of tuples to remove
        for tuplei in self.degreeList:
            ( hi, i, xi, yi, ri, pi, qi ) = tuplei

            #--- remove the tuples for circles already added ---
            if i in self.innerCircles:
                toRemove.append( tuplei )
                continue
            
            #--- find corners with one disk in common with k-th disk---
            #if pi==plast or qi==qlast:
            #    holeDegree = self.holeDegree( xi, yi, ri, pi, qi )
            #    if holeDegree > hi and self.isGood( i, xi, yi, ri ):
            #        toRemove.append( tuplei )
            #        L.append( ( holeDegree, i, xi, yi, ri, pi, qi ) )
            
            if not self.isGood( i, xi, yi, ri ):
                toRemove.append( tuplei )
            else:
                #--- recompute the degree, in case it has changed, and add it to temp list ---
                holeDegree = self.holeDegree( xi, yi, ri, pi, qi )
                if holeDegree > hi:
                    toRemove.append( tuplei )
                    L.append( ( holeDegree, i, xi, yi, ri, pi, qi ) )

        #--- remove what needs to be removed ---
        for tuplei in toRemove:
            self.degreeList.remove( tuplei )

        #--- reorder with highest degree in front of list ---
        self.degreeList = self.degreeList + L                        
        self.degreeList.sort()
        self.degreeList.reverse()

    def getBestHoleDegree( self ):
        ( bestDegree, k, xk, yk, rk, bestP, bestQ ) = bestHoleDegree = self.degreeList[0]
        toRemove = []
        for tuple in self.degreeList:
            ( degree, i, xi, yi, ri, p, q ) = tuple
            if i==k:
                toRemove.append( tuple )
        for tuple in toRemove:
            self.degreeList.remove( tuple )                
        return bestHoleDegree

    def display( self ):
        """display the circles on the window"""
        if self.show == False:
            return
        
        for index in self.circleList:
            self.circleList[index].display()

    def distance( self, x0, y0, x1, y1 ):
        if ( x0, y0, x1, y1 ) in self.distDico:
            return self.distDico[ ( x0, y0, x1, y1 ) ]
        d = sqrt( (x0-x1)**2 + (y0-y1)**2 )
        self.distDico[ ( x0, y0, x1, y1 ) ] = d
        self.distDico[ ( x1, y1, x0, y0 ) ] = d
        return d
    
    def findAdjacentCirclesToC0( self, i, rj ):
       """finds the coordinates of the center of the 2 circles in
       the corner of the ith circle and the containing circle"""
       if   (i,rj) in self.dico1:
           return self.dico1[ (i, rj ) ]
        
       P1 = self.circleList[ i ]
       P0 = self.circleList[ 0 ]

       r1 = P1.getr() + rj
       r0 = P0.getr() - rj

       x0 = P0.getx()
       y0 = P0.gety()

       x1 = P1.getx()
       y1 = P1.gety()
       
       #d = sqrt( (x1-x0)**2 + (y1-y0)**2 )
       d = self.distance( x0, y0, x1, y1 )
           
       a = ( r0**2 - r1**2 +d**2 )/ (2*d)
           
       try:        
           h = sqrt( r0**2 - a**2 )
       except: 
           return ( None, None, None, None )

       if is_nan( h ) or d==0:
           return ( None, None, None, None )
               
       x2 = x0 + a*(x1-x0)/d
       y2 = y0 + a*(y1-y0)/d
       
       x3 = x2 + h*( y1-y0)/d
       y3 = y2 - h*( x1-x0 )/d

       x4 = x2 - h*( y1-y0)/d
       y4 = y2 + h*( x1-x0 )/d

       self.dico1[ (i,rj) ] = ( x3, y3, x4, y4 ) 
       return ( x3, y3, x4, y4)
      
    def findAdjacentCircles( self, i1, i2, r3 ):
       """given the indexes i1 and i2 of the circles to be considered
       and the radius r3 of the 2 circles that will fit in their
       corner, return x3, y3, and x4, y4 of the two circles of radius
       r3 that fits in the 2 corners of the 2 circles.
       Works only for circles that intersect at more than 1 point"""
       if ( i1, i2, r3 ) in self.dico2:
           return self.dico2[ (i1, i2, r3 ) ]
    
       P0 = self.circleList[ i1 ]
       P1 = self.circleList[ i2 ]
       x0 = P0.getx()
       y0 = P0.gety()
       r0 = P0.getr()
       x1 = P1.getx()
       y1 = P1.gety()
       r1 = P1.getr()
       r1 = r1+r3
       r0 = r0+r3
       
       #d = sqrt( (x1-x0)**2 + (y1-y0)**2 )
       d = self.distance( x0, y0, x1, y1 )
       a = ( r0**2 - r1**2 +d**2 )/ (2*d)
           
       try:
           h = sqrt( r0**2 - a**2 )
       except:
           return ( None, None, None, None )

       if d==0 or is_nan( h ):
           return ( None, None, None, None )
        
       x2 = x0 + a*(x1-x0)/d
       y2 = y0 + a*(y1-y0)/d
       x3 = x2 + h*( y1-y0)/d
       y3 = y2 - h*( x1-x0 )/d

       x4 = x2 - h*( y1-y0)/d
       y4 = y2 + h*( x1-x0 )/d
           
       self.dico2[ (i1, i2, r3 ) ] = ( x3, y3, x4, y4 )
       return ( x3, y3, x4, y4)

           
#----------------------------------------------------------------
def waitForClick( win, message ):
   """ waitForClick: stops the GUI and displays a message.  Returns
       when the user has clicked the window, and erases the message."""

   # wait for user to click mouse to start
   #startMsg = Text( Point( win.getWidth()/2, win.getHeight()/2 ), message )
   startMsg = Text( Point( 0, 0 ), message )
   startMsg.draw( win )    # display message
   win.getMouse()          # wait
   startMsg.undraw()       # erase


#----------------------------------------------------------------
def process( conf, containingRadius, win ):
    """ process: tries to find a configuration that fits in the given
    radius"""
    global SAMPLESIZE
    
    # 0) put C0 in configuration
    conf.setContainingRadius( containingRadius )

    #--- start with disks i, j ---
    rangei = range( 1, conf.getNoCircles()-1 )
    if len( rangei )  > SAMPLESIZE:
        rangei = random.sample( rangei, 3 )

    for i in rangei:
        rangej = range( i+1, conf.getNoCircles() )
        if len( rangej ) > SAMPLESIZE:
            rangej = random.sample( rangej, SAMPLESIZE )
        for j in rangej:

            # 1) put Ci in configuration
            # 2) put Cj in configuration
            if not conf.initialize( i, j ):
                print "initialize returns False"
                continue
            
            conf.display()

            # 3) define corner list
            # conf.createCornerList( i, j )  # already done in initialize() 

            # 4) for each in k in [1:N+1] and k not in innerCircles
            #       for each corner (p,q) last added to cornerList
            #            if circle k fits in (p,q)
            #                compute lambda( k, p, q, xk, yk )
            #                degreeList.appen( lambda( k, p, q, xk, yk )

            while True:
                foundK = False

                #--- instead of looking at ALL the remaining circles,
                #--- just pick a sample of them, but only if we are
                #--- in exploration mode (i.e. first is None)
                unusedCircles = conf.listOfUnusedCircles()
                #print "unuseCircles from listOfUnusedCircles()=", unusedCircles
                
                for k in unusedCircles: # conf.listOfUnusedCircles():
                    #print "   k = ", k
                    rk = conf.dist[k]

                    cornerList = conf.getCornerList()
                    for (p,q) in cornerList:
                        #print "Corner (%d, %d)" % ( p, q )
                        if p==0:
                           [x3, y3, x4, y4 ] = conf.findAdjacentCirclesToC0( q, rk )
                        elif q==0:
                           [x3, y3, x4, y4 ] = conf.findAdjacentCirclesToC0( p, rk )
                        else:
                           [ x3, y3, x4, y4 ] = conf.findAdjacentCircles( p, q, rk )

                        if x3 is None or x4 is None:
                            continue

                        if is_nan( x3 ) or is_nan( y3 ) or is_nan( x4 ) \
                           or is_nan( y4 ):
                            print "*** NaN ***", p, q, x3, y3, x4, y4
                            continue
                        
                        if conf.isGood( k, x4, y4, rk ):
                            holeDegree = conf.holeDegree( x4, y4, rk, p, q )
                            conf.addHoleDegree( ( holeDegree, k, x4, y4, rk, p, q ) )
                            foundK = True
                            #print "      Adding new hole"
                            
                        if conf.isGood( k, x3, y3, rk ):
                            holeDegree = conf.holeDegree( x3, y3, rk, p, q )
                            conf.addHoleDegree( ( holeDegree, k, x3, y3, rk, p, q ) )
                            foundK = True
                            #print "      Adding new hole"


                #print "end of for-k loop: foundK = %d len(degreeList)=%d" \
                #      % ( foundK, len( conf.degreeList ) )
                
                # 6) pick best hole-degree circle
                if conf.degreeListEmpty():
                    #print "       *** Degree list empty: break ***"
                    break
                
                tuple = conf.getBestHoleDegree( )
                ( hd, k, xk, yk, rk, p, q ) = tuple
                
                if not conf.isGood( k, xk, yk, rk ):
                    #print "       *** conf is not good: continue ***"
                    continue
                    
                #print "adding best hole-degree circle: ", k,\
                #      "degree=", hd, "between cirlces", p, "and", q
                
                # 7) add circle to configuration
                conf.appendCircle( k, xk, yk, rk ) # add circle, make it
                #print "adding circle", k
                
                conf.addToCornerList( tuple )
                # trim corner list every 10 new corners
                if conf.getLengthCornerList()%10==0:
                    conf.validateCornerList()

                                
                # 7.2) update holeDegree list to reflect new addition
                conf.updateDegreeList( p, q )
                                
                # 8) keep on going as long as we have more circles to add
                if conf.complete():
                    break

            #conf.printStatus( )
            
            if conf.complete():
                #waitForClick( win, "configuration complete" )
                return (True, i, j )
            
            elif conf.degreeListEmpty():
                #waitForClick( win, " S T U C K " )
                continue
    return (False, None, None)

def bestGuestForContainingCircle( L ):
    """ given the list L of inner circles, find their total
    area, and assuming an 85% packing efficiency, guess-timate
    the area of the containing circle, and from that the
    containing radius"""
    totalArea = 0
    for radius in L:
        area = pi * radius**2
        totalArea += area
    conservativeContainingArea = totalArea / 0.80
    approxRadius = sqrt( conservativeContainingArea/ pi )
    return approxRadius

def getDistribution( dataFileName ):
    """get a distribution of circles and their labels from a file"""
    
    text = open( dataFileName, "r" ).readlines()

    Radii = []
    Labels = []
    temp   = []
    N = int( text[0] )
    for i, line in enumerate( text ):
        if i==0:
            N = int( line )
            continue
        words = line[:-1].strip().split()
        temp.append( ( sqrt( float( words[0] ) ), '\n'.join( words[1:] ) ) )
        
    #--- sort circles by decreasing radii ---
    temp.sort()
    temp.reverse()

    #--- pack in two separate lists ---
    for ( radius, caption ) in temp:
        Radii.append( radius )
        Labels.append( caption )

    #for i, label in enumerate( Labels ):
    #   print '%2d: %7.3f "%s"' % ( i+1, Radii[i], label )
    
    return ( Radii, Labels )



"""----------------------------------------------------------------
                     _
     _ __ ___   __ _(_)_ __
    | '_ ` _ \ / _` | | '_ \
    | | | | | | (_| | | | | |
    |_| |_| |_|\__,_|_|_| |_|
   ----------------------------------------------------------------"""

def main():
   global SCREENRADIUSRATIO
   global WIDTH
   global HEIGHT
   win = GraphWin( "DT's Circle Packing", WIDTH, HEIGHT  )
   win.setBackground( "white" )
   
   #--- generate configuration ---
   conf = configurationClass( win )
   conf.setOffsets( 0, 0 )
   conf.setWidthHeight( WIDTH, HEIGHT )

   #--- initialize with list of radii ---
   dataFileName = "data100.dat"
   (radii, labels) = getDistribution( dataFileName )
   conf.setDataFileName( dataFileName )
   conf.setDist( radii )
   conf.setLabels( labels )
   conf.colorCircles()
   conf.setShow( False )
   
   #--- compute best guess for containing circle ---
   (containingRadius ) = bestGuestForContainingCircle( radii )
   print "Starting with radius", containingRadius

   #--- define a new    
   #--- get a configuration for the given radius ---
   winningRadius = None
   #waitForClick( win, "click to start" )
   
   while True:
       conf.setContainingRadius( containingRadius )
       #win.clear()
       conf.setCoords( SCREENRADIUSRATIO * containingRadius, SCREENRADIUSRATIO * containingRadius ) 
       (ok, winningi, winningj) = process( conf, containingRadius, win )
       if ok:
           print "Found good configuration"
           winningRadius = containingRadius
           conf.saveCurrentConfiguration()
           conf.generateSVGFile()
           conf.displayGoodConf()
           conf.printGoodParameters()
           containingRadius = containingRadius * 0.98
           print "Trying radius", containingRadius
       else:
            if winningRadius is None:
                win.clear()
                waitForClick( win, "Could not find a configuration" )
                return
            else:
                print "Displaying final configuration"
                win.clear()
                conf.displayGoodConf()
                conf.generateSVGFile()
                conf.displayGoodConf()
                conf.printGoodParameters()
                break
            
   #--- ready to close window ---
   waitForClick( win, "" )
   win.close()

if __name__=="__main__":
   main()