One size fits all

As an amateur chess player I occasionally watch videos on the subject on youtube. Just as it happens I watched one old video of GM Ben Finegold on isolated queen pawn positions the other day. He started with explaining that most people don't like such positions, because beginners tend to be dogmatic. They are taught a number of rules like ‘isolated pawns are bad’, ‘knights before bishops’ and so on, and they like to stick to the rules without regard to context. So that is when it struck me, we tend do the same in programming, especially true in the case of people with limited experience. That rigid thinking however doesn't do well neither in chess nor in C++.

My favorite object to rant about these days is the beloved by beginners and advanced alike “design patterns” pattern. This was all started by the famous book by the same name written by the so called “gang of four”^[1]. I, being pretty active the last year in the Qt forums, often encounter questions on the topic of design patterns from other users. So I decided to write this very post, where I hope the shed some light on the issue.

As programmers we often encounter similar problems as work progresses along which are solved pretty much the same way, and as humans we love to generalize. So that is the whole premise behind the aforementioned book, to show that a set of similar problems have a corresponding set of “good” solutions. This, of course, is generally true, but many people, especially the ones well conditioned to rigid thinking and/or those with limited experience, are ready to jump overboard and just apply it to every situation, always. So in reality this creates a very profound problem, not only are some patterns inapplicable to certain contexts, but people tend to bend over backwards in their wish to apply them, thus creating a vile mess in the process.

The original ideas were posed for Java, and while it pains me to say it constantly, C++ is not Java. Some of the proposed solutions are inappropriate to be applied directly to C++, or become much more complex if one tries to adapt them. This naturally stems from the fact that C++ is a lower-level language, albeit still object oriented, and as such has a lot of peculiarities that are not relevant to other languages.

A typical example of misunderstanding is the singleton pattern, where only one instance of an object is created, initialized and used. This pattern (or as I lovingly call it: “The singleton antipattern”) has gained widespread use for no good reason. The idea behind it is supposedly to restrict the programmer to using a single instance of a given class. There are a number of pitfalls with this, however:

1. Having a singleton class imposes that you have only one object of that class, conversely needing one object does not sum up to requiring the class to be singleton(ian). If you need one instance of something, just create one!
2. The singleton creates extremely tight coupling between the classes that use it. Think a C global variable. In fact the singleton is just that – a global variable dressed in shiny clothes (i.e. it has methods). The same results one can get by using a number of global functions that make use of a global variable.
3. As a consequence of 1 and 2, the singleton introduces an application global state, which might be hard to manage depending on the context.
4. A C++ specific problem is that most singleton implementations leak memory. Java has no such problem as the memory is managed by its virtual machine, in C++ this may be problematic.
5. A lazy initialization singleton implementation must be explicitly made thread-safe in case the initialization might be done concurrently.

So what the singleton has to show for in return – imposing a single instance.

Other patterns described in the book and on the internet also suffer similarly and the reason is: there is no one-size-fits-all in programming.
Don't be dogmatic, evaluate the context level-headedly and only then settle on a solution.

---

[1] Gamma, E., Helm, R., Johnson, R., Vlissides, J., Design Patterns - Elements of Reusable Object-Oriented Software. Addison-Wesley, 1995.

OpenMPI and qmake

For a project I'm developing I want to use OpenMPI and the Qt toolkit. Naturally, this would mean that I'd like to use QtCreator with the qmake's build chain as well. OpenMPI recommends using their wrapper compiler/linker instead of the default one, so I needed to make aware the qmake system of this. It is possible to change the compiler on a per-project basis, but after reading a bit on qmake, I've decided to solve this a bit more elegantly.

qmake supports features that can be used in different projects by just adding a configuration entry. I've created such a feature and added it as a git repository here: https://bitbucket.org/nye/qmake-openmpi-feature.git. To make qmake aware of a feature one has only to put the .prf file in $QTDIR/mkspecs/features and it's ready to use. To use the OpenMPI compiler/linker just add the following configuration entry to your project file:

CONFIG += openmpi

Count your pluses

This text is a continuation of the “Bubblin' up” post regarding floating point operations, and more specifically here I will explore how to properly sum such numbers. The problem presented is, that because of the specific representation of floating point numbers, every operation done on them is inherently imprecise. Additionally, even if it is assumed that the operations carry no error at all, there is a problem with storing the number. If the number can't be exactly represented in the specified floating point format (as discussed in the previous post), naturally a truncation error will be imposed. What does that mean? Firstly, do not trust floating point numbers. Yes, the statement is a bit exaggerated, but it is a good practice to always keep in mind that floating point numbers are just an approximation. And secondly, special care should be taken when working with them.

