The example below computes sin²(x) + cos²(x) by computing the sine function in R, the cosine function in Python, and summing their squares in Perl. As you’d hope, it returns 1. (Actually, it returns 0.99999999999985 on my machine.)

To execute the code, go to the #+call line and type C-c C-c.

#+name: sin_r(x=1)
#+begin_src R
sin(x)
#+end_src

#+name: cos_p(x=0)
#+begin_src python
import math
return math.cos(x)
#+end_src

#+name: sum_sq(a = 0, b = 0)
#+begin_src perl
$a*$a + $b*$b;
#+end_src

#+call: sum_sq(sin_r(1), cos_p(1))

Apparently each function argument has to have a default value. If that’s documented, I missed it. I gave the sine and cosine functions default values that would cause the call to sum_sq to return more than 1 if the defaults were used.

Here I’ll show some of the most basic possibilities. For much more information, see  orgmode.org. And for the use of org-mode in research, see A Multi-Language Computing Environment for Literate Programming and Reproducible Research.

Source code blocks go between lines of the form

#+begin_src
#+end_src

On the #+begin_src line, specify the programming language. Here I’ll demonstrate Python and R, but org-mode currently supports C++, Java, Perl, etc. for a total of 35 languages.

Suppose we want to compute √42 using R.

#+begin_src R
sqrt(42)
#+end_src

If we put the cursor somewhere in the code block and type C-c C-c, org-mode will add these lines:

#+results:
: 6.48074069840786

Now suppose we do the same with Python:

#+begin_src python
from math import sqrt
sqrt(42)
#+end_src

This time we get disappointing results:

#+results:
: None

What happened? The org-mode manual explains:

… code should be written as if it were the body of such a function. In particular, note that Python does not automatically return a value from a function unless a return statement is present, and so a ‘return’ statement will usually be required in Python.

If we change sqrt(42) to return sqrt(42) then we get the same result that we got when using R.

By default, evaluating a block of code returns a single result. If you want to see the output as if you were interactively using Python from the REPL, you can add :results output :session following the language name.

#+begin_src python :results output :session
print "There are %d hours in a week." % (7*24)
2**10
#+end_src

This produces the lines

#+results:
: There are 168 hours in a week.
: 1024

Without the :session tag, the second line would not appear because there was no print statement.

I had to do a couple things before I could get the examples above to work. First, I had to upgrade org-mode. The version of org-mode that shipped with Emacs 23.3 was quite out of date. Second, the only language you can run by default is Emacs Lisp. You have to turn on support for other languages in your .emacs file. Here’s the code to turn on support for Python and R.

(org-babel-do-load-languages
    'org-babel-load-languages '((python . t) (R . t)))

Update: My next post shows how to call code in written in one language from code written in another language.

Related posts:

Personal organization software
Preventing an unpleasant Sweave surprise

1. Rank by number of questions tagged with that language on Stack Overflow
2. Rank by number of project on GitHub using that language

According to this article, these two measures are well correlated.

I’d be skeptical of either metric by itself. A large number of questions on a language could indicate that it’s poorly documented, for example, rather than popular. And GitHub projects may not representative. But the two measures give similar pictures of the programming language landscape, so together they have more credibility. On the other hand, both measures are probably biased in favor newer languages.

The RedMonk Programming Language Rankings: February 2012

If a matrix A is large and sparse, it may be possible to solve Ax = b but impossible to even store the matrix A^-1 because there isn’t enough memory to hold it. Sometimes it’s sufficient to be able to form matrix-vector products Ax. Notice that this doesn’t mean you have to store the matrix A; you have to produce the product Ax as if you had stored the matrix A and multiplied it by x.

Very often there are physical reasons why the matrix A is sparse, i.e. most of its entries are zero and there is an exploitable pattern to the non-zero entries. There may be plenty of memory to store the non-zero elements of A, even though there would not be enough memory to store the entire matrix. Also, it may be possible to compute Ax much faster than it would be if you were to march along the full matrix, multiplying and adding a lot of zeros.

Iterative methods of solving Ax = b, such as the conjugate gradient method, create a sequence of approximations that converge (in theory) to the exact solution. These methods require forming products Ax and updating x as a result. These methods might be very useful for a couple reasons.

