# Baby steps, continued

I’ve found some more bugs in my code without even looking for them:

I added an entry point for simulating from the Potts model for a fixed value of β, without any external field. This function can be used to estimate the critical value of β by simulation (since β doesn’t depend directly on the data). In doing this, I discovered that my neighbourhood matrix wasn’t a matrix … it was a vector! It appears that the following code silently converts from Rcpp::IntegerMatrix to Rcpp::IntegerVector and thence into an arma::umat with only a single column:

arma::umat neigh = unsign(nR - 1);

This is because Rcpp::Matrix inherits from Rcpp::Vector (shouldn’t it be the other way around, generalization not specialization?) and the binary arithmetic operators are defined to return Rcpp::Vector. Exactly the kind of shenanigans that Felix von Leitner’s polemic warns about…

So, once I fixed bug number 1, this uncovered other problems deeper inside my code:

double total_prob = 0.0; const arma::uvec neighbors = neigh.row(i); const arma::umat labels = z.rows(neighbors); for (unsigned j=0; j < labels.n_cols; j++) { unsigned sum_neigh = arma::sum(labels.col(j)); total_prob += exp(beta*sum_neigh); }

The above code doesn’t even come *close* to performing the vectorized operations that I was expecting. In the end, I gave up and just wrote it as a nested loop:

double total_prob = 0.0; for (unsigned j=0; j < z.n_cols; j++) { unsigned sum_neigh = 0; for (unsigned k=0; k < neigh.n_cols; k++) { sum_neigh += z(neigh(i,k),j); } total_prob += exp(beta*sum_neigh); }

Trying to be too clever just hasn’t paid off with RcppArmadillo.

## Trackbacks & Pingbacks