Truncation, as unseemly as it may appear, is inevitable because of the discreteness of computer memory. Nothing can be done about it, there is no computer with an infinite amount of memory. Here the machine epsilon arises as a natural measurement of the least significant number (in sense of mantissa bits) representable in a specific architecture. Rounding errors are a consequence of truncation and occur because no number can be represented in computers with an infinite amount of significant digits, meaning all floating point operations are approximate and need be rounded, hence the name. By design the machine epsilon is the upper bound on the relative error associated with rounding in floating point arithmetic.

All that said, summation is a classical problem illustrating the accumulation of error. Every addition will carry an error into the result, and even if the error is not large, the accumulation might be significant when adding many terms. Furthermore, because computers normalize the numbers (see wikipedia), if the values' exponents differ much, the smaller numbers' contributions might not be reflected in the sum at all! Although the problem is known, the solution is not as obvious or trivial as it might seem. So lets consider and compare few approaches (error estimations will not be derived, but only presented):

1. The naïve approach consists of just adding up the numbers in an accumulator. The accumulated error will grow linearly on the number of terms – O(n). A snippet in C++ illustrating the method is presented:

double doNaiveSum(double * numbers, int size)
{
    double sum = 0;
    for (int i = 0; i < size; i++)
        sum += numbers[i];
    return sum;
}

2. The naïve approach with presorting in which the numbers are sorted by absolute value before summing them. This case is very similar to the previous one, so no code snippet will be provided. The idea is to sum smaller numbers first so their contribution is not lost. The performance is only marginally better than the former case and the benefit is not guaranteed but depends heavily on the magnitude of the numbers, which might vary to a large degree in different applications.

3. Pairwise summation is a divide and conquer algorithm which consists in dividing each range by half, summing recursively each half and then adding the results. The method offers significant improvement over the naïve approach. The error grows logarithmically O(log(n)) for a minute increase in arithmetic operations. A minor inconvenience is that the natural implementation is recursive. Consider the following example code:

double doPairwiseSum(double * numbers, int size)
{
    if (size < 4)  {  // For 3 numbers or less just perform the regular summation
        double sum = 0;
        for (int i = 0; i < size; i++)
            sum += numbers[i];
        return sum;
    }

    // Calculate the left part size as half the original size, and the right part size as the reminder
    int leftSize = (size >> 1), rightSize = size - leftSize;

    // Perform a recursive sum on each part and return the result
    return doPairwiseSum(numbers, leftSize) + doPairwiseSum(numbers + leftSize, rightSize);
}

4. Kahan summation is a compensated summation technique which uses a running compensation to reduce the numerical error. The idea is to keep an additional variable, where the low-order bits, which would be otherwise lost, can be stored. Although Kahan's algorithm achieves O(1) error growth (a constant not depending on the number of additions), additional operations are required which might not always be appropriate. In such a case pairwise summations is a viable alternative with its logarithmic error growth.

double doCompensatedSum(double * numbers, int size)
{  
    double sum = 0, compensation = 0;
    for (int i = 0; i < size; i++)  {
        // Adjust the input value with the stored compensation value (the low-order bits from previous operations)
        double x = numbers[i] - compensation;
        // Calculate the new value for the sum
        double tempSum = sum + x;
        // Recover the low-order bits lost in the last addition
        compensation = (tempSum - sum) - x;
        // Update the sum variable
        sum = tempSum;
    }
    return sum;
}

When using this algorithm beware of overly aggressive compilers, which could in principle recognize that arithmetically the compensation should always be zero and perform optimizations effectively removing its calculation.

5. Arbitrary precision arithmetic is a brute-force solution where no special care is (implicitly) taken to limit the error growth. For most intents and purposes one of the above (or similar) algorithms would be more appropriate, since arbitrary precision floating points are usually emulated in software (at least to my knowledge) and are carrying a lot of computational and memory overhead.

Bubblin' up

There is a disturbing amount of people willing to compare floating point numbers directly for equality. For most intents and purposes this is harmless, but if one talks about comparing physical quantities it can be a problem. The fact that many people ignore is that real numbers can not be represented accurately in computers. Irrational numbers have infinite number of terms after the decimal points and whatever we do, it is necessary to truncate somewhere to fit a number in the finite memory of our computers.

