Javascript and the Ising Model

I was amazed by Mike Bostock’s algorithms visualization. If you haven’t had a chance to see it, stop reading and go there. The fact is that I felt inspired to learn some Javascript and d3.js. Here are my favorite references: Eloquent Javascript (which is a great book that is open for download), a fun, difficult introduction to d3.js (that’s the name of the tutorial) and mbostock’s blocks (which is a series of examples by Mike Bostock). I think I learned mostly by taking some of Bostock’s examples and trying to modify them. I also learned a lot by just googling and searching StackOverflow.

When someone mentioned Javascript to me, my first thought was: oh, this is going to look like Java ! Actually no ! To my surprise it has nothing (or very little) to do with Java. It is a functional language like Haskell (btw, learning Haskell is a lot of fun. Here is a tutorial that is so much fun you can do just to entertain yourself — and it is free online).

My first impulse was to create an algorithms visualization. When I was trying to find a good one to implement, I saw in my bookshelf the Markov Chains and Mixing Times book that was sitting there for a while and I recalled it had all those nice models from statistical physics that are analyzed algorithmically. The one that is at the same time very simple and interesting is the Ising Model. In the Ising Model, there is a grid (or more generally a graph) in which each cell (or node in the graph) has a spin s \in \{-1,1\}. So a state is an association of a spin for each cell. A transition in the Ising model is as follows: we choose a cell at random, look at its neighbors (the cells adjacent to it) and compute s(N_i) = \sum_{j \in N_i} s_j where N_i is the set of neighbors of i. Then we update the spin of i to \sigma with probability e^{ \sigma s(N_i) u} / (e^{\sigma  s(N_i) u} + e^{- \sigma s(N_i) u}), where u is a parameter known as the inverse temperature parameter. If u = 0 the selected node transition to each spin at random. If u \rightarrow \infty, it transitions to the spin that is more frequent among its neighbors (breaking ties at random).

Here is the visualization: you can choose the inverse temperature (parameter u) as well as the speed. Then in each timestep, one cell is chosen at random, colored in red, and then updated (it is colored in black if the spin is 1 and in white if the spin is -1). We start the simulation from a random configuration (observe that setting the inverse temperature to 0, the spins are scattered, but increasing the inverse temperature, spins start concentrating on certain regions) :

Here is how to do it in Javascript ? First, let me show you a standalone example (i.e. not an example embedded in a blogpost). Follow this link. There you can just look at the source code (in Chrome, View->Developer->View Source) and you get the entire code, which we will discuss in a bit. One other cool thing you can do there is to open the Javascript Console and interact with visualization. In Chrome, go to View->Developer->JavascriptConsole. Then you will open an interface where you can play with the visualization. Try typing something like interv1 and it will return the value of this variable, which is the interval between two updates in the first graph. Now try changing the value by writing interv1=2 and making it a lot faster (you might also want to set the temperature to zero to make updates more evident, either by using the slider or by doing itemp1=0 in the console). Now, to see you can iteract even more, try changing the color of the first cell to green, by doing:


Ok, now to the code. If you look at the html file, the first thing we do is to import the d3.js library, which is the library we use for visualization (in Javascript, like in any programming language we use libraries to get all sorts of cool functionalities). After that we define a function call generateGrid which takes number of horizontal tiles, number of vertical tiles and size of each tile size. The function first creates an SVG canvas, which is a container for graphics. Then adds to it Nx times Ny squares. Notice the for-loop in Javascript is a bit different: we create a range and then call a function for each element of the range, by doing d3.range(n).forEach(function(x) { ... }), which would be equivalent to for(int x = 0; x < n ; ++x) { ... } . The function returns a map from coordinates of type [i,j] to the corresponding square. Now the subsequent code just takes some cell and modifies its color. Here is the part discussed so far:

<script src=""></script>