1. You only have to form products of a sparse matrix and a vector.
2. If don’t need a very accurate solution, you may be able to stop very early.

In Newton’s optimization method, you have to solve a linear system in order to find a search direction. In practice this system is often large and sparse. The ultimate goal of Newton’s method is to minimize a function, not to find perfect search directions. So you can save time by finding only approximately solutions to the problem of finding search directions. Maybe an exact solution would in theory take 100,000 iterations, but you can stop after only 10 iterations! This is the idea behind the Newton-Conjugate Gradient optimization method.

The function scipy.optimize.fmin_ncg can take as an argument a function fhess that computes the Hessian matrix H of the objective function. But more importantly, it lets you provide instead a function fhess_p that computes the product of the H with a vector. You don’t have to supply the actual Hessian matrix because the fmin_ncg method doesn’t need it. It only needs a way to compute matrix-vector products Hx to find approximate Newton search directions.

For more information, see the SciPy documentation for fmin_ncg.

And here’s a parody of it that came out today:

Before the play started, someone told me that the phrase “bidi-bidi-bum” in “If I Were a Rich Man” is a Yiddish term for prayer. I thought “All day long I’d bidi-bidi-bum” was a way of saying “All day long I’d piddle around.” That completely changes the meaning of that part of the song.

When I got home I did a quick search to see whether what I’d heard was correct. According to Wikipedia,

A repeated phrase throughout the song, “all day long I’d bidi-bidi-bum,” is often misunderstood to refer to Tevye’s desire not to have to work. However, the phrase “bidi-bidi-bum” is a reference to the practice of Jewish prayer, in particular davening.

Unfortunately, Wikipedia adds a footnote saying “citation needed,” so I still have some doubt whether this explanation is correct. I searched a little more, but haven’t found anything more authoritative.

Now I wonder whether there’s any significance to other parts of the song that I thought were just a form of Klezmer scat singing, e.g. “yubba dibby dibby dibby dibby dibby dibby dum.” I assumed those were nonsense syllables, but is there some significance to them?

Update: At Jason Fruit’s suggestion in the comments, I asked about this on judaism.stackexchange.com. Isaac Moses replied that the answer is somewhere in between. The specific syllables are not meaningful, but they are intended to be reminiscent of the kind of improvisation a cantor might do in singing a prayer.

The idea behind ceiling of complexity is that every project you complete creates residual responsibilities and expectations. This residual may be small, maybe not even noticeable, but it’s always there. Over time, this residue builds up and adds complexity. Eventually it forms a ceiling and limits further progress until you do something to break through the ceiling and reach a new state of simplicity. The ceiling of complexity is a byproduct of success.

Sullivan’s picture of a ceiling of complexity is sort of existential crisis, something an individual would only face a few times over a career, but I find it useful to use the term for less dramatic situations. It gives a way to talk about the gradual accumulation of small responsibilities that become significant in aggregate.

The idea of a ceiling of complexity can be applied to projects as well as to careers. For example, the entropy of a software code base increases over time. Successful projects may have a faster increase in entropy. The software has to maintain backward compatibility because many people depend on its features. Sometimes even its bugs have to be preserved because people rely on them. It’s much easier to renovate software that nobody uses.

Related posts:

The dark side of linchpins
Abandoning projects
A little simplicity goes a long way

One ought to be ashamed to make use of the wonders of science embodied in a radio set, the while appreciating them as little as a cow appreciates the botanic marvels in the plants she munches.

Source: The Science of Radio by Paul Nahin, first edition

All the books in this series are organized by task. For each task, there is a one-line solution followed by detailed commentary. The explanations frequently offer alternate solutions with varying degrees of concision and clarity. Sections are seldom more than one page long, so the books are easy to read a little at a time.

Programmers who have written a lot of Perl may still learn a few things from Krumins. In particular, those who have primarily written Perl in script files may not be familiar with some of the tricks for writing succinct Perl on the command line.

Other Perl posts:

All languages equally complex?
Periodic table of Perl operators
Three-hour-a-week language

… when you study the evidence, it’s clear that you’re not likely to encounter real interesting opportunities in your life until after you’re really good at something.

If you avoid focus because you want to keep your options open, you’re likely accomplishing the opposite. Getting good is a prerequisite to encountering options worth pursuing.

From Closing your interests opens more interesting opportunities.

