Returning ancillary computation results from a function

Question

Suppose I want to implement a function to solve a very basic Elliptic problem using Finite Element of some degree (It's just an example).

Normally the user of my function is only interested to the numerical solution of the Elliptic problem, i.e. the values of the solution on a proper grid. So, the user call something like

u = EllipticDirichlet[..., range, ...];

Sometimes the caller of my function may be interested to some internal detail of the computation, for example the linear system formulation of the problem, i.e. the stiffness matrix K and the load vector F.

Ideally my function should return these ancillary data tagged in some way, so the user of my function can get these data simply and in a roubust.

I don't want my fuction return everytime the solution and all ancillary data in a list like in the following pseudocode.

{u, K, F} = EllipticDirichlet[..., range, ...];

What are the way to handle this situation?

I can think to the following way.

1. Using Options or Arguments to tell to the function what to return. Not clean in my opinion, not because of the options, but because it would be preferable to keep main result distinct from other ancillary intermediate computation data.
2. Returning an "object" with "properties" like LinearModelFit for example. Probably the best idea but make the function result more difficult to use and to understand for a casual reader.
3. Using Sow. But I'm not sure if this this a use-case for Sow? For example it's legal to Sow from a function if no one Reap? Documentation say nothing about this.

There are other way to accomplish this task?

Update

I'm anything but a MATLAB fan, but what I'm asking it's a very common scenario under MATLAB, and it's very simple to consume and provide functions that returns differently "shaped" answer with different "levels" of details. For example in MATLAB I can call a function with syntax like this

[x, u] = EllipticSolve[...]

if I'm only interested to the x grid used and solution u on grid vertices; or I can call

[x, u, K] = EllipticSolve[...]

if I'm occasionally also interested to the stiffness matrix K of the linear system used to solve the problem; I can also call

[~, ~, K] = EllipticSolve[...]

if I'm interested only to the stiffness matrix K. I think the function EllipticSolve can also do some optimization of the computation recognizing the output requested.

Basically I whish to be able to do something like this.

Answer 1

I do not have a lot of first hand experience with this, as I've never took the time to implement a proper solution for this problem. Also, I don't know a lot about FEM methods. So what I am going to say is mainly based on observing how various Mathematica functions work.

Don't use Sow/Reap for this

I don't think Sow/Reap are designed with this application in mind.

It is recommended to always use tags (seconds argument of Sow) with Sow/Reap. Otherwise things like this may happen:

SetAttributes[f, HoldAll]; 
f[x_] := Reap[x; Sow[1]] 
Reap[f[Sow[2]]]  (* users doing this may not know about the internals of f[] *)

The inner Reap catches everything (incorrectly) and leaves nothing for the outer one.

Having functions which do not Reap Sowed values can be similarly unsafe, even if using tags. Also, having to use tags already decreases the user-friendliness of the approach.

Bundling in lists is simple

The simplest way is to return a list, but of course it's possible to do better when writing a polished package.

In version 10, you can also consider bundling the information in an Association. The result will be very similar to using "properties", as FittedModel does, except that there can't be a default behaviour of the output (similar to how FittedModel can be used as a numerical function).

Objects with properties

I recommend using solution (2): try to imitate basic object oriented programming and return an expression which:

1. Doesn't show its (complex and irrelevant) internal structure by default. Use Format to achieve this.

2. Can be directly used for the main intended purpose, e.g. an InterpolatingFunction or a FittedModel can be directly used as numerical functions.

3. Has a number of properties or methods, which can be accessed as object["someProperty"]

This is very common with built-in Mathematica functions, but usually it is not (well) documented. You can use FittedModel, InterpolatingFunction, SparseArray, MeshRegion, etc. this way. The NDSolve plugin framework also uses this style.

Usually, properties or methods can be queried as object["Properties"] or object["Methods"] (depending the type of the object). If you stick to this convention, then the use of your objects will already be familiar to many users.

To make these objects even more user friendly, you may choose to define functions for retrieving these properties. For example,

StiffnessMatrix[x_FEMSolution] := x["StiffnessMatrix"]

This is what the DifferentialEquations`InterpolatingFunctionAnatomy` package does (take a look at its source).

As I see it, the main disadvantage of having top level functions like this to retrieve properties is that is requires introducing a large number of top-level functions (and choosing a large number of names, each of which has a potential to conflict with other packages). Using the object["StiffnessMatrix"] syntax has the advantage of tying the name "StiffnessMatrix" to this one special object type. The name cannot conflict with anything except other properties or methods of the same type of object.

This technique is really emulating classes from object oriented programming (without things like inheritance), and solves the same type of problem classes do: organization. If a certain function is only going to be used with a certain type of data structre, then it's best to formally define the data structure and tie the name of the function to it.

---

So there are many advantages to using solution (2), but what are the disadvantages?

One possible disadvantage is that if you bundle a lot of rarely used ancillary information into one object, it is going to take up a lot of space. Depending on the particular application, it may not be worth permanently storing all this information when most of the time the user only wants to use the main result only, which might be considerably smaller than the ancillary information.

For this situation I recommend the following:

A task like solving PDEs using a FEM method can usually be broken down into several steps. You may want to implement these steps as separate fuctions. These functions will work with a data structure that contains full information. You can make this (large) data structure accessible through the methods/properties calling style used by FittedModel.

You can also define a higher level function which calls these lower level functions in the correct order and computes the result in one go. This high level function can then return a data structure which contains less ancillary information.

Advanced users of your package can (if they wish to) use the lower level functions and work with the data structure that contains full information (but is more difficult to use).

Others can use the easy-to-use higher level functions and understand that the price to pay for this convenience is that they only have access to some ancillary information.

---

Mathematica 10 contains a FEM framework. You may want to read up on how it works. I did not yet study it in detail, but it seems to work in a similar way to what I described above. See here and related tutorials. You can also study the advanced NDSolve documentation and the NDSolve method plugin framework to get ideas.

Answer 2

I have used Sow/Reap in simple situations, but for what your proposing, I would usually follow your second idea. Something like this toy example:

myfn[___]["Properties"] = {"a", "b", "function"};
myfn[x_, a_, b_, func_]["a"] := a;
myfn[x_, a_, b_, func_]["b"] := b;
myfn[x_, a_, b_, func_]["function"] := func;
myfn[x_, a_, b_, func_][t_?NumericQ] := func[t];
Format[myfn[x_, a_, b_, func_]] := myfn[x];

mycalc[x_] := Module[{a, b, func},
   a = x^2;                             (* sub-computation one *)
   b = x^3;                             (* sub-computation two *)
   func = Evaluate[(1 - #) a + # b] &;  (* "return value" *)
   myfn[x, a, b, func]
   ];

Example:

f = mycalc[4]
(* myfn[4] *)

Using the returned myfn:

f[0]
(* 16 *)

f[1/2]
(* 40 *)

f["b"]
(* 64 *)

f["Properties"]
(* {"a", "b", "function"} *)

A user who only wanted the function and wanted to discard the rest of the data could execute the following:

f = mycalc[4]["function"]
(* 16 (1 - #1) + 64 #1 & *)