Maxima Latex Listings

The listings package is a really nice way to include typeset source-code listings in Latex documents. I wanted to include an appendix in one of my documents to show the Maxima calculations I used to arrive at some results shown in the paper. Unfortunately Maxima is not one of the (many, many) built-in languages recognized by the package. I tried using C, because the comments in Maxima are like C comments, but it doesn't handle single quotes (used to get noun forms like 'diff()) well. So, how to include listings of Maxima code in Latex?

It turns out to be really easy to define your own language in the preamble of the document.

\lstdefinelanguage{Maxima}{
keywords={addrow,addcol,zeromatrix,ident,augcoefmatrix,ratsubst,diff,ev,tex,%
with_stdout,nouns,express,depends,load,submatrix,div,grad,curl,%
rootscontract,solve,part,assume,sqrt,integrate,abs,inf,exp},
sensitive=true,
comment=[n][\itshape]{/*}{*/}
}

All this definition does is italicize comments and bold face the keywords, which is sufficient to make the code more readable and pleasant looking. Here's an example screenshot of the results:

Performance optimization (allocation inside a for loop)

An interesting discussion about an old topic just popped up on the Octave lists. It's not interesting because someone noticed that failing to pre-allocate a vector before entering a loop is slow (doesn't anyone RTFM?), but because of the follow-up discussion on indexing and range objects.

Here's the interesting bit:

octave:1> tic(); n=1e5; retval=1:n; toc()
Elapsed time is 0.000528962 seconds.
octave:2> tic(); n=1e5; retval = (1:n)(1:n); toc
Elapsed time is 0.00593709 seconds.
octave:3> tic();n=1e5;retval=[1:n]; toc
Elapsed time is 0.010952 seconds.

Why the significant difference in performance? According to jwe:

In Octave, an expression like 1:n creates a range object, which contains only the base, limit, and increment as double precision values, so no matter how many elements are in the range, it only takes a few bytes of storage (24 for the data plus some overhead for the internal octave_value object itself).

If you write [1:n], you force a Matrix object with N elements to be created. It will require 8*N bytes of storage, plus the overhead for the internal octave_value object itself.

Similar to the difference between range() and xrange() in Python.

Ares/Orion Slip and Classic Mistake

The latest news about the new NASA Constellation program indicates that they are having cost and schedule problems, this is no surprise, everyone has cost and schedule problems. The surprising thing is that they are deleting tests from the planned program development.

CxP attempted to protect the schedule and budgetary pressures by offsetting these additional strains by deleting test items - notably on the Upper Stage. However, this only proved to cause further disconnects throughout the program.

I suppose it should not be too surprising, this is a classic program management mistake. When the schedule gets tight and the costs get too high above estimates, cut testing. Why? Because it is at the end of the program, and program managers have a bit of hubris that leads them to think that the design is sound because they worked so hard on it up-front. Which is obviously the case because they had to spend more time than they thought in that stage of the program, right? Maybe, but also maybe not, that's why the tests are critically important: to validate the design performance. No slide-show, CAD drawing, CFD or FEA simulation can provide that sort of knowledge.

I can't understand why someone would think that under-testing a new product that is probably delayed and expensive because of its new technology and complexity would produce a good result. Rather than removing unneeded technology or complexity, rather than simplifying or reducing the requirements, they are rushing to field a fragile system. It is unfortunate that the burden of this foolishness is borne by the astronaut corps and not the program managers responsible for these short-sited expediencies. The program manager will probably get a reward, move to a new job, and then astronauts will be at undue risk in an already risky business.

Complex Step Jacobian-Free Newton-Krylov

Newton's method are really useful for solving non-linear systems.

Each iteration of the Newton's method requires solving a linear equation:

where the Jacobian is the matrix of partial derivatives

and then applying the update to the solution vector:

