An Interactive Environment for Advanced Mathematics

Introduction

Systems like Mathematica and Matlab allow users to work interactively with mathematical objects: evaluating expressions, creating and using variables, and operating on them within an interface that immediately executes the commands typed and prints results.

The Meta.Numerics library for scientific programming provides some of the advanced mathematical and statistical functionality of those expensive proprietary systems, but in the form of .NET library APIs that you call from your own programs. There is no "front end" that allows you to access Meta.Numerics functionality interactively, the way Mathematica and Matlab do. Building a "front end" user interface is an expensive and time-consuming proposition. Is there a way we can get some existing .NET interactive parser to expose Meta.Numerics functionality interactively?

IronPython provides a ready-made answer. IronPython is an implementation of the Python programming language using the .NET Dynamic Language Runtime (DLR). Its authors have written a so-called read-evaluate-print-loop (REPL) interface that allows users to interactively write Python programs. Since IronPython is a .NET language, it can use .NET APIs. So even though we are not particularly interested in programming in Python, we can use IronPython's REPL to interact with Meta.Numerics objects. This tutorial shows how.

Getting Started

Install IronPython (available at http://ironpython.codeplex.com/Release/ProjectReleases.aspx) and Meta.Numerics (available at http://metanumerics.codeplex.com/Release/ProjectReleases.aspx).

Unfortunately, the IronPython installer doesn't add IronPython to your Programs menu, but after installation, you will find the executable in an IronPython folder in your Program Files directory (e.g., C:\Program Files\IronPython 2.0.1\ipy.exe). Start the IronPython executable; you will see an interactive interface that looks like this:

To be able to access Meta.Numerics objects from within this interactive interface, you will need to add a reference to the Meta.Numerics assembly and import the objects you want. If you have installed the Meta.Numerics assembly to your computer's GAC (this occurs automatically if you use the MSI installer), you can just type:

>>> import clr
>>> clr.AddReference("Meta.Numerics")

(Just to be extra-clear: you don't type ">>>" -- it's just reproduced here to distinguish what you are typing at the IronPython prompt from what IronPython writes in response.)

If, on the other hand, the Meta.Numerics assembly file is just on your hard drive but not in your GAC (this occurs when you just extract the ZIP file), you can instead type:

>>> import clr
>>> clr.AddReferenceToFileAndPath(r"C:\Program Files\Meta.Numerics\Meta.Numerics.dll")

Of course, you should replace the actual path with the path to the Meta.Numerics assembly on your computer, if it is different.

Now that IronPython has successfully linked to the Meta.Numerics assembly, we just need to tell it which APIs to import into the default namespace. This is basically analogous to issuing the C# using or VB Imports statement. We'll just go ahead and do that now for all the APIs we plan to use.

>>> from Meta.Numerics import *
>>> from Meta.Numerics.Functions.AdvancedMath import *
>>> from Meta.Numerics.Functions.AdvancedIntegerMath import *
>>> from Meta.Numerics.Functions.OrthogonalPolynomials import *
>>> from Meta.Numerics.Functions.FunctionMath import *
>>> from Meta.Numerics.Statistics import *
>>> from Meta.Numerics.Matrices import *

The statement from Meta.Numerics.Statistics import * imports all the classes in the Meta.Numerics.Statistics namespace. The statements from Meta.Numerics.Functions.AdvancedMath import * import all the static methods from the FunctionMath class, so they can be accessed by just typing their names without a class prefix. For example, you will see below that we just type Gamma(0.5), instead of having to type AdvancedMath.Gamma(0.5).

To avoid having to type all that repeatedly, just create a text file containing the Python commands you've typed up to this point. (I can't do it for you because the path strings will depend on your installation.) Call the file "MetaNumerics.py" and put it, say, in your Meta.Numerics directory. The next time you want to use Meta.Numerics interactively within IronPython, you can just fire up IronPython, type...

>>> execfile(r"C:\Program Files\Meta.Numerics\MetaNumerics.py")

...and you'll be ready to go.

Basic and Advanced Mathematics

Since simple arithmetic expressions are just Python expressions, you can just type them directly.

>>> 2 + 3 * 4
14

Notice that Python respects order of operations. You can also assign variables and evaluate expressions containing them.

>>> a = 6 / 2
>>> a + 2
5

Of course, that didn't use Meta.Numerics at all. To do so, let's invoke an advanced function.

