Algorithm.co.il

Math

Fibonacci meals

Posted on by Leave a Comment

Lately I’ve been reminded of the existence of Fibonacci meals. For a long time, I despised those, but now that I’m preparing my own food most of the time, I’m aiming for them, to save time cooking.

A Fibonacci meal is an element of a sequence defined to be:

a(0) = some food
a(1) = some food
a(n) = a(n-1) + (a(n-2) as sauce)

Very simple, and very effective. You can find those in various places, such as camps, school cafeteria, etc…

Generally this sequence does not converge, however there is a special case for cakes and desserts:

c(0) = some cake or dessert
c(1) = some cake or dessert
c(n) = c(n-1) + (c(n-2) as icing)

This sequence does converge, to the chocolate ball. (in Wikipedia, Chokladboll).

Elegant Line Clipping Algorithm

Posted on by Leave a Comment

I actually wrote about that a long time ago, in the wiki that once was algorithm.co.il. Since the wiki is long gone, and I figured this is a very elegant algorithm, I decided I should write about it.

First, some background – while writing Ascii-Plotter, I had to plot lines. To plot lines correctly I had to clip them. Now since this was a fun project, I decided I wanted to figure out how to do that myself, even though the problem is really solved. So I sat down with pen and paper, and did a little thinking, and I came up this nice algorithm.

The idea is simply this: consider a line, that is not parallel to any of the axes. This line intersects all the four lines that make up the limits of our screen. Here is a diagram:

An intersection of a line with a bounding box

Now notice that there are six points of interest: the four intersection points, and the two endpoints of the line. For each point, we find its parameter in the parametric representation of the line, where a and b are the endpoints:

x(t) = a+t(b-a)

To find out the actual endpoints of the clipped line, we sort the six points according to their parameter, and take the two innermost points (with zero based indexing, those would be 2 and 3). What if the line is outside the box you say? We check that the result points’ parameters are between 0 and 1 (which correspond to the original endpoints).

Note that for a line parallel to one of the axes, a similar method applies, but without two of the intersection points. However, it is probably easier to just calculate the clipping directly in that case.

Here is some sample code, without all the surrounding checks (for parallelism with the axes, trivial clipping etc…).
The full version is available in the code of Ascii-Plotter.

points_t = [0.0,1.0]
 
points_t.append( float(rect_bottom_left[0]-line_pt_1[0])/(line_pt_2[0]-line_pt_1[0]) )
points_t.append( float(rect_top_right[0]-line_pt_1[0])/(line_pt_2[0]-line_pt_1[0]) )
points_t.append( float(rect_bottom_left[1]-line_pt_1[1])/(line_pt_2[1]-line_pt_1[1]) )
points_t.append( float(rect_top_right[1]-line_pt_1[1])/(line_pt_2[1]-line_pt_1[1]) )
 
points_t.sort()
if points_t[2] < 0 or points_t[2] >= 1 or points_t[3] < 0 or points_t[2]>= 1:
	return None
result = [(pt_1 + t*(pt_2-pt_1)) for t in (points_t[2],points_t[3]) for (pt_1, pt_2) in zip(line_pt_1, line_pt_2)]
return (result[0],result[1]), (result[2], result[3])

Python module usage statistics – Cont.

Posted on by 3 Comments

Well, very embarrassingly for me, turns out I had a bug in my original post and code. As per Doug’s suggestion, I tried running the script I wrote on the standard library, and got results I didn’t quite believe. So I checked them, opened socket.py, and there was an “import _socket”. However module “_socket” was not listed in my results. After some digging around, (and feeling my face getting red hot…) I found the bug. It was the smallest thing: I forgot to add re.MULTILINE as an argument to re.search, so the start of my regexp, “^[ \t]*import” didn’t match except on the beginning of the file. #@&*!!! Just no words to express how I felt at that moment.

So, I fixed the bug, tested my code a bit, fixed some more (very minor) bugs, and here are the top 25 of the (hopefully correct this time!) results:

  • sys,1426
  • os,1250
  • unittest,566
  • time,446
  • re,383
  • string,321
  • types,298
  • setuptools,264
  • pkg_resources,217
  • cStringIO,184
  • zope.interface,177
  • datetime,173
  • shutil,167
  • os.path,162
  • gtk,143
  • StringIO,143
  • random,136
  • tempfile,132
  • copy,131
  • threading,128
  • distutils.core,127
  • doctest,126
  • md5,125
  • setuptools.command.e,116
  • logging,116

