Centre Blaise Pascal » Bench4GPU

#!/usr/bin/env python

#

# Ising2D model in serial mode

#

# CC BY-NC-SA 2011 : <emmanuel.quemener@ens-lyon.fr> 

#

# Thanks to George Marsaglia : http://en.wikipedia.org/wiki/Multiply-with-carry

import sys

import numpy

#import pylab

from PIL import Image

from math import exp

from random import random

import time

z=0

w=0

def ImageOutput(sigma,prefix):

        Max=sigma.max()

        Min=sigma.min()

        # Normalize value as 8bits Integer

        SigmaInt=(255*(sigma-Min)/(Max-Min)).astype('uint8')

        image = Image.fromarray(SigmaInt)

        image.save("%s.jpg" % prefix)

SIZE=256

def Metropolis(sigma,J,T,iterations,seed_w,seed_z):

        Energy=0

        global z

        global w

        z=seed_z

        w=seed_w

        for p in range(0,iterations):

                # Random access coordonate

                #X,Y=numpy.random.randint(SIZE),numpy.random.randint(SIZE)

                # def znew():

                #         global z

                #         z=numpy.uint32(36969*(z&65535)+(z>>16))

                #         return(z)

                # def wnew():

                #         global w

                #         w=numpy.uint32(18000*(w&65535)+(w>>16))

                #         return(w)

                # def MWC(): return(numpy.uint32((znew()<<16)+wnew()))

                # def MWCfp(): return((MWC() + 1.0) * 2.328306435454494e-10)

                def MWC(): 

                        global w

                        global z

                        z=numpy.uint32(36969*(z&65535)+(z>>16))

                        w=numpy.uint32(18000*(w&65535)+(w>>16))

                        return(numpy.uint32((z<<16)+w))

                def MWCfp(): 

                        global w

                        global z

                        z=numpy.uint32(36969*(z&65535)+(z>>16))

                        w=numpy.uint32(18000*(w&65535)+(w>>16))

                        return(((numpy.uint32((z<<16)+w))+1.0)*

                               2.328306435454494e-10)

                X=numpy.uint32(SIZE*MWCfp())

                Y=numpy.uint32(SIZE*MWCfp())

                #print X,Y

                DeltaE=2.*J*sigma[X,Y]*(sigma[X,(Y+1)%SIZE]+

                                        sigma[X,(Y-1)%SIZE]+

                                        sigma[(X-1)%SIZE,Y]+

                                        sigma[(X+1)%SIZE,Y])

                if DeltaE < 0. or random() < exp(-DeltaE/T):

                        sigma[X,Y]=-sigma[X,Y]

                        Energy+=DeltaE

        return(Energy)

if __name__=='__main__':

    J=1.

    T=0.5

    iterations=numpy.uint32(SIZE*SIZE*SIZE)

    sigma=numpy.where(numpy.random.randn(SIZE,SIZE)>0,1,-1).astype(numpy.int8)

    ImageOutput(sigma,"Ising2D_OpenCL_%i_Initial" % (SIZE))

    seed_w=numpy.uint32(37)

    seed_z=numpy.uint32(91)

    i=0

    step=262144

    while (step*i < iterations):

        start=time.time()

        Metropolis(sigma,J,T,step,seed_w*(i+1),seed_z*((i+1)%SIZE))

        stop=time.time()

        elapsed = (stop-start)

        print "Iteration %i in %f: " % (i,elapsed)

        ImageOutput(sigma,"Ising2D_OpenCL_%i_%.3i_Final" % (SIZE,i))

        i=i+1