Related post:

Demonstrating persistence

The function jinc(x) is defined as J₁(x) / x, so if you have code to compute J₁ then it ought to be a no-brainer. For example, why not use the following C code?

#include <math.h>
double jinc(double x) {
    return j1(x) / x;
}

The problem is that if you pass in 0, the code will divide by 0 and return a NaN. The function jinc(x) is defined to be 1/2 at x = 0 because that’s the limit of J_1(x)(x) / x as x goes to 0. So we try again:

#include <math.h>
double jinc(double x) {
    return (x == 0.0) ? 0.5 : j1(x) / x;
}

Does that work? Technically, it could still fail — we’ll come back to that at the end — but we’ll assume for now that it’s OK.

We could write the analogous Python code, and it would be adequate as long as we’re only calling the function with scalars and not NumPy arrays.

from scipy.special import j1
def jinc(x):
    if x == 0.0:
        return 0.5
    return j1(x) / x

Now suppose you want to plot this function. You create an array of points, say

x = np.linspace(-1, 1, 25)

and plot jinc(x). You’ll get a warning: “ValueError: The truth value of an array with one element is ambiguous. Use a.any() or a.all().” Incidentally, if we called linspace with an even integer in the last argument, our array of points would avoid zero and the naive implementation of jinc would work.

When Python tries to apply jinc to an array, it doesn’t know how to interpret the test x == 0. The warning suggests “Do you mean if any component of x is 0? Or if all components of x are 0?” Neither option is what we want. We want to apply jinc as written to each element of x. We could do this by calling the vectorize function.

jinc = np.vectorize(jinc)

This replaces our original jinc function with one that handles NumPy arrays correctly.

There is an extremely unlikely scenario in which the code above could fail. The value of J₁(x) is approximately x/2 for small values of x. If the floating point value x is so small that 0.5*x returns 0, our function will return 0, even though it should return 0.5. The C code above works for values of x as small as DBL_MIN and even values much smaller. (DBL_MIN is not the smallest value of a double, only the smallest normalized double.) But if you set

x = DBL_MIN / pow(2.0, 52);

then jinc(x) will return 0. If you want to be absolutely safe, you could change the implementation to

#include <math.h>
double jinc(double x) {
    return (fabs(x) < 1e-8) ? 0.5 : j1(x) / x;
}

Why test for whether the absolute value is less than 10^-8 rather than a much smaller number? For small x, the error in approximating jinc(x) with 1/2 is on the order of x²/16. So for x as large as 10^-8, the approximation error is below the resolution of a double. As a bonus, the function jinc(x) will be more efficient for |x| < 10^-8 since it avoids a call to j1.

Related posts:

Jinc function
Sine approximation for small angles
Functions in math.h that seem unnecessary

The sinc function is defined by

sinc(t) = sin(t) / t.

The jinc function is defined analogously by

jinc(t) = J₁(t) / t

where J₁ is a Bessel function. Bessel functions are analogous to sines, so the jinc function is analogous to the sinc function.

Here’s what the sinc and jinc functions look like.

The jinc function is not as common as the sinc function. For example, both Mathematica and SciPy have built-in functions for sinc but not for jinc. [There are actually two definitions of sinc. Mathematica uses the definition above, but SciPy uses sin(πt)/πt. The SciPy convention is more common in digital signal processing.]

As I write this, Wikipedia has an entry for sinc but not for jinc. Someone want to write one?

For small t, jinc(t) is approximately cos(t/2) / 2. This approximation has error O(t⁴), so it’s very good for small t, useless for large t.

For large values of t, jinc(t) is like a damped, shifted cosine. Specifically,

with an error that decreases like O( |t|^-2 ).

Like the sinc function, the jinc function has a simple Fourier transform. Both transforms are zero outside the interval [-1, 1]. Inside this interval, the transform of sinc is a constant, √(π/8). On the same interval, the transform of jinc is √(2/π) √(1 – ω²).

Update: How to compute jinc(x)

Related posts:

How to visualize Bessel functions
Diagram of Bessel function relationships

This reminds me of a standard question that Don Rubin … asks in virtually any situation: “What would you do if you had all the data?” For me, that “what would you do” question is one of the universal solvents of statistics.

Emphasis added.

