Python

Someone sent me email regarding my online calculator for computing the distance between to locations given their longitude and latitude values. He wants to so sort of the opposite. Starting with the longitude and latitude of one location, he wants to find the longitude and latitude of locations moving north/south or east/west of that location. I like to answer reader questions when I can, so here goes. I’ll give a theoretical derivation followed by some Python code.

Longitude and latitude and are usually measured in degrees, but theoretical calculations are cleaner in radians.  Someone using the Python code below could think in terms of degrees; radians will only be used inside function implementations. We’ll use the fact that on a circle of radius r, an arc of angle θ radians has length rθ. We’ll assume the earth is a perfect sphere. See this post for a discussion of how close the earth is to being a sphere.

Moving North/South

I’ll start with moving north/south since that’s simpler. Let R be the radius of the earth. An arc of angle φ radians on the surface of the earth has length M = Rφ, so an arc M miles long corresponds to an angle of φ = M/R radians. Moving due north or due south does not change longitude.

Moving East/West

Moving east/west is a little more complicated. At the equator, the calculation is just like the calculation above, except that longitude changes rather than latitude. But the distance corresponding to one degree of longitude changes with latitude. For example, one degree of longitude along the Arctic Circle doesn’t take you nearly as far as it does at the equator.

Suppose you’re at latitude φ degrees north of the equator. The circumference of a circle at constant latitude φ, a circle parallel to the equator, is cos φ times smaller than the circumference of the equator. So at latitude φ an angle of θ radians describes an arc of length M = R θ cos φ. A distance M miles east or west corresponds to a change in longitude of θ = M/(R cos φ). Moving due east or due west does not change latitude.

Python code

The derivation above works with angles in radians. Python’s cosine function also works in radians. But longitude and latitude are usually expressed in degrees, so function inputs and outputs are in degrees.

import math

# Distances are measured in miles.
# Longitudes and latitudes are measured in degrees.
# Earth is assumed to be perfectly spherical.

earth_radius = 3960.0
degrees_to_radians = math.pi/180.0
radians_to_degrees = 180.0/math.pi

def change_in_latitude(miles):
    "Given a distance north, return the change in latitude."
    return (miles/earth_radius)*radians_to_degrees

def change_in_longitude(latitude, miles):
    "Given a latitude and a distance west, return the change in longitude."
    # Find the radius of a circle around the earth at given latitude.
    r = earth_radius*math.cos(latitude*degrees_to_radians)
    return (miles/r)*radians_to_degrees

Related posts:

What is the shape of the earth?
Spherical trigonometry

According to an article on The Snowtide Blog, MIT has decided to teach beginning CS classes in Python rather than Scheme. (Thanks to procoders.net for the link.) The article paraphrases remarks by Gerdald Sussman at the announcement.The main rationale was that the world is more complicated now than when the course was designed and SICP was written. Now programmers spend more time researching libraries than writing everything from scratch.

Here’s something I expected Sussman to say though apparently he did not: you can use Python as a functional language if you like, you just don’t have to. I don’t know why more people don’t emphasize this. Python is not a functional language, but Python does make it easy to declare functions on the fly, pass functions as arguments, etc. You’re free to write Scheme-like programs in Python if you want to, but you’re also free to use other styles when they are more appropriate.

In a previous post, I discuss my difficulties calling some Python modules from IronPython. In particular I wanted to call SciPy from IronPython and couldn’t. The discussion following that post brought up Ironclad as a possible solution. I wanted to learn more about Ironclad, and so I invited William Reade to write a guest post about the project. I want to thank William for responding to my request with a  very helpful article. — John

Hi! My name’s William Reade, and I’ve spent the last year or so working on Ironclad, an open-source project which helps IronPython to inter-operate better with CPython. Michael Foord recently introduced  me to our host John, who kindly offered me the opportunity to write a bit about  my work and, er, how well it works. So, here I am.

To give you a little bit of context, I’ve been working at Resolver Systems for several years now; our main product, Resolver  One, is a spreadsheet with very tight IronPython integration. We like to describe  it as a “Pythonic spreadsheet”, and that’s clearly a concept that people like.  However, when people think of a “Pythonic spreadsheet”, they apparently expect it  to work with popular Python libraries — such as NumPy and SciPy — and we found that IronPython’s incompatibility put us at a serious disadvantage. And, for some reason, nobody seemed very keen to  solve the problem for us, so we had to do it ourselves.

The purpose of Ironclad is to allow you to use Python C extensions (of which there are many) from inside IronPython without recompiling anything. The secret purpose  has always been to get NumPy working in Resolver One, and in release 1.4 we finally  achieved this goal. Although the integration is still alpha level, you can import  and use NumPy inside the spreadsheet grid and user code: you can see a screencast  about the integration here.