Another more often acknowledged problem is the dynamic range of physical quantities. While the mass of cars, people or other everyday objects are in the comfortable zone of kilograms up to tonnes, the mass of stars are way beyond these magnitudes, and to make matters worse the mass of particles are on the other far end of the scale. It is obvious that some sort of normalization should be employed to bring all that dynamic range into manageable chunks. Naturally this is done by simply scaling the quantities and is not such a problem. Still it is necessary to be able to represent both large and small real numbers, so here the floating point numbers come in play. The structure of a floating point number is pretty simple – there is a fixed point number between zero and one which is multiplied by an exponent (for computers, conveniently, it is a power of 2):

Where m is the mantissa, and p is the power (sl. exponent).

It is pretty obvious that truncating the fixed point number (the mantissa) and the power of the exponent imposes some constraints on the representation of the numbers:

• The dynamic range of the numbers is finite because there is no way to represent any arbitrary integer for the exponent.
• The precision is finite because one can't represent a fixed point number with infinite number of terms after the decimal dot.

Beside disadvantages this representation has its strong sides, mainly convenient implementation and good dynamic range coupled with appropriate precision. It is noteworthy mentioning that while the relative precision is kept, the absolute valued distance between the representable numbers increases with increase of the exponent.

All that said I will go back to the original problem – how to check if floating point numbers are equal?
It should be obvious by now that the imprecise representation excludes just naïvely comparing them with the equality operator and while I am not proponent of the epsilon-delta formalism in mathematics here it is quite appropriate. The idea is to set an epsilon representing the maximum relative difference between the two numbers and use it as a threshold.

Going dynamic

The Fortran people seem to really like the static libraries, why I couldn't really fathom. The idea behind a shared object (or dynamic library in Windows) is quite simple — flexibility. Their explanation is that static libraries (basically an object file) provides all that you need for linking and doesn't cause the problems loaders do. This is not quite true and there are some problems though:

• If you link against a shared binary in the static library code, the linker has no way of incorporating it, thus creating references exactly like these you have in a dynamic library.
• When you want to change something in a statically linked code, you have to rebuild the whole thing.
• (Especially in Linux/Unix) The code is really shared(!), meaning a lot of applications use the same binary, changes in that binary propagate to all and additionally there is only one copy of each one binary unit.

Now dynamic linking is not without cost:

• A person should ensure binary compatibility of the library as to be able to change the implementation without rebuilding the user applications. This could be quite a bit of work, especially for people who don't have a good set of habits for doing it.
• You should ensure all the dependencies are in a place known to the loader.
• Resolving of symbols may be done by the runtime as it is in Linux/UNIX.
• The library interface and exports should be enforced, as in Linux/UNIX by default all symbols are exported, while in Windows none of them are.

All in all dynamic libraries require a bit more consistent and rigid design, while static libraries could be less fuss. In the long rung though, the advantages of dynamic linking are overwhelming and it is the better way, allowing you to be flexible about your code.

A Gaussian is a Gaussian

Interestingly a sum of Gaussian distributions is again Gaussian distribution. It didn't came as a big surprise to me, but my thesis adviser was truly amazed ... strange. Furthermore it holds true for both uncorrelated, and correlated statistics, and the rule for addition is impressively simple.

So let's say that we have two Gaussian distributions:

Then the sum distribution is simply:

Where for the parameters holds:

In the last formula rho is the correlation coefficient, with value defined as:

Exceptions

Before writing the final post on Wiener deconvolution, I decided to share a little code snippet. It occurred to me that from time to time I write an exception class. Mostly they are the same, so why not make it public ... It can be derived of course to provide further specialization, but for the current project I'm working on (which is a console application) it's enough.

So here goes the code:
http://codepad.org/KGvqarqa

Filters and such

This is a continuation of the "Delta spectrum" post, so we start from there. I've realized that if not done incrementally, to obtain a convolution of two signals in the time domain, would be more computationally consuming, than to do it in the frequency domain.

FFT is with complexity, against the complexity of time domain convolution. So I'm going to explore this a little bit.

We start with a signal defined as:

Here H is the frequency response of our system, N is additive noise, X is the real (original signal) and Y is the observed, convoluted signal with noise. We would like to find a function G such as:

Whereas is an estimate of the original signal, which should minimize the mean square error.

Without going into further details, the operation is assumed to be carried out in the frequency domain as follows:

Taking into consideration that G is the Wiener filter defined as:

Now going back and accounting for our original signal and our frequency response, which are correspondingly a delta function and a Gaussian, the last equation is simplified enormously. For the impulse response we have:

And for the original signal one can write:

Where j is the imaginary unit.

The noise could be estimated by doing lowpass filtering and converging to a state where the mean square deviation from the original signal would be minimal. The parameters alpha and tau can be automatically estimated from the signal peaks, by doing nonlinear Gaussian fit in the time domain. The whole point of doing such a complicated process would be to try finding merged energy peaks, low intensity peaks, or estimating Doppler shift in the spectrum.

Read me, write me

Recently I've came up with a timeless problem ... imagine what happens if you try to read from and write to a file simultaneously - junk. Let's say that you are in a single process space, you have the option to synchronize threads and go around the problem, but for more than one process it becomes tricky. From what I can tell there are two options, firstly one can devise a scheme with interprocess communication (IPC) to serialize the read-writes. Unfortunately if you're aiming for platform independence, it won't work, not that it's impossible, but it would take quite an effort to do it. You'd have to implement the IPC for every platform you are going to support.

I've thought about this some time now, and I've came up concluding that there is simpler, easier and portable way to do it - file server. It's not a classical ftp type server per se, but more of a simplistic network oriented file management server. One would create a daemon which listens to some port, and use a simple protocol to communicate with it. The daemon itself would be able to list/create/delete/read/write files on the host machine. Aside from a small overhead, it can work on the local computer (through local port), as well as on a remote machine, which gives it the flexibility I'm looking for. It's not something new in the world of programming, but gives a way to implement file operations serialization primitive which is portable.

The purpose of all this, is a project I have a plan on starting in the near future. And since it would need both network client-server connectivity and local file locking primitive, I've decided to merge the two concepts and use serialization for both cases.

Delta spectrum

Recently we were discussing with Stanimir (another colleague of mine) the skewed Gaussian function. It is pretty obvious that if you want to know what the parameters are, you could perform a nonlinear least squares fit to acquire them. That's okay, but what if you can't spare the computational power, like in my case. Then I started searching for alternatives ...

I was aware that skewed Gaussian you can get by convolving an exponent, Heaviside function and a Gaussian ... so I imagined that probably you would be able to separate them by the reverse process, namely deconvolution. I was thinking about one particular project then, but just now, I realized ... what would happen if I apply the same idea to spectra from nuclear experiments ...

It is pretty well known that scintillation detectors are not so big on resolution, they have nice wide peaks, which accidentally are Gaussian functions. So my mind went on wandering ... if I assume that peaks are in reality delta functions, and they are convolved by a Gaussian, generated from the acquisition system, we would get what we see in the spectrum. But then ... if one applies a deconvolution method in some manner, in theory it is possible to separate every single peak ...

Here it goes a little bit of mathematical motivation, but the idea is sound, at least for now. When I have the time I will implement it, and report the result.

Lets say that you have two functions, their convolution is defined by:

This is just a definition of a functional, and it has some cool properties, but most importantly the following two:

The first is just identity equation, where delta is the delta function (no surprise here). The second one is called the convolution theorem, Where F{} is the Fourier transform.

To keep this short I'll only provide a guideline to the remainer of the method proposed. Let's say we have a system of the type:

Where G is the Gaussian function applied to our spectrum, and n is just an additive noise (uncorrelated to the original signal). I would assume that the noise has Gaussian distribution, which is mostly correct. Knowing that Gaussian transfered to the frequency domain is again Gaussian, means it wouldn't be hard to devise a filter to remove all other components, but the delta functions ... and that is exactly what my target is.

The filter itself may not be in the frequency domain, which is the point of the whole post. It can be created in such a manner, that it's just applied to the spectrum to enhance resolution. The idea can be implemented through the Wiener deconvolution process.

A am far from the idea, that one could get real delta functions, but decrease in the dispersion of existing Gaussian peaks, would be more than enough. So far this is just my imagination going wild, but still it has some potential.

Spectrum of irritation

Usually when I get a nuclear spectrum file, it is just a sequence of numbers separated by newlines. As far as I'm aware, radware knows nothing about such a simple format. So instead of constantly being irritated, I decided to write a simple program to convert such files in radware compatible format. Probably it is not very original, but I'll present the source code anyway. It uses the Qt library, and is pretty much self explanatory.

Here is a link to the program mentioned:
http://codepad.org/kH8Iql2E

A change of heart

