<html><head><meta name="color-scheme" content="light dark"></head><body><pre style="word-wrap: break-word; white-space: pre-wrap;"># -*- coding: utf-8 -*-
"""
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-2015 University of Toronto
    Modifications:  5 May 2015 by David Bailey
                        -   Update to Python 3 style print.
                        -   Tidy up plots.
                    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 &lt;dbailey@physics.utoronto.ca&gt;
                            (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.
Time is measured in steps, which should roughly correspond to vibration
    cycles. We assume 10 Hz oscillations, so each time step is 0.1 seconds.

    """

# Use Python 3 compatible print and literals
from __future__ import print_function
from __future__ import unicode_literals

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 = 20    # Must be an integer

# Time in seconds corresponding to each random walk step
dt = 0.1    # Assume 10 Hz oscillations.

N_measured  = 36  # 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 (&gt; 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.50],m*[0.0],m*[0.50]])
# Check that probabilities add up to 1.0, within 4 * machine tolerance/epsilon
if abs(numpy.sum(p) - m) &gt; 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 &gt;= 0) and numpy.all(x &lt;= 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) &gt;= 0) :
            x = y
            t[i] += dt   # 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 &amp; Fill
#   cumulative=-1, normed=True: gives survival histogram,
#   normed=True: normalizes histogram area to unity, e.g. gives probability
#	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, normed=True)
pyplot.ylabel('Fraction\nsurviving')

# Unknotting Probability Histogram
# 3 rows, 1 column; subplot 2
R_vs_t = fig.add_subplot(312)
# Create &amp; Fill
#     log=True:    log y axis; this helps see single event entries
#     weighted to give unknotting probability per second
contents_R, bins_R, patches = R_vs_t.hist(t, bins=nbins, weights = [dt/N]*N,
                                     range=(0,tmax), log=True)
pyplot.ylabel('Unknotting probability\nper second')

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)
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\nAverages')
pyplot.xlabel('Time (seconds)')

# Note Histogram bins and contents
note(logfile,"") # Insert an empty line
note(logfile,"HISTOGRAMS")
note(logfile,"* Unknotting times:         Sample Average:")
note(logfile,"    Bin Edges, Contents,"+\
             "      Bin Edges, Contents")
bin_string = len(bins_R)*[None]
for i in range(len(bins_R)-1):
    bin_string[i] = "   {0:10.4g}, {1:7.0f},     {2:10.4g}, {3:7.0f}".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,"Time (in seconds) assumed per step      = ", dt)
note(logfile,"Number of samples generated             = ", N_samples)
note(logfile,"      Each sample consisted of            ", N_measured,
                                "unknotting measurements")
note(logfile,"Total number of unknotting measurements = ", 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.
"""
</pre></body></html>