Learning to walk again
Wow, I’ve found a lot of bugs in my code this week.
To be fair, I haven’t written any serious C++ code since 2003. A lot of idioms from Java and R have infected my programming style, so I’ve been making more than my usual number of typos, but these generally get picked up at compile time.
My Gibbs sampler is less than 400 lines of code, so it barely even qualifies as “serious” software development. Nevertheless, the flaws that I’ve uncovered are embarrassing newbie mistakes, such as the following:
- setting the precision to the inverse square root of the standard deviation (power of -0.5 instead of -2)
- using the sum of squares instead of the sum of squared differences (SSD)
- forgetting to take the square root of the inverse gamma when returning the standard deviation
- forgetting to set the old labels z to zero before sampling the new labels
- calling abs(..) instead of fabs(..) and thereby rounding to the nearest integer
- forgetting that z has one more row than y in my algorithm (this simplifies the definition of the neighbourhood matrix)
- using n_cols instead of n_elem on an arma::rowvec didn’t behave as I expected
- calling size(..) on an arma::rowvec definitely didn’t behave as I expected!
The only way that I tracked these down was by going back to basics. A hidden Potts model can be viewed as a specialization of an independent mixture model, which in turn is a specialization of a single distribution (i.e. a mixture model with only one component, k=1). With the assumption of additive Gaussian noise, we end up with the following inheritance diagram:
By adding these simpler classes, it enabled me to test methods like gibbsStdDev(…) in isolation. I also ran the code against simulated data, instead of my enormous cone-beam CT scans. The moral of the story is: baby steps, or you’ll trip over yourself!