For some time now I use a nice transformation for double integrals. Okay in my case they are not exactly double, but 6th order integrals, nonetheless it is the same. The problem I'm trying to solve is a double integral over the radius vectors of two particles, so it is 3rd order by the coordinates of a single particle (I have two particles). You would ask why I need to do that ... well I have to simplify otherwise very heavy calculation. This double volume integral is calculated for permutations (which are quite a lot) of the quantum numbers of the particle system. For convenience, because of the problem symmetries, I work in cylindrical coordinate system. Up until now I've managed to separate the linear part, which is pretty straightforward, and was able to factorize the double integration through the transformation of coordinates mentioned.

Here it goes the simple but elegant method:

Say we have a double volume integral over two variables x and y. And they are used as sums and differences throughout the function we are integrating. It is obvious we would want to separate them in some manner, so we would calculate two one dimensional integrals, instead of one two dimensional. It's not hard to imagine that:

would do the trick. But there is also a kick in all that. I've realized yesterday that this is nothing more than rotating the integration area ... fun huh?

This transformation being a orthogonal one (rotation at fixed angle) means it has a constant Jacobian. It even gets better. If we normalize the new variables (divide by square root of 2) the change preserves volume:

This last change of variables is called Moshinsky transformation. There is a little bit more on the subject but goes out of the scope of this post, maybe someday I'll write down the Gogny separation method as well.

I had a discussion with Maya

Recently I had a discussion with Maya (a colleague of mine from the university). I tried to convince her, and I hope I did, about the sweetness and bitterness of software. Okay we have some software for nuclear spectra analysis, but mostly in any lab I've been there is a different set of programs used. It really would be great if there was a unified software. Furthermore I really believe it is possible to create such a tool. If I was not busy I'd even start to create it (opensource of course) ... but for the time being it's just dreams.

I'd still describe my idea though. The proposed software would consist of server and client (terminal). The terminal would be a graphical program where scientists would be able to see what's going on over the server. The server would be just a daemon connected to the detectors and all that complicated stuff we use to collect data for analysis. The data would be collected in a SINGLE unified file format on the server, and would be accessible through the terminal program (including exporting/importing). More or less it should work as the X server. The application would be very modular, organized in plugins which are binary compatible, not some monolithic stuff which wouldn't compile at all, if something is not right with a single function (which is probably not used anyway). To be honest, there are some programs doing similar things, but I don't like them ... and I'll even give an argument why.

Let's take Midas - it has only two problems ... it is ugly, which is manageable, and it is slow, which is not. In fact it is so slow that I'd rather kill myself (not literary) or destroy the computer I'm working on (quite literary), instead of using it. But sometimes one doesn't have a choice ... unfortunately. I don't know if it can be run on both Windows and Linux, but with that speed it doesn't really matter.

Okay next on the block is radware - it is quite good program, with one exception ... getting it to build and run is a pain, and it only works on Linux (naturally?!). Aside from these two technical difficulties, it is really good for processing spectra. Ah, and I almost forgot ... it has terrible learning curve ... you have a lot of different files for different things and bad support for external file formats. No good option to export your spectra in graphic files is also a drawback.

In our labs we have a local program - ANL - which is just an antique, it runs only under DOS and in practice cannot do much.

In HIL, Poland (where I was an Erasmus student) there is another local program widely used - SMAN. It has quite a lot of capabilities, including scripting and drawing two dimensional spectra. It's not bad really, but I saw a little problem in it's monolithic nature. Also it runs under DOS (come on, who uses DOS anyway?!). I doubt it runs under Linux at all, although I'm not sure, it may. Another thing I find inappropriate is that it's written in pure C ...
I really don't get it. Why people still use pure C, it is old stuff, we have classes, and interfaces, and tools to create new software which is bug free (for the most part) and easy to maintain. If I didn't create that program how would I ever find my way in its millions lines of source code?! The guy who created that software even tried (unsuccessfully) to convince me that C++ is no good ... well it is way better than C, but hey, he is entitled to have an opinion.

I won't talk about root, since I haven't used it. But from what I hear ... it's not very convenient, or easy to learn.

For stuff written in Fortran I won't even start. This language is just stuck in the 60's. I can relate that physicists are not programmers, but Fortran ... really?! It is just archaic, with bad syntax, a lot of reserved words, a lot of implicit conventions (like variable types depending on name?!) ... do I have to go on ...

In conclusion, you can see my frustration - just to perform some experiment you'd have to know different programs, use different scripting languages. This is the reality, but I'm pretty sure that the time is right to create something good for the community, which has steep learning curve, and is fast enough to be used in real experiments.