Posts Tagged ‘C++’

Don’t do this!

Monday, September 21st, 2009

Don’t ever do this:

try {
   ...
} catch (std::exception e) {
   throw e;
}

It doesn’t do what you think it does!

Unless you think it construct and then deconstructs an exception and throws away the type information about the original exception, in which case it does exactly what you think it does.

Of course, that might be the whole point of this piece of code, but why anyone would want to throw away information about which exception was thrown is beyond me.

Catching by value just means that you are creating a new object, but an exception is already copied when thrown, so there is no point.  You save the copying by catching by reference.

Throwing the explicit copy means you slice away the dynamic type, so the (re-)thrown exception has a different type than the original exception, so you have no way of knowing which kind of exception was originally thrown.

Catching by reference doesn’t help you here, since “throw e” when “e” has a type means you throw with the static type of “e”, not the dynamic type, so you will still slice away the type.

Consider this small example:

#include <iostream>

struct Exception {};
struct SpecialException : public Exception {};

void foo() { throw SpecialException(); }

void bar() {
  try { foo(); }
  catch (Exception e) { throw e; }
}

void baz() {
  try { foo(); }
  catch (Exception &e) { throw; }
}

int main() {
  try { foo(); }
  catch (SpecialException &e) {
    std::cout << "foo threw special" << std::endl;
  } catch (Exception &e) {
    std::cout << "foo threw plain" << std::endl;
  }

  try { bar(); }
  catch (SpecialException &e) {
    std::cout << "bar threw special" << std::endl;
  } catch (Exception &e) {
    std::cout << "bar threw plain" << std::endl;
  }

  try { baz(); }
  catch (SpecialException &e) {
    std::cout << "baz threw special" << std::endl;
  } catch (Exception &e) {
    std::cout << "baz threw plain" << std::endl;
  }

  return 0;
}

Here “bar” will throw away the type information, while “baz” will not.

The result is this:

$ ./foo
foo threw special
bar threw plain
baz threw special

In the original example, I’m also a bit puzzled as to the reasoning behind it.  Not catching the exception in the first place seems to me a better approach than catching it and just rethrowing it, without any further processing.  The code really only contributes by adding a bit of overhead and removing useful type information.

PS. Yes, I stumbled over code exactly like this in a library I’m working with, which prompted this post…

264-297=-33

Problems with boost::interval

Saturday, May 23rd, 2009

I need to work with exponentiation and logarithms with boost interval arithmetic, but nothing seems to work for me.

Here’s a simple little test program:

#include <iostream>
#include <boost/numeric/interval.hpp>
using namespace boost::numeric;
using namespace interval_lib;

int main()
{
  interval<double> a = 1, b = 1;
  std::cout << width(exp(a + b)) << std::endl;
  return 0;
}

If I try to compile this,I get the errors:

/usr/include/boost/numeric/interval/transc.hpp: In function ‘boost::numeric::interval<T, Policies> boost::numeric::exp(const boost::numeric::interval<T, Policies>&) [with T = double, Policies = boost::numeric::interval_lib::policies<boost::numeric::interval_lib::rounded_math<double>, boost::numeric::interval_lib::checking_strict<double> >]‘:
foo.cc:9:   instantiated from here
/usr/include/boost/numeric/interval/transc.hpp:34: error: ‘struct boost::numeric::interval_lib::rounded_math<double>’ has no member named ‘exp_down’
/usr/include/boost/numeric/interval/transc.hpp:34: error: ‘struct boost::numeric::interval_lib::rounded_math<double>’ has no member named ‘exp_up’

Why aren’t exp_down and exp_up defined?  According to the documentation, they should be defined for the built in types:

By default, it is interval_lib::rounded_math<T>. The class interval_lib::rounded_math is already specialized for the standard floating types (float , double and long double). So if the base type of your intervals is not one of these, a good solution would probably be to provide a specialization of this class. But if the default specialization of rounded_math<T> for float, double, or long double is not what you seek, or you do not want to specialize interval_lib::rounded_math<T> (say because you prefer to work in your own namespace) you can also define your own rounding policy and pass it directly to interval_lib::policies.

This is really annoying me, ’cause I have no idea what I am doing wrong or how to fix it.  No idea whatsoever, and I cannot continue with my project before figuring it out.

Sometimes, the boost library is just a little too complex for me, I guess…

142-152=-10

Linear algebra in C++ is no fun

Tuesday, April 28th, 2009

One of the reasons I’m in the UK this week is to work with David Balding on HapCluster, a method developed in his group that I have extended in various directions.

There’s a few extensions we want to add, and maybe (finally) get a paper written about all the extensions.  The latter probably won’t happen … even with lots of extensions is is hard getting a paper accepted that just adds features here and there, and it is always necessary to compare with a lot of other tools (when really no single tool can beat all the others unless you find just the right “sweet spot” in parameter space where you are lucky, and that takes months of fiddling with parameters and I just can’t be bothered).

Anyway, one of the extensions we are looking at is mapping of quantitative traits.  Ed Waldron already did this in his thesis, but I never got around to implementing it in my tool, ’cause I didn’t fully understand the mathematics.  I still don’t, but I think that I get enough to implement the formulas now at least.

The only problem is that it is a lot of matrix expressions that I need to implement, and that is just hell to do in C++.  At least if you, as I do, use GSL and CBLAS.

A simple expression such as

\mathbf{V}^* := \left(\mathbf{V}^{-1} + \mathbf{B}^T\mathbf{B}\right)^{-1}

becomes

gsl_matrix *
calc_Vstar(const gsl_matrix *B, const gsl_matrix *invV)
{
	gsl_matrix *copy = gsl_matrix_alloc(invV->size1, invV->size2);
	gsl_matrix_memcpy(copy, invV);
	gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1.0, B, B, 1.0, copy);
	gsl_matrix *Vstar = invert_matrix(copy);
	gsl_matrix_free(copy);
	return Vstar;
}

with a pre-computed \mathbf{V}^{-1}.

Computing

\mathbf{m}^* := \left(\mathbf{V}^{-1} + \mathbf{B}^T\mathbf{B}\right)^{-1}\left(\mathbf{V}^{-1}\mathbf{m}+\mathbf{B}^T\mathbf{Y}\right)

we can re-use the computation from before, but it still ends up as this long function:

gsl_matrix *
calc_mstar(const gsl_matrix *Vstar, const gsl_matrix *invV,
                const gsl_matrix *B, const gsl_matrix *m,
                const gsl_matrix *y)
{
	assert(Vstar->size1 == Vstar->size2);
	int p = Vstar->size1;

	// x := V^-1 m
	gsl_matrix *x = gsl_matrix_alloc(p, 1);
	gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, invV, m, 0.0, x);

	// z := B'Y
	gsl_matrix *z = gsl_matrix_alloc(p, 1);
	gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1.0, B, y, 0.0, z);

	// x := V^-1 m + B'Y = x + z
	for (int i = 0; i < p; ++i)
		gsl_matrix_set(x, i, 0, gsl_matrix_get(x, i, 0) + gsl_matrix_get(z, i, 0));

	// m* := (V^-1 + B'B)^-1 (V^-1 m + B'Y) where V* = (V^-1 + B'B)^-1
	gsl_matrix *mstar = gsl_matrix_alloc(p, 1);
	gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, Vstar, x, 0.0, mstar);

	gsl_matrix_free(x);
	gsl_matrix_free(z);

	return mstar;
}

You have to break down the expressions into tiny little chunks, like you are writing assembly programs for arithmetic.  It is very frustrating.

It is low-level, of course, because the direct translation of matrix expressions to code is pretty inefficient, with lots of temporary matrices being created and deleted.  Still, that is pretty much how my code looks anyway, I am just manually doing the direct translation, ’cause I don’t want to spend too much time trying to optimise the code before I know that it actually works (and is too slow to be worth optimising).

Template expressions, like in Blitz++, should be able to handle this, but I have never managed to learn that library.  I don’t really deal with matrix expressions that often, so I haven’t bothered.  If I end up having to do much more matrix hacking, I will.  Life is just too short to manually coding it up…


118-134=-16

This just isn’t right

Friday, April 10th, 2009

So, I’m working on some genome scale analysis and am a bit worried about numerical precision issues.  After all, even relatively stable operations might run into problems after milions of additions and multiplications.

So I figured I could try using interval arithmetic to get a feeling for the uncertainty in the numerical values.

Boost, as so often, is helpful here and has a library for interval arithmetic.

Great, I thought, I’ll just chuck in the interval type in my algorithm (that is parameterised with the numerical type anyway so I can use either float, double or long double) and see how it goes.

Well, no luck.  All my intervals ends up being a single value, and that just isn’t right.

A simple example like this, that that adds 1000 random floats, results in an empty interval:

#include <iostream>

#include <boost/numeric/interval.hpp>
using namespace boost::numeric;
using namespace interval_lib;

#include <ctime>
#include <boost/random.hpp>

int main (int argc, char * const argv[])
{
	typedef interval<double> I;

	static boost::mt19937 rng(std::time(0));
	static boost::uniform_real<> uni_dist(0,1);
	static boost::variate_generator<boost::mt19937&, boost::uniform_real<> > random_01(rng, uni_dist);

	I x = 0.0;
	for (int i = 0; i < 1000; ++i) {
		x += random_01();
	}
	std::cout << '[' << x.lower() << ',' << x.upper()
                  << "](" << (x.upper() - x.lower()) << ')'
                  << std::endl;

    return 0;
}

There is no way there is no precision lost in adding 1000 values!

Clearly I am doing something wrong, but what?

100-122=-22

Longest ordered subsequence

Friday, March 13th, 2009

If you are into string algorithms (and typical bioinformatics applications) you will love these two posts:

Yeah, I know it doesn’t look like bioinformatics at first sight, but it is very similar to the kind of string algorithms we use there.

It is great posts with a nice mix of coding and algorithmcs.

72-90=-18