However, while Resolver One is a great tool, you aren’t required to use it to get the benefits: Ironclad has been developed completely separately, has no external  dependencies, and is available under an open source license. If you consider  yourself adequately teased, keep reading for a discussion of what Ironclad actually  does, what it enables you to do, and where it’s headed.

[click to continue...]

There’s a new math blog, Numerical Recipes. The initial post says the blog will be understandable to anyone with a high school math background, will solve all problems in Python, and will not use any external libraries. The blog will also be available in Spanish.

It’s difficult to use SciPy from IronPython because much of SciPy is implemented in C rather than in Python itself. I wrote an article on CodeProject summarizing some things I’d learned about using Python modules with IronPython. (Many thanks to folks who left helpful comments here and answered my questions on StackOverflow.) The article gives stand-alone code for computing normal probabilities and explains why I wrote my own code rather than using SciPy.

Here’s the article: Computing Normal Probabilities in IronPython

Here’s some more stand-alone Python code for computing mathematical functions and generating random numbers. The code should work in Python 2.5, Python 3.0, and IronPython.

Related post: IronPython is a one-way gate

I often find it strange how some programming books arbitrarily classify content as either “beginner” or “advanced.” The advanced material may not be advanced at all. Maybe by “advanced” the author means “the stuff I didn’t learn right away.” However, Expert Python Programming by Tarek Ziadé objectively is an advanced book. I had heard good things about the book from the Python 411 podcast and had it on my list of books to buy when the book’s publisher gave me a copy to review.

Expert Python Programming assumes the reader has a solid knowledge of Python. In some ways, it is the Python counterpart to Scott Meyers’ Effective C++ series, giving advice about recommended practice. But Ziadé’s Python book is broader than Meyer’s C++ books because Expert Python Programming goes beyond best practices for language use. The latter chapters discuss how to package and distribute Python applications, software life cycle management, documentation, test-driven development, optimization, and design patterns. I’d like to see more books follow this pattern, covering these important topics in the context of a particular language or tool set. I was particularly pleased to see 27 pages devoted to documentation.

True to its name, Expert Python is not for beginners; I’d recommend Wesley Chun’s Core Python for a first book. But I’d heartily recommend Expert Python Programming as a second Python book. The book has an impressive range of practical material including numerous links to even more resources.

Update: Thomas Guest also reviewed this book.

Sergey Fomel left a valuable comment on my post about computing the error function erf(x). He gave some sample code for exposing C functions to Python via SWIG on Linux. I haven’t used SWIG, but I’ve heard good things about it. I believe Google uses SWIG extensively to make C++ code callable from Python.

Here’s Sergey’s code.

bash$ cat erf.i

%module erf
#include
double erf(double);

bash$ swig -o erf_wrap.c -python erf.i
bash$ gcc -o erf_wrap.os -c -fPIC -I/usr/include/python2.4 erf_wrap.c
bash$ gcc -o _erf.so -shared erf_wrap.os
bash$ python
>>> from erf import erf
>>> erf(1)
0.84270079294971489

I’ve seen several people ask lately how to compute the distribution (CDF) function for a standard normal random variable, often denoted Φ(x). They want to know how to compute it in Java, or Python, or C++, etc. Every language has its own standard libraries, and in general I recommend using standard libraries. However, sometimes you want to minimize dependencies. Or maybe you want more transparency than your library allows. The code given here is in Python, but it is so compact that it could easily be ported to any other language.

I just posted Python code for computing the error function, erf(x). The normal density Φ(x) is a simple transformation of erf(x). Given code for erf(x), here’s code for Φ(x).

def phi(x):
    return 0.5*( 1.0 + erf(x/math.sqrt(2)) )

After deriving the transformations between erf(x) and Φ(x) several times, including their complements and inverses, I wrote them down to save. See the PDF file Relating Φ and erf.

See also stand alone code for computing the inverse of the standard normal CDF.

The question came up on StackOverflow this morning how to compute the error function erf(x) in Python. The standard answer for how to compute anything numerical in Python is “Look in SciPy.” However, this person didn’t want to take on the dependence on SciPy. I’ve seen variations on this question come up in several different contexts lately, including questions about computing the normal distribution function, so I thought I’d write up a solution.

Here’s a Python implementation of erf(x) based on formula 7.1.26 from A&S. The maximum error is below 1.5 × 10^-7.

import math

def erf(x):
    # constants
    a1 =  0.254829592
    a2 = -0.284496736
    a3 =  1.421413741
    a4 = -1.453152027
    a5 =  1.061405429
    p  =  0.3275911

    # Save the sign of x
    sign = 1
    if x < 0:
        sign = -1
    x = abs(x)

    # A & S 7.1.26
    t = 1.0/(1.0 + p*x)
    y = 1.0 - (((((a5*t + a4)*t) + a3)*t + a2)*t + a1)*t*math.exp(-x*x)

    return sign*y

