Math.NET Project News

Revised Interpolation Toolkit in Iridium

Christoph Rüegg — Sat, 10 May 2008 15:31:00 GMT

The next release of Math.NT Iridium, Iteration 16, comes with a revised interpolation architecture and implementation.

Up to now, the interpolation classes have been a bit akward to use, partially because of the SampleList collection class you had to use, and because of the design in general. It also provided only two interpolation algorithms, which are both somewhat outdated these days.

The new implementation provides some newer more stable algorithms (like Floater and Hormann's algorithm for pole-free rational interpolation) together with a cleaner design. Additionally, some of the algoriths can also provide the first and second derivative and in the case of splines even a definite integration. There is also a new facade/portal class that reduces building an interpolation to one simple method call. Usually you would just use this facade class to build/precompute an interpolation, but all the algorithms are also publicly available in the Algorithms-namespace, so if you know what you're doing you can use them directly.

Sample code, which uses a pole-free rational barycentric interpolation (the default algorithm) with 5 given sample pairs (t, x(t)):

double[] t = new double[] { -2, -1, 0, 1, 2};
double[] x = new double[] { 1, 2, -1, 0, 1};

IInterpolationMethod method =
     Interpolation.Create(t, x);

double res = method.Interpolate(0.5);

Simple, isn't it?

Unfortunately it was not possible to fit the new design into the old classes, so the new classes and interfaces replace the old classes completely. These old classes are still there for now and continue to work, but they're marked as obsolete and we recommend strongly to upgrade your code base to the new architecture.

The new architecture has already been checked in to the source repository. If you're interested, please have a look at it and provide feedback - it's not released yet so we can still change it completely :).

Iridium 2008 April Release (v2008.4.14.425)

Christoph Rüegg — Sun, 06 Apr 2008 11:50:00 GMT

Iteration 14 of Math.NET Iridium (Numerics) is now available, as 2008 April Release with Version 2008.4.14.425. Grab it here.

Sorry for the very short release cycle (just one week after iteration 12). The reason is that I won't be able to work on Math.NET for the next three weeks and that some of the fixes and changes are important enough to not let you wait three weeks for no reason.

Please continue reporting issues and bugs you find, it's very useful and helps making the whole project better. There's also a big chance that the issue will actually be fixed: in the last few releases we always managed to fix all bugs we were aware of at that point. Thanks!

Team: Christoph Rüegg, Joannès Vermorel, Matthew Kitchin

Summary:

Bugs: All known 3 bugs have been fixed.
Better special function precision (Gamma, Beta, Erf, Distributions etc): now up to 12 - 14 digits.
New direct/real gamma function, new harmonic number function.
Interpolation: usability enhancements (better double-array support, less user code)

New Features:

IRID-122: Core - New direct Gamma function (additional to GammaLn) with negative value support
IRID-123: Core - New Special Function: Harmonic Numbers

Enhancements:

IRID-121: Core - Better numerical precision for Gamma function
IRID-125: Interpolation - Additional interpolation and sample list constructors for double arrays.
IRID-126: Interpolation - Better interpolation order access and defaults

Fixed Bugs:

IRID-119: Interpolation - Polynomial Extrapolation in positive direction throws IndexOutOfRangeException
IRID-120: Linear Algebra - Infinite recursion
IRID-124: Linear Algebra - Matrix.CopyToArray - wrong indexer in inner loop condition.

For more details have a look at the Iridium tracker change log.

Iridium 2008 March Release (v2008.3.12.405)

Christoph Rüegg — Mon, 31 Mar 2008 09:41:00 GMT

Iteration 12 of Math.NET Iridium (Numerics) is now available, as 2008 March Release with Version 2008.3.12.405. Grab it here.

Team: Christoph Rüegg, Joannès Vermorel, Matthew Kitchin
New Contributions: Mike Shugai

Summary:

Bugs: All know 3 bugs have been fixed.
Linear Algera: Matrix Kronecker product, new Vector class, performance work (caching).
Core: New Sinc function, Neper/Decibel ratio routines.
Probability Distributions: New distributions: Student's T-Distribution, F-Distribution, Skew-Alpha Stable Distribution.
Assemblies no longer signed with a certificate (because verification caused network access).

New Features:

IRID-113: Core - Neper and Decibel Helper and Factors
IRID-116: Core - Sinc Function
IRID-111: Linear Algebra - Matrix Kronecker Tensor Product
RID-59: Linear Algebra - New Vector class (related to the Matrix class)
IRID-99: Probability Distributions - New Distribution: Skew Alpha Stable Distribution
IRID-108: Probability Distributions - New Distribution: F-Distribution
IRID-109: Probability Distributions - New Distribution: Student's-T

Enhancements:

IRID-106: Linear Algebra - Cache for on-demand computations (like decompositions)

Fixed Bugs:

IRID-107: Core - Complex: Unexpected power behavior at zero
IRID-97: Linear Algebra - Matrix.Identity allocation bug in non-square cases
IRID-98: Probability Distributions - ArbitraryDistribution NextInt32 does not consider offset.

Other:

IRID-110: Remove Certificate Signing (cert validation causes network access)

For more details have a look at the Iridium tracker change log.

Math.NET Tightened Security - working fine for you?

Christoph Rüegg — Sat, 09 Feb 2008 11:57:00 GMT

Since the February 2008 Release, Math.NET Iridium and Neodym both adapted a new security model:

Strong Name: The official assemblies are now strong named. This means that you can put them into the global assembly cache, configure .NET code access policies, and prevent anyone from injecting bad code into the assembly or replace it completely.
Delayed and Test Signing: Other than other opensource projects we don't make the private key of our strong name signing key publicly available. Instead we use delay and test signing for everyone who develops directly from the source code repository: the assemblies are delay signed with the official public key, but then are test-signed with individual keys you can easily create and register yourself (by running a batch file). If you have concerns about relaying on strong named assemblies where you don't have access to the private key, simply re-sign the assembly with your own key pair of choice...
Code Access Security: The assemblies don't run as full-trust anymore but instead are locked down and refuse nearly all system-provided permissions including any file I/O, network access, reflection and unmanaged code. With this we guarantee that we don't do anything bad on your system and that nobody can misuse the library to do bad things through this library. That's why we now also allow partially trusted callers - i.e. you can use the library even from a network share or directly from the internet, or from assemblies the user doesn't need trust fully.
Signed with a certificate: All the official assemblies are now also signed with a certificate to certify authenticy. However, this has no effect on using the library in your own project, and the root certificate used to create the code signing certificate is not generally known (and of course not distributed with Microsoft Windows).

This tightened security aspects might also have effects on your code. Please let me know if you have an issue with these changes, if it breaks your code, makes it much slower or whatever. Thanks.

Iridium 2008 February Release (v2008.2.10.364)

Christoph Rüegg — Sun, 03 Feb 2008 11:39:00 GMT

Iteration 10 of Math.NET Iridium (Numeric) is now available, as 2008 February Release with Version 2008.2.10.364. Grab it here.

This is mostly a service release:

Bugs: All known 8 bugs have been fixed.
Performance: The linear algebra implementation has ben optimized, resulting in nearly 50% perf gain.
Api References: Inline Xml Documentation has been improved (but is still far from where we'd like it to be).
Security: The released binaries now have a strong name, are locked down with code access security and allow partial trusted callers. We now use test-signing internally. Also, the official assemblies are now signed with a certificate.
Build/Release Integration: We now moved completely to custom msbuild targets, releases are now fully automated (incl. documentation generation) and the continuous integration system has been upgraded. Since releasing is now much easier, you can expect new releases more often.

For the complete list of changes have a look at the Iridium tracker change log.

At the same time I also released a first version of Math.NET Neodym (Signal Processing), Iteration 2 what makes it Version 2008.2.2.364. Grab it here. Hopefully I'll have more time in the future to work on Neodym...

Matrix data structure optimization: Nearly 50% perf gain

Christoph Rüegg — Sun, 14 Oct 2007 19:28:00 GMT

Some time ago I did some very basic performance analysis of solving linear equation systems with Iridium. Since then we decided to rewrite the Matrix class to use jagged arrays instead of rectangular ones, see this forum discussion. There already was a discussion about that issue some long time ago, but at that time we decided to go for the more clean and safe way of rectangular arrays. Unfortunately the C# compiler today still can't optimize loops on rectangular arrays as good as loops on jagged arrays. So we finally moved forward to jagged arrays, and indeed, we got a performance improvement (solving a linear equation system) by nearly 50%.

Unfortunately, the change of the data structure comes at a cost: The semantics of the two following members changes, as they now do deep-copies instead of using the data structure directly as internal data structure. Have a look at the mentioned forum discussion on why this might be an issue.

public Matrix(double[,] A)
public static implicit operator double[,] (Matrix m)

If you want to avoid deep-copying, e.g. for performance reasons, then use double[][] instead of double[,] to fill the matrix.

The changes are already submitted to the subversion repository, and will be included in the next iridium release.

Math.NET Iridium on Linux and Mono

Christoph Rüegg — Sun, 19 Aug 2007 16:43:00 GMT

One of the principles of the Math.NET Iridium numerics library is to not depend on special hardware or some external super-optimized library, but just on the core .Net Framework 2.0. Thanks to this, Iridium also works with Mono on Linux with no special treatment - just add a reference to the assembly in MonoDevelop and start coding...

Apparently some users are not used to reference .Net assemblies on Mono yet, so I've created some screen shots showing how to add such a reference step-by-step in MonoDevelop on Ubuntu Linux:

Math.NET on Mono (and Linux)

Iridium Performance Analysis

Christoph Rüegg — Sun, 15 Apr 2007 00:30:00 GMT

I started performing some very basic performance analysis of interesting parts of the Math.NET Iridium (Numerics) Library, motivated by a forum post.

Machine Configuration

I tested it on the following computers:

A: Desktop, Intel P4, 3.4 GHz, XP SP2
B: Notebook, Intel Core 2 Duo, 2 GHz, XP MCE SP2
C: Desktop, Intel Core 2 Duo, 2.?? GHz, XP SP2

On both machines there are lots of other processes running, hence I start the app with high priority.

Test Strategy

Quite simple (for a start): I run different problem sets. For each set I first run it two times (warm up, to let the JIT compiler do its work), then start a timer/stopwatch, run it five times, stop the timer and divide the resulting time span through five, to get an average.

Case 1: Fourier Transform

I tested the forward fourier transform on a real sample set of doubles with seven different sample set lengths:

Set Length:	A:	B:	C:
1024	0 ms	0 ms	0 ms
4096	0 ms	0 ms	0 ms
16384	1 ms	1 ms	1 ms
65536	6 ms	5 ms	5 ms
262144	33 ms	28 ms	30 ms
1048576	503 ms	578 ms	507 ms
2097152	1339 ms	1375 ms	1032 ms

Note that on B and C the algorithm only used one of the two cores. Looks like we need to tweak that a bit, since multi-core CPUs are standard these days.

Case 2: Solving a linear equation system

I tested solving a square real linear equation system. Internally the LU decomposition algorithm is used for such square systems. A Size of 1000x1000 means that we solve for 1000 unknowns:

Unknowns:	Matrix Elements:	A:	B:	C:
100	10000	2 ms	2 ms	2 ms
200	40000	18 ms	18 ms	17 ms
400	160000	129 ms	129 ms	121 ms
600	360000	415 ms	418 ms	396 ms
800	640000	1056 ms	968 ms	981 ms
1000	1000000	2041 ms	1911 ms	1901 ms

Again only one of the two cores on machine B and C was used.

Conclusion

I think the performance is quite acceptable for now. Of course it's far away from high performance libraries like the Intel MKL (blas, lapack etc. which take advantage of all specialized features of Intel CPUs), but in the end we target another kind of developer/application anyway. If you really want to get all out of your machine, you'll hardly do that in a managed programming environment like .Net. What we offer instead is a very easy to use infrastructure that is completely managed (no unsafe wrapper), and that runs on any platform that runs .Net (including your PDA) and that is still fast enough for most cases.

After all, only 0.006 seconds for a Fourier transform of 65536 samples is not that bad, neither is 0.018 seconds to solve a linear equation system with 200 unknowns.

Update: Added Machine C. The differences between the measured values on A,B and C don't make much sense, I think I need a better testbed and a larger sample size. Or there are simply other factors that play an important role (e.g. memory alignment -> cache behavior).

Iridium Class Reference

Christoph Rüegg — Fri, 06 Apr 2007 13:05:00 GMT

I've finally tried Sandcastle to automatically generate class references from the inline xml comments, as a replacement for the discontinued NDoc. The result for Math.NET Iridium is now available for download (in both HtmlHelp CHM and Ms Help 2 HxS format). Please note that because NDoc wasn't working anymore, inline XML documentation wasn't a top priority for the last years. I didn't have time to fix it yet, but thought it probably still helps more than none and thus made it available for download. I plan to provide such references for all future Iridium releases and expect its quality to improve clearly in the next release.

Yttrium Diagramming

Christoph Rüegg — Fri, 30 Mar 2007 13:44:00 GMT

There was quite a huge check-in lately on the yttrium source code repository. I've been refactoring the yttrium code base to a new microkernel based architecture (motivation: getting rid of cyclic dependencies and the high degree of complexity). The most striking change is the number of assemblies: Instead of a single dll the project is now split in multiple modules & services.

I'm also working on the diagramming frontend Whiteboard again. As you may have noticed, François is back working on Diagramming. Although he's finally moved forward towards WPF/.NET3 with his new Unfold project, he also released the older GDI+/.NET2 based Netron 3 code to the public (GPL; thanks, François!). I'm currently porting Whiteboard to the newly refactored yttrium codebase and extending the Netron code to better match the yttrium model. It looks great so far, keep in touch. Btw, I also include adjusted shapes for the new PetriNet demo package.

Iridium Release 8: 2007 March (v2007.3.8)

Christoph Rüegg — Wed, 14 Mar 2007 21:47:00 GMT

I've just released a new version of the C# Math.NET Numerics toolkit Iridium, v2007.3.8. Grab it here.

What's new in this release?

Completely revised FFT implementation
Integration of the Torschutz Random Library, revised and extended with pdf/cdf functions.
Various minor enhancements and fixes

I allowed myself to include some breaking changes; you'll have to update your code if you're working with either FFT transformations or with probability distributions and random generators. Nevertheless, I still recommend to upgrade, as I don't expect breaking changes in this extent (if any) in the next few releases.

Revised FFT Core for Iridium

Christoph Rüegg — Tue, 27 Feb 2007 23:41:00 GMT

I'm currently working on the Iridium (Math.NET Numerics) Fast Fourier Transformation Implementation, to finally fix known bugs and improve its usability (make it easier to work with). I expect a new Iridium release this week.

What's new:

Unit Tests: The code base is now tested (including MATLAB consistency and symmetry preservation).
Completely redesigned backend (e.g. removed the previous inheritance scheme).
Better support for real transformations. New member design, also for complex transformations.
Optional Conventions parameter. Controls scaling and integrand exponent behavior, e.g. to behave like Maple (symmetric scaling), MATLAB (asymmetric scaling) or Numeric Recipes (asymmetric scaling & inverse exponent sign).
Stateless* instead of stateful classes. You don't have to pass over the data on the constructor anymore, but as method parameters instead.
*except the precomputation cache. However, the classes are not static so you can get rid of the cache if you don't need it anymore.

The bad news: You'll have to slightly adapt your code if you work with Iridium's FFT implementation. However, the code will only become simpler...

PowerShell for Math.NET

Christoph Rüegg — Sat, 23 Sep 2006 10:40:00 GMT

I think it would be really interesting to propagate the new MS PowerShell as a new interface for Math.NET, especially for the Yttrium project. Yttrium already provides very simple script support, but PowerShell has so much more to offer. I think about writing some PowerShell cmdlets and functions, and maybe a provider for yttrium systems, and make this the foundation of any future interactive interfaces.

Let me know what you think about this approach!

Handling Floating Point Numbers

Christoph Rüegg — Wed, 16 Aug 2006 21:23:00 GMT

Did you know, that:

a double can represent 0xFFDFFFFFFFFFFFFF = 18437736874454810623 distinct numbers?
there are 450359962736 doubles between 1.0001 and 1.0002?
54975581 between 10'000.0001 and 10'000.0002?
52 between 10'000'000'000.0001 and 10'000'000'000.0002?
the smallest difference between neighbor numbers is 4.941E-324 for 0, and the biggest difference 1.996E+292 for double.MaxValue?
the difference between 1.0 and the biggest number smaller than 1.0 is 1.110E-16
(this value is often called negative epsilon (*2 = the (positive) epsilon) , but not in .Net! In the framework, double.Epsilon is defined as 4.941E-324, the number from above!)

It's common knowledge that floating point numbers are somewhat tricky to work with, especially if rounding-off errors and exact values become important. Although this is (hopefully) discussed in any computer science and engineering study course, there has been quite some blogging about it lately:

Floating-Point Numbers (Wesner Moise)
.999... = 1 (Wesner Moise)
The trouble with rounding floating point numbers (Dan Clarke)
Quiz of the Month: Double Trouble (Kathy Kam)
Why you can't Double.Parse(double.MaxValue.ToStr... (Scott Hanselman)
Floating Point Arithmetic (Wesner Moise)
Numbers in .NET (Wesner Moise)
Binary floating point and .NET
Comparing floating point numbers (Bruce Dawson)
"What Every Computer Scientist Should Know About Floating-Point Arithmetic" (PDF)

One of the conclusions is that it's dangerous to compare floating point numbers directly. Of course this is a problem in Math.NET Iridium, so I tried to add some helper functions to Iridium for handling double numbers. You'll find them in the new static Number class directly in the MathNet.Numerics namespace. The most important function, AlmostEquals, is inspired by Bruce Dawson's article (see above). They're propositions, so please let me know if you think they're flawed...

Integer Representation of Doubles

In order to work with doubles in an exact way, we have to treat them as integers. We can convert them either by doing some unsafe pointer casting, by emulating a union (StructLayout and FieldOffset attributes), or simply by using the BitConverter (which uses the first pointer approach internally; System.Runtime.InteropServices namespace).

Unfortunately, floating point numbers are stored as signed magnitudes, while integers use the two's complement to represent negative numbers. We thus have to convert negative numbers if we need them to be lexicographically ordered (such that the distance from 0 to 1 equals to the distance from -1 to 0, and that -0 = +0):

public static long SignedMagnitudeToTwosComplementInt64(long value)
{
  return (value >= 0) ? value : (long)(0x8000000000000000 - (ulong)value);
}

public static ulong ToLexicographicalOrderedUInt64(double value)
{
  long signed64 = BitConverter.DoubleToInt64Bits(value);
  ulong unsigned64 = (ulong)signed64;
  return (signed64 >= 0) ? unsigned64 : 0x8000000000000000 - unsigned64;
}

The expected behavior (excerpt from the unit tests):

ToLexicographicalOrderedUInt64(2 * double.Epsilon) == 2
ToLexicographicalOrderedUInt64(1 * double.Epsilon) == 1
ToLexicographicalOrderedUInt64(0.0) == 0
ToLexicographicalOrderedUInt64(-1 * double.Epsilon) == 0xFFFFFFFFFFFFFFFF

SignedMagnitudeToTwosComplementInt64(1) == 1
SignedMagnitudeToTwosComplementInt64(0) == 0
SignedMagnitudeToTwosComplementInt64(-9223372036854775808) == 0
SignedMagnitudeToTwosComplementInt64(-9223372036854775808 + 1) == -1

Incrementing Doubles

Thanks to this two's complement integer representation we immediately have the power to increment or decrement floating point numbers, and iterate through all numbers a double can represent in a given range:

public static double Increment(double value)
{
  if(double.IsInfinity(value) || double.IsNaN(value))
    return value;

  long signed64 = BitConverter.DoubleToInt64Bits(value);
  if(signed64 < 0)
    signed64--;
  else
    signed64++;
  if(signed64 == -9223372036854775808) // = "-0", make it "+0"
    return 0;

  value = BitConverter.Int64BitsToDouble(signed64);
  return double.IsNaN(value) ? double.NaN : value;
}

Again an excerpt of the expected behavior:

double x = double.Epsilon;
x = Decrement(x);
--> x == 0
x = Decrement(x);
--> x == -double.Epsilon
x = Increment(x);
--> x == 0
x = Increment(x);
x = Increment(x);
--> x == 2 * double.Epsilon
x = double.MaxValue;
x = Decrement(x);
--> x < double.MaxValue
x = Increment(x);
--> x == double.MaxValue

Evaluating the Precision: Epsilon

The incrementation step depends on the exponent, so it may differ between different numbers. It is important for example in iterative approximation algorithms to know when to stop. The following function returns the decrementation step near a given number, often referred to as a scaled negative epsilon:

public static double EpsilonOf(double value)
{
  if(double.IsInfinity(value) || double.IsNaN(value))
   return double.NaN;

  long signed64 = BitConverter.DoubleToInt64Bits(value);
  if(signed64 == 0)
  {
    signed64++;
    return BitConverter.Int64BitsToDouble(signed64) - value;
  }
  else if(signed64-- < 0)
    return BitConverter.Int64BitsToDouble(signed64) - value;
  else
    return value - BitConverter.Int64BitsToDouble(signed64);
}

And some examples:

EpsilonOf(1.0).ToString() == "1.11022302462516E-16"
EpsilonOf(0.0).ToString() == "4.94065645841247E-324"
EpsilonOf(-1.0e+100).ToString() == "1.94266889222573E+84"
EpsilonOf(1.0e-100).ToString() == "1.26897091865782E-116"
EpsilonOf(double.MinValue).ToString() == "1.99584030953472E+292"

Numbers between two doubles

It may be interesting to know how many numbers can be represented by a double variable, which are all between two given numbers. Using the methods introduced above this is simple:

public static ulong NumbersBetween(double a, double b)
{
  // left out (see repository): handling NaN and Infinity 

  ulong ua = ToLexicographicalOrderedUInt64(a);
  ulong ub = ToLexicographicalOrderedUInt64(b);

  return (a >= b) ? ua - ub : ub - ua;
}

Examples:

NumbersBetween(1.0, 1.0) == 0
NumbersBetween(0, double.Epsilon) == 1
NumbersBetween(-double.Epsilon, 2 * double.Epsilon) == 3
double test = Math.PI * 1e+150;
NumbersBetween(test, test + 10 * Number.EpsilonOf(test) == 10
NumbersBetween(1.0001, 1.0002) == 450359962737
NumbersBetween(10000000000.0001, 10000000000.0002) == 53
NumbersBetween(double.MinValue, double.MaxValue) == 18437736874454810622

Compare doubles for (almost) equality

Now we come to the final but most important method, for a safe and more sensible way of comparing two floating point numbers for almost equality. Call this method with two doubles and a parameter specifying how many numbers may be between the two numbers at most (+1):

public static bool AlmostEqual(double a, double b, int maxNumbersBetween)
{
  // left out (see repository): argument checks

  // NaN's should never equal to anything
  if(double.IsNaN(a) || double.IsNaN(b)) //(a != a || b != b)
    return false;
  if(a == b)
    return true;

  // false, if only one of them is infinity or they differ on the infinity sign
  if(double.IsInfinity(a) || double.IsInfinity(b))
    return false;

  ulong between = NumbersBetween(a, b);
  return between <= (uint)maxNumbersBetween;
}

Examples:

Convert.ToDouble("3.170404") == 3.170404 --> TRUE
Convert.ToDouble("4.170404") == 4.170404 --> FALSE
Number.AlmostEqual(Convert.ToDouble("3.170404"), 3.170404, 0) --> TRUE
Number.AlmostEqual(Convert.ToDouble("4.170404"), 4.170404, 0) --> FALSE
Number.AlmostEqual(Convert.ToDouble("4.170404"), 4.170404, 1) --> TRUE
Number.AlmostEqual(0.0, 0.0 + double.Epsilon, 0) --> FALSE
Number.AlmostEqual(0.0, 0.0 + double.Epsilon, 1) --> TRUE
double max = double.MaxValue;
Number.AlmostEqual(max, max - 2 * Number.EpsilonOf(max), 0) --> FALSE
Number.AlmostEqual(max, max - 2 * Number.EpsilonOf(max), 1) --> FALSE
Number.AlmostEqual(max, max - 2 * Number.EpsilonOf(max), 2) --> TRUE

New Math.NET Website Organization

Christoph Rüegg — Wed, 16 Aug 2006 11:38:00 GMT

In the last few years the Math.NET project websites became scattered around web, including SourceForge and my personal web site. That's why I've been reorganizing the websites the last weeks (and continue to do so the next days).

Our new central website, replacing SourceForge and my personal website as entry point for the project:

http://mathnet.opensourcedotnet.info/

Every other project pages on the web are (or will be) linked from this website, mainly:

All the information of the old websites (including http://nmath.sourceforge.net and http://www.cdrnet.net/projects/nmath/) will be integrated into the new website or one of the resources listed above. Please consider updating your links.