I had not heard Don Rubin’s question before, but I think I’ll be asking it often. It reminds me of Alice’s famous dialog with the Cheshire Cat:

“Would you tell me, please, which way I ought to go from here?”

“That depends a good deal on where you want to get to,” said the Cat.

“I don’t much care where–” said Alice.

“Then it doesn’t matter which way you go,” said the Cat.

Related post: Irrelevant uncertainty

I’ll tell you what I see. I see some kind of vague showy, wiggling lines  — here and there an E and a B written on them somehow, and perhaps some of the lines have arrows on them — an arrow here or there which disappears when I look too closely at it. When I talk about the fields swishing through space, I have a terrible confusion between the symbols I use to describe the objects and the objects themselves. I cannot really make a picture that is even nearly like the true waves. So if you have difficulty making such a picture, you should not be worried that your difficulty is unusual.

From The Feynman Lectures on Physics, volume II.

Other Feynman posts:

Richard Feynman and Captain Picard try to prove Fermat’s Last Theorem
Most important event of the 19th century
God is in the details

Thank you everyone for following and for providing feedback.

Related posts:

How to subscribe to a Twitter post via RSS
Daily tip accounts broader than their names imply

This isn’t much of a paradox; the resolution is quite simple. A graph may look wiggly because the scale is wrong. If the graph is flat, graphing software may automatically choose narrow vertical range, accentuating noise in the graph. I haven’t heard a name for this, though I imagine someone has given it a name.

Here’s an extreme example. The following graph was produced by the Mathematica command Plot[Gamma[x+1] - x Gamma[x], {x, 0, 1}].

This is unsettling the first time you run into it, until you notice the vertical scale. In theory, Γ(x + 1) and x Γ(x) are exactly equal. In practice, a computer returns slightly different values for the two functions for some arguments. The differences are on the order of 10^-15, the limit of floating point precision. Mathematica looks at the range of the function being plotted and picks the default vertical scaling accordingly.

In the example above, the vertical scale is 15 orders of magnitude smaller than the horizontal scale. The line is smooth as glass. Actually, it’s much smoother than glass. An angstrom is only 10 orders of magnitude smaller than a meter, so you wouldn’t have to look at glass under nearly as much magnification before you see individual atoms. At a much grosser scale you’d see imperfections in the glass.

The graph above is so jagged that it demands our attention. When the horizontal axis is closer to the proper scale, say off by a factor of 5 or 10, the problem can be more subtle. Here’s an example that I ran across yesterday.

The curves look awfully jagged, but this is just simulation noise. The function values are probabilities, and when viewed on a scale of probabilities the curves look unremarkable.

Now listen to the rule of the last inch. The realm of the last inch. The job is almost finished, the goal almost attained, everything possible seems to have been achieved, every difficulty overcome — and yet the quality is just not there. The work needs more finish, perhaps further research. In that moment of weariness and self-satisfaction, the temptation is greatest to give up, not to strive for the peak of quality. That’s the realm of the last inch — here, the work is very, very complex, but it’s also particularly valuable because it’s done with the most perfect means. The rule of the last inch is simply this — not to leave it undone. And not to put it off — because otherwise your mind loses touch with that realm. And not to mind how much time you spend on it, because the aim is not to finish the job quickly, but to reach perfection.

Via Still I Am One

It can be hard to know when something deserves the kind of polish Solzhenitsyn talks about. Sometimes you’re in the realm of rapidly diminishing return and it’s time to move on. Other times, the finishing touches are everything.

Related post:

Scripting and the last mile problem

The most expensive function in the simulation has four integer arguments. Each time the function is called, one of the arguments is incremented by 1. A random process determines which argument is incremented each time. This function has a set of recurrence relations that make it faster to calculate a new value if you have a previously computed value for adjacent arguments.

The recurrence relation approach is about 100x faster than the memoization approach for small runs. However, the time per run with memoization decreases with each run. The longer it runs, the more values it caches, and so the greater the chances that a value has already been computed. In the end, the recurrence relation is about 2x faster than the memoization approach. With an even larger run, I suspect the time for the two approaches would be roughly equal.

The new approach is faster, but it was more work to implement, and it depends critically on the assumption that consecutive arguments differ only by one increment of one argument.

Related post:

20 weeks down to 20 minutes