Except the larger numbers, the arrangement pretty much stayed the same.

(This seems to at least confirm my claim that the results were representative)
Here are the results for the standard library (/usr/lib/python2.5/):

  • sys,1309
  • os,1065
  • ctypes,588
  • re,496
  • string,493
  • types,435
  • time,374
  • numpy,297
  • warnings,254
  • os.path,204
  • cStringIO,196
  • common,185
  • math,159
  • traceback,158
  • gettext,152
  • codecs,147
  • StringIO,147
  • copy,133
  • __future__,128
  • tempfile,126
  • random,119
  • threading,117
  • unittest,108
  • numpy.testing,105
  • errno,100

These results seem different. ctypes seems to be the biggest change. Note that these results might be slanted, as some stdlib modules are implemented in C, and I didn’t test a ‘clean’ Python installation (but rather one with all my non-default modules installed).

Here is a link to the new results, the new stdlib results, and the new fixed graph. I created it using the first 325 module. The fixed R squared value is now approx. 0.99. Here is the code.

I would like to apologize to anyone who was misled by the original results, and again state that it is quite possible that there are still more bugs there, but I still claim the results to be represntative.

Python module usage statistics

Posted on by 5 Comments

IMPORTANT UPDATE: The code I used to create these statistics had some bugs. The fixed statistics are available here.

After reading Doug Hellman’s post about python stdlib modules he needs the documentation to, I commented there that I need the documentation for logging because I don’t use it too frequently. Later, I thought a little bit more about it, and I wanted to check which modules are actually used often, and which rarely.

My first step was to “grep import *.py” on my python scripts directory. Later I progressed to writing a simple script which basically does the same thing, but also rates each imported module according to frequency (it looks for “import <modules>” and “from <module> import bla”). The results weren’t interesting enough. Before going any further, I sat down, and wrote my expected list of frequent stdlib modules:

  • sys
  • os
  • re
  • math
  • thread
  • urllib2
  • time
  • random
  • struct
  • socket
  • itertools
  • operator

My next step was to try Google Codesearch. Some time ago I wrote a little script to harvest results from Google Codesearch, but enough time has passed, and the script doesn’t work anymore. Parsing their output is messy enough, and their results don’t seem that helpful anyway. (They don’t return all the results in a project). So I thought about using koders.com, when I had a flash. I can use Cheeseshop! easy_install must have a way to just download the source of a project! (For those who don’t know: you can do easy_install SQLObject, and easy_install will look for the module in the PyPI, download it and then install it for you). Indeed, easy_install had this option (-eb), and I was in luck. Quickly, I got a list of all the modules in PyPI, shuffled it, and picked 300 randomly. (Actually, downloading them wasn’t that quick, and about 20 didn’t download right :).

I ran my statistics script, but non-stdlib modules were also listed. I fixed some bugs, improved the script a little, and (partially) disallowed inner project imports using a little heuristic. These are the improved results:

  • sys,113
  • os,107
  • setuptools,57
  • common,50
  • unittest,50
  • __future__,25
  • distutils.core,24
  • zope.interface,22
  • re,19
  • pygame,18
  • time,13
  • datetime,11
  • string,10
  • zope,10
  • pyglet.gl,9
  • types,8
  • random,7
  • pkg_resources,7
  • gtk,6
  • struct,6
  • ez_setup,6
  • zope.component,5
  • math,5
  • logging,5
  • sqlalchemy,5

(Please note that I do not argue that these results are completely true, as my code might still contain problems, and it probably filters too much. However, I do claim that these results are representative).

Well, some surprises… I didn’t expect setuptools and types to be there. Especially types… who needs it anymore (and anyway)? Seeing unittest so high in the list gives me a pretty good feeling about opensource. __future__ being there is amusing… Seems that we just can’t wait for the next version. And indeed, logging is pretty low, but not as low as I thought it will be.

Another curiosity: the module frequencies seem (unsurprisingly) to obey Zipf’s law. Here is a small graph showing the results: module imports on a log-log graph

(plotted with Google Docs). The R squared value is 0.92. For the original results (where inner dependencies were allowed), it seemed that the law was even stronger.

For those interested in the full results, here is a csv file containing them.

One final note: while these results are true for projects in Python, they don’t represent usage in interactive mode. Since I use interactive mode a lot, and during interactive work I use os and re quite a lot, I bet the results would be changed if we somehow could count interactive mode usage.


Pagerank recommends movies

Posted on by 1 Comment