>>> Gamma(0.5)
1.77245385091
>>> Gamma(5) 
24.0

Note that Γ(1/2) is π^1/2 and Γ(5) is 4! You can also compute incomplete Gamma functions, Beta functions, and Psi functions. How about a Riemann Zeta function?

>>> RiemannZeta(3.0)
1.20205690316

No problem. You can also compute Dirichlet's Eta function. For a complete list of available special functions, consult the Meta.Numerics documentation at the project website.

Now let's get a little fancier. Let's find the first root of the zero^th order Bessel function J₀. Meta.Numerics's FindZero function expects to be given a function that takes one parameter: a double. But Meta.Numerics's BesselJ function takes two parameters: an integer (the order) and a double (the argument). To get over this small hurdle, we just define a Python function that returns BesselJ(0,x).

>>> def J0(x): return(BesselJ(0,x))
>>> x0 = FindZero(J0, 3.0)
>>> x0
2.40482555770

(The second argument of the FindZero function is just a starting point in the search for a root.) Now let's find the integral of J₀ from the origin out to that first root.

>>> Integrate(J0, Interval.FromEndpoints(0.0,x0))
1.47030004338

Pretty cool.

Statistics

Now let's move on to some data analysis. Meta.Numerics defines different classes for different kinds of experimental data. Each class offers the appropriate descriptive statistics, statistical tests, and fitting procedures. We will give several examples.

Sample

Suppose we have some time-to-failure measurements for a component. Let's put them in an appropriate data container.

>>> s = Sample()
>>> s.Add(1.0)
>>> s.Add(1.3)
>>> s.Add(1.5)
>>> s.Add(1.7)
>>> s.Add(2.9)

Oops, that last data point was supposed to be 1.9. Good thing we are working in an interactive framework.

>>> s.Remove(2.9)
True
>>> s.Add(1.9)

So, what's the mean time-to-failure?

>>> s.PopulationMean
1.48 ± 0.156204993518133

Notice that Meta.Numerics gave you an uncertainty on its estimate of a population parameter. It will also give estimates, with uncertainty, for arbitrary higher moments.

Suppose we want to guarantee a mean time-to-failure above 1.25. Do our data give us 95% confidence that we are doing so? To find out, do a t test.

>>> t = s.StudentTTest(1.25)
>>> t.Statistic
1.47242411923
>>> t.LeftProbability
0.892554887728

Not quite. Our data makes us 89% confident that we have met the requirement, but not 95% confident.

DataSet

Next let's work with some data with error bars. Suppose we have taken the following measurements of, say, the water level in a reservoir as a function of time.

Time (h)	Water Level (m)
2.0	2.5 ± 0.5
3.0	2.0 ± 1.0
5.0	-1.0 ± 0.5
8.0	-4.0 ± 0.5
10.0	-1.5 ± 0.5
12.0	0.0 ± 1.0

Create a DataSet object to work with this data.

>>> D = DataSet()
>>> D.Add(2.0,2.5,0.5)
>>> D.Add(3.0,2.0,1.0)
>>> D.Add(5.0,-1.0,0.5)
>>> D.Add(8.0,-4.0,0.5)
>>> D.Add(10.0,-1.5,0.5)
>>> D.Add(12.0,0.0,1.0)
>>> D.Count6

Let's try to model the data with a steady rate of change.

>>> f1 = D.FitToLine()
>>> f1.GoodnessOfFit.LeftProbability
0.999999999654

That model is pretty well ruled out: there is only a 0.0000000001 chance of getting such a bad fit under the null hypothesis that the model is correct. How about a sinusoidal oscillation?

>>> def ff(a,x): return(a[0] * Math.Sin(2.0 * PI * x / a[1] + a[2]))
>>> f2 = D.FitToFunction(ff,(2,10,0))
>>> f2.GoodnessOfFit.LeftProbability
0.690194092616

That fits pretty well. (Note that (2,10,0) was just an initial guess for the amplitude, period, and phase parameters. As long as your initial guess isn't too crazy, its precise value isn't important and Meta.Numerics will find the best-fit point.) So, what's the best-fit oscillation period?

>>> f2.Parameter(1)
14.5075811109964 ± 0.989537307149617

About 14 hours.

BinaryContingencyTable

Let's work with one more kind of experimental data. Suppose we randomly assign patients to receive or not receive an experimental treatment and obtain the following results.

	Recovered	Died	Total