For small systems it is tractable to analytically calculate the Jacobian and invert it directly. This rapidly becomes intractable for realistic engineering problems because the number of unknowns can be on the order of 1e6 (if you have a good workstation) to 1e9 (if you have a big cluster). Rather than explicitly creating a 1e6-squared matrix and inverting it at every iteration, an iterative method that only requires the action of the matrix on a vector can be used (such as conjugate gradient or GMRES). Then we can approximate this matrix vector product by

which is a good improvement in storage over the matrix approach, but we can improve it still further by using a complex step

This means we only need one system evaluation to approximate the Jacobian, and much more accurate estimates can be made. As shown in a previous post it is second order accurate, and not as limited by subtractive cancellation as the standard difference approach is.

Here's a simple example of the idea in Maxima

/* trying out complex step finite differences for jacobian free newton
krylov methods */

/* define a simple system of equations */
F : zeromatrix(2,1);
F[1,1] : x[1]*x[2] + 2*x[2];
F[2,1] : 2*x[1] + x[1]*x[2];

/* analytic Jacobian */
J : zeromatrix(2,2);
J[1,1] : diff(F[1,1],x[1]);
J[1,2] : diff(F[1,1],x[2]);
J[2,1] : diff(F[2,1],x[1]);
J[2,2] : diff(F[2,1],x[2]);

u : zeromatrix(2,1);
v : zeromatrix(2,1);
u[1,1] : 1.5;
u[2,1] : 1.0;
v[1,1] : 0.2;
v[2,1] : 0.1;

c_step : 0.0 + 1.0e-200*%i; /* complex step */

/* approximate Jacobian using complex step */
Jv_approx : (
ev(F, x[1] = u[1,1] + c_step*v[1,1], x[2] = u[2,1] + c_step*v[2,1])
) / imagpart(c_step);
imagpart(Jv_approx);

/* evaluate the analytic Jacobian to compare */
Jv_analyt : ev(J.v,x[1]=u[1,1],x[2]=u[2,1]);

A good survey article introducing the JFNK method and describing some applications:
Jacobian-free Newton-Krylov methods: a survey of approaches and applications

Nice short paper about using complex steps in finite differences:
Using Complex Variables to Estimate Derivatives of Real Functions

Three Dimensional Poisson's Equation

This is a follow up to the post I wrote a while back about implementing numerical methods. The goal is to demonstrate a work-flow using tools that come included in the Fedora Linux (and many other) distributions.

1. In Maxima
   
   • Define governing equations
     
   • Substitute finite difference expressions for differential expressions in the governing equations
     
   • Output appropriate difference expressions to Fortran using f90()
     
2. In Fortran: wrap the expression generated by Maxima in appropriate functions, subroutines and modules
   
   
3. In Python
   
   • Compile modules with f2py, which comes packaged along with numpy
     
   • Write the high-level application code, parsing input decks, performing optimization, grid convergence studies, or visualizations that use the compiled Fortran modules
     

Maxima

The governing equation is the three-dimensional Poisson's equation. In Cartesian coordinates The Maxima code to define this equation is straightforward:

depends(u,x); depends(u,y); depends(u,z); /* Poisson's equation: */ p : 'diff(u,x,2) + 'diff(u,y,2) + 'diff(u,z,2) = -f;

Now that we have the continuous differential equation defined we can decide on what sort of discrete approximation we are going to solve. The simple thing to do with second derivatives is to use central differences of second order accuracy, which require only information from closest neighbours in the finite difference grid. Replacing the differential expressions with difference expressions is accomplished in Maxima with the ratsubst() function.

