Ulam spiral

So, bored in a math talk Ulam starts writing integers in a spiral and circling primes. He soon notices something weird:


20 layers of Ulam sprial. The red dot is the origin and it rotates counter clockwise starting with 1 to the right of the red dot

Here is a bigger picture:


200 layers of Ulam spiral. That is a 400×400 grid.

Here is the sage code that I wrote to plot these:

def ulam_spiral(n):
    # n is the number of layers
    # output is list of coordinates of points corresponding to prime numbers in the spiral
    layers = n
    P = []
    count = 0
    for k in xsrange(layers):
        for s in range(4):
            for i in range(2*k):
                count += 1
                if s == 0:
                    if count.is_prime():
                if s == 1:
                    if count.is_prime():
                if s == 2:
                    if count.is_prime():
                if s == 3:
                    if count.is_prime():

## plot them
n = 20 #number of layers
p = points(ulam_spiral(n),size=100000/n^2,axes=False)
q = point((0,0),color='red')

You can read more detailed version of the story here. Then there came the Sacks spiral, and other variations. Here is a nice article on number spirals and some python code.

Here are all the primes less than 50000:


Primes less than 50000 in Sacks spiral.

And here is the sage code to generate the plot:

def sakcs_spiral(n):
    # all primes less than n
    # output is list of coordinates of points corresponding to prime numbers in the spiral
    layers = n
    P = []
    for i in xsrange(n):
        if i.is_prime():
            x = -cos(sqrt(i)*2*pi)*sqrt(i)
            y = sin(sqrt(i)*2*pi)*sqrt(i) 

## plot them
n = 50000 #number of layers
p = points(sakcs_spiral(n),size=8,axes=False)
q = point((0,0),color='red')

And of course here is a numberphile video on the topic.


2017 Fields Medal Symposium in honour of Martin Hairer

I recently participated at the Fields Medal Symposium held at the Fields Institute in Toronto, in honour of Martin Hairer, one of the recipients of the medal in 2014. His work revolves around stochastic partial differential equations.


I am amazed by how good the speakers could talk about their deep deep ideas that even a non-specialist could understand their main ideas, something that doesn’t usually happen in math conferences! The videos will be available here. Here are a few snapshots from some of the talks:

Self-avoiding walk, spin systems, and renormalisation
Gordon Slade, Univeristy of British Columbia

2017-10-16 09.40.47

Gordon was interested in self-avoiding walks and their critical exponents. He started his talk by mentioning the results in enumerating such walks.

2017-10-16 09.34.47

One interesting fact is that for exponents less than the critical ones, the walk tends to follow a geodesic curve, but for exponents less than the critical ones, the walk becomes like a Brownian motion. This came up in a couple of other talks. Here is a list of some experimentally found critical values:2017-10-16 09.40.42

A nice side thing that he showed was the size of a self-avoiding walk. In the following picture the left is 10^8 steps of a self-avoiding walk and on the right (the tiny thing) is 10^8 steps of a simple random walk!

2017-10-16 09.42.17

Then he gave some theorems and talk about continuous-time weakly self-avoiding walks, and mentioned some problems with physicist’ approach.

2017-10-16 09.50.14.jpg

Then he talked about long-range self-avoiding walks (from 1972), and mentioned that they don’t converge to a Brownian motion for some values of their parameters, they rather approach a stable process. The he connected it to the Ising model.

2017-10-16 09.57.15.jpg

2017-10-16 09.59.03.jpg

2017-10-16 10.27.33.jpg

Statistical motion of a convex body in a rarified gas close to equilibrium
Laure Saint-Raymond, École Normale Supérieure


Laure started by describing how (sphere looking) atoms of the same size interact with each other:

2017-10-16 11.03.04.jpg

And then she put a big (non-spherical) molecule in there, which suddenly made the dynamics much much more complex:

2017-10-16 11.05.43.jpg

One of the audiences asked if the big molecule (if too big) can be considered as a plane, and hence as an approximation of a huge huge sphere, which didn’t really get a response. Then she talked about the weak coupling limit, and showed the limit equations, and talked about the collision trees.

2017-10-16 11.21.40.jpg