So how was your Yom Kippur? After Yom Kippur ended, I sat down to write something that was nagging me for quite some time. I wanted to see what pagerank had to say about movies. I’ve always liked to picture things as graphs, and looking at movies and actors as nodes in a graph is very intuitive. I also think that the idea has some merit. I like most of the movies that Edward Norton appears in, so it makes sense to give him a high score, and see which other movies Edward Norton “recommends”.

So I fired up my IDE, and started to work.

Well, at first I was a little bit disappointed, because IMDbPY doesn’t have code to read the IMDB charts. Well, I was going to extend it a little bit, but then I remembered that I already downloaded all the data I needed some time ago (more than half a year I believe), and that was good enough. So I opened up some old pickle file I had lying around, and Voila! after some playing around with the format, I had a weighted graph at my disposal. I weighted the graph according to the IMDB listing – if an actor appeared higher on the cast, his (or hers) weight to and from that movie will be higher.

I implemented a quick and dirty pagerank, and then added some code to score my recommendations higher. I did that by adding each iteration my recommendation’s weight to the movie’s score.

Without recommendations, the top spots were:

  • 24
  • The Sopranos
  • The Wizard of Oz
  • JFK
  • Cidade de Deus
  • Citizen Kane
  • King Kong
  • Voyna i mir
  • Dolce vita, La
  • The Patriot

After adding “Edward Norton” and “Tom Hanks” as recommendations, I got the following results:

  • American History X
  • Fight Club
  • 24
  • Red Dragon
  • The Sopranos
  • Catch Me If You Can
  • The Wizard of Oz
  • JFK
  • Forrest Gump
  • Cidade de Deus

Well, I was surprised a little to see 24 and the Sopranos rank so high, and also a little bit disappointed. The recommendations didn’t really add enough information, and more work will have to be done for this to work properly. It has the right ‘smell’ however. It looks like it has potential to be useful – I can just see the program, having a small table with a ‘seen/unseen’ checkbox, and a place to write your favorites, and the program will tell which movies you should watch next. This reminds me a little bit of the recommendations on soulseek and Amarok. It would be really fun to write a total recommendation engine, with the best of all worlds. Someone probably already implemented it.

In any case, I’ll work on it some more. This is interesting.

If anyone wants the code, tell me, and I’ll make it ‘online worthy’…

A new fold on the origami page: Qua Qua Egg

Posted on by Leave a Comment

In Israel, we have a fold we call ‘qua-qua-de-la-oma” or for short, “qua-qua”. I believe in english it is called a ‘color changer’, but I’m not sure of that. In any case, this fold yields a 3D shape similar to two connected tetrahedrons. (It is only similar because the shape is not completely symmetric). I couldn’t find a name for this shape, and since the same fold yields a qua-qua, I called it a qua-qua-egg. This fold is very useful as a makeshift ball to throw around the office.

Origami Qua-Qua-Egg

Understanding the Nyquist Limit

Posted on by 2 Comments

The Nyquist Limit of half the sampling rate, often seemed to me as a ‘magical limit’. After I read the proof I understood it better, and the graphical representation of the overlapping spectra (also seen on the wikipedia page of the proof) is quite convincing. Yet now I found a very interesting and ‘easy’ explanation, which I’m going to relate. The explanation may not be complete, but it certainly gave me a good feel of the limit.

It all started when some late night I was driving with my girlfriend from Haifa to Tel-Aviv. (In Israel, many mathematical school problems start that way :). As it sometimes happens, part of the road was under repairs, and there were ‘danger lights’ alongside the road. These danger lights were many beacons, spaced a few meters apart along the road, each blinking in its own time. Driving alongside them, you could see it as if the ‘light’ was going through them ( see illustration, 0 is unlighted, x is lighted):

… x 0 0 x 0 …

… 0 x 0 0 x …

… 0 0 x 0 0 …

and so on. To me it seemed as if the light was advancing. Let us assume now a light pattern of one lighted, two off, and so on, going forward, one beacon at a time. This is the pattern in the illustration. This pattern however, is identical to the pattern of one lighted, two off, going backwards, two beacons at a time. Just like watching a rotor – you could parse it either way, and both ways are correct. This could be interpreted as aliasing.

At this point – this will already start to ‘smell’ like something is going on with the sampling theorem… Come to think of it, watching a rotor is exactly an example of the application of the sampling theorem! Cool!

(Warning: this is very reminiscent of group theory. Every now and then I will use it a little. To the non-mathematically-inclined: bear with me :)

