pkg-config

A reasonably brief post about pkg-config, which we used to help build matplotlib the other day (here).

The web page is here and there is a guide here.

As it says on the first page:

pkg-config is a helper tool used when compiling applications and libraries. It helps you insert the correct compiler options on the command line so an application can use gcc -o test test.c `pkg-config --libs --cflags glib-2.0` for instance, rather than hard-coding values on where to find glib (or other libraries). It is language-agnostic, so it can be used for defining the location of documentation tools, for instance.

And from the man page:

DESCRIPTION
The pkg-config program is used to retrieve
information about installed libraries in the
system. It is typically used to compile and
link against one or more libraries. Here is a
typical usage scenario in a Makefile:

program: program.c
cc program.c $(pkg-config --cflags --libs gnomeui)

pkg-config retrieves information about pack-
ages from special metadata files. These files
are named after the package, and has a .pc
extension. On most systems, pkg-config looks
in and
for these files. It will additionally look
in the colon-separated (on Windows, semicolon-
separated) list of directories specified by
the PKG_CONFIG_PATH environment variable.

The package name specified on the pkg-config
command line is defined to be the name of the
metadata file, minus the .pc extension. If a
library can install multiple versions simulta-
neously, it must give each version its own
name (for example, GTK 1.2 might have the
package name "gtk+" while GTK 2.0 has
"gtk+-2.0").

What's happening is that the matplotlib setup.py script uses pkg-config if it's present. The binary is in /usr/local/bin/pkg-config. It looks for a subdirectory of /usr/X11/lib (and /usr/local/lib) called pkgconfig that constains .pc files for packages that specify metadata about themselves to take advantage of pkgconfig functionality:

> cat /usr/X11/lib/pkgconfig/libpng.pc prefix=/usr/X11 exec_prefix=${prefix} libdir=${exec_prefix}/lib includedir=${prefix}/include/libpng12 Name: libpng Description: Loads and saves PNG files Version: 1.2.44 Libs: -L${libdir} -lpng12 Libs.private: -lz Cflags: -I${includedir}

---

We used Homebrew only to grab pkg-config. You could get it yourself: 0.26 is the newest, but requires glib-2.0; 0.25 bundles glib with pkg-config. 0.25 builds with no dependencies using the standard ./configure; make; make install. Download page for version 0.25 and 0.26 is here. Or of course, you can use Homebrew like we did.

Not clear why it doesn't ship with OS X. The license is GPLv2. But then, so is BASH.

Some minor issues:

• when I tried using pkg-config on my iMac (Snow Leopard, littered with MacPorts and other stuff), matplotlib seemed to build and install but died with:

> python x.py Traceback (most recent call last): File "x.py", line 1, in import matplotlib.pyplot as plt File "/Library/Python/2.6/site-packages/matplotlib/pyplot.py", line 95, in new_figure_manager, draw_if_interactive, _show = pylab_setup() File "/Library/Python/2.6/site-packages/matplotlib/backends/__init__.py", line 25, in pylab_setup globals(),locals(),[backend_name]) File "/Library/Python/2.6/site-packages/matplotlib/backends/backend_macosx.py", line 243, in class TimerMac(_macosx.Timer, TimerBase): AttributeError: 'module' object has no attribute 'Timer'

For reasons that I don't understand, the version of the file in the build directory

/Users/telliott/Software/matplotlib-1.0.1/build/lib.macosx-10.6-universal-2.6/matplotlib/backends/backend_macosx.py

is not the same as the one in site-packages

/Library/Python/2.6/site-packages/matplotlib/backends/backend_macosx.py

Simply replacing the site-packages version fixed the problem:

> sudo cp /Users/telliott/Software/matplotlib-1.0.1/build/lib.macosx-10.6-universal-2.6/matplotlib/backends/backend_macosx.py /Library/Python/2.6/site-packages/matplotlib/backends/backend_macosx.py > python -c "import matplotlib.pyplot as plt" >

• the Homebrew "Formula" in /usr/local/Library/Formula/pkg-config.rb tells the installed pkg-config which paths to search for lib.pc files:

/usr/local/lib/pkgconfig /usr/lib/pkgconfig /usr/X11/lib/pkgconfig

I'm not quite sure how that would work for a manual install.

• the libpng in /usr/local/lib that I installed shadowed the one in /usr/X11/lib so to use the X11 version, I removed the files from /usr/local

• The notation libpng.pcfileLibs.private: -lz indicates that libpng does depend on zlib, but it's already been linked against it (hence private). I'm not sure where it is, though.

matplotlib on OS X Lion--revised

In the previous post, I repeated an installation of matplotlib (now for OS X Lion) following the same method that I've linked to a number of times on the blog. It worked this time as well, but lead me to consideration of other possible methods. I found one on the web here. So that's what this post is about.

BTW, whatever you do, do not follow the instructions that the matplotlib developers provide. You do not need another Python, including MacPython or the Endthought distribution or anything else..