/* substitue difference expressions for differential ones: */
p : ratsubst((u[i-1,j,k] - 2*u[i,j,k] + u[i+1,j,k])/dx**2 , 'diff(u,x,2), p);
p : ratsubst((u[i,j-1,k] - 2*u[i,j,k] + u[i,j+1,k])/dy**2 , 'diff(u,y,2), p);
p : ratsubst((u[i,j,k-1] - 2*u[i,j,k] + u[i,j,k+1])/dz**2 , 'diff(u,z,2), p);
/* substitute the correct array value of f:*/
p : ratsubst(f[i,j,k], f, p);

This assumes that the solution u and the forcing function f are stored in three dimensional arrays, if not then the indexing in the difference expressions becomes more complicated. If we want to apply a Gauss-Seidel method to solve our elliptic problem we just need to solve p for u[i,j,k], and then dump that expression to Fortran format.

gs : solve(p, u[i,j,k]);
del : part(gs[1],2) - u[i,j,k];
/* output for fortran: */
with_stdout("gs.f90", f90(d = expand(del)));

This gives the expression for the update to the solution at each iteration of the solver.

Fortran

The Maxima expressions developed above need control flow added so they get applied properly to our solution array.

do k = 1, nz
do j = 1, ny
do i = 1, nx
d = (dy**2*dz**2*u(i+1,j,k)+dx**2*dz**2*u(i,j+1,k)+dx**2*dy**2* & u(i,j,k+1)+dx**2*dy**2*dz**2*f(i,j,k)+dx**2*dy**2* & u(i,j,k-1)+dx**2*dz**2*u(i,j-1,k)+dy**2*dz**2*u(i-1,j,k))/ & ((2*dy**2+2*dx**2)*dz**2+2*dx**2*dy**2)-u(i,j,k) u(i,j,k) = u(i,j,k) + w*d
end do
end do
end do

The expression for the update was generated from our governing equations in Maxima, one of the things to notice is the w that multiplies our update, this is so we can do successive over relaxation. Also, we'll go ahead and package several of these iterative schemes into a Fortran module

module iter_schemes
implicit none
contains

...bunches of related subroutines and functions ...

end module iter_schemes

Python

Now we want to be able to call our Fortran solution schemes from Python. Using F2py makes this quite simple:

[command-line]$ f2py -m iter_schemes -c iter_schemes.f90

This should result in a python module iter_schemes.so that can be imported just as any other module.

import numpy
from iter_schemes import *

... do some cool stuff with our new module in Python ...

Conclusion

This approach might seem like over-kill for the fairly simple scalar equation we used, but think about the complicated update expressions that can be generated for large vector partial differential equations like the Navier-Stokes or Maxwell's equations. Having these expressions automatically generated and ready to be wrapped in loops can save lots of initial development and subsequent debugging time. It also makes generating unit tests easier, and so encourages good development practice. But does it work? Of course, here's the convergence of the Gauss-Seidel and SOR schemes discussed above along with a Jacobi scheme and a conjugate gradient scheme applied to a 50 by 50 by 50 grid with a manufactured solution. The plot was generated with the Python module matplotlib.

Vectorize Matrix Slice Extraction

There was an interesting thread on the Octave mailing lists recently about vectorizing matrix slice extraction. That is avoiding the dreaded for-loop in something like:

%% create a random sparse matrix
n = 10000; %% matrix dimension
nz = 300; %% number of non-zero elements

ia = floor(rand(nz,1)*n); %% dimension 1 indices
ja = floor(rand(nz,1)*n); %% dimension 2 indices
A(:,:) = sparse(ia, ja, rand(nz,1), n, n, n);

%% taking slices with a for loop, slow, poor scaling to large matrices:
for jj = nz:-1:1
slow_columns(:,jj) = A(:,ja(jj));
endfor

%% example using built-in sub2ind, fast, good scaling to large matrices:
jj = sub2ind(size(A)(1),ja);
fast_columns = A(:,jj);

This is a simplified example, the one in the thread shows how you can access a 5-dimensional matrix with 4 indexes. Knowing the built-in functions that come with your high-level language is well worth it, even if it means you have to do weird things to stay in the high-level language and still get reasonable speed. This allows you to avoid slow looping and use the quick compiled functions that someone else has already spent time and effort developing.

Unblurring Images

Tikhonov regularization is a pretty cool concept when applied to image processing. It lets you 'de-blur' an image, which is an ill-posed problem. The results won't be as spectacular as CSI, but still cool nonetheless.

Original

The original image, CSI fans pretend it is a license plate or a criminal's face in a low quality surveillance video, was made by taking a screenshot of some words in a word processor. The image doesn't have to be character based, but characters make good practice problems.

Blurring

We then blur the image with a Toeplitz matrix. In Octave that's pretty easy:

im = imread("foo-fiddle.png");
[m, n] = size(im);

sigma_true = 10;
p = exp(-(0:n-1).^2 / (2 * sigma_true^2))';
p = p + [0; p(n:-1:2)];
p = p / sum(p);

G_true = toeplitz(p);

im = G_true*double(im)*G_true';

Generally after doing some floating point operations on an image you won't end up with something still scaled between 0 and 255, so the little function shown below rescales a vector's values to the proper range for a grayscale image.

function im = image_rescale(double_vector)
%% scale values linearly to lie between 0 and 255
m = max(double_vector);
n = min(double_vector);
scale = max(double_vector) - min(double_vector);
im = double_vector - (min(double_vector) + scale/2);
im = im ./ scale;
im = (255 .* im) + 127.5;
endfunction

De-Blurring

This is where the regularization comes in. Since the problem is poorly conditioned we can't just recover the image by straightforward inversion. We'll use the singular value decomposition of our blurring matrix to create a regularized inverse, this is similar to the pseudoinverse. If the singular value decomposition of the blurring matrix is


Then the regularized inverse is


Where the entries in the diagonal matrix are given by


So as alpha approaches zero we approach the un-regularized problem. Chapter 4 of Watkins has good coverage of doing SVDs in Matlab. The Octave to perform this process is shown below in the function tik_reg().

function im_tik = tik_reg(im,sigma,alpha)
[m,n] = size(im);
p = exp(-(0:n-1).^2 / (2 * sigma^2))';
p = p + [0; p(n:-1:2)];
p = p / sum(p);

G = toeplitz(p);

[U, S, V] = svd(G,0);

Sinv = S / (S^2 + alpha*eye(length(S),length(S))^2);

im_tik = V*Sinv*U'*double(im)*V*Sinv*U';
endfunction

This gives us two parameters (sigma and alpha) we have to guess in a real image processing problem (because we won't be the ones blurring the original). Another weakness of the method is that we assumed we knew the model that produced the blur, if the real physical process that caused the blur is not similar to Gaussian noise, then we won't be able to extract the original image.

Here's a couple for-loops in Octave that try several combinations of alpha and sigma.

sigma = [9,10,11];
alpha = [1e-6,1e-7,1e-8];

for j = 1:length(sigma)
for i = 1:length(alpha)
unblur = tik_reg(im, sigma(j), alpha(i));
unblur = reshape(image_rescale(reshape(unblur,m*n,1)),m,n);
fname = sprintf("unblur_%d%d.png",i,j);
imwrite(fname, uint8(unblur));
endfor
endfor

Optical character recognition probably wouldn't work on the output, but a person can probably decipher what the original image said.

CFD Integrity

Quantifying uncertainty in simulation predictions is something at which most modelling communities fall short, and the Computational Fluid Dynamics community is no exception. There is plenty published on uncertainty estimation, but the state of the practice lags the state of the art considerably in this area.

It does not help matters that it is very hard and expensive to do a thorough job of quantifying the uncertainty of the result of a complex calculation. If it takes a week to get one data point, doing a Monte Carlo sensitivity analysis is probably out of the question, and differentiating your code could take years. It also doesn't help that many of the decision makers who are the consumers of the CFD product don't really want to hear how uncertain the results are, often they want pretty pictures and a warm fuzzy feeling (colorful fluid dynamics is great for marketing too).

I highly recommend the article on quantifying uncertainty by P.J. Roache. It is well written, and it is a special treat to read a CFD paper with a sentence like,

Clearly, we are not interested in such worthless semantics and effete philosophizing, but in practical definitions, applied in the context of engineering and science accuracy.

Roache is addressing the hand wringing made popular in the "soft" academic communities about the ability to define "truth".

Other than being a good basic engineering practice, making sure users of your results are well aware of the uncertainty is a matter of integrity too. There's an obligation that comes with expertise. Boole's words on the subject still ring true. More recently Feynman said much the same thing about integrity as an expert advising others:

I'm talking about a specific, extra type of integrity that is not lying, but bending over backwards to show how you are maybe wrong, that you ought to have when acting as a scientist. And this is our responsibility as scientists, certainly to other scientists, and I think to laymen.

Convective Flux Jacobians

The Euler equations are an important system of equations in computational fluid dynamics (CFD). They are basically the full continuum equations of motion without the viscosity. Adding in the viscosity gives you the more realistic Navier-Stokes equations.

The conserved variable vector is

And the flux vector is

The flux Jacobian is simply the partial derivative of the flux with respect to the conserved variables (mass, momentum, energy).

The maxima to calculate the Jacobian is straightforward:

U : matrix([u1], [u2], [u3]);

f1 : u2;
f2 : u2^2 / u1 + p_exp;
f3 : (u2 / u1) * (u3 + p_exp);
F : matrix([f1], [f2], [f3]);

A : zeromatrix(nd, nd);

for j : 1 thru nd do
(
for i : 1 thru nd do
(
A[i,j] : diff(F[i,1], U[j,1])
)
);

Diagonalization

Using the chain rule lets you write the Euler equations in a quasi-linear form in terms of the flux Jacobian, now called A.

The flux Jacobian can be diagonalized into a matrix of right eigenvectors, diagonal matrix of eigenvalues, and a matrix of left eigenvectors.

The code to accomplish this in maxima is a little more involved, because of some substitutions needed to get everything in terms of u,gamma,a.

load("eigen");
VLV : eigenvectors(simp_A);

a_exp : (g * (g - 1) / u1) * (u3 - u2^2 / (2 * u1));
a_exp : sqrt(a_exp);
simp_a : ratsimp(ev(a_exp, u1 = rho, u2 = rho * u, u3 = rho * E));
a_eqn : a^2 = simp_a^2;
c_exp : g*(2*E-u^2);
a_eqn : facsum(a_eqn, c_exp);

/* put the eigenvectors into a matrix: */
V : columnvector( VLV[2] );
V : addcol( V,VLV[3],VLV[4] );
V : ratsimp( V );

/* put the eigenvalues into a diagonal matrix: */
load("diag");
L : diag( VLV[1][1] );

radsubstflag : true;
/* substitute in the speed of sound: */
V : ratsimp(ratsubst(a, simp_a, ratsimp(V)));
V : factor(ratsimp(ratsubst(a, simp_a, rootscontract(V))));
/* factor out gammas */
V : facsum(V,g);
c_a : part(solve(ratsubst(c,c_exp,a_eqn),c)[1],2);
V : ratsubst(c,c_exp,V);
V : ratsubst(ev(c_a),c,V);
V : facsum(V,u);

/* substitute in the speed of sound: */
L : ratsimp(ratsubst(a, simp_a, ratsimp(L)));

/* left eigenvectors: */
Vinv : factor( ratsimp( invert( V ) ) );
/* factor out gammas */
Vinv : facsum(Vinv,g);
d_exp : (g-1)/a^2;
Vinv : ratsubst(d,d_exp,Vinv);

Conclusion

The thing to note about this result is that there are three distinct wave speeds in the one dimensional Euler equations, u-a, u+a, and u, where a is the speed of sound. There are lots of possible ways to write the eigenvectors, in terms of the enthalpy, or the energy, or the pressure, but the form above with only the flow velocity, the speed of sound, and the ratio of specific heats is pretty elegant. It really illustrates what a fundamental quantity the speed of sound is for this set of equations.

Chapter 3 of Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction has a good work-up of this development and shows some alternative formulations for the flux Jacobian.