The most expensive part of the simulation is a function that takes in four integer parameters and does a numerical integration. Each parameter could be any value from 0 to 120. So in principle there are 121^4 = 214,358,881 possible parameter combinations. I knew that some of these combinations were extremely unlikely or impossible, so I thought I could reduce the execution time by caching the most commonly called parameters. Maybe a few thousand arguments would account for most of the function calls.

I saved the function arguments and computed values in map, a.k.a. dictionary, associative array, etc. (The buzz word for this is memoization.) I set a maximum size on the map because I thought it could grow very large. Even if, say, a thousand arguments accounted for most of the function calls, maybe there were many more that only appeared once.

It turned out that the function was only called with 408 distinct parameter combinations. That means the same integrals were being evaluated millions of times. Caching these values reduced the execution time by a factor of about 800. In other words, the software could retrieve a previously computed value 800x faster than compute it from scratch.

The simulation that would have taken 300 hours took only 22 minutes, and so there was no need to make it run in parallel.

Update: I found a bug in my hashing algorithm. After fixing the bug, the time to run my simulations increased to 38 minutes. So my speed-up was a about a factor of 420 rather than 800, but still enough to run on one thread over lunch. The number of distinct arguments to the numerical integration was much more than 408 in the corrected program, but still small enough to cache.

Many man pages are hard to read, but I think that the grand prize for difficulty has to go to the man page for bash. As I was doing my research for this book, I gave it a careful review to ensure that I was covering most of its topics. When printed, it’s over 80 pages long and extremely dense, and its structure makes absolutely no sense to a new user.

On the other hand, it is very accurate and concise, as well as being extremely complete. So check it out if you dare, and look forward to the day when you can read it and it all makes sense.

* If you’re not familiar with Unix lingo, “man” stands for “manual”. Man pages are online documentation.

Related post:

Review of The Linux Command Line

True to its name, the book is about using Linux from command line. It’s not an encyclopedia of Linux. It doesn’t explain how to install Linux, doesn’t go into system APIs, and says little about how to administer Linux. At the same time, the book is broader than just a book on bash. It’s about how to “live” at the command line.

The introduction explains the intended audience.

This book is for Linux users who have migrated from other platforms. Most likely you are a “power user” of some version of Microsoft Windows.

The book has a conversational style, explaining the motivation behind ways of working as well as providing technical detail. It includes small but very useful suggestions along the way, the kinds of tips you’d pick up from a friend but might not find in a book.

The book has four parts

1. Learning the shell
2. Configuration and the environment
3. Common tasks and essential tools
4. Writing shell scripts

The book could have just included the first three sections; the forth part is a bit more specialized than the others. If you’d prefer, think of the book has having three parts, plus a lengthy appendix on shell scripting.

The Linux Command Line is pleasant to read. It has a light tone, while also getting down to business.

Related post:

Perverse hipster desire for retro-computing

From TED

With ordinary differential equations (ODEs), the initial conditions are often an afterthought. First you find a full set of solutions, then you plug in initial conditions to get a specific solution.

Partial differential equations (PDEs) have boundary conditions (and maybe initial conditions too). Since people typically learn ODEs first, they come to PDEs expecting boundary values to play a role analogous to ODEs. In a very limited sense they do, but in general boundary values are quite different.

The hard part about PDEs is not the PDEs themselves; the hard part is the boundary conditions. Finding solutions to differential equations in the interior of a domain is easy compared to making the equations have the specified behavior on the boundary.

No model can take everything into account. You have to draw some box around that part of the world that you’re going to model and specify what happens when your imaginary box meets the rest of the universe. That’s the hard part.

Related posts:

Three views of differential equations
Nonlinear is not a hypothesis
Approximating a solution that does not exist

Some people are asking why go to college when you can download a college education via iTunes U.

First, you would have no credentials when you’re done.

Second, you almost certainly would not put in the same amount of work as a college student without someone to pace you through the material and to provide external motivation. You’d be less likely to struggle through anything you found difficult or uninteresting.

Third, you’d have no feedback to know whether you’re really learning what you think you’re learning.

The people that I hear gush about online education opportunities are well-educated, successful, and ambitious. They may be less concerned about credentials either because they are intrinsically motivated or because they already have enough credentials. And because of their ambition, they need less accountability. They may need less feedback or are resourceful enough to seek out alternative channels for feedback, such as online forums. Resources such iTunes U and The Teaching Company are a godsend to such people. But that doesn’t mean that a typical teenager would make as much of the same opportunities.