To recap, matplotlib lists dependencies (in the make.osx file) of: libpng, libfreetype, and zlib. I read somewhere on the matplotlib site today (but can't find the link now) that zlib is not a required dependency. These days, the other two are actually provided by Apple:

> ls -al /usr/X11/lib/libpng* -rwxr-xr-x 1 root wheel 296864 Jul 29 16:11 /usr/X11/lib/libpng.3.dylib lrwxr-xr-x 1 root wheel 14 Jul 29 16:11 /usr/X11/lib/libpng.dylib -> libpng15.dylib -rwxr-xr-x 1 root wheel 294160 Jul 29 16:11 /usr/X11/lib/libpng12.0.dylib -rwxr-xr-x 1 root wheel 305008 Jul 29 16:11 /usr/X11/lib/libpng15.15.dylib lrwxr-xr-x 1 root wheel 17 Jul 29 16:11 /usr/X11/lib/libpng15.dylib -> libpng15.15.dylib > ls -al /usr/X11/lib/libfreetype* -rwxr-xr-x 1 root wheel 1060416 Jul 29 16:11 /usr/X11/lib/libfreetype.6.dylib lrwxr-xr-x 1 root wheel 19 Jul 29 16:11 /usr/X11/lib/libfreetype.dylib -> libfreetype.6.dylib

I'm not sure at the moment whether X11 came with Xcode (as it used to) or was present in the original Lion install.

In any case, I spent an hour or two trying to figure out how to use the build commands from the makefile that comes with matplotlib, but point at these libraries. In the process I found what appears to be a bug:

/usr/X11/include/png.h:666: error: forward declaration of ‘struct png_info_def’

(and see Prashant's answer here), but we're going to use a different approach so it doesn't matter. I'm not sure how our method solved this in the end. As mentioned, the approach is to use Homebrew, in particular, something called pkgconfig. I used a Ruby script that downloads and installs Homebrew, as described here. The install gave this warning:

Warning: The following *evil* dylibs exist in /usr/local/lib They may break builds or worse. You should consider deleting them: /usr/local/lib/libfreetype.6.dylib /usr/local/lib/libpng.3.dylib /usr/local/lib/libpng12.0.dylib /usr/local/lib/libz.1.2.5.dylib

These are, of course, the libraries we just installed in order to get matplotlib to build.

That's what got me working on the other method more seriously. I don't think the danger is all that great, but clearly it is better to use libraries that are (i) still available from the maintainers and (ii) have been vetted at least to some extent by other folks including Apple. AFAIK there aren't any known security issues with the versions we installed previously. I ripped them out anyway (they seemed to confuse pkgconfig).

If we ask Homebrew for the libraries, it just points us to the ones we saw in /usr/X11:

> brew search libpng Apple distributes libpng with OS X, you can find it in /usr/X11/lib. However not all build scripts look here, so you may need to call ENV.x11 in your formula's install function. > brew search freetype Apple distributes freetype with OS X, you can find it in /usr/X11/lib. However not all build scripts look here, so you may need to call ENV.x11 in your formula's install function.

Following the instructions in the blog post I first installed pkgconfig:

> sudo brew install pkgconfig Cowardly refusing to `sudo brew' > brew install pkgconfig

!!

That's all the Homebrew we need. This is followed by:

cd ~/Desktop git clone git://github.com/matplotlib/matplotlib.git cd matplotlib/ python setup.py build sudo python setup.py install

And we can see that I really did overwrite the first matplotlib install, and that we're actually using the libraries from /usr/X11, by first doing

export DYLD_PRINT_LIBRARIES=1

Then when we run a script that imports matplotlib.pyplot, the Terminal shows (among much else) this:

dyld: loaded: /usr/X11/lib/libfreetype.6.dylib .. dyld: loaded: /usr/X11/lib/libpng15.15.dylib

So that's what I'd recommend and it seems to be working fine. This simple script works exactly as you'd expect.

import matplotlib.pyplot as plt Y = [1,4,9,16] plt.scatter(range(len(Y)),Y,s=250,color='r') plt.savefig('example.png')

So the next thing to do is to figure out pkgconfig and Homebrew work their magic!

P.S. You will still need to make and edit ~/.matplotlibas discussed last time.

matplotlib on OS X Lion--old

[ UPDATE: I'm going to leave this post up, for the record, but I found a better way, and that's in the following post. ]

Just a note to say that I got matplotlib installed with only a few issues to solve. As before, I relied on Gavin Huttley's instructions (here). Although the matplotlib instructions say you should install a different Python, I've found that just leads to confusion down the road, and I want to try to use the system version only if I can get away with it. YMMV.

The sticking points for the install always have to do with the basic dependencies, which are listed at the top of the file make.osx in the matplotlib source (from here). The source is the last link at the bottom.

ZLIBVERSION=1.2.3 PNGVERSION=1.2.39 FREETYPEVERSION=2.3.11

As silly as it seems, although OS X provides excellent functionality for dealing with compression, png images and fonts, we need to install these libraries for matplotlib to use. Unfortunately, I'm not expert enough to alter matplotlib to get rid of these dependencies, so we're stuck with it.

---

libpng

I got libpng-1.2.39 from here. I used this version b/c that's what matplotlib asks for, despite the fact that the current release is 1.5.4 (here). Standard magic incantations:

./configure make sudo make install

---

freetype

curl -O http://ftp.twaren.net/Unix/NonGNU/freetype/freetype-2.3.7.tar.bz2

According to Gavin's notes, we need these compiler flags for the other two (but they interfere with the libpng build):

export MACOSX_DEPLOYMENT_TARGET=10.7 export CFLAGS="-arch i386 -arch x86_64" export FFLAGS="-arch i386 -arch x86_64" ./configure make sudo make install

---

zlib

zlib is a compression library. The current version is 1.2.5. I tried this:

curl -O http://www.zlib.net/zlib-1.2.3.tar.gz

but as I've seen before, although you get a (small) file back that has the right icon, it is really a 404 Page Not Found, which doesn't do me any good. It seems the older versions have been disappeared off the web and I couldn't find them poking around on the site. What I should've done at this point is try the version they do offer, but instead I went to my other machine and grabbed 1.2.3. I guess that since I'd used it before to build zlib, I found I needed to clean it first:

make clean make sudo make install

---

matplotlib

Finally, I got the matplotlib source. I altered the file make.osx:

PREFIX=/usr/local PYVERSION=2.7 PYTHON=python${PYVERSION} ZLIBVERSION=1.2.3 PNGVERSION=1.2.39 FREETYPEVERSION=2.3.7 MACOSX_DEPLOYMENT_TARGET=10.7 OSX_SDK_VER=10.7 ARCH_FLAGS="-arch i386-arch x86_64"

make -f make.osx mpl_build sudo python setup.py install

mkdir ~/.matplotlib cp matplotlibrc.template ~/.matplotlib/matplotlibrc

Edit the above file to have:

backend : MacOSX

> python Python 2.7.1 (r271:86832, Jun 16 2011, 16:59:05) [GCC 4.2.1 (Based on Apple Inc. build 5658) (LLVM build 2335.15.00)] on darwin Type "help", "copyright", "credits" or "license" for more information. >>> from matplotlib import * >>>

Yes!

[ UPDATE: I see I used freetype 2.3.7 rather than 2.3.11, by mistake.

And if you don't like my approach, you could use Homebrew and follow this post. I don't have any experience with this yet. ]

OS X Lion

OS X Lion has been released for a week now, so I thought I'd try it out on one of my machines. A definitive review is here. I'll be too busy to really play with it for a while, or to do much programming, but I thought I should at least try. An additional motivation is that my laptop (an aluminum MacBook from about December 08), had developed an issue. It completely drains the battery (starting from full charge) in sleep mode, even when no apps are running. The battery hasn't been particularly stressed:

Health Information:
Cycle Count: 568
Condition: Normal
Battery Installed: Yes
Amperage (mA): 0
Voltage (mV): 12546

And, it doesn't drain when powered off. This problem appeared recently, and I am interested to see if installing Lion might fix it, but that hasn't been tested yet. I'll update when I have those results.

This post is just to document that I got Lion onto a USB flash drive and used it to do the install and it worked fine. Instructions from here. I did this because it's the only way to do a clean install. One odd thing: if the machine on which you are attempting to do the download can't run Lion, you won't be allowed to do the download at all! I had to bring my laptop to work (and get the IT guys to reauthorize it for our network).

Also, I have two Apple IDs---an old one and one that goes with my MobileMe account. When I set up Mail, the ID for the computer was automatically set to be the latter, but I originally did the Lion download with the former. We'll see if and when that causes a problem.

It seems clear that the Mac OS will (in time) become more iOS-like (locked down, user not allowed to risk hurting himself, even if he wants to) and less easy to hack around with. But for now, it's still a win over Linux, for me.

Python is 2.7:

> python Python 2.7.1 (r271:86832, Jun 16 2011, 16:59:05) [GCC 4.2.1 (Based on Apple Inc. build 5658) (LLVM build 2335.15.00)] on darwin Type "help", "copyright", "credits" or "license" for more information.

The default scrolling is like that on the iPhone---as if you're moving the content, rather than a scrollbar. After 2 or 3 hours, I changed the setting back to the old method. Maybe I will try to retrain my brain at a later date.

I also downloaded Xcode 4. It quit in the middle for no reason, but completed on the second try. Weirdly, in the middle of the install I got this alert panel:



I say "weirdly" because iTunes is not running. Starting and quitting iTunes had no effect, either. Finally I went in Terminal and looked for an iTunes-related process:

163 ?? 0:00.03 /Applications/iTunes.app/Contents/MacOS/iTunesHelper.app/Contents/MacOS/iTunesHelper -psn_0_57358

after I killed it, the install finished fine. There is a mention in the big Ars Technica review that apps don't actually quit (or perhaps not always), and may not show the little dot for an active app, but there was no process for iTunes so I don't think so. It's weird that a bug like that should still be present at this stage of the game.

And I like the new looks in Terminal:



[ UPDATE: Yep, the battery does not drain on sleep anymore. (Well, 2% in 2 hr). That's the good news. The bad news is that the extension I had to help autofill passwords when forms contain autocomplete="off" doesn't work with Safari any more (post here, extension here). I sure hope Apple didn't do that on purpose. ]

UPDATE 2: I just needed to check to the right box in the Prefs.

Verification

I posted two derivations for Euler's famous equation (here and here):

eix = cos x + i sin x

This can be verified in a particularly simple way.
The series representation of the exponential function:

ex = 1 + x + x2/2! + x3/3! + x4/4! + ..

is especially neat because each term in the series is the derivative of the term following, and the result of that is:

d/dx ex = ex

Which is, indeed, one definition of this function.
Substitution of ix leads to a simple shift in the pattern:

eix = 1 + ix - x2/2! - ix3/3! + x4/4! + ..

repeating with a period of 4. Powers of x like:

4n + 1 are multiplied by i 4n + 2 -1 4n + 3 -i 4n 1

But remembering the series for sine and cosine, and multiplying the former by i:

sin x = x - x3/3! + x5/5! - x7/7! + .. i sin x = ix - ix3/3! + .. cos x = 1 - x2/2! + x4/4! - x6/6! + ..

Adding:

cos x + i sin x = 1 + ix - x2/2! - ix3/3!+ x4/4! + .. = eix

Note on trig substitution

A fact we needed in a recent post deriving Euler's formula is the value of a particular integral:

∫ dx / √(x2 + 1) = ln (x + √(x2 + 1))

This gave me a lot more trouble than expected, so I thought I'd work through it here for reference. The easiest way is to differentiate the answer. Let:

u = x + √(x2 + 1) du = [1 + 2x / 2 √(x2 + 1)] dx = [1 + x / √(x2 + 1)] dx y = ln(u) dy = 1/u du = [1 / (x + √(x2 + 1)) ] [ 1 + x / (√(x2 + 1)] dx

The trick is to notice that when we find the common denominator for the part at the far right and then add the terms, we generate the denominator of the part to the left.

dy = [1 / (x + √(x2 + 1) ] [ (√(x2 + 1) + x) / (√(x2 + 1)] dx = 1 / (√(x2 + 1)

The integrand is normally written more generally as:

∫ 1 / √(x2 + a2)

and that doesn't change anything.

The forward version starts with a trigonometric substitution.

x = a tan y y = tan-1 (x/a)

Remembering that

(u/v)' = (vu' - uv')/v2 tan' = (cos2 + sin2)/cos2 = 1/cos2 = sec2

So

dx = a sec2 y dy

The other part of the integral is:

1 / √(x2 + a2) x = a tan y x2 = a2 tan2 y x2 + a2 = a2 tan2 y + a2 = a2 ( 1 + tan2 y)

If we start with the standard identity:

sin2 y + cos2 y = 1

and divide by cos2 y:

tan2 y + 1 = sec2 y x2 = a2 ( 1 + tan2 y) = a2 sec2 y

Thus the integral reduces to:

∫ (a sec2 y) / (a sec y) dy = ∫ sec y dy

This small integral is itself a bit tricky. The answer, which I found on a really nice site (here), is to multiply top and bottom by:

sec y + tan y sec y ------------- sec y + tan y Let u = sec y + tan y

If you work through it you'll see that

du = (sec y tan y + sec2 y) dy

Thus, this is

∫ du/u = ln(u) = ln (sec y + tan y) = ln [ ((√(x2 + a2) + x ) / a ]

As I said, not so easy as other trig substitutions.

Euler's Gem 2

Here is a sketch of a second derivation of Euler's famous formula:

eiθ = cosθ + i sinθ

as presented by William Dunham in his book Euler, The Master of Us All. First post here.

The first step is to recall a standard trigonometric substitution in calculus:

y = sin x x = sin-1 y √(1 - y2) = cos x

We're interested in the integral:

∫ dy / √(1 - y2)

Substituting with x we see that:

dy = cos x dx

And the integral is

∫ (1/cos x) cos x dx = ∫ dx = x x = ∫ dy / √(1 - y2)

Now Euler makes a complex change of variable:

y = iz x = ∫ dy / √(1 - y2) = ∫ i dz / √(1-(iz)2) = i ∫ dz / √(1 + z2) = i ln [√(1 + z2) + z]

The last step is another standard result from calculus which I will assume for the time being (more here).

Undo the substitution:

z = y/i = sin x / i z2 = -sin2 x √(1 + z2) = √(1 - sin2 x) = cos x x = i ln (cos x + sin x / i)

We will use two identities involving i:

u / i = - i u 1 / (cos u - i sin u) = (cos u + i sin u)

(For the second one, see the previous post). Now:

x = i ln (cos x + sin x / i) x = i ln (cos x - i sin x) ix = - ln (cos x - i sin x) = ln [ 1 / (cos x - i sin x) ] = ln (cos x + i sin x)

Just eponentiate:

eix = cos x + i sin x

Wow, again!

Euler's gem

Here is a sketch of the derivation of Euler's famous formula:

eiθ = cosθ + i sinθ

as presented by William Dunham in his book Euler, The Master of Us All.

The first part of the proof is similar to when we used Euler's formula to derive other formulas for trig functions of sums and differences of angles (post), only backward. Start from the definition of i:

i = √-1

To begin with, having i allows us to factor new expressions:

1 = cos2 s + sin2 s = (cos s + i sin s)(cos s - i sin s)

(I'm going to use s and t, as before, rather than θ and φ).

This shows where the original idea of cos + i sin comes from. (Of course, we could just as well do sin + i cos, that would result in a different convention for the orientation of the complex plane).

Suppose we have two angles s and t, we can multiply and then use the formulas from before (obtained by the geometric proof):

(cos s + i sin s)(cos t + i sin t) = = (cos s cos t - sin s sin t) + i(sin s cos t + cos s sin t) = cos(s + t) + i sin(s + t)

Set s = t:

(cos s + i sin s)2 = cos(2s) + i sin(2s)

In fact Euler showed it works for fractional n but I'll assume that part:

[1] (cos s + i sin s)n = cos(ns) + i sin(ns) n >= 1

If we multiply the difference rather than the sum:

(cos s - i sin s)(cos t - i sin t) = = (cos s cos t - sin s sin t) - i (sin s cos t + cos s sin t) = cos(s + t) - i sin(s + t)

Again, with s = t we have:

(cos s - i sin s)2 = cos(2s) - i sin(2s) [2] (cos s - i sin s)n = cos(ns) - i sin(ns)

Adding [1] and [2] we have:

2 cos(ns) = (cos s + i sin s)n + (cos s - i sin s)n

---

The middle part of the proof is where the magic happens. Let:

s = x/n

As

n -> ∞ s -> 0 cos s -> 1 sin s -> s

So..

cos x = cos ns = 1/2 [(cos s + i sin s)n + (cos s - i sin s)n] cos x = 1/2 [(1 + is)n + (1 - is)n] cos x = 1/2 [(1 + ix/n)n + (1 - ix/n)n]

But..

eix = (1 + ix/n)n as n -> ∞

So

cos x = 1/2 [eix + e-ix]

---

By very similar manipulation to what's in the first part we can also handle the sine:

2i sin(ns) = (cos s + i sin s)n - (cos s - i sin s)n

We will obtain:

sin x = 1/(2i) [eix - e-ix]

Now it's just a matter of addition:

cos x + i sin x = 1/2 [eix + e-ix + eix - e-ix] = 1/2 [eix + eix] = eix

Wow!

Rotating a hyperbola: general case

Just to finish up quickly with the analysis of hyperbolas from the other day (here and here), following Stewart (pdf), we derived the following equations for x and y in terms of a coordinate system based on u and v (orthogonal, rotated through angle θ):

x = u cosθ - v sinθ y = u sinθ + v cosθ

The general formula for a quadratic is:

Ax2 + Bxy + Cy2 + Dx + Ey + F = 0

When we transform to the new coordinates, we get:

Ax2 = A[u2 cos2θ - 2uv sinθ cosθ + v2 sin2θ] Bxy = B[u2 cosθ sinθ + uv cos2θ - uv sin2θ - v2 sinθ cosθ Cy2 = C[u2 sin2θ + 2uv sinθ cosθ + v2 cos2θ

Gathering the terms in uv we obtain the coefficients:

(C-A)(2sinθ cosθ) + B(cos2θ - sin2θ)

Remember the double angle formulas:

sin(s+t) = sin s cos t + cos s sin t sin(2θ) = 2 sinθ cosθ cos(s+t = cos s cos t - sin s sin t cos(2θ) = cos2θ - sin2θ

So we obtain:

(C-A) sin(2θ) + B cos(2θ)

for the coefficients of xy. These must equal zero for all the xy terms to disappear. Thus:

tan(2θ) = B / (C-A)

This approach runs into a problem if C = A, as it does for our example:

xy = 1

But we can just invert the step at the end:

cot(2θ) = (C-A)/B = 0

The cotangent is zero when the cosine is zero, e.g. for

2θ = π/2

Thus, if

θ = π/4

then all the xy terms vanish, as we found before.

Mistakes

A few months ago I posted briefly about Sal Khan's academy (here, his link here). He has put up an extraordinary number of talks (1000?) where you see video of what he's "writing" on the screen and Sal talks you through explanations of various topics. It's mostly math, but also some biology, even finance. If you don't know about it, check it out.

There's a lot to like about this project. For starters, Sal knows his linear algebra and calculus, not to mention algebra. Also, it's free. Also, the lectures come as bite-sized chunks which is just right for the modern student with an attention span measured in single-digit minutes if not seconds. One negative aspect is that the system for projecting what he writes is pretty marginal---some kind of a pen would be a lot better.

But the reason I'm re-posting about this is I found a rather egregious error in one of the videos, and it reminded me of my own experience. The error is central to the example (in probability and statistics), in which he explains that "the standard deviation is just the average distance from the average," and proceeds to calculate just as you would if that were true. Of course, it's not true.

Why it's done the way it is, that's a whole 'nother topic.

But that's not my point. This error reminded me of my first teaching experience, long ago. I was a first-year graduate student, and I had barely noticed that I was scheduled to run a "recitation section" for a eukaryotic genetics class. I left my DNA samples on ice and ran to the first meeting of class. Unfortunately, the students expected me to actually explain stuff to them, like what was in the reading assignment on meiosis.

Hmm... meiosis.

I thought I remembered that. Unfortunately, I remembered wrong. Homologs segregate at division 1, not division 2. Oops.

Almost none of the students in my section turned up for week 2! In fact, the other TA accused me of screwing up on purpose, so she would have more students in her section. (Cindy, I'm sorry, I really am).

Anyway, I did learn my lesson. You must know the material when you walk in the door. And if you don't, be honest (and show up next time with the answer).

I think a similar thing happened to Sal. He has so many videos that he failed to prepare. As James Baker said in his biography: 'prior preparation prevents poor performance.' The five P's. A motto to live by.

Rotational transformation: geometrical view

A bit more about rotational transformations (see the previous post for the setup to this). We derived equations describing the coordinates of a point at x,y with respect to a new coordinate system in u,v that is rotated counter-clockwise by an angle θ:

x = u cos(θ) - v sin(θ) y = u sin(θ) + v cos(θ)

We can also solve these equations for u and v. (It helps to know the answer). Notice that if we multiply the first equation by cos(θ) and the second one by sin(θ) and add, we'll end up with u by itself on the right-hand side:

x cos(θ) = u cos2(θ) - v sin(θ) cos(θ) y sin(θ) = u sin2(θ) + v cos(θ) sin(θ) x cos(θ) + y sin(θ) = u cos2(θ) + u sin2(θ) x cos(θ) + y sin(θ) = u

Similarly, if we multiply the first equation by sin(θ) and the second one by cos(θ) and subtract the first from the second, we end up with v by itself:

x sin(θ) = u sin(θ) cos(θ) - v sin2(θ) y cos(θ) = u sin(θ) cos(θ) + v cos2(θ) - x sin(θ) + y cos(&theta) = v cos2(θ) + v sin2(θ) - x sin(θ) + y cos(&theta) = v

I came up with geometric explanations for these relationships, the one for x and y in terms of u and v is at the top of the post, and below is the same drawing but explaining u and v in terms of x and y.



A quicker way to the same place is to start by considering u and v to be the original coordinate system, and rotate clockwise through the angle θ. Then θ is negative, sin(θ) = -sin(-θ) and the cosine stays the same. We end up switching the signs of the sine terms, as we found.

Different views, same hyperbola

In the algebra study guide that my son is using (ALEKS) they set up a hyperbola in two "standard" forms as:

(x-h)2/a2 - (y-k)2/b2 = 1 (y-k)2/b2 - (x-h)2/a2 = 1

or more generally:

(x-h)2/a2 - (y-k)2/b2 = +/- 1

Let's simplify our lives and consider an example centered at the origin:

x2/a2 - y2/b2 = 1

Then they divide the world into those hyperbolas that open up and down versus those that open up left and right.

By looking at the equation, I think we can agree that if x² equals 0 we've got a problem, since no value of y can satisfy the equation, and in fact x must not be less than 1. Hence, we can already predict that this system opens left and right, as the plot shows (a2 = b2 = 2):



The ALEKS review goes on to explain that for this type of equation, the vertices of the hyperbola (points of closest approach to the origin at h,k) are equal to

h +/- a, k

and the asymptotes have slopes equal to

+/- b/a

This all works great. The problem I had was that the simplest parabola I can think of is:

xy = 1

which doesn't fit the system.



The answer to my confusion is that any parabola may be rotated around its origin in the xy-plane. In 2 selected orientations its equation will have only x2 and y2 terms (when the vertices are on the x or the y-axis), whereas in 2 other selected orientations it may have only xy terms (when the vertices are on y= +/- x. The rest of the time it contains both.

This is explained in an excerpt from Stewart's Calculus I found here.

Consider a point P in the plane with coordinates x,y and distance from the origin r. Now, let's establish a new coordinate system u,v which is rotated counter-clockwise by the angle θ. A vector from the origin through P is rotated an angle φ with respect to the u axis and φ + θ with respect to the x-axis. Here's a screenshot from the pdf:



(Note, he uses capital X and Y for the second set of coordinates).

In the u,v system, the point P has coordinates

u = r cos(φ) v = r sin(φ)

In the x,y system the coordinates are:

x = r cos(φ + θ) y = r sin(φ + θ)

We remember (my post here) that:

sin(s+t) = sin s cos t + cos s sin t cos(s+t) = cos s cos t - sin s sin t

Hence:

x = r cos(φ) cos(θ) - r sin(φ) sin(θ) = u cos(θ) - v sin(θ) y = r sin(φ) cos(θ) + r cos(φ) sin(θ) = v cos(&theta) + u sin(θ) = u sin(&theta) + v cos(θ)

Now, consider the hyperbola:

xy = 1

Rewrite this in terms of u,v and multiply, giving 4 terms

+ u cos(θ) u sin(θ) + u cos(θ) v cos(θ) - u sin(θ) v sin(θ) - v sin(θ) v cos(θ)

If θ equals π/4, then

cos(θ) = sin(θ) = 1/√2

The two middle terms drop out and leave us with:

u2/2 - v2/2 = 1

If the middle terms don't cancel, we're left with a mixture including some fraction of
u²/2, v²/2, and u times v.

Note: we rotated the coordinate system counter-clockwise, which has the effect of rotating the plot clockwise, when the coordinate system is viewed in standard orientation.

One more thing about 3x3 determinants

While we're on the topic of determinants (see yesterday's post), rather than pounding on arithmetic drills, the other thing I would show algebra students is that you can use any row or column to set up the computation. That can be useful if one choice contains zeros.

The standard approach uses the top row:

a b c d e f g h i

And the determinant is:

a(ei-hf) - b(di-fg) + c(dh-eg)

Where the sign of the second term is negative by the checkerboard rule:

+ - + - + - + - +

Multiplying out we obtain:

aei - afh - bdi + bfg + cdh - ceg

The 6 terms contain three components, each taken from a different row and column. For example, the components of bdi are from:

b = (1,2) d = (2,1) i = (3,3)

The checkerboard rule makes the sign come out correctly.

If we're working with the top row or the middle column and so processing b x (di - fg) or d x (bi - ch), we'll need the minus sign; whereas if we're obtaining this term from i x (ae - bd) we already have a minus sign.

Let's try using the last row. We have:

a b c d e f g h i g(bf - ce) - h(af - cd) + i(ae - bd) bfg - cdg - afh + cdh + aei - bdi

Compare with the first example to see that all the terms are present.

Can we do it by the diagonal? Try ceg:

c(dh - eg) - e(ai - cg) + g(bf - ce)

Nope. Some terms are correct, but some are duplicates.

Cramer's rule calculation

We have algebra homework that involves using Cramer's rule to solve not only 2 x 2 but also 3 x 3 systems. It seems kind of silly since this method is overkill for 2 x 2, and would never be used for 4 x 4 or larger.

(Note on the wikipedia article, start about halfway down, where it says "Explicit formulas for small systems".)

Also, and this gets closer to the point, drilling by solving 3 x 3 matrices is not really about the rule, which is pretty simple. It's about making an easy problem harder by stuffing a lot of arithmetic into it. And to me, that is symptomatic of a big difficulty with math education as I'm encountering it through my son. Computers are much better at computing sums than humans. It's just silly to drill students on arithmetic. If you want to do something complicated, why not derive Cramer's rule?

So, I decided to write a solver for 3 x 3 systems in Python. I wouldn't say it's thoroughly debugged yet, so let me know if you run into a problem. With the example shown, I did get the same answer as this online calculator.

The first code segment contains the equations explicitly entered as an array. I'm sure you know how to modify it to read input from a file.

test.py

import numpy as np import Cramer def test_Cramer(): L = [2, 3, 0, 5, 1, 1, 1, 3, 2,-1, 3, 7] A = np.array(L) A.shape = (3,4) result = Cramer.solve(A) if result: x,y,z = result print 'solution' print 'x =', x print 'y =', y print 'z =', z, '\n' Cramer.check(A,x,y,z) test_Cramer()

The output looks like this:

> python test.py solve [[ 2 3 0 5] [ 1 1 1 3] [ 2 -1 3 7]] compute 3 x 3 det of [[ 2 3 0] [ 1 1 1] [ 2 -1 3]] D = 5 compute 3 x 3 det of [[ 5 3 0] [ 3 1 1] [ 7 -1 3]] Dx = 14 compute 3 x 3 det of [[2 5 0] [1 3 1] [2 7 3]] Dy = -1 compute 3 x 3 det of [[ 2 3 5] [ 1 1 3] [ 2 -1 7]] Dz = 2 solution x = 2.8 y = -0.2 z = 0.4 check row 0 = [2 3 0 5] 2.0*2.8 + 3.0*-0.2 + 0.0*0.4 = 5.0 row 1 = [1 1 1 3] 1.0*2.8 + 1.0*-0.2 + 1.0*0.4 = 3.0 row 2 = [ 2 -1 3 7] 2.0*2.8 + -1.0*-0.2 + 3.0*0.4 = 7.0

Cramer.py

import numpy as np def det2x2(A, v=False): if v: print 'compute 2 x 2 det of' if v: print A assert A.shape == (2,2) return A[0][0]*A[1][1] - A[0][1]*A[1][0] def det3x3(A): print 'compute 3 x 3 det of' print A assert A.shape == (3,3) a,b,c = A[0] c1 = a * det2x2(A[1:3,[1,2]]) c2 = b * det2x2(A[1:3,[0,2]]) c3 = c * det2x2(A[1:3,[0,1]]) return c1 - c2 + c3 def solve(A): print 'solve' print A, '\n' assert A.shape == (3,4) D = det3x3(A[:,:3]) print 'D = ', D, '\n' if D == 0: print 'no solution' return Dx = det3x3(A[:,[3,1,2]]) print 'Dx = ', Dx, '\n' Dy = det3x3(A[:,[0,3,2]]) print 'Dy = ', Dy, '\n' Dz = det3x3(A[:,[0,1,3]]) print 'Dz = ', Dz, '\n' return Dx*1.0/D, Dy*1.0/D, Dz*1.0/D def check(A,x,y,z): print 'check' for i,r in enumerate(A): print 'row', i, '=', r pL = list() for coeff,var in zip(r[:3],(x,y,z)): c = str(round(coeff,2)) v = str(round(var,2)) pL.append(c + '*' + v) print ' + '.join(pL), print ' =', r[0]*x + r[1]*y + r[2]*z, '\n'

Note about the sum of cosines formula

I was showing someone the derivation of the formulas for sums and differences of sines and cosines (my post here). Unfortunately, I have some trouble remembering these. The trick I used was to try to recall that the derivation started by analyzing the difference cos(s-t) and it's a particularly easy form:

cos(s-t) = cos s cos t + sin s sin t

Then it occurred to me that there is a fairly obvious point about this that should make it even clearer. Just remember that pattern is sine sine, cosine cosine, both terms positive.
Then suppose s = t, we have

cos2(s) + sin2(s) = 1

So, which function and for what combination of s with itself would we always get 1? Well, it's obviously the difference, which always equals zero (the sum, 2s, could be any angle). And which function always gives 1 with an argument of 0? The cosine, of course.

cos(s-t) = cos s cos t + sin s sin t

Getting to the formula for cos(s+t) just involves realizing that if we plug in u = -t we have

cos(s+u) = cos s cos(-u) + sin s sin(-u)

but

cos(-u) = cos(u) sin(-u) = - sin(u)

So it's the sine term in the formula that changes sign when we add.

cos(s+u) = cos s cos u - sin s sin u

As for the other one, perhaps the easiest is Euler:

eis = cos s + i sin s ei(s+t) = cos(s+t) + i sin(s+t) ei(s+t) = eis eit = (cos s + i sin s) (cos t + i sin t) = cos s cos t - sin s sin t + i (sin s cos t) + i (cos s sin t)

The real part gives us what we had before,

cos(s+t) = cos s cos t - sin s sin t

and the imaginary part is equal to the imaginary part of the sum from the previous line:

i sin(s+t) = i (sin s cos t) + i (cos s sin t) sin(s+t) = sin s cos t + cos s sin t

In fact, maybe this is enough by itself. :)

Law of sines, and cosines

Continuing with some homework, we're going to use vector algebra to prove two geometric theorems.

[This is a similar diagram to the one from last time, but with labels switched around--sorry for any confusion].

We have vectors (a, b and c) and the lengths of the corresponding sides (|a| = a, etc.); also, the angle opposite side a is labeled A and so on.

The law of sines states that the ratio of the length of each side to the sine of the angle opposite is the same:

a/sin(A) = b/sin(B) = c/sin(C)

Recall that the area of the parallelogram formed by a and b is given by the absolute value of the cross-product:

|a X b| = a b sin(C)

And the area of the triangle is one-half that. But we must obtain the same area no matter which two vectors we use to compute the cross-product, and no matter which orientation. Thus:

|c X -a| = a c sin(B)
|-b X c| = b c sin(A)

also have the same area.

a b sin(C) = a c sin(B) = b c sin(A)

This leads directly to the law of sines. Now let's relabel the triangle slightly.


If b - a looks a little funny, just consider that

a + b - a = b

The law of cosines states that:

c2 = a2 + b2 - 2 a b cos(C)

where c = |b - a|.
We can obtain this simply by expanding the dot product:

(b - a) • (b - a) =
= b • b - b • a - a • b + a • a
= b2 + a2 - 2 a b cos(C)

But

(b - a) • (b - a) = c2

So finally:

c2 = b2 + a2 - 2 a b cos(C)

Using vectors makes it easy.

Ceva using vectors--special case

Some time ago we looked at Ceva's theorem (post). I'm starting on a book about Vector Calculus, and saw this question early in Chapter 1.

Using vectors, prove that the lines from each vertex of a triangle to the midpoint of the opposite side cross at a single point. The picture is as shown below:


We solve this by constructing parametric equations for the midpoint lines. The first one starts at O and travels along the vector a + (b-a)/2, which is the diagonal of the parallelogram formed by a and b. The equation is:

u[a + (b-a)/2]

where u is the parameter. Similarly, construct the line extending from A to its opposing side. It starts from A and travels along the vector b/2 - a. The equation is:

a + v(b/2 - a)

At the point where the vectors cross, these are equal:

u[a + (b-a)/2] = a + v(b/2 - a)

I had a little trouble at this point, and I must confess I peeked at the answer. We can solve this by considering that both the a and the b terms must balance. (Even though a isn't perpendicular to b, the part of b which is perpendicular to a is a constant fraction k of the whole).

This leads to:

½ k ub = ½ k vb
u = v = ⅔

ua - ½ ua = a - va
½ u + v = 1
u = v = ⅔

Given this suggested solution u = v = ⅔ , we can easily verify that

u[a + (b-a)/2]
⅔[a + (b-a)/2]
= ⅔ a + ⅓ b - ⅓ a = ⅓ (a + b)

a + v(b/2 - a)
a + ⅔(b/2 - a)
= a + ⅓ b -⅔ a = ⅓ (a + b)

Now consider the third side. The midpoint line starts from B and travels along the vector a/2 - b. The equation is:

b + w(a/2 - b)

We observe that for w = ⅔ this becomes:

b + w(a/2 - b)
b + ⅔(a/2 - b)
b + ⅓ a - ⅔ b = ⅓ (a + b)

Thus, all three lines intersect at the same point.

I haven't got an extension to the general case just yet (i.e. any point, not just the midpoint).

Flag update

AM Armenia BD Bangladesh CU Cuba DO Dominican Republic JO Jordan KH Cambodia LA Laos LB Lebanon LV Latvia MK Macedonia MW Malawi NP Nepal SY Syria UG Uganda UZ Uzbekistan VG Virgin Islands (British)

Thanks for reading!

Why is time?

Watch and listen.

Intro to ANOVA

This is an introductory post on ANOVA (analysis of variance). We ask the question: given three (or more) groups of observations, is one or more of the group means significantly different from the others. We will compute an F-statistic, and compare that with an F-distribution (carry out an F-test). If the statistic exceeds the 95% quantile, we will reject the null hypothesis that the means are the same.

According to wikipedia:

the two-group case can be covered by a t-test (Gosset, 1908). When there are only two means to compare, the t-test and the ANOVA F-test are equivalent; the relation between ANOVA and t is given by F = t².

ANOVA is a versatile (and complex) set of methods. This is just an elementary application, where we'll use the R implementation on three simple groups of data, and then compute the result ourselves in Python to see how it works internally.

To begin with, we follow the simple example from Dalgaard. You will need R and the ISwR package (or just construct the "data frame" yourself). R code:

> library(package=ISwR) > data(red.cell.folate) > summary(red.cell.folate) folate ventilation Min. :206.0 N2O+O2,24h:8 1st Qu.:249.5 N2O+O2,op :9 Median :274.0 O2,24h :5 Mean :283.2 3rd Qu.:305.5 Max. :392.0 > red.cell.folate folate ventilation 1 243 N2O+O2,24h 2 251 N2O+O2,24h 3 275 N2O+O2,24h 4 291 N2O+O2,24h 5 347 N2O+O2,24h 6 354 N2O+O2,24h 7 380 N2O+O2,24h 8 392 N2O+O2,24h 9 206 N2O+O2,op 10 210 N2O+O2,op 11 226 N2O+O2,op 12 249 N2O+O2,op 13 255 N2O+O2,op 14 273 N2O+O2,op 15 285 N2O+O2,op 16 295 N2O+O2,op 17 309 N2O+O2,op 18 241 O2,24h 19 258 O2,24h 20 270 O2,24h 21 293 O2,24h 22 328 O2,24h

We have 22 values in 3 groups.

> attach(red.cell.folate) > anova(lm(folate~ventilation)) Analysis of Variance Table Response: folate Df Sum Sq Mean Sq F value Pr(>F) ventilation 2 15516 7757.9 3.7113 0.04359 * Residuals 19 39716 2090.3 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 >

The attach gives us access to the names of the columns of values (folate) and factors (ventilation). We make a plot (shown at the top of the post):

> plot( folate ~ ventilation, data = red.cell.folate ) > stripchart(x=folate~ventilation, pch=16,vertical=T,add=T,col='blue')

The value 0.04359 indicates that we have P < 0.05.

We write the data to disk, and remember that the groups have (respectively) 8, 9 and 5 values.

> setwd('Desktop') > write.table(folate,'data.txt',row.names=F,col.names=F)

We use the Python script below to compute the F-statistic:

python script.py .. MS_W 2090.32 MS_X 7757.88 F_stat 3.71

If you look in the R output above, you'll see the same values as given here for MS_W and MS_X and the F-statistic.

To get the right F-distribution, we need to know that the degrees of freedom are:

k-1 = 2 # k = number of groups N-k = 19 # N = total observations

Since 3.71 is just higher than the 95% quantile of this F-distribution we can reject the null hypothesis H₀.

I found a calculator online. You can see the results in the screenshot.



The underlying calculation is pretty simple. We compute sumsq, the sum of the squares of the differences from the mean for several sets of values and the relevant means.

For the within groups comparisons, using mathematical notation this is (i groups with j observations in each group):

Σ Σ (xij - mi)2 (over i) (over j)

In Python, for each group we sumsq for the group compared with the group mean, and add the results for all three groups to give SSD_W.

To carry out the between groups comparisons, we first compute the grand mean m (of all of the samples). Then for each group we compute:

Σ Σ (mi - m)2 (over i) (over j)

Since the squared value is the same within each group, this is equivalent to:

Σ ni (mi - m)2 (over i)

In the Python code this becomes:

len(g)*(mean(g)-m)**2

and sum these over all the groups to give SSD_X. This is a sum of squares of the group means.

Finally, we compute:

MS_W = SSD_W/(N-k) MS_X = SSD_X/(k-1) F_stat = MS_X/MS_W

As to why we do this, for now you will have to go read the article. The explanation in Dalgaard is particularly clear, indeed, the whole book is excellent.

There is a SciPy function to carry out ANOVA (stats.f_oneway), but I don't have SciPy installed right now at home, and this post is long enough already. That's for another day.

Python code:

fn = 'data.txt' FH = open(fn,'r') data = FH.read().strip().split() FH.close() data = [int(n) for n in data] A = data[:8] B = data[8:17] C = data[17:] #A = [243,251,275,291,347,354,380,392] #B = [206,210,226,249,255,273,285,295,309] #C = [241,258,270,293,328] def mean(L): return sum(L)*1.0/len(L) def sumsq(L): m = mean(L) print 'sumsq' rL = [(x-m)**2 for x in L] for n,p in zip(L,rL): print n, round(p,1) S = sum(rL) print 'total', round(S,1), '\n' return S def ANOVA(G): # variation within groups SSD_W = 0 for g in G: SSD_W += sumsq(g) # a bit awkward, just flattening the list of lists # to get the mean and N T = list() for g in G: T.extend(g) m = mean(T) # variation between groups (X for 'cross') SSD_X = 0 for g in G: SSD_X += len(g)*(mean(g)-m)**2 N = len(T) # 22 k = len(G) # 3 MS_W = SSD_W*1.0/(N-k) MS_X = SSD_X*1.0/(k-1) F_stat = MS_X/MS_W return MS_W, MS_X, F_stat MS_W, MS_X, F_stat = ANOVA([A,B,C]) print 'MS_W', round(MS_W,2) print 'MS_X', round(MS_X,2) print 'F_stat', round(F_stat,2)