# Import modules # First, the ones we wrote... from Automata import * from Basic_Flow import * from Wrapped2DArray import * from analysis import * #General python modules import time from random import randint from numpy import * from pylab import * from Scientific.Statistics.Histogram import Histogram from Scientific.Functions.LeastSquares import leastSquaresFit try: import Numeric as N import pygame import pygame.surfarray as surfarray except ImportError: raise ImportError, "Numeric and PyGame required." # SET PARAMETERS FOR ITERATING AND DISPLAYING ############################################# #Display parameters CellSize = 9 # how many pixels each cell takes up (so we can actually see it) nSteps = 25 # number of steps between calculating speed and printing result hgStep = 20 numBins = 10 disp_hg = False # **NOT WORKING YET** disp_settle_time = True disp_fit =True disp_multi_run_avg =False# True disp_all_runs = False #True frame_rate = 1 # refresh display every *frame_rate* iterations make_movie = False #Iteration parameters nIter = 2 # number of iterations total before simulation stops nSites =10 # size (one dimension) of array to be iterated nRuns = 1 # number of times to run program (new seed each time) dorandom = 11 # 0 - load file from array # 1 - load random IC # n - load the nth seed #Flow Rule Parameters #NOTE: all of these params may be superseded # by main program below (see last two dozen lines or so) c0 = 0.3 # parameter determining how much mass moves # valid range is [0.0,1.0] k0 = 0.1 # sine wave number w0 = 0.01 # sine frequency mode = 0 # 0 - original averaging mode # 1 - parameter varying mode # 2 - parameter and seed varying mode (averaging done right) ############################################# def Iterate(x,nSites): x.step() return x.cur.ThisArray # Spatial configuration display def SpDisplayState(state,nSites,surface,CellSize): ConfigRect = pygame.Rect(0,0,nSites*CellSize,nSites*CellSize) surface.fill( (0,0,255) ,ConfigRect) for i in xrange(nSites): for j in xrange(nSites): fc = int((1-state[i,j])*255) # grayscale based on state value if fc>255: # make sure we're within bounds of RGB vals fc = 255 if fc<0: fc = 0 fill_Color = fc,fc,fc surface.fill(fill_Color,[i*CellSize,j*CellSize,CellSize,CellSize]) def One_Run(nSites,flow_type,params): #create instance of class if flow_type == 1: x = Adjacent_Flow() #nearest neighbor if flow_type == 2: x = C8_Flow_Half() #eight cell neighborhood if flow_type == 3: x = C5_Flow_Half() #5 cell neighborhood if flow_type == 4: x = Rev_Flow() #reverse 1 flow if flow_type == 5: x = C8_Flow_Full() #eight cell neighborhood for both if flow_type == 6: x = C5_Flow_Full() #five cell neighborhood for both if flow_type == 7: x = Sine_Simple_Flow_Full() #1 driven by sine wave if flow_type == 8: x = Sine_C8_Flow_Full() #8 driven by sine wave if flow_type == 9: x = Sine_C5_Flow_Full() #5 driven by sine wave #Set parameters x.c = params[0]/4.0 #where our parameter c is a number between 0 and 1.0 x.k = params[1] x.w = params[2] # Set initial configuration x.seed(dorandom,nSites) nSites = x.curSeed.arraysize # since if dorandom == 0 the file loaded from # array may have changed size from number above #Set up display size = CellSize*nSites pygame.init() OffColor = 0, 0, 255 # Get display surface screen = pygame.display.set_mode( (size,size) ) pygame.display.set_caption('2D Cellular Automaton Simulator') # Clear display pygame.display.flip() # Create RGB array whose elements refer to screen pixels # Strangeness: sptmdiag is no longer used, but if this line # is removed, then display updates slow down considerably! sptmdiag = surfarray.pixels3d(screen) #iterate state, updating display each time t=0 t0 = time.clock() while (t <= nIter): for event in pygame.event.get(): if event.type == pygame.QUIT: sys.exit() # Iterate state state = Iterate(x,nSites) add_total.append( add.reduce(add.reduce( abs( x.add.ThisArray ) ) ) ) # Uncomment next couple lines to check whether flow rule is mass preserving # (Should always be because of how Automata calls flow rules) # mass_total = add.reduce( add.reduce( x.cur.ThisArray ) ) # print mass_total # Display state if( t % frame_rate == 0 or t == nIter) or t == 0: SpDisplayState(state,nSites,screen,CellSize) if t == nIter: pygame.image.save(screen, "Final_State.bmp") # diff = x.cur.ThisArray - x.curSeed.ThisArray # SpDisplayState(diff,nSites,screen,CellSize) if t == 0: pygame.image.save(screen, "Initial_State.bmp") pygame.display.flip() t += 1 if (t % nSteps) == 0: #Calculate how many iterations per second we are managing print "iteration number: " + str(t) t1 = time.clock() print "Iterations per second: ", float(nSteps) / (t1 - t0) t0 = t1 if (t % hgStep) == 0 and disp_hg == True: # save the state for later histogram analysis saved_states.append( zeros( (nSites,nSites),int ) ) for i in range(0,nSites): for j in range(0,nSites): saved_states[ int(t/hgStep)-1 ][i][j] = state[i][j] if make_movie: fname = "mass_movie.txt" x.cur.save(fname) #appends the state to a text file return state #ACTUALLY CALL SOME FUNCTIONS #Ask user which type of flow we are implementing: flow_type = 0 num_flow_types = 9 while flow_type == 0: print "1 Nearest Neighbor(Basic)" print "2 Radius = 2 Neighborhood (Experimental)" print "3 Radius = sqrt(2) Neighbor (Experimental)" print "4 Reverse Nearest Neighbor (Expel mass)" print "5 Radius = 2 Neighborhood (Both Sides)" print "6 Radius = sqrt(2) Neighborhood (Both Sides)" print "7 Nearest Neighbor with Wave" print "8 Radius = 2 Neighborhood (Both Sides) with Wave" print "9 Radius = sqrt(2) Neighborhood (Both Sides) with Wave" try: flow_type = int(raw_input("Which type of flow?(Enter number)...")) except: print "Invalid entry. please enter an integer" flow_type = 0 if (flow_type <= 0) or (flow_type > num_flow_types): print "Invalid number. try again." flow_type = 0 data_run_list = [] #c0 = 0.3 c_list =[] fit_list = [] c_multi_list = [] fit_multi_list = [] if mode == 2: nTimes = 20 for h in range(0,nTimes): for g in range(0, nRuns): c = c0*(h+1) k = k0 w = w0 params = (c,k,w) dorandom = 42 + g saved_states=[] add_total = [] print "RUN NUMBER " + str(g+1) state = One_Run(nSites,flow_type,params) data_run_list.append(add_total) # Save information about how much mass was transferred # each iteration during this run # CALL ANALYSIS TOOLS massHigh = 0.95 massLow = 0.05 A_highlow(nSites,massHigh,massLow,state) if disp_hg: A_hist() # Mass density vs time -- NOT WORKING CURRENTLY if ((disp_settle_time) or (disp_fit)) or (disp_multi_run_avg): xlist, settle_time = A_make_xlist(nIter,nSites,add_total) if disp_settle_time: A_settle(xlist,add_total,settle_time) # Prints settle time number and plots mass transfer vs time if disp_fit: fit = A_fit_exponential(xlist,add_total) # Fits an exponential to a mass transfer vs time plot if disp_multi_run_avg: multi_data = A_multi_run_avg(data_run_list,disp_all_runs) # averages mass transfer vs time plots over many runs # can choose to display plots for all runs being averaged over # by setting disp_all_runs=True in parameter settings at top of file fit_multi = A_fit_exponential(xlist,multi_data) fit_multi_list.append(fit_multi) c_multi_list.append(c) wait_for_input = raw_input("type anything to continue...may desire to close window?...") f1 = [] f0 = [] for i in range(0,len(fit_multi_list)): f1.append(fit_multi_list[i][0][1]) f0.append(fit_multi_list[i][0][0]) plot(c_multi_list,f1,'ob') plot(c_multi_list,f0,'or') show() wait_for_input = raw_input("type anything to close windows...") else: # older version without averaging over multiple runs at each value of c for g in range(0, nRuns): c = c0 + g*(0.99-c0)/nRuns if mode != 1: c = c0 k = k0 w = w0 params = (c,k,w) saved_states=[] add_total = [] print "RUN NUMBER " + str(g+1) state = One_Run(nSites,flow_type,params) data_run_list.append(add_total) # Save information about how much mass was transferred # each iteration during this run # CALL ANALYSIS TOOLS massHigh = 0.95 massLow = 0.05 A_highlow(nSites,massHigh,massLow,state) if disp_hg: A_hist() # Mass density vs time -- NOT WORKING CURRENTLY if ((disp_settle_time) or (disp_fit)) or (disp_multi_run_avg): xlist, settle_time = A_make_xlist(nIter,nSites,add_total) if disp_settle_time: A_settle(xlist,add_total,settle_time) # Prints settle time number and plots mass transfer vs time if disp_fit: fit = A_fit_exponential(xlist,add_total) # Fits an exponential to a mass transfer vs time plot fit_list.append(fit) c_list.append(c) if disp_multi_run_avg: A_multi_run_avg(data_run_list,disp_all_runs) # averages mass transfer vs time plots over many runs # can choose to display plots for all runs being averaged over # by setting disp_all_runs=True in parameter settings at top of file A_fit_exponential(xlist,data_run_list) if mode == 1: wait_for_input = raw_input("type anything to continue...may desire to close window?...") f1 = [] f0 = [] for i in range(0,len(fit_list)): f1.append(fit_list[i][0][1]) f0.append(fit_list[i][0][0]) plot(c_list,f1,'ob') plot(c_list,f0,'or') show() wait_for_input = raw_input("type anything to close windows...")