Introduction

The Measurement Equation Analysis Library (MEAL) is a set of classes, implemented in C++, that may be used to perform non-linear least-squares and error propagation through matrix equations. Arbitrary functions are built by modular construction of elementary functions, which may be parameterized by any number or type of arguments.

It was written by Willem van Straten (@straten) to support the research and development of

Measurement Equation Modelling (MEM) Radio Astronomical Polarimetry and Point-Source Calibration, The Astrophysical Journal Supplement Series, 152:129 (2004)
Matrix Template Matching (MTM) Radio Astronomical Polarimetry and High-Precision Pulsar Timing, The Astrophysical Journal, 642:1004 (2006)
Measurement Equation Template Matching (METM) High-fidelity Radio Astronomical Polarimetry Using a Millisecond Pulsar as a Polarized Reference Source, The Astrophysical Journal Supplement Series, 204:13 (2013)

Function Characteristics

A mathematical function or expression is characterized by its return type and input variables; a distinction is made between function parameters, \(\boldsymbol{a}=(a_0, a_1, ... , a_N)\) and independent variables, or arguments, \(\boldsymbol{x}=(x_0, x_1, ... , x_M)\). Parameters are double-precision floating point values that can vary freely during modeling. Arguments have arbitrary type and cannot be treated as free variables during modeling.

The parameter and argument access interface is defined by the MEAL::Function base class, and the return type is defined by evaluation base classes that inherit Function. The evaluation base classes define an evaluate method that returns a result, \(f(\boldsymbol{x},\boldsymbol{a})\), and the partial derivatives of this result with respect to each of the parameters, \(\partial f/\partial a_i\).

The behaviour of function expression classes is organized into three categories: MEAL::Parameter, MEAL::Argument, and MEAL::Evaluation. These behaviours are implemented by the children of three base classes, each of which inherits the MEAL::FunctionPolicy base class. The behaviour policies are mutually exclusive, and derived classes can incorporate pre-defined behaviours by setting the appropriate policy attribute.

Parameter Policy

A function expression may have an arbitary number of scalar parameters, \(\boldsymbol{a}=(a_0, a_1, ... , a_N)\). Each parameter has an associated name, estimated variance, and a flag that indicates if the parameter is free or fixed. Parameter management and access is implemented by children of the MEAL::ParameterPolicy base class.

Argument Policy

A function expression may be further parameterized by an arbitrary number of independent variables, or arguments, \(\boldsymbol{x}=(x_0, x_1, ... , x_M)\). Unlike parameters, arguments have no specified type, no estimated variance, and can never be free. Because they have no specified type, the interface between a Function and its arguments is mediated through the MEAL::Argument and MEAL::Argument::Value abstract base classes. Argument management and behaviour is implemented by children of the MEAL::ArgumentPolicy base class.

Evaluation Policy

The return type of a function expression is unspecified in the MEAL::Function base class definition. Therefore, derived classes inherit MEAL::Function via the MEAL::Evaluable template, which has a template argument that determines the return type, Return, and defines the following method.

  virtual Result evaluate (std::vector<Result>* gradient = 0) const;

This method returns a result \(M\) and, if gradient is not null, the partial derivative of the result with respect to each of the model parameters, \(\partial M/\partial a_i\).

For convenience of notation, there are intermediate classes that inherit the MEAL::Evaluable template. For example,

MEAL::Scalar - inherits MEAL::Evaluable<double> and represents functions that return a scalar value, \(f({\bf a}; {\bf x})\)
MEAL::Complex2 - inherits MEAL::Evaluable<Jones<double>> and represents functions that return a \(2\times2\) complex-valued matrix, \(J({\bf a}; {\bf x})\)

Model Composition

A number of template classes may be used to simplify the modular construction of more complicated functions. These templates implement the following basic rules of differentiation.

MEAL::ChainRule - an arbitrary function in which one or more parameters is set equal to the ordinate of a MEAL::Scalar function
MEAL::BinaryRule - an associative binary operation, such as the MEAL::SumRule or MEAL::ProductRule.

For example, use MEAL::ChainRule to create a feed in which the receptors have equal ellipticities.

  auto feed = new MEAL::Feed;
  auto feed_chain = new MEAL::ChainRule<MEAL::Complex2>;
  feed_chain->set_model (feed);
  auto ellipticity = new MEAL::ScalarParameter;
  feed_chain->set_constraint (0, ellipticity);
  feed_chain->set_constraint (2, ellipticity);

Use MEAL::ProductRule create multiply a boost along the Stokes V axis and a rotation about the Stokes U axis (in a basis defined by linearly polarized receptors).

  auto product = new MEAL::ProductRule<MEAL::Complex2>;
  auto boost = new MEAL::Boost1 (Vector<3, double>::basis(2));
  product->add_model (boost);
  auto rotation = new MEAL::Rotation1 (Vector<3, double>::basis(1));
  product->add_model (boost);

Modular Construction

TO DO: Show how new functions can be built up from more basic elements.

Example Usage

Non-linear Least-Squares Estimation

TO DO: Document lmfit

Error Propagation

The Estimate template class is very useful for storing a value and its estimated variance. There are also operators and functions which enable the propagation of error estimates to derived quantities; for example:

  Estimate<float> y (0.596,0.0034);
  Estimate<float> x (-0.83,0.0072);
  Estimate<float> z = pow (y,x);

automatically computes the variance of the new variable, z. However, the Estimate template class fails when a variable appears more than once in an expression; e.g.

Estimate<float> ratio = x/x;

should yield \(1\pm0\); however, the Estimate template class does not recognize that the numerator and denominator are the same variable, and incorrectly sums the weighted variances.

The problem of correctly computing the partial derivatives of an expression with respect to its variables makes use of the exact same functionality used to generate the gradient and Hessian matrix in non-linear least squares fitting.

A simplified interface to this functionality is implemented by the MEAL::ScalarMath class. MEAL::ScalarMath objects may be conveniently initialized as a single parameter and its estimated variance using the Estimate template class. As with float and double types, ScalarMath objects may be combined using normal arithmetic operations and basic mathematical functions, creating Scalar functions of any number of parameters. For example:

  MEAL::ScalarMath x (Estimate<double> (0.9587, 0.00058));
  MEAL::ScalarMath y (Estimate<double> (-0.283, 0.00034));
  cerr << "Polar angle = " << atan2 (y, x) << endl;

yields the output

  Polar angle = (-0.287039 +/- 0.0189612)

As with any native type, the ScalarMath class can be used as a template argument, e.g.

  complex<MEAL::ScalarMath> z (Estimate<double> (0.87, 0.0041),
                               Estimate<double> (2.38, 0.0095));

  complex<MEAL::ScalarMath> w (Estimate<double> (1.74, 0.0081),
                               Estimate<double> (-.63, 0.0043));

  Jones<MEAL::ScalarMath> jones (z, conj(z),
                                 conj(w), w);

enabling error propagation through increasingly complex expressions.