Ready for review. Russell, could you take a look? It's pretty big, but there are some things I'd particularly like your opinion on, mostly minor details of the interface.

I've created GitHub pull requests for the three affected packages (shapelet, meas_extensions_multiShapelet, and meas_multifiit), which are linked from the JIRA issue page.

DM-1181 review

Based on the code in meas_multifit, making an array of builders seems rather complex. One needs a vector of factories, and then to use that to compute the workspace size. I see very similar code to do this in psf.cc and UnitTransformedLikelihood.cc.

Jim's response in email:

Yeah, I thought about unifying that in a single interface in shapelet (or meas_multifit) as well, and it's not so much that I don't want to do it as I didn't want to do it on this ticket; it was already getting too large, and I figured it'd be best to defer this to another one, as the design of such a unification will take a bit of thought.

shapelet MatrixBuilder.h: how does the shared workspace get passed to the MatrixBuilder? I don't see any suitable arguments in MatrixBuilder constructors, nor any suitable member functions. So far I've worked out that MatrixBuilderFactory calls _impl::apply(workspace) but I'm not yet sure what that does.

Jim response in email:

MatrixBuilder also has a private constructor that just takes an Impl. That's what's used by MatrixBuilderFactory, which is what you have to use if you want to pass an external workspace. In fact, all of the other constructors actually delegate to that one too, because they just create a MatrixBuilderFactory, invoke that, and assign the result to this.

My response to that:

Thank you. Please do document how the factory builds the MatrixBuilder; I think it would help others understand the code. Also, I am not very comfortable with a factory function using a private interface (I assume it has something to do with the use of shared memory, but surely you could just pass in a workspace that is explicitly shared and document it as such?)

There are many different implementation classes in MatrixBuilder.cc and I'm finding it tough to figure out which one is for which class. It might help to split these things up (the C++, at least). I'm pretty sure the right one is:

class< Factory : public MatrixBuilderFactory<T>::Impl

but its apply simply calls apply(workspace) for each component, so the real work is occurring even further down the rabbit hole.

Jim response:

Note that there are two Impl hierarchies: one inheriting from MatrixBuilder::Impl, and one inheriting from MatrixBuilderFactory::Impl. The derived classes for the latter are defined within the derived classes for the former, as nested classes

A MatrixBuilderFactory::Impl's job is to create a MatrixBuilder::Impl (what its apply() method does), and provide a way to compute the amount of workspace needed before the workspace is allocated (what computeWorkspace() does).

A MatrixBuilder::Impl's job is to actually evaluate the matrix (what its apply() method does).

It sounds like one early takeaway from the review is that I should rename both of these apply methods (and possibly one or both of the public operator() methods that delegate to them) to keep them from being confused with each other.

My response:

Anything to clear it up would help. There's a lot of different similar-sounding objects being defined and declared in one pair of .h/.cc files (especially since the implementations have base classes). I think it would help to split the .cc (and probably the .h as well) into 3 files: MatrixBuilder, MatrixBuilderFactory and MatrixBuilderWorkspace.

Why is this useful: workspaceSize = std::max(workspaceSize, factories.back().computeWorkspace()); in particular, what is special about the last factory?

Jim's response:

It's the one we just constructed in the line above. It might be more clear to have instead to put the new factory in a local variable, call computeWorkspace() on it, and only then add it to the factories vector, but that'd also be less efficient (an unnecessary copy of the new factory), at least until we have C++11 move.

My response:

I understand that the last factory is the one used to call computeWorkspace; what I was asking is why you could get away with that and ignore the other factories that came before it. Maybe I should have said "why does this work" instead of "why is this useful". My guess is that in some way each factory allocates space or bumps a counter, so that the final workspace has the final value of the allocated space or the fully incremented counter. But if so...it seems a bit odd not to ask the object managing the space or the counter (a workpace-like object) how much room is required. Also note that the documentation for CompoundImpl explicitly says "the workspace needed by the compound implementation is actually the maximum needed by any of its components", which suggests that every factory needs to be taken into account when computing the workspace.

MatrixBuilderWorkspace seems like a generic memory pool for matrices and vectors. I don't see anything specific to shapelets or Matrix

Jim's reponse:

Probably, someday, but I didn't want to move it to a lower package until we had a concrete use for it elsewhere, and I didn't want to take the time now to ensure the design was sufficiently general to be of use for other tasks (though I think it probably is anyhow).

shapelet MatrixBuilder.cc I am a bit nervous seeing so much implementation go into anonymous namespaces. I realize you can get it it by calling the right constructors of the public classes for which the implementation is defined, but it would be nice to be able test the implementations directly in unit tests and see their documentation in Doxygen.

Minor issues and questions:

• shapelet MatrixBuilder.h typos in otherwise very helpful doc strings (warning: I removed line wraps):
• "it may be useful to use one workspace for several a sequence of throwaway"
• "share the same memory, but main different "current" pointers, and hence create arrays"
• shapelet MatrixBuilder.cc typo in doc string:
• "It olds the original coordinate arrays" (olds instead of holds?).
• "remapped to an different"
• shapelet MatrixBuilder.cc request for clarification of doc string:
• "In the code, we call that the "lhs" matrix": please clarify what "that" refers to. Also, in the next paragraph you refer to a "higher-order lhs basis" shortly afterwards, without ever quite defining it. Maybe they are the same thing. In any case, those two paragraphs could use a bit of clarification.