Yes, it is a crazy idea to assume that the collisions form a tree. But that was the main idea, where multiple collisions happen, things go wrong in the calculations, they eventually explode into infinity, the art is to show enough things cancel each other on average, hence the error is small!

2017-10-16 11.24.00.jpg

2017-10-16 11.28.28.jpg

[I’ll add rest of the talks here later]

Listen to your pc music through your phone

So I have this bunch of music on my external hard drive that is connected to my PC at home and I need to get over there to play it while I’m at home, but that’s usually a hassle, well, to some extent. So, here is a partial solution to it. I ssh into my PC from my phone and play the music, and listen to it from my bluetooth loud speakers or headphone, while doing my other things at home (specially in the morning having breakfast).

Starting an ssh server on linux:

service sshd start

This should open the port 22. Check it by:

netstat -ant | grep 22
tcp     0     0*           LISTEN 
tcp     0     0 *           LISTEN 
tcp     0     0     ESTABLISHED
tcp     0   468  ESTABLISHED
tcp     0     0     ESTABLISHED
tcp6    0     0     :::22                :::*                LISTEN

If I want the server to start on boot I can do this:

chkconfig sshd on

Then I need to figure out what is my PC’s local IP address so that I can ssh into it. A quick way to do it is by looking at the output at ifconfig and going through to find the corresponding IP. Right now mine looks like this:

enp0s25: flags=4099<UP,BROADCAST,MULTICAST>  mtu 1500
        ether **:85:**:0b:**:6c  txqueuelen 1000  (Ethernet)
        RX packets 0  bytes 0 (0.0 B)
        RX errors 0  dropped 0  overruns 0  frame 0
        TX packets 0  bytes 0 (0.0 B)
        TX errors 0  dropped 0 overruns 0  carrier 0  collisions 0
        device interrupt 19  memory 0xf0500000-f0520000  

lo: flags=73<UP,LOOPBACK,RUNNING>  mtu 65536
        inet  netmask
        inet6 ::1  prefixlen 128  scopeid 0x10
        loop  txqueuelen 1  (Local Loopback)
        RX packets 468  bytes 36452 (35.5 KiB)
        RX errors 0  dropped 0  overruns 0  frame 0
        TX packets 468  bytes 36452 (35.5 KiB)
        TX errors 0  dropped 0 overruns 0  carrier 0  collisions 0

virbr0: flags=4099<UP,BROADCAST,MULTICAST>  mtu 1500
        inet  netmask  broadcast
        ether 00:00:00:00:00:00  txqueuelen 1000  (Ethernet)
        RX packets 0  bytes 0 (0.0 B)
        RX errors 0  dropped 0  overruns 0  frame 0
        TX packets 0  bytes 0 (0.0 B)
        TX errors 0  dropped 0 overruns 0  carrier 0  collisions 0

wlp48s0: flags=4163<UP,BROADCAST,RUNNING,MULTICAST>  mtu 1500
        inet  netmask  broadcast
        inet6 ****:6477:****:73e2:****:d8ff:****:e8b2  prefixlen 64  scopeid 0x0
        inet6 fe80::****:d8ff:****:e8b2  prefixlen 64  scopeid 0x20
        ether **:16:**:36:**:b2  txqueuelen 1000  (Ethernet)
        RX packets 970380  bytes 1163277882 (1.0 GiB)
        RX errors 0  dropped 0  overruns 0  frame 0
        TX packets 469268  bytes 56819419 (54.1 MiB)
        TX errors 0  dropped 0 overruns 0  carrier 0  collisions

and that second line under wlp48s0 will be my ip. Since I’m using a dynamic ip assignment, this might change every time my PC connects to network, so I have written a script that runs at startup and looks this up (and my public ip) every one hour and writes it to a file in my dropbox so that I can access it from my phone! Here it is:

while true; do
    ifconfig wlp48s0 | grep cast > ~/Dropbox/opt/myip
    wget http://ipecho.net/plain -O - -q >> ~/Dropbox/opt/myip
    sleep 1h

ssh clients for android and some shortcuts