Treated	12	6	18
Untreated	8	15	23
Total	20	21	41

This is a 2 X 2 contingency table, a kind of experiment for which Meta.Numerics offers a tailored class.

>>> B = BinaryContingencyTable()
>>> B[0,0] = 12
>>> B[0,1] = 6
>>> B[1,0] = 8
>>> B[1,1] = 15

The treatment appears to help, but the sample size is small: how sure can we be that this isn't just a statistical fluke? Some entries are too small for the assumptions of Pearson's X² test to be justified. Fortunately, Meta.Numerics implements Fisher's exact test, which is appropriate even for small samples.

>>> B.FisherExactTest().LeftProbability
0.061603593338

There's only a 6% chance of getting such a strong result under the null hypothesis that the treatment does nothing. So we can be pretty confident, but not quite 95% confident as the convention for medical acceptance usually demands.

So, about how much better are my chances with the treatment than without it?

>>> B.OddsRatio
3.75 ± 2.49217525467211

About 4 times better, but with a rather large uncertainty.

Consult the Meta.Numerics documentation for information on additional capabilities of these data analysis classes, and to learn about a few more data analysis classes.

Matrix Operations

Meta.Numerics defines specific classes for different kinds of matrices: e.g., SquareMatrix, SymmetricMatrix, TridiagonalMatrix. Each class has methods appropriate to that kind of matrix, with implementations optimized to take advantage of the matrix structure.

For example, let's create a 3 X 3 symmetric matrix.

>>> M = SymmetricMatrix(3)
>>> M[0,0] = 6
>>> M[0,1] = 4
>>> M[0,2] = 1
>>> M[1,1] = 5
>>> M[1,2] = 2
>>> M[2,2] = 3
>>> M
{  6  4  1  }
{  4  5  2  }
{  1  2  3  }

Notice that, because Meta.Numerics knew that our matrix was symmetric, we only had to tell it half the matrix elements. It filled in the others by symmetry, and, underneath the hood, is using only half the storage space that SquareMatrix would require to store the same entries.

A Cholesky Decomposition is a decomposition M = S S^T, where S is lower-triangular. (S is often called the "square root" of the matrix M.) It exists only for positive definite symmetric matrices, but if you have such a matrix, Cholesky Decomposition is faster and more stable than LU decomposition. The SymmetricMatrix class offers this decomposition.

>>> CD = M.CholeskyDecomposition()
>>> S = CD.SquareRootMatrix()
>>> S
{   2.44948974278               0              0  }
{   1.63299316186   1.52752523165              0  }
{  0.408248290464  0.872871560944  1.43924583426  }

Let's verify that the decomposition is correct.

>>> S * S.Transpose()
{  6  4  1  }
{  4  5  2  }
{  1  2  3  }

Note that the evaluation of the expression requires matrix multiplication. You can use the Cholesky Decomposition to solve systems of equations...

>>> CD.Solve((1.0,2.0,3.0))
{  0  }
{  0  }
{  1  }

...or to invert the matrix.

>>> MI = CD.Inverse()

Meta.Numerics can also find eigenvalues and eigenvectors. For a fun application of Markov matrix diagonalization to analyze the board game Monopoly, see the CodeProject article "Markov Monopoly".

In addition to giving us the ability to dynamically interact with our library, the IronPython interface lets us reflect over our library in some useful ways. Want to know the type of that last result?

>>> type(MI)
<type 'SymmetricMatrix'>

Want to know what properties and methods that type makes available?

>>> dir(MI)
['CholeskyDecomposition', 'Clone', 'ColumnCount', 'Dimension', 'Eigensystem', 'E
igenvalues', 'Equals', 'GetHashCode', 'GetType', 'Inverse', 'Item', 'MemberwiseC
lone', 'ReferenceEquals', 'RowCount', 'ToString', 'Trace', 'Transpose', '__add__
', '__class__', '__delattr__', '__div__', '__doc__', '__eq__', '__getattribute__
', '__getitem__', '__hash__', '__init__', '__mul__', '__ne__', '__new__', '__rad
d__', '__reduce__', '__reduce_ex__', '__repr__', '__rmul__', '__rsub__', '__seta
ttr__', '__setitem__', '__str__', '__sub__']

Conclusion

You have seen how to use the Meta.Numerics library interactively within the IronPython interface. This technique has added a lot of value with very little work. The idea can easily be applied to allow the interactive use of any .NET library.