This problem is typical in two ways: A&S has a solution, and you’ve got to know a little background before you can use it.

The formula given in A&S is only good for x ≥ 0. That’s no problem if you know that the error function is an odd function, i.e. erf(-x) = -erf(x). But if you’re an engineer who has never heard of the error function but needs to use it, it may take a while to figure out how to handle negative inputs.

One other thing that someone just picking up A&S might not know is the best way to evaluate polynomials. The formula appears as 1 - (a₁t¹ + a₂t² + a₃t³ + a₄t⁴ + a₅t⁵)exp(-x²), which is absolutely correct. But directly evaluating an nth order polynomial takes O(n²) operations, while the factorization used in the code above uses O(n) operations. This technique is known as Horner’s method.

My previous post showed how it is possible to do partial function application in C++ by using function objects. Here I’ll show how much simpler this can be done in Python.

As before, we want to evaluate a function of one variable and three parameters: f(x; a, b, c) = 1000a + 100b + 10c + x. Here’s an example how we could do this in Python using the functools module.

from functools import partial

def f(a, b, c, x):
    return 1000*a + 100*b + 10*c + x

g = partial(f, 3, 1, 4)

print g(5)

The code will print 3145.

The function f above is so simple that it may be hard to imagine why you would want to do such a thing. The earlier post gives a more realistic application, using partial function application in numerical integration and root-finding.

Fonts, translations, and programming languages have one thing in common: they work best when you don’t notice them.

If someone says “Hey, look at this cool font I just found!” you probably wouldn’t want to read a book set in that font. At least to an untrained eye, a great font will not stand out in a list of small samples. You have to see large blocks of text set in a font to appreciate it. Even then, most people will not consciously appreciate a very readable font.

Translations are similar. If you find yourself saying “What an interesting translation!” then the translator has probably fallen down on the job. A good translation is neither archaic nor trendy. It does not draw attention to itself but allows you to focus on the original content. I believe the English Standard Version achieves that with Bible translation.

Python is like a good font or a good translation. For years I’d look into Python briefly when someone would recommend it. I’d thumb through a Python book, but it all looked rather plain. Only later did I come to appreciate that the beauty of Python is that it is rather plain. It doesn’t call attention to itself. It just gets out of your way and lets you write programs. It seems to me that compared to other programming language communities, the Python community brags less about their language per se and more about what they’re able to do with it.

I needed to find a spell checker I could call from Python, so I did a Google search and ran across GNU aspell. I tried installing it but got contradictory warning messages: aspell not installed, aspell already installed, etc. Then I remembered what an awful time I’d had before when I’d tried to use aspell and gave up.

Next I tried Ryan Kelly’s PyEnchant and it worked like a charm. I downloaded the installer for Windows and ran it. Then I opened up a Python console and typed an example following the online tutorial.

>>> import enchant
>>> d = enchant.Dict("en_US")
>>> d.check("Potatoe")
False
>>> d.check("Potato")
True

It just works.

Sometimes code is easier to understand than abstract math. A few days ago I was having a hard time conveying bias, consistency, and efficiency in a statistics class. Writing some pseudo-code on the board seemed to help clear things up. Loops and calls to random number generation routines are more tangible than discussions about random samples.

Later I turned the pseudo-code into Python code (after all, Python is supposed to be “executable pseudo-code”) and fancied it up a bit. The following page gives some explanation, some plots of the output, and the source code.

The difference between an unbiased estimator and a consistent estimator

Symmetric APIs are easier to use. I was reminded of this when doing some regular expression programming in Python and comparing it to Perl. Perl’s regular expression operators for search and replace are symmetric in a way that their Python counterparts are not.

Perl uses m/pattern/ for matching and s/pattern/replacement/ for substitution. Both apply to the first instance of a pattern in a string by default. The g option following a match or substitute operator causes the command to apply to all instances of the pattern. The i option after either a match or substitute command causes the pattern to apply in a case-insensitive manner. Matching and substitution are symmetric.

Python uses re.search() for matching and re.sub() for substitution. The search function can only apply to the first instance of a pattern; to match all instances of a pattern, use re.findall(). The function re.sub() applies to all instances by default, but it has a max parameter that can be set to limit the number of instances it applies to. To make a search pattern case-insensitive, pass in re.IGNORECASE flag. To make a substitution case-insensitive, modify the regular expression itself by adding (?i).

In general, I find Python syntax much cleaner than Perl, but regular expressions are implemented more elegantly in Perl.