The next step is to run a ssh client on my phone. Currently I’m using JuiceSSH, but I think I’ll be moving to ConnectBot or some other free apps soon, or probably I should just run a terminal on my phone and ssh from there. After connecting to on port 22 I can browse through the files on my PC and run them. Since my external hard drive has a long name and is mounted somewhere that requires me to remember a lot of things I just added an alias to the .bashrc file so that I can directly go to my music folder:

alias cdmusic="cd 'path/to/my music/folder/'"

Then I find the folder I wanna play and use mplayer to play it. Since I had this other script to personalize my mplayer a little bit, I keep using that to play all the files in a folder in a fashinable way:

alias play="mplayer -msgcolor -cache 500 20 *.*"

You can find the full list in the mplayer man page: man mplayer (See the bottom of the post). And if that’s too much, here is short list of  essential ones from mplayer --help:


You could also turn your phone into a remote control for your video playback this way! I should still find good ways to do this over internet, and stream the music to my phone rather than to a bluetooth device! I mean, I don’t want an entire personal cloude running on my PC, but why not!


bluetooth and audio settings

In order to have a bluetooth device automatically connect to your PC you can do this:


This will show a list of all the bluetooth devices, their MAC addresses, and their names. Then I can “trust” them by typing

trust de:vi:ce:ma:ca:dd:re:ss

This makes my device to connect to the PC automatically when it’s on (read the help for complete set of commands and options). Now to change the sound output to each device I can do this, assuming pulse audio is in use:

pacmd list-cards

I get a list of audio devices and lots of descriptions. Somewhere there you should be able to find your bluetooth device and an “index” related to it. mine is “3” right now (but it changes each time the device disconnects and connects) for my sound box, and “0” for the PC speakers. Then all I need to do is:

echo "set-default-sink 3" | pacmd


echo "set-default-sink 0" | pacmd

to go back to PC.

Partial Update:

Apparently there are ways to assign a static IP to your computer, for example see here, but there are some problems that needs to be worked out afterwards, e.g. disconnecting from internet etc. Fortunately, my modem can set a static IP address for each MAC address that I could manage through its settings. That is, I don’t need to read the IP of my PC each time.


  • figure out how to find the correct local ip automatically
  • add the above to my script!
  • figure out how to find the correct index for the sound output device automatically.

to run a script so that you can change the current directory run it like this:

. myscript



Here is a full list of keyboard shortcuts for mplayer that one can use (even via ssh):