Let’s move on. Imagine now a finite number of beacons. Let’s say five. Now consider a single light traveling through the beacons. What is the maximum speed in which it can travel, where we still can discern the direction? (this is equivalent to an infinite number of beacons, with a pattern of 1 on, 4 off. Just like Z5 ~= Z/[5])

Let’s check it out. (note: since 5 is prime each speed will generate one pattern only and not more.) We will obviously ignore the zero speed, and move on to the speed of 1. At this speed, the pattern generated is:

x 0 0 0 0
0 x 0 0 0
0 0 x 0 0
0 0 0 x 0
0 0 0 0 x

or in reverse.

At the speed of 2, the pattern generated is:

x 0 0 0 0
0 0 x 0 0
0 0 0 0 x
0 x 0 0 0
0 0 0 x 0

or in reverse.

Notice that we can discern both cases which direction the pattern is going. But a speed of 3, will be exactly the reverse of the speed of two! Equivalently, a speed of four will not be discernible from a reverse of speed 1! If you don’t believe me, try to see it in the patterns above. This is equivalent of saying that over Z5, -1 == 4, and -2 == 3. We can carry on with this, and show it is the same for any n. (In the additive group Zn, it is easy to show that the inverses lie in the different halves of the group). To summarize – given a finite number N of beacons, we will be able to discern the directions of speeds up to floor((N-1)/2). Cool.

I know that for the mathematically inclined this might seem obvious, and for the rest of the world, this could seem a little cryptic. But I’m really happy with finding this connection, and such a simple example of the Nyquist Limit.

Mathematical note: I bet I got a few mathematical things wrong here. I didn’t write a rigorous proof (but I did proof read it). Tell me your corrections, and what you think. I’ll be glad to hear.

Origami crease-patterns are NP-complete

Posted on by 1 Comment

After talking about it with someone I met, I looked it up, and it does seem to be so. There is of-course a lot of research on the subject that I didn’t yet get the chance to read, but I think this is the general picture.

In any case, this shows that my suggestion of using origami as a hash function is valid! Really cool stuff… I wonder if it could actually be useful somewhere.

Computation over Zp in Python

Posted on by 2 Comments

Lately I’ve been working a lot on my Algebric Structures homework. One of the reasons I don’t blog as much as I should. While working on my homework, I had to factor some polynomials over Z5 – the field containing the numbers {0,1,2,3,4}. The trick was that given a 4th degree polynomial without any roots, you could only factor it to two 2nd degree polynomials.

That turns out not to be that hard, but a bit of work – especially if you dislike just computing values. So I wrote a simple script to solve this little problem for me. While writing that script, I saw that it could be really fun to write a general Zp class. So naturally, I fired Google up, and tried looking for existing implementations. I didn’t find any (except some python implementations of many number thoery algorithms). Being quite a small amount of work, I hacked up something. After finishing with it, I wrote a polynomial class that can accept any given “number class”. There is actually a public implementation of real valued polynomials in SciPy, but not for integer valued or Zp-valued polynomials. After some more work my polynomial class was also finished. I actually used that class for some computations later on.

I added polynomial division, just for the heck of it, even though I didn’t use it. I hope this little script will be useful to someone.

Riddle of the Week – LStrings

Posted on by Leave a Comment

So I wrote about lstrings. I intend to write about them again in a short while – I already finished the basic script a few days ago, but I’m waiting, until I will be satisfied with it.

In the meantime, here is a curious riddle (That I came up with):

Version 1 (not final)

Assume you are given some iteration of an lstring – S. How do you discover the original lstring used to create S?

Discussion

Well, those who are quick to answer will first try to discover the length of the original string – it could be computed directly from S. At this point you should come to the conclusion that without being given the exact number of iterations – the riddle is not that interesting… Iterating on an lstring is an associative operation. If we denote the substitution of a string A into a string B as multiplication, when both are iterations on an original string T , we will get:

A*B =T^n * T^k = T^(n+k) = T^k* T^n

This will result in the obvious answer to the original riddle – “Why, the lstring S is the original string, with zero iterations!”. This answer is correct – and obviously useless :)

So, let us rephrase the riddle:

Version 2 (final)

Assume you are given some iteration of an lstring – S. How do you discover the minimum possible lstring that could be used to create S?

Discussion – But not a solution

This is the proper riddle – have at it! I will be glad to read your thoughts…

note: The acceptable solution should be either an algorithm, or a constructive mathematical proof. I don’t like reading only statements of existence :)