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.

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.
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.
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.

There is no one "right" way to implement what you describe. There are many ways is could done and which is best is both a matter of opinion and would require detailed knowledge of the context you are working in. Therefore, I am going to vote close this question as too broad. As regards to the question Reap and Sow you tacked on at the end, you should be able answer that yourself by performing a simple experiment or two with Mathematica. — m_goldberg, Jul 09 '14 at 14:06
I edited the question so that it's no more opionion-based. About the Sow/Reap, because at present, i'm just planning, I don't have a context to experiment with. Basically I'm asking if there are know issues using this approach. — unlikely, Jul 09 '14 at 14:13
@m_goldberg I've only read the updated version of this question, but I do think it's a relevant question. There is indeed a subjective element to the decision, but it's more complicated than that. The choice of the initial design will have implications which may not be easy to see/predict. Once the design is fixed, it may not be easy to change it. So choosing the best design at the beginning is important. (This situation is familiar to me from C++: I see several ways to design something, but I only notice the problems with my chosen design when it's too late to turn back...) — Szabolcs, Jul 09 '14 at 14:42
I don't have time to write an answer right now (maybe later today), but I'd go with no. 2, and would definitely avoid 3. — Szabolcs, Jul 09 '14 at 14:52
@Szabolcs. In my opinion the OP has not given enough information concerning the design issues specific to his problem. Therefore, it seems to me that answers could only discuss general design patterns, which I feel is too broad a topic for this site. If the community feels differently, as it might, then the question will simply stay open. — m_goldberg, Jul 09 '14 at 15:41
@m_goldberg What information you want I add to my question? I wrote an example of a specific problem, and I think it's enough to explain what I want to do, but please see also my next edit to the question. — unlikely, Jul 09 '14 at 16:02
@Szabolcs Why would you avoid 3? (I'd also feel uncomfortable doing that, but am not sure why) — acl, Jul 09 '14 at 16:25
What I would like to know why returning the full list {u, K, F} is not satisfactory. The user can always ignore the quantities he/she is not interested in. I would also like to know if recalculating EllipticDirichlet to get get the ancillary values is an issue. Like Szabolcs I am inclined to think implementing a simple option would be the quick and easy way; the option, if False (default), means the function returns u, but if True, the function returns {u, K, F}. But this is only one of six possibilities I can think of. An answer that discussed them all would be too long for this site. — m_goldberg, Jul 09 '14 at 17:00
@m_goldberg I think returning always all ancillary results is confusing for the consumer of the my function: taking parts of returned expression every time we call the function it's annoying and... confusing... because the real main output of the function is the solution. Retuning a list from the function with index-based accessing prevent to change the implementation of the function, to add new ancillary data in the future. If you want, you can just outline the other 4 possibilities you think of. From what I can understand propose another approach... — unlikely, Jul 09 '14 at 17:29
@unlikely I hadn't seen this question and your update earlier, but I had asked a similar before (re: multiple outputs matlab style) http://mathematica.stackexchange.com/q/314/5 I would not recommend it though... the question was purely for curiosity sake :) — rm -rf, Aug 01 '15 at 16:13

score 7 · Accepted Answer · edited Jun 16 '20 at 09:23

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:

Doesn't show its (complex and irrelevant) internal structure by default. Use Format to achieve this.
Can be directly used for the main intended purpose, e.g. an InterpolatingFunction or a FittedModel can be directly used as numerical functions.
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.

I wrote this up in one go ... should go through it and make it more coherent when I have time. — Szabolcs, Jul 09 '14 at 17:36

score 5 · Answer 2 · answered Jul 09 '14 at 18:04

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 & *)

This seems basically what @Szabolcs is saying.... – Michael E2 Jul 09 '14 at 18:05 — Michael E2, Jul 09 '14 at 18:05

Returning ancillary computation results from a function

Update

2 Answers2

Don't use Sow/Reap for this

Bundling in lists is simple

Objects with properties