keyboard control
              LEFT and RIGHT
                   Seek backward/forward 10 seconds.
              UP and DOWN
                   Seek forward/backward 1 minute.
              PGUP and PGDWN
                   Seek forward/backward 10 minutes.
              [ and ]
                   Decrease/increase current playback speed by 10%.
              { and }
                   Halve/double current playback speed.
                   Reset playback speed to normal.
              < and >
                   Go backward/forward in the playlist.
                   Go forward in the playlist, even over the end.
              HOME and END
                   next/previous playtree entry in the parent list
              INS and DEL (ASX playlist only)
                   next/previous alternative source.
              p / SPACE
                   Pause (pressing again unpauses).
                   Step forward.  Pressing once will pause movie, every consecutive press will play one frame and then go into pause mode again (any other key unpauses).
              q / ESC
                   Stop playing and quit.
                   Stop playing (and quit if -idle is not used).
              + and -
                   Adjust audio delay by +/- 0.1 seconds.
              / and *
                   Decrease/increase volume.
              9 and 0
                   Decrease/increase volume.
              ( and )
                   Adjust audio balance in favor of left/right channel.
                   Mute sound.
              _ (MPEG-TS, AVI and libavformat only)
                   Cycle through the available video tracks.
              # (DVD, Blu-ray, MPEG, Matroska, AVI and libavformat only)
                   Cycle through the available audio tracks.
              TAB (MPEG-TS and libavformat only)
                   Cycle through the available programs.
                   Toggle fullscreen (also see -fs).
                   Toggle stay-on-top (also see -ontop).
              w and e
                   Decrease/increase pan-and-scan range.
                   Toggle OSD states: none / seek / seek + timer / seek + timer + total time.
                   Toggle frame dropping states: none / skip display / skip decoding (see -framedrop and -hardframedrop).
                   Toggle subtitle visibility.
              j and J
                   Cycle through the available subtitles.
              y and g
                   Step forward/backward in the subtitle list.
                   Toggle displaying "forced subtitles".
                   Toggle subtitle alignment: top / middle / bottom.
              x and z
                   Adjust subtitle delay by +/- 0.1 seconds.
              c (-capture only)
                   Start/stop capturing the primary stream.
              r and t
                   Move subtitles up/down.
              i (-edlout mode only)
                   Set start or end of an EDL skip and write it out to the given file.
              s (-vf screenshot only)
                   Take a screenshot.
              S (-vf screenshot only)
                   Start/stop taking screenshots.
                   Show filename on the OSD.
                   Show progression bar, elapsed time and total duration on the OSD.
              ! and @
                   Seek to the beginning of the previous/next chapter.
              D (-vo xvmc, -vo vdpau, -vf yadif, -vf kerndeint only)
                   Activate/deactivate deinterlacer.
              A    Cycle through the available DVD angles.

              (The following keys are valid only when using a hardware accelerated video output (xv, (x)vidix, (x)mga, etc), the software equalizer (-vf eq or -vf eq2)  or
              hue filter (-vf hue).)

              1 and 2
                   Adjust contrast.
              3 and 4
                   Adjust brightness.
              5 and 6
                   Adjust hue.
              7 and 8
                   Adjust saturation.

              (The following keys are valid only when using the quartz or corevideo video output driver.)

              command + 0
                   Resize movie window to half its original size.
              command + 1
                   Resize movie window to its original size.
              command + 2
                   Resize movie window to double its original size.
              command + f
                   Toggle fullscreen (also see -fs).
              command + [ and command + ]
                   Set movie window alpha.

              (The following keys are valid only when using the sdl video output driver.)

                   Cycle through available fullscreen modes.
                   Restore original mode.

              (The following keys are valid if you have a keyboard with multimedia keys.)

                   Stop playing and quit.
              PREVIOUS and NEXT
                   Seek backward/forward 1 minute.

              (The following keys are only valid if you compiled with TV or DVB input support and will take precedence over the keys defined above.)

              h and k
                   Select previous/next channel.
                   Change norm.
                   Change channel list.

              (The following keys are only valid if you compiled with dvdnav support: They are used to navigate the menus.)

              keypad 8
                   Select button up.
              keypad 2
                   Select button down.
              keypad 4
                   Select button left.
              keypad 6
                   Select button right.
              keypad 5
                   Return to main menu.
              keypad 7
                   Return to nearest menu (the order of preference is: chapter->title->root).
              keypad ENTER
                   Confirm choice.

              (The following keys are used for controlling TV teletext. The data may come from either an analog TV source or an MPEG transport stream.)

                   Switch teletext on/off.
              Q and W
                   Go to next/prev teletext page.

       mouse control
              button 3 and button 4
                   Seek backward/forward 1 minute.
              button 5 and button 6
                   Decrease/increase volume.

       joystick control
              left and right
                   Seek backward/forward 10 seconds.
              up and down
                   Seek forward/backward 1 minute.
              button 1
              button 2
                   Toggle OSD states: none / seek / seek + timer / seek + timer + total time.
              button 3 and button 4
                   Decrease/increase volume.

Non-teaching aspects of teaching at University of Calgary

I’ve been a math postdoc at University of Calgary for two years now, and I’ll be here for at least one more year. The teaching environment is very pleasant, and the resources are abundant. Nonetheless, there are a few non-teaching aspects of it that bothers me. The way that the teaching is handled in the mathematics department for postdocs is that we are hired as a sessional instructors, with a no strings attached for the teaching. A course is assigned to us, with a fixed number for the pay. There are no benefits coming with it, and no pension associated to it. However, after 1 year of being a sessional instructor we can apply for $200 reimbursement in professional development activities.

