""" knot_random_walk.py Simulates a knot in a beaded chain on a vibrating plate, assuming the movement of the crossing points in the knot are constrained random walks. Copyright (c) 2011,2012 University of Toronto Modifications: 2 June 2012 by David Bailey - Fixed asymmetry in starting position for even number of beads. Original Version: 9 July 2011 by David Bailey Contact: David Bailey (http://www.physics.utoronto.ca/~dbailey) License: Released under the MIT License; the full terms are this license are appended to the end of this file, and are also available at http://www.opensource.org/licenses/mit-license.php In this simulation, random walk positions are not allowed to cross, but the microscopic (and hence computational) details of this constraint are not clear. For now we'll only consider the simplest (computational) possibility: the probabilities for moving left/centre/right don't change when random walk positions are near to each other, but if two random walk positions try to cross then the step is rejected and ignored as if it never happened. This may not be an adequate model for the effective "forces" between the crossing points.""" import time ticks = [] # timer for speed optimization ticks.append(time.clock()) import sys, datetime from matplotlib import pyplot import numpy import scipy.stats # number of dimensions of random walk # i.e. number of initial crossing = points in knot m = 3 # Total allowed range of steps, e.g. the length of the chain in steps. # This length is the actual chain length minus the mininum knot size, N0. L = 1 # Must be an integer N_measured = 30 # Number of iterations per sample, i.e. walks generated # This is effectively the number of unknotting times # measured for this chain length. N_samples = 100 # Number of samples generated # Total number of walks to generate # To get good statistics requires large N_samples (> 1000), which will # take a few minutes to run. The default settings (m=3, L=25, N_measured=30, # N_samples=100), takes 19 seconds on a 2010 low-end MacBook Pro. N = N_measured*N_samples # Probability of moving left, not moving at all, or moving right p = numpy.array([m*[0.5],m*[0.0],m*[0.5]]) # Check that probabilities add up to 1.0, within 4 * machine tolerance/epsilon if abs(numpy.sum(p) - m) > 4*numpy.finfo(numpy.float).eps : print "Probabilities (p) don't add up to 100%." print p sys.exit("End Script") t = [] # time counter end = [] # tracks which end unknots ticks.append(time.clock()) for i in range(N) : # Starting position for the random walk is set to middle of chain for all # dimensions (i.e. all crossing points). If the chain has an even number # of beads, i.e. a non-integer midpoint, the starting point alternates # between the beads above and below the midpoint (see "-i%2")/ x = numpy.array(m*[(L-i%2)/2]) t.append(0) end.append(0) # Loop as long as all random walks on chain. while numpy.all(x >= 0) and numpy.all(x <= L-1) : ran = numpy.random.random((1,m)) step = numpy.floor(ran-p[0]) + numpy.floor((ran+p[2])) y = step+x # Crossing points are not allowed to cross, so they must remain # in order; if not, try again. if numpy.all(numpy.diff(y) >= 0) : x = y t[i] += 1 # Increment time counter # End is 1 if crossing 1 unknots, 2 if crossing 2 unknots, and 3 if both # unknot simultaneously. Note that if p[1]=0, then the first and last # crossing points will always be separated by an even number of beads, # so if L is even, the chain can never unknot simultaneously at both # ends. end[i] = 2*int(y[0][2]==L) + int(y[0][0]==-1) ticks.append(time.clock()) # Monitor how long simulation iterations took # Plot fig = pyplot.figure( facecolor='none') # initialize figure, no background colour # Number of bins in histograms (fewer bins if less iterations) nbins = min(40, N/4) # Maximum time in histograms (rounded to nearest 10) tmax = 10*int(1+numpy.max(t)/10) # Survival Probability Plot # 3 rows, 1 column; subplot 1 S_vs_t = fig.add_subplot(311) # Create & Fill # cumulative=-1, normed=True: gives survival histogram, # i.e. each bin gives the number of knots surviving longer than the bin time S_vs_t.hist(t, bins=nbins, range=(0,tmax), cumulative=-1) pyplot.ylabel('Fraction surviving') # Unknotting Probability Histogram # 3 rows, 1 column; subplot 2 R_vs_t = fig.add_subplot(312) # Create & Fill # normed=True: normalizes histogram area to unity, e.g. gives probability # log=True: log y axis; this helps see single event entries contents_R, bins_R, patches = R_vs_t.hist(t, bins=nbins, range=(0,tmax), log=True) pyplot.ylabel('Unknotting probability') ticks.append(time.clock()) # Monitor how long plotting took # A little function that writes output to both the terminal and a log file def note(logfile, *text): output_string = "" for item in text: output_string += str(item)+" " print(output_string) logfile.write(output_string+"\n") # Output file names are created from current date and time date_time = datetime.datetime.now().strftime("%Y_%m_%d_%H_%M_%S") logfile = open('krw_log__'+ date_time +'.txt',"w") plotfile = 'krw_plot_'+ date_time +'.png' # accepted formats are eps, ps, png, pdf, svg # Calculate average, rms, and most probable values for unknotting times tau = numpy.average(t) std = numpy.std(t)/tau peak = (bins_R[contents_R.argmax()]+bins_R[contents_R.argmax()+1])/2 # Uncertainty in the averages of the samples is now calculated sample_taus = numpy.average(numpy.split(numpy.array(t),N_samples),axis=1) sample_std = numpy.std(numpy.split(numpy.array(t),N_samples),axis=1) taus_vs_t = fig.add_subplot(313) contents_T, bins_T, patches = taus_vs_t.hist(sample_taus, bins=nbins, range=(0,tmax), log=True) pyplot.ylabel('Sample Averages') pyplot.xlabel('Time (in steps, i.e. arbitrary units)') # Note Histogram bins and contents note(logfile,"") # Insert an empty line note(logfile,"HISTOGRAMS") note(logfile,"* Unknotting times: Bin Edges, Bin Contents,"+\ "Sample Average: Bin Edges, Contents") bin_string = len(bins_R)*[None] for i in range(len(bins_R)-1): bin_string[i] = "{0:10.4g}, {1:7d}, {2:10.4g}, {3:7d}".format( bins_R[i],contents_R[i],bins_T[i],contents_T[i]) note(logfile,bin_string[i]) # Print timing information note(logfile,"") # Insert an empty line note(logfile,"TIMING INFORMATION") note(logfile,"Initialization took : ", ticks[1]-ticks[0],"seconds") note(logfile,"Loops took : ", ticks[2]-ticks[1], "seconds") note(logfile,"Plotting took : ", ticks[3]-ticks[2], "seconds") # Print settings note(logfile,"") # Insert an empty line note(logfile,"PARAMETER SETTINGS") note(logfile,"Chain Length = ", L) note(logfile,"Crossing Points in Knot = ", m) note(logfile,"Number of samples generated = ", N_samples) note(logfile,"Number of iterations per sample = ", N_measured) note(logfile,"Total number of iterations = ", N) note(logfile,"Probability of moving left, not moving, or moving right : ", p[0][0],p[1][1],p[2][2]) # Print Statistical Information note(logfile,"") # Insert an empty line note(logfile,"STATISTICS") note(logfile,"* All Generated Unkotting Times") note(logfile," Average = ", tau) note(logfile," Standard Deviation / Average = ", std) note(logfile," Most Probable / Average = ", peak/tau) note(logfile," Median / Average = ", numpy.median(t)/tau) note(logfile," Unknotting End [L, R, Both] = ", numpy.bincount(end)[1:4] ) note(logfile,"* Samples") note(logfile," Median / Average = ", numpy.median(sample_taus)/tau) note(logfile," Standard Deviation / Average = ", numpy.std(sample_taus)/tau) note(logfile," Naive Expectation = ", std/numpy.sqrt(N_measured-1)) note(logfile," -1 sigma (16% cdf) / Average = ", scipy.stats.scoreatpercentile(sample_taus, 15.865)/tau-1) note(logfile," +1 sigma (84% cdf) / Average = ", scipy.stats.scoreatpercentile(sample_taus, 84.135)/tau-1) logfile.close() pyplot.savefig(plotfile, dpi=(300)) # save the plot pyplot.show() # show the plot # End of knot_random_walk.py """ Full text of MIT License: Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. """