I’m starting by reading the Feynman lectures on E&M. After that I plan to read something on electronics. If you have resources you recommend, please let me know.

I’ve started new Twitter account, @GrokEM. I figure that tweeting about E&M will help me stick to my goal. My other Twitter accounts post on a regular schedule (plus a few extras) and are scheduled weeks in advance. GrokEM will be more erratic, at least for now. (In case you’re not familiar with grok, it’s a slang for knowing something thoroughly and intuitively.)

Here’s what Feynman said about mathematicians learning physics, particularly E&M.

Mathematicians, or people who have very mathematical minds, are often led astray when “studying” physics because they loose sight of the physics. They say: “Look, these differential equations — the Maxwell equations — are all there is to electrodynamics … if I understand them mathematically inside out, I will understand the physics inside out.” Only it doesn’t work that way. … They fail because the actual physical situations in the real world are so complicated that it is necessary to have a much broader understanding of the equations. … A physical understanding is a completely unmathematical, imprecise, and inexact thing, but absolutely necessary for a physicist.

Heinlein coined grok around the same time that Feynman made the above remarks. Otherwise, Feynman might have said that only studying differential equations is not the way to grok electrodynamics.

Related posts:

What it means to understand an equation
Most important event of the 19th century

Rather than write a long post, I’ll restrain myself and simply say that I’d like to hear more people talk about “educational monoculture.”

Related post:

Don’t standardize education, personalize it

We can prove that the ratio oscillates by starting with the formula

where φ = (1 + √5)/2 is the golden ratio.

From there we can work out that

This shows that when n is odd, F(n+1) / F(n) is below the golden ratio and when n is even it is above. It also shows that the absolute error in approximating the golden ratio by F(n+1) / F(n) goes down by a factor of about φ² each time n increases by 1.

Related posts:

Honeybee genealogy
Fibonacci numbers at work
Breastfeeding and the golden ratio

Beginning algebra students rush to enter numbers into a calculator. They find it comforting to reduce expressions to floating point numbers as frequently as possible. This is understandable, but it’s a habit they need to break for numerous reasons: it’s more work, harder for a grader to follow, less accurate, etc. Better to do as much calculation as possible with symbols, then stick in numbers at the end.

A similar phenomena happens in scientific programming. We’re anxious to see numbers, so we print out numbers as soon as we produce them. There’s a tendency to sprinkle code with printf statements, then write scripts to scrape text output to gather results.

It’s better to keep numbers in data structures as long as possible, then dump the data to text as the last step. Why? For one thing, the output format might change: instead of a text dump, maybe you’ll want to put the data in a database or format it as an HTML table. More importantly, the “last step” will often change. What was the last step now becomes the next-to-last step because you think of something else to do. So you do more with the data in memory before producing output, rather than scraping the text output.

I quickly learned to delay plugging numbers into algebraic expressions. It took me longer to learn not to output numeric results until the last moment.

Mr. McCabe thinks that I am not serious but only funny, because Mr. McCabe thinks that funny is the opposite of serious. Funny is the opposite of not funny, and of nothing else. … Whether a man chooses to tell the truth in long sentences or short jokes is a problem analogous to whether he chooses to tell the truth in French or German. Whether a man preaches his gospel grotesquely or gravely is merely like the question of whether he preaches it in prose or verse. … The truth is, as I have said, that in this sense the two qualities of fun and seriousness have nothing whatever to do with each other, they are no more comparable than black and triangular.

Emphasis added. From Heretics. Text available online from Project Gutenberg.

Given their likely acquaintance such a cognitively rich world of work, it is hardly surprising that when Henry Ford introduced the assembly line in 1913, workers simply walked out. One of Ford’s biographers wrote, “So great was labor’s distaste for the new machine system that toward the close of 1913 every time the company wanted to add 100 men to its factory personnel, it was necessary to hire 963.”

A dozen years ago people would talk of building “software factories” to crank out software projects. Back then someone tried to get me excited about joining an effort to create such a factory. I told him I did not want to work in a factory. He tried to back-peddle, saying that it’s not what it sounds like. But I’m sure it was exactly what it sounded like.