During the first semester that I was here, I didn’t have a teaching duty, but I was offered to teach a course for extra pay, which I took it. However, by the end of semester I realized that I haven’t been paid for my teaching at all. It turned out that it was an error of the accounting personnel of the department which hadn’t put me on payroll (teaching salary is separate from my postdoc salary). Eventually, I got a lump sum for my Fall semester teaching in the following February. A year later when filing for taxes, I realized that getting paid all of it in February wasn’t a good idea, since it resulted in my paying about an extra $1000 in taxes for that year. The case isn’t resolved yet. Well, to be honest, I still don’t know how to even go after fixing this!

Let me tell you a little about how our department counts teaching. Each “regular” course (whatever that means) in a semester counts as one “half course equivalent”, which means a sessional instructor gets paid about $6250 for teaching it, and it counts as “one” teaching duty for a full time faculty member. It turns out that the department counts teaching, for example, a Calculus 1 course (Math 249 or Math 265) with up to 180 students as 1 “half course equivalent” for the faculty, and up to 240 students as 1.5, and up to 360 students as 2 “half course equivalents”. That is, if you are a full time faculty here with teaching load of 2 for the Fall semester, all you need to do is to teach a Calculus 1 course with 360 students in it. The weird part of it is that you don’t get to know any of these things by default, and I couldn’t find any documents in the department talking about this, maybe for lack of trial!

What I have experienced here as a postdoc (sessional instructor) though, is that all of the courses that I’ve taught with 240 students still counts as 1 “half course equivalent”, that is I get paid about $6250 for it, and also I am supposed to run one of the labs for whatever course I teach. I haven’t had any concerns regarding this so far, until there came up a situation about my teaching duties next year. It all started with a simple confusion, in my postdoc contract I was supposed to be teaching one course as part of my responsibilities. So, as I said above, the way that this is handled is that they reduce my negotiated pay by around $6250 and then hire me as a sessional instructor and then pay that amount to me via a separate series of pay checks. The confusion was that the associate head of teaching of the department thought I am supposed to teach two courses as part of my duties, he then scheduled me for one large section (360 students) of Calculus 1 (Math 265), which my new postdoc supervisor objected, as it will consume too much time from me, hence I won’t be able to spend enough time on my research. After a few emails back and forth between people (I wasn’t involved here) I was told that I will be teaching the large section and get paid for two half course equivalents, and even as an “exception” they’ll drop the lab for me! It sounded like an awesome deal. So I agreed.

Until last week that they sent me the contract to be signed, which said I’m getting paid for one half course equivalent. Upon inquir”es” it turned out that the department has decided that “it is pretty fair” if I teach a large section instead of a small section + lab. I’ll probably take the offer, but it leaved me with this question that why there are such double standards segregating the faculty and sessionals, and why the policies in the department are not transparent.

So you wanna move to Canada?

I’ve recently moved to Canada and applied for permanent residency through the Express Entry system. If your situation is similar to me (you have a higher education degree, a few years of part time — yes TA and RA count — or full time experience, you can speak English fair enough) then you might be eligible to enter Canada the same way. To do so, you don’t need to be present in Canada, and you don’t have to remain in Canada the entire time if you don’t think that’s something you’ll want to keep.

There are several ways for moving to Canada, but I know of only one, and I only know a few things about it which is obviously about my case! The way that I applied for immigration to Canada was through the Express Entry process, a somewhat new method which started some time in 2014 I guess. The idea is based on your education, language skills, work experience, age, marital status etc you get a score and enter a pool. Roughly every two weeks the government draws the top people from the pool and invites them to apply for the Express Entry program. That basically means if you have claimed everything correctly, and if you provide all the required documents, and if you have a clean background, then almost surely you’ll get to become a Canadian permanent resident (a five year non-bounding status to live in Canada).

What are your chances to be eligible for Express Entry?

First get a rough estimate of your rating here: http://www.cic.gc.ca/english/immigrate/skilled/crs-tool.asp

Then head over here http://www.cic.gc.ca/english/express-entry/rounds.asp to see the minimum ratings drawn from the pool in the last draw and how others are doing in the pool!

For example, as of my writing, if your score is above 451, andif  in the next round they draw 800 people, you’ll be selected.

What do you need to enter the pool?

