Sunday, January 17, 2016

On calling Ceres from Scala - 4.1: Hello World explained

This is the fifth installment of a series documenting my progress on porting the Ceres solver API to Scala. All my sources quoted below are available on github.

Hello World Explained


In this post I'll go through the nontrivial parts of my first SWIG wrap of the ceres-solver library. The code is in the sandbox section of the skeres github repo, and includes a port to pure Scala of Ceres's helloworld.cc example.

Build

The library and example are built using a simple Bourne shell script (we'll tackle sbt configuration later). The ceres-library headers of interest are swigged into a com.google.ceres Java package, and the JNI auto-generated .cxx source file is produced in the toplevel directory, as is the shared library (.dylib file in Mac OSX) they are eventually linked into. The ported example itself belongs to a separate org.somelightprojections.skeres package.

Code

The example's code closely follows the C++ original, with one basic difference: the cost function derivative is computed analytically, rather than through automatic differentiation, as the latter will be the subject of the next step in this port. Therefore the residual term code looks like this:

// Concrete cost functor: cost(x) = 10.0 - x
// Base class CostFunction is the SWIG-wrapped ceres one.
class CostFunctor extends CostFunction {
// This functor adds one par. block.
mutableParameterBlockSizes.add(1)
// ... with just one residual.
setNumResiduals(1)
// Implements the CostFunction::Evaluate API.
override def evaluate(
parameters: DoublePointerPointer, // C++: double const* const*parameters
residuals: DoublePointer, // C++: double* residuals
jacobians: DoublePointerPointer // C++: double** jacobians
): Boolean = {
val x = parameters.get(0, 0)
residuals.set(0, 10.0 - x)
if (!jacobians.isNull) {
jacobians.set(0, 0, -1.0)
}
true
}
}
As you can see, it follows closely the C++ implementation (i.e. what the C++ codes would look like if it didn't use autodiff): it describes the size of its parameter and residual blocks, and then overrides the Evaluate method to actually compute the cost. Note the use of camelCase-style identifiers, in keeping with the usual Scala style conventions.

Probably the only part deserving an explanation is the manipulation of the residual array and of the parameter and jacobian matrices. This takes advantage of a very simple C++ helper declared inline in the ceres.i SWIG configuration file:

// Convenient access to matrices represented as double** (e.g.
// the Jacobian passed to the cost function).
// See implicits in skeres.scala to see how these are
// used in practice.
%include "carrays.i"
%array_class(double, DoubleArray);
%inline %{
struct DoubleMatrix {
static bool isNull(double** matrix) {
return matrix == NULL;
}
static double* row(double** matrix, int i) {
return matrix[i];
};
};
%}
Remember the discussion in a previous post about the use of the carrays.i SWIG library. The %array_class stanza declares DoubleArray to be a class wrapping an array of double's. The inline DoubleMatrix helper exposes a null test and access to the matrix by rows. Thanks to the magic of implicit conversions, this is all we need to make these low-level classes behave in a more idiomatic Scala way.

First, for one dimensional arrays of doubles we make use of the enrichment pattern to add the usual get/set accessors, as well as conversions of DoubleArray objects to pointers and Scala arrays:

case class RichDoubleArray(a: DoubleArray) {
def get(i: Int): Double = a.getitem(i)
def set(i: Int, x: Double): Unit = a.setitem(i, x)
def toPointer: DoublePointer = a.cast
def toArray(length: Int): Array[Double] = {
val aScala = Array.ofDim[Double](length)
aScala.indices.foreach(i => aScala(i) = a.getitem(i))
aScala
}
}

We use a similar enrichment for matrices:

case class RichDoubleMatrix(data: DoublePointerPointer) {
def isNull: Boolean = DoubleMatrix.isNull(data)
def doubleArrayRow(i: Int) = DoubleArray.frompointer {
DoubleMatrix.row(data, i)
}
def get(i: Int, j: Int) = doubleArrayRow(i).getitem(j)
def set(i: Int, j: Int, x: Double) = doubleArrayRow(i).setitem(j, x)
}
Note the call to the (swig-wrapped) DoubleMatrix helper to access a matrix row. Note also that these calls involve only tiny JVM object creations, with no copies of the data arrays themselves.

Finally, we define implicit conversions to instantiate these enrichments from the swig-wrapped low-level pointers:

package object skeres {
// Readable renaming of SWIG-generated types.
type DoublePointerPointer = SWIGTYPE_p_p_double
type DoublePointer = SWIGTYPE_p_double
// Implicit conversions betwee pointer objects and
// arrays or matrices.
implicit def doublePointerPointerToRichDoubleMatrix(p: DoublePointerPointer) =
RichDoubleMatrix(p)
implicit def doublePointerToRichDoubleArray(p: DoublePointer) =
RichDoubleArray(DoubleArray.frompointer(p))
implicit def doubleArraytoRichDoubleArray(a: DoubleArray) =
RichDoubleArray(a)
}

Odds and ends

The ceres.i configuration file only includes as many ceres headers as are needed to compile the example. In addition, however, it also exports all the ceres-predefined "loss function" classes that are used to robustify the residuals. I mention this here to highlight another feature provided by SWIG, namely ownership transfer (i.e. memory leak avoidance). The relevant stanza looks like this:

%newobject PredefinedLossFunctions::trivialLoss;
%newobject PredefinedLossFunctions::huberLoss;
// ... more like the above
%inline %{
struct PredefinedLossFunctions {
static ceres::LossFunction* trivialLoss() { return new ceres::TrivialLoss; }
static ceres::LossFunction* huberLoss(double a) { return new ceres::HuberLoss(a); }
// ... more
};
%}
The inline-declared struct allows us to create loss function instances from Scala using calls like PredefinedLossFunctions.huberLoss(2.0).  The %newobject declarations cause the destructors for the newly created-and-wrapped C++ objects to be called whenever the wrapping Java object is garbage-collected.

A similar issue is addressed in the HelloWorld.scala source where the Problem object is declared:

val problemOptions = new Problem.Options
problemOptions.setCostFunctionOwnership(Ownership.DO_NOT_TAKE_OWNERSHIP)
problemOptions.setLossFunctionOwnership(Ownership.DO_NOT_TAKE_OWNERSHIP)
val problem = new Problem(problemOptions)
view raw ceres ownership hosted with ❤ by GitHub

Here we use the (wrapped) ProblemOptions to instruct Ceres not to take ownership of the cost and loss function objects. This is necessary because they are JVM-owned objects, and the default behavior in which the ceres Problem instance takes ownership and deletes them at in its destructor would likely cause a double-free heap corruption.

No comments:

Post a Comment