• meas_multifit UnitTransformedLikelihood.cc: please document the arguments in the various functions (they were not in the old code, either).

• meas_multifit psf.cc: Where does template parameter "Pixel" come from and what is it? This is not a change; I'm just curious.

• meas_multifit MatrixBuilder.h: why did useApproximateExp go away? I assume you had a good reason, I'm just curious.

From Russell:

Thank you. Please do document how the factory builds the MatrixBuilder; I think it would help others understand the code. Also, I am not very comfortable with a factory function using a private interface (I assume it has something to do with the use of shared memory, but surely you could just pass in a workspace that is explicitly shared and document it as such?)

I'll be sure to add that documentation. But I have to disagree with your comment about factory functions and private interfaces; I think one of the main reasons one uses factory functions and classes is to provide a public construction interface that is more complex than a simple constructor would allow. Whether the constructor the factory ultimately calls should be public should be based on whether it is useful and safe for other public users of the class, and in this case, it simply isn't (because it takes an Impl object the user can't construct).

Anything to clear it up would help. There's a lot of different similar-sounding objects being defined and declared in one pair of .h/.cc files (especially since the implementations have base classes). I think it would help to split the .cc (and probably the .h as well) into 3 files: MatrixBuilder, MatrixBuilderFactory and MatrixBuilderWorkspace.

I'm torn about this; I agree that the file is large enough that it's getting a bit unwieldy, but it'd be impossible to split it up along these lines, as the implementation of each builder/factory pair is much to closely intertwined (each pair definitely has to be in since source file). I could split up just the .h files, and not the .cc files, but I tend to think that's more confusing that just have a large file for both. I also considered splitting up the factory/builder pairs into separate header and source files (i.e. one for the base Impls, one for Shapelet, one for ConvolvedShapelet, etc.). But that puts a lot of purely private implementations in header files, and my feeling is that shouldn't be done unless you need to do it for testing purposes or public class definitions, as it pollutes the user-visible public interfaces and can cause too much stuff to be included downstream (in e.g. Swig wrappers, where our tendency is to include everything when something goes wrong, because it's so hard to determine exactly what's needed), making the overall library API more fragile.

I understand that the last factory is the one used to call computeWorkspace; what I was asking is why you could get away with that and ignore the other factories that came before it. Maybe I should have said "why does this work" instead of "why is this useful". My guess is that in some way each factory allocates space or bumps a counter, so that the final workspace has the final value of the allocated space or the fully incremented counter. But if so...it seems a bit odd not to ask the object managing the space or the counter (a workpace-like object) how much room is required. Also note that the documentation for CompoundImpl explicitly says "the workspace needed by the compound implementation is actually the maximum needed by any of its components", which suggests that every factory needs to be taken into account when computing the workspace.

I'm not ignoring the factories that came before; they set the previous value of the workspace size, so at each iteration, I'm updating the workspace size to be equal to the maximum of all the workspace sizes I've encountered so far. When the iteration is done, that workspace size variable is indeed the maximum workspace size needed by any individual factory. Once I've done that, I can allocate a workspace object that thus has enough space for any individual builder, and iterate over the factories, using the same workspace pointer to construct each builder. The builder(s) that needed the most workspace will be the only one(s) that will use all of it. Calling computeWorkspace() on a factory doesn't increment any counters or allocate any space; it just reports how much workspace that factory would need to construct a builder. It's only when constructing a builder that the workspace counter is incremented (and in this case, we circumvent that by copy-constructing the workspace object, to make it so the builders all reuse the same workspace).

shapelet MatrixBuilder.cc I am a bit nervous seeing so much implementation go into anonymous namespaces. I realize you can get it it by calling the right constructors of the public classes for which the implementation is defined, but it would be nice to be able test the implementations directly in unit tests and see their documentation in Doxygen.

My goal here is very much to keep these out of the public API and IMO that means they should stay out of the public documentation (why I used standard comments in the .cc file instead of Doxygen comments). In the future I can imagine adding other derived classes for more specific cases that may need to be optimized (e.g. Gaussians, which are just zeroth-order shapelets) and I don't that want to produce an API change. It's also just easier, because I can avoid dealing with the pain of Swig for things that don't translate well to Python, like Eigen::Map objects. I agree that this approach does make testing trickier - the tests need to know about some implementation details, if only to know when the coverage is complete - but I think I've managed to test everything without exposing any of the internals, and I hope I'll be able to continue that going forward.

I'll make the various doc fixes you requested; thanks for the detailed inspection.

meas_multifit psf.cc: Where does template parameter "Pixel" come from and what is it? This is not a change; I'm just curious.

It's actually a global typedef in meas_multifit (in constants.h), much like the one in Kernel/Psf in afw.

meas_multifit MatrixBuilder.h: why did useApproximateExp go away? I assume you had a good reason, I'm just curious.

It turned out not to work in the most important cases, and it turned out not to be necessary:

• The approximate exponent implementation we had was jagged if you zoomed in enough, and that played havoc with a nonlinear optimizer that used finite differences for derivatives.
• Switching from double precision to single precision moved exp from the obvious bottleneck in the code to somewhere around 10% of the total processing.

Ok, I've made the changes I promised above (method renaming, doc fixes and improvements). Will merge master later today, but I'll give you a chance to reply before then if you like.