Here are the first things you need to do:

  • Take a general IELTS exam at your earliest time. When you are signing up, it asks you what do you need it for, choose for immigration purposes.
  • Evaluate your degrees (all degrees which are not from a Canadian institute). There are a few companies that the Canadian government work with, out of which I only remember WES. Create an account there and fill the forms to tell them what degrees you want to evaluate. Then it will tell you how much you need to pay. After paying it, it will give you a code that you need to keep, and then it will tell you what documents shall you send to them and in what condition. Make sure that the code is mentioned on all the documents.

These are all you need to get into the pool. Do these first and enter the pool as soon as possible. Basically, you claim all the rest and when you get drawn you provide documents to support your claims. So, then start working on the rest of the documents:

  • Work experience. Gather letters from your bosses (preferably your direct supervisors for all the previous jobs, but HR also works) to show your work history. For some countries like Iran you might need to provide your insurance history too. Here is what you need to be in the letter:
    • written on company letterhead,
    • signed by the responsible supervisor,
    • have the printed name and title of the responsible supervisor beneath the signature,
    • show the company’s full address, telephone and fax numbers, e-mail and website addresses,
    • stamped with the company’s official seal (if applicable).
    • the specific period of your employment with the company,
    • the positions you have held during the period of employment and the time spent in each position, and the status of the current positoin
    • main responsibilities and duties in each position,
    • total annual salary plus benefits,
    • the number of hours worked per week

      Your TA and RA also count.  For TA your department head probably is the supervisor. For RA your advisor. Ask them to be very specific about dates, and number of hours, etc. As specific as they can be! Nonetheless all this information is in your contracts.

  • Police records: If you live in US you need one from FBI. At the minimum you want police reports for the countries that you lived in in the past 10 years for at least 6 months. If you need additional reports, they will let you know when you are drawn from the pool.
  • Proof of financial assets: you need about 13000 CND per person in the application. Though per application it might be a different situation. But if that’s the case you need to provide an official letter indicating your financial profile, including
    • list of all your bank (chequing and savings) accounts
    • the account numbers
    • dates each account was opened, and
    • the balance of each account over the past six months,
    • list all outstanding debts, i.e. loans, credit cards, etc,
    • be printed on the letterhead of the financial institution, and include your name and the contact information of the financial institution (address, telephone number and e-mail address).
    • proof of assets, like cars, houses, etc.


That’s all for now. I’ll keep adding details here.

Good luck.


Miroslav and the separating wand

In a huge graph, it is hard to determine a bottleneck by just looking at it, or doing any type of search. And by a bottleneck, I mean part of graph that is not very much connected to the other part. In different setting this could mean different things, for example, if your graph is a communication network, then a bottleneck mean the information doesn’t get to all of the nodes very easily, there will be traffic around the bottleneck.

images.duckduckgo.comOne of the celebrated results in algebraic graph theory is by Miroslav Fiedler, a Czech mathematician [1, 2], about the algebraic connectivity of a graph. It simply says the second smallest eigenvalue of the Laplacian matrix of a graph determines how connected the graph is. In particular, if it is 0, the graph is disconnected. Furthermore, it shows the eigenvector corresponding to this eigenvalue distinguishes the vertices in the two connected components. In the case that the graph is connected, this eigenvector will find bottlenecks of the graph. Below I first generate a random graph that has a visible bottleneck, then I’ll use Fiedler’s result to find it computationally.

First, let’s generate a random graph with a bottleneck. In order to do so, I’ll generate two random graphs, then I look at their disjoint union and put few random edges between them.

def random_bottleneck_graph(n1,n2,p1,p2,p):
    #n1 = number of vertices of first part
    #n2 = number of vertices of second part
    #p1 = probability of edges of first part
    #p2 = probability of edges of second part
    #p  = probability of edges between two parts
    import random
    #generate two parts randomly and put them together
    H1 = graphs.RandomGNP(n1,p1) #generate a random GNP graph for first part
    H2 = graphs.RandomGNP(n2,p2) #generate a random GNP graph for second part
    G = H1.disjoint_union(H2) #form the disjoint union of the two parts
    #add edges between two parts
    for i in range(n1): #pick each vertex in H1
        for j in range(n2): #pick each vertex in H2
            a = random.random() #generate a random number between 0 and 1
            if a > 1-p: #with probability p add an edge between the two vertices of H1 and H2