// Generates and Nx by Ny grid with pointers to each tile in the grid
// return a map from [i,j] to the corresponding tile
function generateGrid(Nx, Ny, side_length) {
  var svg ="body").append("svg")
    .attr("width", side_length*(Nx+2))
    .attr("height", side_length*(Ny+2));
  var gTile = svg.append("g")
    .attr("stroke", "#000000")
  var gPointmap = {};

  d3.range(Nx).forEach(function(x) {
    d3.range(Ny).forEach(function(y) {
      gPointmap[[x,y]] = gTile.append("rect")
        .attr("x", side_length*(x+1))
        .attr("y", side_length*(y+1))
        .attr("width", side_length)
        .attr("height", side_length)
        .attr("fill", "#FFFFFF");
  return gPointmap;

var spincolor = {"1":"#000000", "-1":"#FFFFFF"};

So far we have only a helper function but we don't have yet anything. So if you just add the line pm = generateGrid(10,10,20); you create a chessboard. Actually not quite because it is all white. In fact, if you go to this link, open the Javascript console and type:

pm = generateGrid(10,10,20);

you will see your 10 by 10 white grid appearing in the end. If you want to color some of its nodes like a chessboard, you can do something like:

d3.range(10).forEach(function(x) {
  d3.range(10).forEach(function(y) {
    if ((x+y)%2 ==0) pm[[x,y]].attr("fill","red");

then you will get a chessboard. You can create others and keep changing their colors. But now, let's go for the Ising Model simulation: first we create a grid, then we pick a spin at random for each of the nodes. Finally, we define a main loop iteration. We do it as a callback: i.e. a function that returns another function, which is the mainloop execution. We do that since d3.timer(...) takes a function and executes it after a certain interval, so the main loop callback returns a function

  // create a white grid 20 by 20 grid
  var n2 = 20;
  var i2_pointmap = generateGrid(n2,n2,20);
  itemp2 = 10;
  interv2 = 200;
  directions = [[-1,0],[1,0],[0,-1],[0,1]];
  spin2 = {};

  // random spin for each cell and color it accordingly
  d3.range(n2).forEach(function(x) {
    d3.range(n2).forEach(function(y) {                                                                
      spin2[[x,y]] = 1-2*(~~(2*Math.random()));                                                       
      i2_pointmap[[x,y]].attr("fill", spincolor[spin2[[x,y]]]); 

  // main loop (in fact a callback, returns a function that executes
  // one iteration of the main loop)
  function square_ising_callback() {
    return function() {
      // samples a random cell
      var x = ~~(n2*Math.random());                                                                   
      var y = ~~(n2*Math.random());
      // colors this red to display with cell is being updated
      i2_pointmap[[x,y]].transition().attr("fill", "#FF0000");
      // calculate the sum of neighbors spins
      var neighbors_spin = 0;
      d3.range(4).forEach(function(i) {                                                               
        var x1 = x+directions[i][0], y1 = y+directions[i][1];
        if (x1 >= 0 & x1 < n2 & y >= 0 & y < n2) {
          neighbors_spin += spin2[[x1,y1]];
      // compute the probability of changing to spin 1 and samples
      // the new spin according to this probability
      var prob_pos_spin = Math.exp(itemp2 * neighbors_spin) / ( Math.exp(itemp2 * neighbors_spin)+ Math.exp(-itemp2 * neighbors_spin) );
      if (Math.random() < prob_pos_spin) {
        spin2[[x,y]] = 1;
      } else {
        spin2[[x,y]] = -1;
      // colors the cell with the color of its new spin. We add
      // a small transition with delay interv2/2 such that the transition
      // is visible
      i2_pointmap[[x,y]].transition().delay(interv2/2).attr("fill", spincolor[spin2[[x,y]]]);
      // after interv2 time, calls back the main loop
      return true;
  // call the main loop for the first time
  d3.timer(square_ising_callback(), interv2);

And that's it ! I was surprised how little effort it took. And I guess it would have taken a lot less if I took time to write the appropriate functions to simulate the Ising model instead of just writing everything in the main loop, which admittedly is not a good practice.


PS: You might wonder how to include Javascript in WordPress. The solution I got was to create an html file and add to WordPress as an iframe. To do that, simply add something as follows in the textual editor (and of course, adding this to any webpage will display the diagram in this blog post):

<iframe src="" style="border:0; width: 700px; height:550px"></iframe>

Bakeries in West Village / Soho

Say you are in the West Village or Soho (around the NYU area) and feel the need for French pastries. This is what this blog post is about. I tried a few of them in my last year and this is what I got:

  • Francois Payard Bakery (FPB): this is my favorite one. I judge French bakeries by the quality of the chocolate eclair. The one in FPB is hard to beat. I usually go there and get the same things every time: 1 chocolate eclair, 1 gateau basque and 1 lemon tart (because my wife loves lemon tarts). Sometimes I run very early, then I stop there and get fresh bread — it never disappoints us.
  • Lafayette Bakery: This is a small bakery inside a very chic French restaurant. I’ve never been to the restaurant (although I am constantly thinking about going there in a special occasion), but I stop by bakery almost every weekend. I get the chocolate eclair (as usual), get caneles for my wife (as she is crazy about caneles and this is the best I’ve found in NYC) and sometimes the praline. The praline is really good. I usually get the sweets and bring them home, but they also have 2 or 3 tables for the bakery and this is a great place to stop for a coffee with some sweets. Their bread is also very good and sometimes I buy their baguette. If you go there, you might also want to spend some time in the neighborhood. Across the street, there is Astor Wines — the best wine store in town. The staff is really knowledgeable and they have wines in all price ranges — and this is usually safe to buy cheaper wines there, as their selection is amazing. So probably you won’t be disappointed. Also around is La Colombe, which is a very nice coffee shop. Around is also a very high end butcher, where I haven’t bought anything — since this is a bit outside my possibilities, but sometimes I like to go there to admire the meat on display. This is also a place surrounded in mistery.
  • Claude Patisserie: Claude is a hidden treasure — it is located in an ugly street and the store front is not welcoming. But once you are there, the sweets are just incredible. Their eclaires are amazing, so is their mille-feuille and their cakes and their petit feurs. We get their cakes for our birthdays.
  • Macarons & Cookies is a small stand in Prince Street — when we are out of the subway we like to stop there for caneles (did I mention we love caneles?) and their alfajores.
  • La Maison du Macaron: a bit further on the 23rd st, there is La Maison du Macaron. First, the chocolate eclaire is great and so are the macarons — they are very fresh and one of the few places in NYC I found the macarons worthy (even Laduree in NYC is not great, as they fly their macarons from France and they don’t get to NYC very fresh). La Maison is the opposite: you can see the kitchen where they are being made. They are so fresh, delicious and with creative flavors. Everything they do tastes good: their coffee, their croissant, their bread, …

You might want to know my opinion about coffee shops since I am a person that freely gives opinions about bakeries. I might tell you that coffee shops deserve their own post, but nevertheless I’d mention some that I really like: Joe Coffee (they are all over the place — the one near Columbia is a great place to work), Stumptown (where I usually get beans from), Gimme Coffeee (the most honorable outpost of rural New York in this big metropole), Ninth Street Espresso (I usually go to the one in Chelsea Market), Blue Bottle Coffee (the SF chain has now several stores in Manhattan. The one in Chelsea has a mezzanine with a Japanese coffeee bar, which is entertaining to go if you haven’t been to one of those) and Grumpy (which has great coffee, but they ask you not to work there).

List of Japanese Restaurants I like in NYC

Here are some of the Japanese places I like in Manhattan (this is more or less equivalent to the list of restaurants I like in Manhattan, since 90% of the time I go out to restaurants, I end up going to Japanese places). Also, I must warn you that I already did post list of restaurants I like in the past and upon looking at those, I regret many choices. This happens for various reasons, a once great restaurant can become mediocre, my taste can change or I can be introduced to much better versions of what I liked in the past and cease to like previous ones. That being said, the list below is quite accurate for 09/27/2014. I might update or post a new list when this substantially changes.

The first thing to do when you go to NYC with the goal of getting Japanese food is that you should visit East 9-th St, between 2nd and 3rd avenue. There is an enormous concentration of Japanese places: there Sunrise Market, a Japanese supermarket,  there is are two stores offering Japanese street food:

  • Otafuku: serves mainly takoyaki, which are fried balls filled with octopus. Among all fried things ever done in this world, takoyaki is my favorite thing. And this is coming from someone that avoids eating fried things, but takoyaki is not one you can (or should) resist. They have medetai, a pastry in the form of a fish filled with either beans or chocolate and banana. Both are delicious and our first choice when we go for a simples and cheap snack.
  • Yonekichi: it is a beautiful rice burger, i.e., instead of a bun, there is rice, something more or less similar to an onigiri, but really in the form of a burger. There is a variety of choices for what comes inside: my favorite is the eel and avocado (unagi burger). The one with Salmon is also delicious. It is also quite light and healthy, much unlike your usual burger — in fact, with a burger it shares mostly the name and the form, but other then that, it is much more similar to onigiri.

Besides street food, there is also desert:

  • Cha-An: a Japanese tea-house with a variety of deserts — I am a big fan of their mochi and their black-sesame brulee. I also recommend the desert sampler in the first time you visit the place. While their specialty are deserts, toasts and tea, they also have one dinner set in the menu. I had that once and that was quite satisfying. But mostly we go somewhere else for dinner and come to Cha-An for deserts. Just beware that they are cash-only, so come prepared.

There is a great sake store — Sakaya — also nearby. This is the best sake store I found in Manhattan so far. There is a great yakitori place — Yakitori Taisho — that is super cheap one block away. Ok, now you might be thinking that I pointed you to street food, yakitori, sake, tea house, supermarkets, … But if you came to this post for sushi, sashimi and izakayas, I should also say something. Let me start (all the options in the following block are more of less equivalent in price, being $40 per person plus $10 with sake):

  • Robataya: also in East 9th there is Robataya (just across the street from Cha-An) that has a great menu — both izakaya food (Japanese tapas) and some sushi and sashimi. The kamameshi rice they have is particularly delicious — I got the one with eel (when I am faced with the possibility of eating eel, I don’t think twice). We also got the mille-feuille and the saba sashimi and everything was delicious.
  • Izakaya Ten (23rd st and 10th ave): probably my favorite Izakaya restaurant in Manhattan. I love the octopus they have, the eel is fantastic, good yakitori, grilled fish, … Sometimes they have sashimi and even uni once in a while. Whenever this happens, it always makes me happy.
  • Shuka Dining Bar: this used to be the home of our absolute favorite Japanese place in Manhattan — Jinya. Jinya closed and it became Shuka, with a similar crew, but a more upscale menu. While I miss the previous restaurant, I must say I love the new menu, specially the Uni toast. The yakitori is quite good and we got various small dishes and all of them were great. I recently noticed that they are serving brunch on Sunday (for $22) which is really good.
  • Sakagura: Sakagura is hidden in the basement of an office building and it is completely impossible to find if you don’t know the exact address. It is a hidden gem. They have an izakaya menu (and some sashimi, but no sushi). They are also known to have the largest collection of sake in NYC. During the week, they also have lunch menus which are really great and considerably cheaper than dinner.

Now sushi and sashimi, but beware I didn’t explore sushi and sashimi as much as I did Izakaya type restaurants. I am writing about some of the ones I’ve been. I like all of them and came back many times to each:

  • Neta: Neta is the best sushi I had so far in NYC. The attention to details is incredible. If you sit in the sushi counter you can observe the amount of work that goes into every piece of sushi. They are also open for lunch everyday. A sushi set during lunch is around $33. During dinner, Omakase menus are around $100. It is worth every cent of it. For desert, they manage to have the best ice-cream I had in the city.
  • Ootoya: Ootoya is a Japanese chain and the main thing in the menu are rice bowls with sashimi on top, in other words, Chirashi ! I get usually the Kaisen Don, Hanabi Don, Maguro Don or the Chirashi. They have other things, such as soba, tonkatsu, … but I haven’t had a chance to try them. I constantly come back there for the rice bowls. Prices are around $20 per person during lunch and $30 per person for dinner. Add a bit more for sake.
  • Marumi: Marumi is our neighborhood sushi spot. I come there often and I love it. Every time I like it more then the previous time. I always order the soft-shell crab as a starter (actually, sometimes I order two of them!). It is just too delicious. They have also some potatoes with uni for starter that are great. Sushi is of great quality and I also really like the unagi don and some of the other dishes.
  • Lobster place: inside the fish market in the Chelsea Market building there is a sushi counter. Their sushi is quite expensive, but their chirashi, chu-toro don and unagi don are reasonable and delicious. There is also charm in eating in the middle of a fish market. At least all the raw fish around makes me super hungry. I come there often, mostly for the chirashi. They also have sushi to go — even though I usually don’t like taking food to go, I’ve once or twice taken food from this place to eat in the Highline Park nearby. The sushi to go is also freshly made and a bit cheaper then in the counter. It comes from the same place and the quality is the same — except you get in a plastic box instead of the pretty wooden box.

Oh, at this point you are thinking of asking me about ramen and soba ? Here we go:

  • Cocoron: I learned about Cocoron next week and it is delicious.I got the Mera Mera dip soba, which seems to be their speciality. Amazing. all appetizers were also great.
  • Ippudo: I came here a couple of times and it is always great — it is also a Japanese chain, but don’t go there on Fridays or in the weekends. It is such a huge line and impossible to get a table. Once they told me it was a 3 hour wait. And also they don’t take reservations. But if you go there in the middle of the week, or very early or late, it is likely you will leave happy.
  • Momofuku Noodle bar: also great noodles and prices. Tasty appetizers. The place is small and they ask you to share a table with other people. Most likely you get the counter. but again, hard to regret — noodles are great.

I might end this post by taking about Bohemain, but I can’t, since Bohemian is a secret. Instead I leave you with some pictures (not of Bohemian, of course) in the form of a slideshow:


This slideshow requires JavaScript.

Multiplying polynomials in nearly-linear time without leaving the ring of integers

The problem is really simple: given two polynomials: p(x) = p_0 + p_1 x + p_2 x^2 + \hdots + p_{a} x^{a} and q(x) = q_0 + q_1 x + q_2 x^2 + \hdots + q_{b} x^{b} with integer coefficients, how to compute the product r(x) := p(x) \cdot q(x). The product is a polynomial of degree a+b where the coefficient r_k = \sum_{i,j \vert i+j=k} p_i \cdot q_j. The definition suggests an algorithm of complexity O(a\cdot b) = O(n^2) where n = a+b+1.

Surprisingly, one can do it a lot faster: in fact, one can calculate the coefficients of p(x) \cdot q(x) in time O(n \log n) if we are allowed to use complex numbers. This is one of those “a-ha” algorithms: the idea is very simple but incredibly ingenious — and a proof that surprising algorithm do exist. The algorithm is based on the following idea: if we know the value of a polynomial of degree d on d+1 points, then we can reconstruct it exactly. So, if one can evaluate polynomials p and q on ab+1=n points, say \{ x_0 x_1, \hdots, x_{n-1}\} then we can calculate r(x_i) = p(x_i) \cdot q(x_i) and reconstruct r. Notice that once we calculate p(x_i), q(x_i), obtaining r(x_i) is takes only O(n) operations.

So far, we only shifter the complexity around: while the actual multiplication takes O(n) we still need to evaluate the polynomials n times and then reconstruct r from the values r(x_0), r(x_1), \hdots, r(x_{n-1}). The idea is to choose values x_0, \hdots, x_{n-1} carefully in such a way that we can evaluate those polynomials really fast.

Enter complex numbers ! A choice that works is to take x_i = \omega_n^i where \omega_n = e^{2 \pi i / n} \in \mathbb{C}, i.e., n-th root of the unit. And here is the reason why we can multiply polynomials so fast:

\begin{bmatrix} p(\omega_n^0) \\ p(\omega_n^1) \\ \vdots \\ p(\omega_n^{n-1})    \end{bmatrix} =   \begin{bmatrix} 1 & 1 & \hdots & 1 \\ 1 & \omega_n^1 & \hdots & \omega_n^{n-1} \\  \vdots & \vdots & & \vdots \\ 1 & \omega_n^{n-1} & \hdots & \omega_n^{(n-1)^2}   \end{bmatrix} \cdot  \begin{bmatrix} p_0 \\ p_1 \\ \vdots \\ p_n    \end{bmatrix}

In other words, if we see p and p(\omega^i_n) = [p(\omega_n^i)]_i as vectors of dimension n, then we can write them as p(\omega^i_n) = A_n \cdot p where A is a matrix n \times n with entries \omega_n^{ij} for i,j=0,\hdots,n-1. Now, we use the power of recursion. Assuming n is a power of 2 to simplify things, we first note that:

\begin{aligned} p(x_i) & = \sum_{j=0}^{n-1} p_j \omega_n^{ji} = \sum_{j=0}^{n/2 - 1} p_{2j} \omega_n^{2ji} +  \sum_{j=0}^{n/2 - 1} p_{2j+1} \omega_n^{(2j+1)i}  = \\ & = \sum_{j=0}^{n/2 - 1} p_{2j} \omega_{n/2}^{ji} +  \omega_n^i \cdot \sum_{j=0}^{n/2 - 1} p_{2j+1} \omega_{n/2}^{ji}  \end{aligned}

since \omega_n^2 = (e^{2 \pi i / n})^2 = e^{2 \pi i / (n/2)} = \omega_{n/2}. So, if we define p_{\text{even}} = [p_{2j}]_{j=0..n/2 - 1} and p_{\text{odd}} = [p_{2j+1}]_{j=0..n/2 - 1}, then we see that we can compute in linear time the coefficients of p(\omega^i_n) from the coefficients of A_{n/2} p_{\text{even}} and A_{n/2} p_{\text{odd}}.

So, the problem of computing the coefficients of p(\omega^i_n) corresponds to computing two problems of half the size plus some O(n) computation. This gives a total time of O(n \log n). The algorithm works as follows:

% assumes length of p is a power of 2
function w = fft(p)
  n = length(p);  omega = exp(2 * pi * i * sgn / n);
  if (n == 1) 
    w = p;
  p_even = p(1:2:n);  w_even = fft(p_even); 
  p_odd  = p(2:2:n);  w_odd  = fft(p_odd);
  for j = 0 : n-1
    w(j+1) = w_even(mod(j,n/2)+1) + omega^j * w_odd(mod(j,n/2)+1); 

We call the method fft since this is also known as the Fast Fourier Transform. We are almost there: we know how to compute p(\omega_n^i), q(\omega_n^i) and therefore r(\omega_n^i). We only need to get the coefficients of r. Notice that r(\omega_n^i) = A_n r, so we need only to solve this linear system.

How to do it ? Another “a-ha”: the matrix A_n^{-1} = \frac{1}{n}[\omega_{n}^{-ij}] for i,j = 0, \hdots, n-1. The reason for that is that: \sum_{k=0}^{n-1} \omega_n^{ik} \omega_n^{-kj} = \sum_{k=0}^{n-1} \omega_n^{k(i-j)}. If i=j, then this is clearly n. If i \neq j, this zero, since this is the sum of all roots of the unit. Therefore, we can calculate A_n and its inverse in a very easy way:

function w = fft_aux(p, sgn)
  n = length(p);  omega = exp(2 * pi * i * sgn / n);
  if (n == 1) 
    w = p;
  p_even = p(1:2:n);  w_even = fft_aux(p_even, sgn); 
  p_odd  = p(2:2:n);  w_odd  = fft_aux(p_odd, sgn);
  for j = 0 : n-1
    w(j+1) = w_even(mod(j,n/2)+1) + omega^j * w_odd(mod(j,n/2)+1); 

function w = fft(p)
  w = fft_aux(p, +1);

function w = fft_inv(p)
  w = fft_aux(p,-1) ./ length(p);

Now, we have all the pieces to multiply to polynomials:

function r = multiply_poly(p, q)
  n = length(p) + length(q) - 1;
  p = expand_to_power_of_2(p, n);
  q = expand_to_power_of_2(q, n);
  r = fft_inv(fft(p) .* fft(q));
  r = resize(r, 1, n);

% pads p with zeros until it has length equal to 
% the first power of 2 greater of equal n
function p = expand_to_power_of_2(p, n)
  plen = 1;
  while (plen < n) plen = plen * 2; endwhile
  p = resize(p, 1, plen);

And now we are done ! We multiplied two polynomials in time O(n \log n). Even though our initial objective was to multiply two polynomials with integer coefficients, we had to leave the ring of integers, go to the field of complex numbers, multiply things over there and then come back to the ring of integers. Those things happen all the time in algebra, that the answer for a question involves leaving the space in which a question is asked, solving the question in a higher space and coming back in the end. But here we ask: can we multiply two polynomials doing only operations over the ring of integers ?

Schonhage and Strassen give an algorithm that runs in time O(n \log n \log \log n) that does exactly that (this later improved here and here). The idea is to try to mimic the FFT-algorithm without using complex numbers. What we needed was an element \omega with three properties: (i) \omega^n = 1, (ii) \sum_{i=0}^{n-1} \omega^i = 0 and (iii) \{\omega^0, \omega^1, \hdots, \omega^{n-1}\} are n distinct elements.

The idea of the algorithm is to work with some sort of symbolic \omega that has the properties listed above. Let's begin by fixing some notation, let \mathbb{Z}[x] denote the ring of polynomials over integers, which is just another name for the set of polynomials with integer coefficients. Now, fixed a certain polynomial q(x) \in \mathbb{Z}[x], we define the class of residues module \mathbb{Z}[x] / (q(x)) as the residual we get from dividing a polynomial in \mathbb{Z}[x] by q(x). More formally, associate with each p \in \mathbb{Z}[x] the set \bar{p} := \{p(x) + q(x) \cdot a(x); a(x) \in \mathbb{Z}[x]\}, then: \mathbb{Z}[x] / (q(x)) = \{ \bar{p}; p \in \mathbb{Z}[x]\}.

This is just some formalism. here we are only concerned with two such rings: (i) \mathbb{Z}[x] / (x^n - 1). This is the ring in which x^n = 1. This means that for every p(x) = \sum_{i=0}^K p_i x^i \in \mathbb{Z}[x] we associate the polynomial \bar{p} =   \sum_{i=0}^K p_i x^{i \text{ mod } n} \in \mathbb{Z}[x] / (x^n - 1). (ii) \mathbb{Z}[x] / (x^n + 1). This is the ring in which x^n = -1. This means that for every p(x) = \sum_{i=0}^K p_i x^i \in \mathbb{Z}[x] we associate the polynomial \bar{p} =   \sum_{i=0}^K p_i (-1)^{\lfloor i / n \rfloor} x^{i \text{ mod } n} \in \mathbb{Z}[x] / (x^n - 1).

The great thing is that if instead of \mathbb{C} we have \mathbb{Z}[x] / (x^n + 1), then we have an element like \omega, which is x^2. Notice that (x^2)^n = x^{2n} = (-1)^2 = 1, and \sum_{i=0}^n x^{2i} = 0 since x^{2(i + n/2)} = x^{2i} \cdot x^n = -x^{2i}. This is essentially the main trick: we will do arithmetic on \mathbb{Z}[x] / (x^n + 1), doing FFT and all... Now, comes the last piece of the puzzle: notice that if \text{deg}(p) \cdot \text{deg}(q) < n[/latex], then multiplying [latex]p \cdot q \in \mathbb{Z}[x][/latex] is equivalent to [latex]\bar{p} \cdot \bar{q} \in \mathbb{Z}[x] / (x^n + 1)[/latex], so we only need to multiply in the ring [latex]\mathbb{Z}[x] / (x^n + 1)[/latex].  The idea to multiply in this ring is as follows: write [latex]\bar{p} = \sum_{i=0}^n p_i x^i = \sum_{i=0}^{\sqrt{n}} (\sum_{j=0}^{\sqrt{n}} p_{j + \sqrt{n}i} x^j) x^{i \sqrt{n}}[/latex]. Notice that this maps to a polynomial with degree [latex]\sqrt{n}[/latex] taking values in [latex]\mathbb{Z}[x] / (x^{\sqrt{n}} + 1)[/latex]. Now, the coefficients (which are in [latex]\mathbb{Z}[x] / (x^{\sqrt{n}} + 1)[/latex]) have a suitable [latex]\omega[/latex] we can use and multiplying them can be done recursively by calling calling the multiplication for  [latex]\mathbb{Z}[x] / (x^{\sqrt{n}} + 1)[/latex]. Following this plan, a multiplication in  [latex]\mathbb{Z}[x] / (x^{n} + 1)[/latex], generates [latex]\sqrt{n}[/latex] multiplications in [latex]\mathbb{Z}[x] / (x^{\sqrt{n}} + 1)[/latex], each of which generates [latex]\sqrt{n}[/latex] multiplications in [latex]\mathbb{Z}[x] / (x^{\sqrt{\sqrt{n}}} + 1)[/latex], and so on. Since taking square roots [latex]\log \log (n)[/latex] gets to [latex]O(1)[/latex], this gives an algorithm of complexity [latex]O(n \log n \log \log n)[/latex].   In order to make this plan work, we need to care about various implementation details. The first step is to reduce the problem to multiplying in the ring [latex]\mathbb{Z}[x] / (x^{n} + 1)[/latex]. To facilitate recursion we will make [latex]n[/latex] be a power of [latex]2[/latex] in the following way:   <pre class="prettyprint">  function r = multiply_poly(p,q)    n = length(p) + length(q) - 1;    [p,k] = expand_to_power_of_2(p,n);    [q,k] = expand_to_power_of_2(q,n);    r = multiply_poly_in_Zxn(p,q,k);    r = resize(r, n, 1);  endfunction    function [p,k] = expand_to_power_of_2(p, n)    p = reshape(p, length(p), 1);    plen = 1; k = 0;    while (plen < n)      plen = plen * 2;      k = k + 1;    endwhile    p = resize(p, plen, 1);  endfunction  </pre>  Ok, now we need to provide an implementation for <code>multiply_poly_in_Zxn(p,q,k)(</code> which multiplies two polynomials in [latex]\mathbb{Z}[x] / (x^{n} + 1) where n = 2^k. The first step is to write each p and q as polynomials with coefficients in \mathbb{Z}[x] / (x^{2K} + 1) for some K = O(\sqrt{n}). We will do as follows: let L = 2^{\lfloor (k-1)/2\rfloor} and K = 2^{\lceil (k-1)/2\rceil} and write: p(x) = \sum_{i=0}^{2L-1} (\sum_{j=0}^{K-1} p_{Ki+j} x^j) x^{Ki} and let's think of $later \bar{p}_i = \sum_{j=0}^{K-1} p_{Ki+j} x^j$ as a polynomial in \mathbb{Z}[x] / (x^{2K} + 1). Now, we need an element in \mathbb{Z}[x] / (x^{2K} + 1) that will serve as a 2L-root of the unit (since this is the degree of the polynomials we are trying to multiply). Notice that x^{4K} = 1 If K = L, then \omega = x^2 is a 2L-root of the unit. If K = 2L, then \omega = x^4 is a 2L-root of the unit. This makes us almost ready to apply FFT. There is only one catch left and this is very subtle.

A brief intermission: given a \in \mathbb{Z}, there is the evaluation map E_a : \mathbb{Z}[x] \rightarrow \mathbb{Z} that maps p(x) \mapsto p(a) \in \mathbb{Z}, simply by evaluating the polynomial. Now, say we want to extend this map to \mathbb{Z}[x]/(q(x)). The elements of \mathbb{Z}[x]/(q(x)) are (infinite) sets of polynomial which are equivalent modulo q(x). If they all the polynomials in the equivalence class \bar{q} evaluate to the same, then this defines a very evaluation map. It is simple to see that this happens iff a is a root of q(x), i.e., q(a) = 0. In such a case, the evaluation map extends naturally to \bar{E}_a : \mathbb{Z}[x]/(q(x)) \rightarrow \mathbb{Z}. If a is not a root, then there is no evaluation map \bar{E}_a such that E_a = \bar{E}_a \circ \pi, where \pi is the projection \pi: \mathbb{Z}[x] \rightarrow \mathbb{Z}[x]/(q(x)) .

We are looking at polynomials as elements in A[y] / (y^{2K} + 1) where A = \mathbb{Z}[x] / (x^{2K} + 1). So each element is in fact a set of polynomials, i.e., all polynomials equivalent modulo y^{2K} + 1. If we want to evaluate them at a certain point \omega we must make sure all polynomials evaluate the same. Here is how to do it: define \psi = \sqrt{\omega}, i.e., \psi = x if K = L and \psi = x^2 if K = 2L. Then the map \bar{p}(y) \mapsto \bar{p}(\psi y) maps from A[y] / (y^{2K} + 1) to A[y] / (y^{2K} - 1), in which \omega is a root. So, now we are fine. Notice that we can also map back using the inverse \bar{p}(y) \mapsto \bar{p}(\psi^{-1} y) = \bar{p}((x^{4K} /\psi) y). This point is really subtle and it was the one that took me more time to understand. After we understand why this needs to be done, things get really easy: we just need to do FFT as before, in other words, evaluate in 1, \omega, \hdots, \omega^{2L-1} and then reconstruct using the inverse FFT. The code in Octave is below:

function r = multiply_poly_in_Zxn(p,q,k)
  if k<=2
  L = 2^floor((k-1)/2); K = 2^(k-1) / L; psi = 2;
  if K==L psi = 1; endif

  % map Z[x] --> (Z_{2K})_{2L}
  p = [ reshape(p, K, 2*L); zeros(K, 2*L) ];
  q = [ reshape(q, K, 2*L); zeros(K, 2*L) ];
  for i = 0 : 2*L - 1
    % map Z_{2L}[y] / (y^k+1)  -->  Z_{2L}[y] / (y^k - 1)
    % implemented via   p(y)  |-->  p(psi*y)
    p(:,i+1) = multiply_mono(p(:,i+1), i*psi);
    q(:,i+1) = multiply_mono(q(:,i+1), i*psi);
  w_p = fft(p, 2*psi);
  w_q = fft(q, 2*psi);
  w_r = [];
  for i = 1:2*L
    w_r(:,i) = multiply_poly_in_Zxn(w_p(:,i), w_q(:,i), 
  r = fft_inv(w_r, 2*psi);
  rr = [];
  for i = 0 : 2*L - 1
    % back to Z_{2L}[y] / (y^{2L} +1)
    rr(:,i+1) = multiply_mono(r(:,i+1), i*(4*K-psi));
  r = rr;
  r1 = reshape(r(:,1:2:2*L), 2*K*L, 1);
  r2 = reshape(r(:,2:2:2*L), 2*K*L, 1);
  r2 = [-r2(2*K*L-K+1:2*K*L); r2(1:2*K*L-K)];
  r = r1 + r2;

The auxiliary methods are below:

function r = naive_multiply_poly(p,q)
   n = length(p); r = zeros(n, 1);
   for i=0:n-1 for j =0:n-1
      r(1+mod(i+j,n)) = r(1+mod(i+j,n))
                        + (-1)^(floor((i+j)/n)) * p(1+i) * q(1+j);
   endfor endfor

function r = multiply_mono(p, t)
  n = length(p);
  r = zeros(n, 1);
  for j = 0:n-1
    r(1+mod(j+t,n)) += p(j+1) * (-1)^(floor((j+t)/n));

function w = fft_aux(p, omega)
  n = length(p(1,:));
  if (n == 1)
    w = p;
  p_even = p(:, 1:2:n); w_even = fft_aux(p_even, 2*omega);
  p_odd =  p(:, 2:2:n); w_odd  = fft_aux(p_odd,  2*omega);
  for j = 0 : n-1
    w(:, j+1) = w_even(:, mod(j,n/2)+1)
                + multiply_mono(w_odd(:,mod(j,n/2)+1), j*omega);

function w = fft(p, omega)
  w = fft_aux(p, omega);

function w = fft_inv(p, omega)
  n = length(p(1,:));
  w = fft_aux(p, (n-1)*omega) ./ n;

Notice however, that the implementation provided above is not super efficient since it does data copies, memory allocation and recursive calls. The implementation above was optimized not for speed but for making the algorithm clear.

Solving Linear Systems and Inverting a Matrix is equivalent to Matrix Multiplication

Matrix multiplication consists on given two matrices A and B, say for now both of dimensions n \times n, compute C = A \cdot B with C_{ij} = \sum_k A_{ik} \cdot B_{kj}. The definition implies a naive algorithm of running-time O(n^3) for multiplying two matrices. When I first heard about that, I was somewhat surprised that one can multiply matrices faster then that. The simpler algorithm beating the naive O(n^3) was given by Strassen (1969) and has running-time O(n^{\log_2 7}) and is based on a simple observation. Say n = 2^k and break each matrix in 4 blocks of dimensions 2^{k-1} \times 2^{k-1}. So we get something like that:

 \begin{bmatrix} C_{11} & C_{12} \\ C_{21} & C_{22} \end{bmatrix}  = \begin{bmatrix} A_{11} & A_{12} \\ A_{21} & A_{22} \end{bmatrix} \cdot  \begin{bmatrix} B_{11} & B_{12} \\ B_{21} & B_{22} \end{bmatrix} = \begin{bmatrix} A_{11} B_{11} + A_{12} B_{21}  & A_{11} B_{12} + A_{12} B_{22}  \\ A_{21} B_{11} + A_{22} B_{21}  & A_{21} B_{12} + A_{22} B_{22}  \end{bmatrix}
which are 8 multiplications of smaller matrices. Strassen found a recursion that computes the C_{ij} factors using 7 matrix multiplication, which leads to O(n^{\log_2 7}). This number has been improved various times by Coppersmith and Winograd, Stothers, Vassilevska-Williams and Le Gall. The current best algorithm has running-time O(n^{2.3728}). No lower bound is known for this problem other then O(n^2) which is the time needed to read the matrices. The constant \omega is defined as the minimum constant for which it is possible to multiply matrices in time O(n^\omega). So far it is known that 2 \leq \omega \leq 2.3728....

From multiplying matrices to solving linear systems

What I learned this week is that solving linear systems, inverting a matrix and computing determinants can be done also in time O(n^\omega) with an ingenious idea by Hopcroft and Bunch (1974). Next we show how to do it following the presentation in Mucha and Sankowski (2004). In one line, the idea consists in performing a lazy LU decomposition.

The LU decomposition is another way to refer to the Gaussian elimination procedure. We decompose a matrix A in the form A = L U where L is a lower-diagonal matrix with ones in the diagonal and U is an upper-diagonal matrix with ones in the diagonal. Here are some great lecture notes on the topic, but here is a concise description. The following code is in Octave:

function [L,A] = naive_lu(A)
  n = length(A);
  L = eye(n);
  for i = 1 : n
    L(i+1:n,i) = A(i+1:n,i) / A(i,i);
    A(i+1:n,i:n) = A(i+1:n,i:n) - L(i+1:n,i)*A(i,i:n);

This implementation is very naive in the sense that it doesn’t do pivoting, it doesn’t work for non-singular matrices and so on, but it should be fine for a random matrix, for example. This captures the essence of Gaussian elimination: one makes the matrix triangular by eliminating rows one-by-one. The complexity of the naive implementation is O(n^3) since in each step we update a matrix of size n-i \times n-i+1.

The great thing of having an LU decomposition is that now solving Ax = b is very easy, since solving LUx = b involves solving L y = b and then U x = y. Since L and U are both triangular, this can be solved in O(n^2) time. Also, we can compute determinants easily, since det(A) = det(L) \cdot \det(U) = \prod_{i=1}^n u_{ii}. The inverse of A is given by A^{-1} = U^{-1} \cdot L^{-1} which is twice the time to invert a triangular matrix plus the time to multiply to matrices. The following argument shows that inverting a triangular matrix can be done in the same time as multiplying two matrices: consider an upper triangular matrix U, for example, then:

 \begin{bmatrix} U_{11} & U_{12} \\  & U_{22} \end{bmatrix}  \cdot \begin{bmatrix} Y_{11} & Y_{12} \\  & Y_{22} \end{bmatrix} = \begin{bmatrix} I & 0 \\ 0 & I \end{bmatrix}

So, we can solve this system by solving Y_{11} = U_{11}^{-1}, Y_{22} = U_{22}^{-1} and since U_{11} Y_{12} + U_{12} Y_{22} = 0, by pre-multiplying this system by Y_{11}, we get Y_{12} = - Y_{11} U_{12} Y_{22}. Therefore, the time I(n) to invert a triangular n \times n corresponds to two inversions of \frac{n}{2} matrices plus two matrix multiplications, i.e., I(n) = 2 I(\frac{n}{2}) + O(n^\omega). Solving the recursion we get I(n) = O(n^\omega).

The previous arguments reduce the problem of computing determinants, solving linear systems and computing inverses to the problem of computing an LU decomposition in O(n^\omega) time.

Lazy LU decomposition

Note that the expensive step of the LU decomposition consists in updating the matrix A by subtracting a rank one matrix. The idea behind the lazy LU decomposition is to postpone the updates such that they are done many row updates are done together in an exponential fashion, i.e., for each j \geq 1 a 2^{-j} fraction of the updates modifies only 2^{j-1} rows. The following picture show the traditional LU decomposition:


The k-th iteration, uses entry A_{ii} in order to zero the entries in the blue block marked with k. The Hopcroft-Bunch version of the LU factorization zeros the entries in a different fashion, shown in the following picture:


The first few steps work like that: the first step will be almost like the traditional decomposition except that we will use A_{11} to make entry A_{21} zero (instead of the entire 1-column): this is done by setting B \leftarrow B - A C^{-1} D and A \leftarrow 0. Notice that this is the same as doing [A,B] \leftarrow [A,B] - A C^{-1} [C,D].


Next we do the same thing for block 2, we perform the same matrix operation: B \leftarrow B - A C^{-1} D. Notice that at this point C will be an upper triangular non-singular matrix. We continue the same procedure until we the entire matrix become upper-diagonal. We store the modifications we performed (i.e., the matrices A C^{-1} in the lower diagonal matrix L in the same way we did or the usual LU decomposition. The few next steps are:


Generalizing this for an matrix of size 2^k \times 2^k we notice that there are 2^{k-j} operations with matrices of size 2^j for j = 1 to k-1. Each iteration consists of inverting a O(2^{j\omega}) and a multiplication of a 2^j \times 2^j matrix by a 2^j \times n matrix, which takes time \frac{n}{2^j} O(2^{j\omega}), since this can be decomposed in \frac{n}{2^j} multiplications of matrices of size 2^j \times 2^j. The overall time is then \sum_{j=1}^{\log(n)} \frac{n}{2^j} \cdot \frac{n}{2^j} O(2^{j\omega}) = n^2 \sum_j O(2^{j(\omega - 2)}) = n^2 \cdot 2^{(\omega-2)\log n} = O(n^\omega) .

The code for the lazy LU decomposition is surprisingly simple and similar to the code for the traditional LU decomposition. Essentially in each iteration i, we compute the size of the block g as the maximum value of 2^j such that 2^j divides i (in the code this is done using a bitwise operation).

function [L,A] = lazy_lu(A)
 n = length(A);
 L = eye(n);
 for i = 1 : n-1
   g = bitxor(i, bitand(i, i-1));
   L(i+1:i+g,i-g+1:i) =  
     A(i+1:i+g,i-g+1:i) *  inv(A(i-g+1:i,i-g+1:i)) ;
   A(i+1:i+g,i-g+1:n) =  
     A(i+1:i+g,i-g+1:n) - L(i+1:i+g,i-g+1:i) * A(i-g+1:i,i-g+1:n);

The Octave code has running time O(n^\omega) where \omega is the exponent associated with the matrix multiplication algorithm used. Notice that above we are still using Octave’s native matrix multiplication and inversion, which are O(n^3). In order to improve on that, one would need to substitute the inversion by a custom inversion (following the receipt given earlier in this post) and the multiplication operations A(i+1:i+g,i-g+1:i) * inv(A(i-g+1:i,i-g+1:i)) and L(i+1:i+g,i-g+1:i) * A(i-g+1:i,i-g+1:n) by custom matrix multiplication.

Also, the code below should produce a correct decomposition for a random matrix, which is guaranteed to be well behaved. For specific matrices (like the matrix with all zeros and one in the secondary diagonal), the LU will fail since it will try to invert a singular matrix (the same way that a naive LU decomposition — our first algorithm in this page — will fail). The solution is to do pivoting and write an LUP decomposition, where L and U are as before and P is a permutation matrix. Pivoting here can be done as usual, here is a code snippet that is guaranteed to work for any non-singular matrix:

function [L,A,P] = lazy_lu_with_pivoting(A)
 n = length(A);
 L = eye(n);
 perm = [1:n];
 for i = 1 : n-1
  g = grow(i);
  f = min(i+g, n);

  # pivoting step
  [_, piv] = max(A(i,i:end)); piv = piv + i - 1;
  if (i != piv)
   temp = A(:,i); A(:,i) = A(:,piv); A(:,piv) = temp;
   temp = perm(i); perm(i) = perm(piv); perm(piv) = temp;

  L(i+1:f,i-g+1:i) =  
    A(i+1:f,i-g+1:i) *  inv(A(i-g+1:i,i-g+1:i)) ;
  A(i+1:f,i-g+1:n) = 
    A(i+1:f,i-g+1:n) - L(i+1:f,i-g+1:i) * A(i-g+1:i, i-g+1:n);
 P = eye(n)(perm, :);

Lastly, one might wonder if solving a linear system or inverting a matrix is perhaps easier then multiplying two matrices. Hopcroft and Bunch point out the following fact that shows that inverting a matrix in O(n^\omega) implies an algorithm for multiplying matrices in the same time:

 \begin{bmatrix} I & A & 0 \\ 0 & I & B \\ 0 & 0 & I \end{bmatrix}^{-1}  = \begin{bmatrix} I & -A & AB \\ 0 & I & -B \\ 0 & 0 & I  \end{bmatrix}

Tokyo Food Diary

Tokyo is the city that lives in my heart. Last week Carol and I took vacations and went to Japan ! Tokyo is like a paradise: first of all, food is awesome ! Very unlike NYC, everywhere you go, food is good. There is not much need to choose much, you are choosing between great and super awesome food. It may be that if you compare the best restaurants in Tokyo with the best restaurants in NYC or Paris, they might be similar. But comparing a random place, there is probably nowhere that can beat Tokyo. That being said, a guidebook is not so necessary, but I got two anyway: the first was the Rough Guide to Tokyo, which turned out not very useful and we ended up leaving it in the hotel most of the days. We had Google Maps (T-mobile now has free worldwide roaming) and the main reason to have a Rough Guide was to orient in the city, was mostly gone. We still checked it every morning to read about the neighborhood we were visiting, but I believe I could get a similar thing from Wikipedia or Wikitravel. The second guide we bought was called Food Sake Tokyo and it was a guidebook only about food. This one was super useful and we carried it most of the time. Guide was nice to navigate the options Tokyo offered and the variety of cuisines in Japan. Half is restaurant recommendations, and half is just a book about Japanese food, specially things outside the usual or hidden restaurants that were very popular. In the end we used this guide half the time and half the time we picked places at random. Another book that I read was Pretty Good Number One. It is a travel diary of an American family (a couple and their 8-year old daughter) spending one month in Tokyo and falling in love with Japanese food. This is not a guidebook but I found a really nice inspiration for a trip.

We arrived in Japan in the night of April 12. Sort of anxious, since we wanted to see cherry blossoms. A week before, a friend of mine living in Tokyo wrote saying that trees were blossoming, I should come quickly. I thought I’d arrive and they would be over. The first thing upon we arrived was to go to the gardens in Shinjuku Gyoen to see if we could still spot some cherry blossoms. We were super happy when saw the park was still full of cherry blossoms. Some of the trees are late bloomers and the flowers appear one week after the first cherry blossoms. We then went to the park, bought some sweets with cherry blossom themes, of course, and filled with cherry boossom (sakura) paste and stayed there for appreciating the flowers. Here are us and the sweets:


As the title of this blog post is saying, this was intended to be a food diary in which I describe what we ate in Japan. So let’s go for it. Forgive me if I sound like a snob, but food in Tokyo is just too awesome not to write about it). Here we go:

Sunday: After spending the morning in the park, we start to get hungry (even after all those cherry blossom sweets). So we go to Nihombashi to eat eel. I am fascinated by eel and the first thing I felt like doing in Japan was to try how eel tasted like in Japan. Most of the time I eat eel in the US is unagi, i.e. fresh water eel. So we went to try anago, the leaner sea-water eel in a restaurant in Nihombashi (Anago Tamai) which served a single dish, which was anago grilled over rice. They served together with tea and recommended to eat the final 1/4 of the dish as ochazuke.


After lunch, we went to the depachika in a nearby department store in Nihombashi. A depachika is a floor devoted to food. We spent at least one hour browsing eating all kinds of sweets. The stores are usually divided in a section with Japanese sweets (Wagashi), sweets made of beans, green tea and fruits and a section with Western sweets (Yogashi), mostly French.

For dinner we met Leandro Hideo, who took us to an Izakaya specialized in food from Hokkaido, the northern region from Japan. Besides getting sashimi (in particular, we got uni, which I was super excited to eat) and a crab (which I believe it was a snow-crab). Here is me, Hideo and the crab.


Monday: After visiting the Zozo-ji in the morning, we walked to the Tsukiji Fish Market to have an omakase sushi lunch at Sushi Bun. We waited in line for a while, but it was definitely worth it. Fish was super fresh and nice. In the evening, we walked in Ginza and went to a place at random. It was a yakitori izakaya specialized in chicken yakitori. We got a tasting menu, which brought us skewers with some less conventional parts of the chicken. As a Brazilian, I am a big fan of chicken hearts. I also learned in this meal that chicken necks can be super tasty.

Tuesday: We got some really good street food near the meiji shrine. First, we got takoyaki, little octopus balls — which is my all-time-favorite fried thing. Second we got some pastries in the shape of a fish filled with sweet red beans. Finally, we sampled a bunch of sweet shops. The simple lunch sort of counter-balanced the dinner, in which we indulged outselves in a kaiseki dinner . We went to Tsukiji Tamura, just outside the fish market. Kaiseki meals are multi-course menus, known for the careful presentation and seasonal ingredients. In total we had eight courses, which were more of less like that: (1) the first course resembled an appetizer. It was a round wooden dish with a piece of raw fish, a small cooked fish, a small bowl with two small squids, a skewer with vegetables, a piece of beef, some cooked roots, … (2) second dish was a tofu dish served in a bowl with a lid. (3) third dish was a bowl filled with ice and sashimi on top (3) a piece of broiled fish and snail sashimi, eaten directly from the shell. Here I was lucky to eat my own snail and the one Carol’s  dish. (4) a vegetable course; (5) a second vegetable course; (6) a hearty rice dish. (7) a light fruity desert (8) a sweet bean soup; (9) a cup of matcha.

Wednesday: On Wednesday we took the train to Kamakura, a city 1 hour from Tokyo known for its many temples and shrines. Food-wise the city is well-known for the Zen cuisine, very light and based on vegetables and rice. Carol got for lunch rice with curry and vegetables. I got tamago gohan, which is a dish of rice mixed with a raw egg, spices, soy souce and sesame oil. It was quite nice, specially because I got a bowl of rice, spices, oil, sauce separate, and an unbroken egg and assembled the dish myself. Later we further got some deserts in the city:


In the evening, we got okonomiyaki at the train station — which is a pancake with pork (it actually can contain pretty much anything, but the one we got had pork inside) and a bunch of things on top, like mayonnaise, mustard, bonito flakes, …

Thursday: For lunch we got Chirashi in Kagurazawa (Futaba) — which is a bowl of rice with small pieces of fish on top. The restaurant again had only one dish and one price (which was around $15 dollars). Having only one option is great, specially in Japan. For two reasons: they really focus on preparing this option perfectly: the ingredients are fresh and seasonal. Besides, you don’t need to decide — which is hard since all the menus and explanations are in Japanese.  They asked us to pay to a cashier, showed a table and brought our bowl of chirashi:DSC_0564Our final dinner in Japan we chose to walk around the neighborhood we were staying (Nishi-shinjuku) and sit in the first sushi place we found. As it was Japan, it was great !

Friday: In the last day, we had to catch our flight back. We had a little time, so we walked around Roppongi and Azabu. We finish with a fantastic view of the city we got from the Roppongi Hills:

IMG_20140418_122100Here are our travel pictures.




Thermometer (using Arduino and RaspberryPi)



Here is a small project I did to play with my RaspberryPi and my Arduino. I bought both a while ago but didn’t get much chance to play with them. The plan was as follows: to capture the temperature in my living room using a temperature sensor using the arduino, make the arduino write to the Serial port to the RaspberryPi, have a Python program running in the RaspberryPi that reads from the Serial port and then stores the info in an AppEngine Server. In the end, you can point your browser to a certain address and get a plot of the temperature in my living room:




  • arduino and TO-92 (the sensor): here is a guide how to physically connect the sensor to the arduino.
  • raspberrypi and arduino are connected via an USB cable
  • raspberrypi is powered using a power cord, connected to the internet using an Ethernet cable and has an SD card with Raspbian pre-loaded (if the card came with the Raspberry, it might have already the Raspbian or similar distribution pre-loaded).

Logging in the RaspberryPi

If you have a monitor and a keyboard, you can plug them to the Raspberry. If like me, you don’t, the easiest thing to do is to ssh to it. It is really simple. Make sure you are in the same local network and follow the instructions here. Here is quick summary of what is there: do

ifconfig -a

to find out the address from your wlan. If you have a linksys router, then it will be something like 192.168.1.x. If so, then do the following to figure out the IP of the devices connected to your local network:

sudo nmap -sP

Now, you just try to ssh to all of them. For one of them, the following will work:

ssh pi@

If it works, it will ask for a password. Typically, the initial password for Raspbian is ‘raspberry’. Then, you will get something like:

pi@raspberrypi ~ $

Ok, now we are logged in the Raspberry. It will be useful to have both python and ino installed. Ino is a tool used to compile and upload it to arduino.


Code for arduino

The first step is to use ino init to create a temperature project. You can follow the Ino Quickstart here. Here is the sketch.ino that I wrote to read from the temperature sensor and write to the serial port:

void setup() {

void loop() {
  int valT = analogRead(0);
  float mV = (valT / 1024.0) * 5000;
  float temperature = (mV - 500)  / 10;

  Serial.print("T: ");

  delay(300000); //5 minutes

Note that I am using the analog pin A0 to read from the sensor. Then we do ino build and ino upload (as described in the quickstart). And at that point the arduino should start sending messages with temperature to the Serial port. Notice I added a delay(300000); so that the arduino sends a temperature every 5 min (300000 miliseconds). If you are following this, you can initially start with a smaller interval to make sure everything is working (say a couple of seconds).


Reading from the serial port

Code now should be up and running in the arduino and it should be writing to the serial port. Now, let’s read from it. By typing

pi@raspberrypi ~ $ino serial

you can see the output of the serial port and see if it is doing the right thing. You will also see what the serial port actually is. In my case it is /dev/ttyACM0. Now, let’s see how to read from it in python:

from serial import Serial
import re

serial_pattern = r"T: (\d+\.\d*)\n";
serial_port = '/dev/ttyACM0';
serial_bauds = 9600;

def open_serial_port() :
  s = Serial(serial_port, serial_bauds);
  line = s.readline();
  return s

def read_temperature(s):
  line = s.readline();
  m = re.match(serial_pattern, line);
  return float(


AppEngine Server

Now we are ready to create a server in the Google AppEngine and start uploading the temperature we read to it. This tutorial will walk you through the main steps in setting up a server. The server we will set up is actually simpler the the one in the example. I did mine in Go, as I was trying to learn the language (check Golang Tour, it is real fun), but you can do the tutorial in Python as well. Let’s stick with Go. You need essentially three files in your server: app.yaml (which is pretty much standard — see the tutorial), index.yaml (which is autogenerated, so no need to worry about this one either) and the go file that actually specifies the server. I called mine hello.go. The file is longer, but you can find it in this link.

I might dicuss it in more details in a further post (for example, how I am using the Google visualization api), but for now, I will just post the python code I am running in the RaspberryPi that is communicating with the appengine server.

from serial import Serial
import re
import urllib

serial_pattern = r"T: (\d+\.\d*)\n";
serial_port = '/dev/ttyACM0';
serial_bauds = 9600;
url = r''

if __name__ == "__main__":
  print "Opening serial port"
  s = Serial(serial_port, serial_bauds);
  print "Reading first line from port"
  line = s.readline();
  print "Initializing communication"
  while 1:
    line = s.readline()
    m = re.match(serial_pattern, line);
    params0 = {"temperature", "source":"sensor_living_room"}
    params1 = urllib.urlencode(params0)
    print params1
    data = urllib.urlopen(url, params1).read()


Temperature plot

And finally, you can check the temperature here or the tables here.


Final comments

One last comment: there are probably easier (and more elegant ways) of doing those things — for example, it is possible to read sensor data directly from the RaspberryPi without using the Arduino. This small project came out of me trying to learn all those things (Arduino, RaspberryPi and Go) at the same time. One advice that I got once was that the best way to learn something is to come up with a small project that uses all those things and implement it.

renatoppl . stl

Renato from

Yesterday I went to the Makerbot store in Soho and took a picture in their 3D photo booth. Store is awesome and definitely worth a visit next time you come to NYC. They have a bunch of 3D printers printing all sorts of cool stuff (Manhattan skylines, a tiny model of NYC inside an apple, house maquettes, …). They print a bunch of bracelets that they give away to visitors — Carol got one. After us, a small girl walked in the store and they were out of them: so they said, would you wait 5 min, we are almost finishing to print another one: in 5 min, she took is fresh out of one of the printers in the display and gave the bracelet to the girl.

Now that I have a 3d model of myself, I am getting encouraged to go over the blender tutorials to learn a bit of 3d modeling so that I can do things like insert my face in a monkey or make a myself-themed coffee cup. Here is the full gallery of all the pictures me and Carol took there ($5 give you 3 pictures).

Some games you might want to play (even if you don’t usually play games)

I used to be a big fan of games when I was a child — come back home and spend the whole time playing. After high school I didn’t play much, but once in a while you learn about a game that is so good that you want to spend some hours playing again. Let me tell you about two games I learned about recently and well worth playing (even if games are not your thing). Those two are the best games I ever played — they are both short (one of them you can finish in one hour and the other one you can probably finish in a total of four hours or so). Both have great stories and deep philosophy behind. They are really art masterpieces that were incredibly well though,  where things have hidden meanings, challenges your perception of reality.

  • Braid: the game goes consists of various phases with titles such as “Time and Forgiveness”, “Time and Mystery”, “Time and Place”, … Each involves solving puzzles that have to do with controlling the flow of time in different ways. The puzzles don’t require you to be able with games or quickly (it never asks you to press this button in the right time), but require thinking the possibilities very precisely. The story of the game (Tim trying to rescue a princess) is by itself quite interesting and mysterious and involves various layers with different meanings. I also can say the ‘last’ world of the game is the most ingenious thing I ever played.
  • Journey:when I went to buy Journey in Amazon, I saw tons of incredibly enthusiastic reviews. one of them said that Journey was more then a game and more like a spiritual journey. Other review asked me not to read reviews any further, since the less you know about it before you start the better. I agree with the comment. So, let me only describe the first minute playing.  You start in the middle of a desert with no instructions. You can move, jump and press one button that makes a sound, but not clear what is this for. There is a city or a mountain in the horizon and going in that direction makes sense. You finish the entire game in one hour or so but you start understanding it after you play it a couple of times. I also recommend to play it online — this makes a big difference for this game.

Also, if you take airplanes and are bored in airports, here is a game you might consider playing: Don’t Starve. I found it really addictive (maybe too addictive, I’d say).

Engagement ring

In my last post, I talked about my recent adventures in 3d printing and I one idea was to illustrate it using my first project: designing the engagement ring with which I proposed to my (by-then girlfriend and now) wife. A 3d printed ring was perfect for both us because it married my interest in computer science and her in material science (her PhD is on polymer science and essentially about plastic). So it had all to be a great success. Here we go with the project:

  • ring mesh: first I wrote a python script that would generate the mesh around the ring. the mesh is composed by various hexagons rotates with respect to each other. (in the process I learned — or was forced to remind myself — that a consistent triangulation must have all the normals pointed in the same direction. enforcing this was tricky, but quite fun to figure out how to do it).



  • the bunny: first I played around with diamond-like shapes, but I decided to go with a bunny ! So I downloaded the Stanford bunny from Thingverse.
  • the letters: and at least, I wanted to inscribe our names in the ring, so I chose to use write the first two letters of our names: R (Renato) and C (Carolina). As you might have already realized, our first letters are also the names of two of the most important rings:
    R  and C . (I also told Carol that she is my algebraic completion). I thought my reader would enjoy the pun of two rings  in a ring. So I chose the letters in mathbb to write in the ring just opposite to the bunny.

Then I brought all those together in then OpenSCad file and compiled the ring:



Then I generated the STL file. (See my previous blog post to see how to do it.) And also a box, since a ring needs a box — that also with our mathbb letters.

makerbot_screenshotand voilà: