#
# Hardcore model on grid
# Bryan Clair, June 2012
#
import random

class HardcoreGrid:
    def __init__(self,sizex, sizey, fugacity=1.0):
        """Planar square grid of sizex X sizey vertices.
           Fugacity is like pressure and controls the density of sites"""
        self._sizex = sizex
        self._sizey = sizey
        self._contents = []
        for x in range(sizex):
            self._contents.append([False]*sizey)
        self._fugacity = fugacity

    def size(self):
        """Returns a pair (sizex, sizey)"""
        return (self._sizex,self._sizey)

    def occupied(self,x,y):
        """True if site (x,y) is occupied."""
        return self._contents[x][y]

    def density(self):
        """Returns density of occupied sites"""
        p = 0
        for col in self._contents:
            for site in col:
                if site:
                    p += 1
        return p/float(self._sizex * self._sizey)

    def polarity(self):
        """Returns fraction of occupied sites which are even polarity"""
        r = 0
        b = 0
        for x in range(self._sizex):
            for y in range(self._sizey):
                if self.occupied(x,y):
                    # exactly one of r, b will increase
                    r += (x+y) % 2
                    b += (x+y+1) % 2
        if r+b == 0:
            return 0
        return r/float(r+b)

    def set(self,x,y):
        """Put a particle at site (x,y)"""
        self._contents[x][y] = True

    def clear(self,x,y):
        """Clear site (x,y)"""
        self._contents[x][y] = False

    def toggle(self,x,y):
        """Toggle state of site (x,y)"""
        self._contents[x][y] = not self._contents[x][y]

    def __str__(self):
        """Convert to string, . for empty site, * for occupied."""
        out = ''
        for y in range(self._sizey):
            for x in range(self._sizex):
                out += '*' if self.occupied(x,y) else '.'
            out += '\n'
        return out[:-1]

    def neighbors(self,x,y):
        """Return number of occupied neighbors."""
        n = 0
        if x > 0 and self.occupied(x-1,y):
            n += 1
        if x < self._sizex-1 and self.occupied(x+1,y):
            n += 1
        if y > 0 and self.occupied(x,y-1):
            n += 1
        if y < self._sizey-1 and self.occupied(x,y+1):
            n += 1
        return n
        
    def oneStepMetropolis(self):
        """Take one step of the Metropolis  Markov chain.
           Return False if state was unchanged."""
        # pick location at random
        x = random.randint(0,self._sizex-1)
        y = random.randint(0,self._sizey-1)
        if self.neighbors(x,y) == 0:
            # Possibly change state
            # compute probability of switching state
            if self.occupied(x,y):
                switchprob = min(1.0,1.0/self._fugacity)/2.0
            else:
                switchprob = min(1.0,self._fugacity)/2.0
            # see if state switches
            if random.random() <= switchprob:
                self.toggle(x,y)
                return (x,y)
        return False

    def oneStepGlauber(self):
        """Take one step of the Glauber dynamics Markov chain.
           Return False if state was unchanged."""
        # pick location at random
        x = random.randint(0,self._sizex-1)
        y = random.randint(0,self._sizey-1)
        if self.neighbors(x,y) == 0:
            # Possibly change state
            if random.random() <= 1/(1.0+self._fugacity):
                # clear site
                if self.occupied(x,y):
                    self.clear(x,y)
                    return (x,y)
            else:
                # place particle
                if not self.occupied(x,y):
                    self.set(x,y)
                    return (x,y)
        return False

if __name__ == '__main__':
    print 'Neighbor count testing'
    g = HardcoreGrid(5,5)
    g.set(0,0)
    g.set(3,4)

    def test(x,y):
        g.set(x,y)
        print g
        g.clear(x,y)
        print (x,y), 'has', g.neighbors(x,y), 'neighbors'
        print

    for (x,y) in [(0,1), (1,0), (1,1), (4,3), (2,3), (3,2)]:
        test(x,y)

    test(3,3)
    g.set(4,3)
    test(3,3)
    g.set(2,3)
    test(3,3)
    g.set(3,2)
    test(3,3)
    
    print 'Metropolis chain'
    raw_input()
    print 'starting state:'
    print g
    for steps in range(10):
        k = 0
        while not g.oneStepMetropolis():
            k+=1
        print 'after',k,'tries, moved to:'
        print g

    print 'Fugacity testing, 10000 steps'
    raw_input()
    for fugacity in [0.1,0.5,2.0,10.0]:
        g = HardcoreGrid(10,10,fugacity)
        for steps in range(10000):
            g.oneStepGlauber()
        print 'fugacity',fugacity
        print g
        print 'density',g.density()
        print 'polarity',g.polarity()