And here is a sample run:

n1 = 40 #number of vertices of first part
n2 = 60 #number of vertices of second part
p1 = .7 #probability of edges of first part
p2 = .6 #probability of edges of second part
p = .01 #probability of edges between two parts

G = random_bottleneck_graph(n1,n2,p1,p2,p)

Of course when I look at the graph, sage who has generated it knows that the graph has two main parts, and even the vertices are numbered in a way that anyone can recognize it. The left picture is generated by the above ‘show’ command, and the right picture is using the following code.


def show_random_bottleneck_graph(G,layout=None):
    EdgeColors={'red':[], 'blue':[], 'black':[]}

    for e in G.edges():
        if not e[0][0] == e[1][0]:
        elif e[0][0] == 0:

So, in order to really get a random graph with a bottleneck, I take the above graph and try to loose the information that I had from it.

def relable_graph(G,part=None):
    # input: a graph G
    # output: if part is not given, it gives a random labeling of nodes of graph with integers
    #         if part is given, then it partitions the vertices of G labeling each vertex with pairs (i,j) where i in the number of partition the vertex is in part, and j is an integer.
    V = G.vertices()
    E = G.edges()
    H = Graph({})
    map = {}
    if part == None:
        new_part = None
        import random
        for i in range(len(V)):
            map[V[i]] = i
        new_part = [ [] for p in part]
        count = 0
        for i in range(len(part)):
            for w in part[i]:
                map[w] = (i,w)
                count += 1
    for v,w,l in E:


H = relable_graph(G)[0]
H.show(layout='circular',vertex_labels=False, figsize=[10,10])

Now if I look at the graph I won’t be able to recognize the bottleneck:


Now that I truly have random graph, which I know has a visible bottleneck, let’s see if we can find it. Below I use ‘networkx’ to compute the Fiedler vector, an eigenvector corresponding to the second smallest eigenvalue of the Laplacian matrix. Then I’ll partition the vertices according to the sign of the corresponding entry on this vector:

def partition_with_fiedler(G):
    # input: a graph G
    # algorithm: using networkx finds the fiedler vector of G
    # a partitioning of vertices into two pieces according to positive and negative entries of the Fiedler vector
    import networkx
    H = G.networkx_graph()
    f = networkx.fiedler_vector(H, method='lanczos')
    V = H.nodes()
    P = [[],[]]
    for i in range(len(V)):
        if f[i] >= 0:        

Let’s see how it is partitioned:

P = partition_with_fiedler(H)
H.show(partition=P, layout='circular', vertex_labels=True, figsize=[10,10])


Well, not much can be seen yet. Now let me put the red vertices in one side and the blue vertices in the other side:

def position_bottleneck(G):
    P = partition_with_fiedler(G)
    V = G.vertices()
    new_pos = {}
    for v in V:
        if v in P[0]:
            new_pos[v] = [6*(1+random.random()),6*(2*random.random()-1)]
            new_pos[v] = [-6*(1+random.random()),6*(2*random.random()-1)]

pos = position_bottleneck(H)
H.show(partition=P, pos=pos, vertex_labels=True, figsize=[10,10])


Now it’s easy to see the bottleneck of the graph! We can even recover the original picture as follows:

(K,newP) = relable_graph(H,part=P)
K.show(partition=newP, layout='circular', vertex_labels=False, figsize=[10,10])


And here is the representation of it using the colored edges:

show_random_bottleneck_graph(K,layout='circular') #suggested layouts are "circular" and "spring"



Since I posted this, I’ve been thinking about multiway Fiedler clustering. Here is what I got so far. First I define a list N of number of vertices in each part, and a matrix A which is basically a probability matrix where entry (i,j) is the probability of an edge between part i and part j, hence it is symmetric. If you choose diagonal entries close to 1 and off-diagonal entries close to zero, then you’ll get a graph with multiple communities. I’ll shuffle the vertices around after creating the multi-community graph.

def multi_bottleneck_graph(N,A):
    #N = list of size n of number of vertices in each part
    #A = an nxn matrix where entry (i,j) shows the probability of edges between part i and part j of vertices, and n is the number of parts.
    import random
    n = len(N)
    G = Graph({})
    #add vertices for each part
    for i in range(n):
        for j in range(N[i]):
    V = G.vertices()
    #add edges
    for v in V:
        for w in V:
            if not v == w:
                a = random.random() #generate a random number between 0 and 1
                p = A[v[0],w[0]]
                if a > 1-p: #with probability p add an edge between the two vertices of H1 and H2

def relable_graph(G,part=None):
    # input: a graph G
    # output: if part is not given, it gives a random labeling of nodes of graph with integers
    #         if part is given, then it partitions the vertices of G labeling each vertex with pairs (i,j) where i in the number of partition the vertex is in part, and j is an integer.
    V = G.vertices()
    E = G.edges()
    H = Graph({})
    map = {}
    if part == None:
        new_part = None
        import random
        for i in range(len(V)):
            map[V[i]] = i
        new_part = [ [] for p in part]
        count = 0
        for i in range(len(part)):
            for w in part[i]:
                map[w] = (i,w)
                count += 1
    for v,w,l in E:


def shuffled_multi_bottleneck_graph(N,A):
    G = multi_bottleneck_graph(N,A)
    H = relable_graph(G)[0]

# run
N = [10,15,20]
A = matrix([[.9,.01,.01],[.01,.9,.01],[.01,.01,.9]])
G = shuffled_multi_bottleneck_graph(N,A)
G.show(layout='circular',vertex_labels=False, vertex_size=30, figsize=[10,10])


Then I break the graph into two pieces using Fiedler vector, compare the two pieces according their algebraic connectivities (second smallest eigenvalue of the Laplacian) and break the smaller one into two pieces and repeat this process, until I get n communities.

def partition_with_fiedler(G):
    # input: a graph G
    # algorithm: using networkx finds the fiedler vector of G
    # a partitioning of vertices into two pieces according to positive and negative entries of the Fiedler vector
    import networkx
    H = G.networkx_graph()
    f = networkx.fiedler_vector(H, method='lanczos')
    V = H.nodes()
    P = [[],[]]
    for i in range(len(V)):
        if f[i] >= 0:        

def multi_fiedler(G,n=2):
    if n == 2:
        P = partition_with_fiedler(G)
        import networkx

        P = partition_with_fiedler(G)
        H = [ G.subgraph(vertices=P[i]) for i in range(2) ]
        C = [ (networkx.algebraic_connectivity(H[i].networkx_graph(), method='lanczos')/H[i].order(),i) for i in range(2) ]
        l = min(C)[1]
        while n > 2:
            n = n-1
            temP = partition_with_fiedler(H[l])
            nul = P.pop(l)

            temH = [ G.subgraph(vertices=temP[i]) for i in range(2) ]
            nul = H.pop(l)
            temC = [ (networkx.algebraic_connectivity(temH[i].networkx_graph(), method='lanczos')/temH[i].order(),i) for i in range(2) ]
            nul = C.pop(l)
            l = min(C)[1]
# run
G.show(partition=K, layout='circular', figsize=[10,10])


Now I can separate same colour vertices to see the communities in the graph:

def show_communities(G,n):
    #G: a graph
    #n: number of communities > 1
    if n == 1:
        P = [G.vertices()]
    elif n > 1:
        P = multi_fiedler(G,n)
        raise ValueError('n > 0 is the number of communities.')
    L,newP = relable_graph(G,part=P)
    L.show(partition=newP, layout='circular', vertex_labels=False, figsize=[10,10])

# run


There is only one important thing in using networkx. It seems like the default method (tracemin) for algebraic_connectivity and for fiedler_vector has a bug. So, I used the method=’lanczos’ option there.

You might be interested to see what happens if you ask for more or less than 3 communities. Here are a couple of them:


The next question is how to decide about number of communities.

  1. M. Fiedler, Algebraic connectivity of graphs. Czechoslovak Math. J. 23(98):298–305 (1973).
  2. M. Fiedler, A property of eigenvectors of nonnegative symmetric matrices and its application to graph theory, Czech. Math. J. 25:619–637, (1975).