yomm2

Abstract

It is possible to use YOMM2 methods in combination with C++ templates. This tutorial starts with an overview of YOMM2’s public interface, then shows how it can be used to specialize ordinary or templatized methods with templatized definitions.

In any case, the user of a library that uses open methods with templates is required to initialize the library. However, it is possible for the author of the library to arrange for the initialization code to be minimal.

The problems with macros

YOMM2’s macros provide a convenient way of declaring methods, and adding definitions to them. It would be nice if the matrix example could be re-implemented for templatized matrices like this:

template<typename T>
declare_method(string, to_json, (virtual_<matrix<T>&>));

template<typename T>
define_method(string, to_json, (dense<T>& m)) {
    return "json for dense matrix...";
}

There are several reasons why this is not possible.

  1. YOMM2 macros create more than one definition behind the scenes. For example, declare_method declares a struct, defines an inline function that serves as the entry point for the method, and creates a static object that registers the method, so update is aware of it. The template-id (i.e. the template<...> specifier) would need to be applied to each of the entities generated by the macro. The template syntax is complex. For example, template arguments need to be repeated when a template is specialized (as in template<typename T> struct X<vector<T>>). Conceivably, an extra set of macros could be crafted, that would take both the template-id and the template arguments. These macros would be cumbersome (arguably, not much more than the C++ syntax). They would eventually fall behind the constant evolution of the language’s syntax (think of C++20’s requires clauses). Finally, they would run afoul of the next problem.

  2. Macros do not interact well with templates. For starters, angle brackets are ordinary characters, as far as the C preprocessor is concerned. As a consequence, a macro call like MACRO(std::pair<int, double>) will cause two arguments to be passed to MACRO, namely: std::pair<int, and double>. This is probably not what is intended.

  3. While C++’s template instantiation mechanism would cope with templatized method declarations gracefully - if it’s used, instantiate it -, it would see no reason to instantiate templatized method definitions, as they are not referenced directly by user code. Instead, YOMM2 wires the definitions into the appropriate methods by means of statically constructed objects. In presence of method definition templates, the static objects would not be instantiated at all. The situation is not completely specific to YOMM2: C++ has an explicit template instantiation mechanism to make it possible to cope with situations where the automatic instantiation falls short (fox example, in certain situations involving dynamically loaded shared object). Likewise, YOMM2 templatized definitions need to rely on a similar construct.

Using the YOMM2 API

YOMM2 can be used without macros, as explained in the API tutorial. If you have not read it yet, I suggest you do now, before going any further in this tutorial.

We are going to implement two libraries, both inspired by linear algebra. The first is a fully type-erased vector class that supports mixed type vector operations. The second is a more complex example involving matrix of different categories.

The libraries are just skeletons: they don’t contain any real math.

In both cases we will use a simple letter-envelope approach to manage object lifetime. Here is is:

class ref_count {
    mutable std::size_t refs{0};
  public:
    virtual ~ref_count() {}
    void add_ref() { ++refs; }
    void remove_ref() { if (--refs == 0) delete this; }
};

template<typename Representation>
class handle {
    Representation* rep;
  public:
    using body_type = Representation;

    handle() : rep(nullptr) {}
    explicit handle(Representation* r) : rep(r) { rep->add_ref(); }
    ~handle() { if (rep) { rep->remove_ref(); }}
    // etc
};

A vector library

In the first example, we implement a library to manipulate vectors of any numeric type. Mixed operations are supported, but the vector class itself is type erased. Here is is:

namespace vector_algebra {

struct abstract_vector : public ref_count {
    virtual ~abstract_vector() {}
};

template<typename T>
struct concrete_vector : abstract_vector {
    concrete_vector(std::initializer_list<T> values) : entries{values} {}
    template<typename U>
    concrete_vector(const std::vector<U>& values) : entries{values} {}
    template<typename U>
    concrete_vector(std::vector<U>&& values) : entries{values} {}
    concrete_vector(const concrete_vector&) = default;
    concrete_vector(concrete_vector&&) = default;

    const std::vector<T> entries;
};

using vector = handle<abstract_vector>;

vector ints(new concrete_vector<int>{1, 2, 3});
vector reals(new concrete_vector<double>{4., 5., 6.});

Operations on vectors.

We are going to implement three operations on vectors: addition, subtraction, and comparison. The first two are similar, in that they both take two vectors and yield a vector. The third takes two vectors and yields a boolean.

#include <yorel/yomm2/core.hpp>
#include <yorel/yomm2/symbols.hpp>

using namespace yorel::yomm2;

struct YOMM2_SYMBOL(addition);

using addition = method<
    YOMM2_SYMBOL(addition),
    vector(
        virtual_<abstract_vector&>,
        virtual_<abstract_vector&>)>;

struct YOMM2_SYMBOL(subtraction);

using subtraction = method<
    YOMM2_SYMBOL(subtraction),
    vector(
        virtual_<abstract_vector&>,
        virtual_<abstract_vector&>)>;

struct YOMM2_SYMBOL(comparison);

using comparison = method<
    YOMM2_SYMBOL(comparison),
    bool(
        virtual_<abstract_vector&>,
        virtual_<abstract_vector&>)>;

For the user’s convenience, we wrap the method call in an operator:

inline vector operator+(const vector& a, const vector& b) {
    return addition::fn(*a.get(), *b.get());
}

inline vector operator-(const vector& a, const vector& b) {
    return subtraction::fn(*a.get(), *b.get());
}

inline bool operator==(const vector& a, const vector& b) {
    return comparison::fn(*a.get(), *b.get());
}

Now we need to provide definitions for these methods. But which ones?

We could decide for the user that only (say) vectors of doubles and ints are supported. We could use add_function or add_definition to define four specializations, covering all the possible combinations (i.e. the Cartesian product):

concrete_vector<{int, double}> x concrete_vector<{int, double}>

This “closed” approach is lazy, and defies the purpose of open methods. It is not for the author of the library to decide which types to support. The user may have no interest in integer vectors, but may need complex vectors. Or vectors of decimal floating-point numbers. Or, types from another library, or perhaps created by the user himself. An advanced user may also want to add new operations or vector subtypes (like large sparse vectors). All this must be possible, and not unduly complex. Extensibility is what open methods are all about. Templates are pretty open too, because they can be specialized post-hoc.

After a lot of experimenting, I came up with a framework. It consists of a design pattern, a small meta-programming library, and, for library designers, a set of goals.

The pattern is:

  1. Put the definitions in one or several templates, where the method being defined is the first template argument.
  2. Create a meta-function that generates interesting combinations of methods and types.
  3. Instantiate the container(s) for these combinations. The resulting classes must satisfy the requirements of a definition container, or be derived from class not_defined. Add the definition containers to the method, extracted from the first template argument.

The meta-programming library comprises:

These constructs will be introduced as needed in the following examples.

The goals for library designers are:

  1. Make it very easy for a new user to use the library with all the defaults.
  2. Make it easy for a “normal” user to specify which parts of the library to use.
  3. Make it possible for an advanced user to extend the library with new types and new operations.

Let’s apply these ideas to the vector methods. First we declare a container template:

template<typename Method, typename...>
struct definition;

Let’s create partial specializations that defines addition, subtraction and comparison for vectors or two underlying types:

template<typename T, typename U>
struct definition<addition, T, U> {
    static vector fn(concrete_vector<T>& a, concrete_vector<U>& b) {
        using numeric_type = std::common_type_t<T, U>;
        std::vector<numeric_type> result(a.entries.size());
        std::transform(
            a.entries.begin(), a.entries.end(), b.entries.begin(), result.begin(),
            [](auto a, auto b) { return a + b; });
        return vector(new concrete_vector<numeric_type>(std::move(result)));
    }
};

template<typename T, typename U>
struct definition<subtraction, T, U> {
    static vector fn(concrete_vector<T>& a, concrete_vector<U>& b) {
        using numeric_type = std::common_type_t<T, U>;
        std::vector<numeric_type> result(a.entries.size());
        std::transform(
            a.entries.begin(), a.entries.end(), b.entries.begin(), result.begin(),
            [](auto a, auto b) { return a - b; });
        return vector(new concrete_vector<numeric_type>(std::move(result)));
    }
};

template<typename T, typename U>
struct definition<comparison, T, U> {
    static bool fn(concrete_vector<T>& a, concrete_vector<U>& b) {
        auto other = b.entries.begin();
        for (auto& value : a.entries) {
            if (value != *other++) {
                return false;
            }
        }
        return true;
    }
};

Let’s suppose that we want to perform those three operations for vectors of integers and vectors of doubles. We need to add definition<addition, int, int>, definition<addition, int, double>, etc, to method addition. And the same for subtraction and comparison. In other words, we want to apply template definition to the Cartesian product:

{addition, subtraction, comparison} x {int, double} x {int, double}

…then add the resulting definitions to their respective method.

For the Cartesian product, we can use two constructs provided by YOMM2: the types container, and the product meta-function.

types<typename...> is a simple list of types. It is similar to Boost.Mp11’s mp_list. Unlike mp_list, types does not have a body. This helps detect meta-programming bugs.

product<typename... TypeLists> takes any number of types lists, and returns a types list containing all the combinations of one element taken from each input list, each combination being itself wrapped in a types. For example:

static_assert(
    std::is_same_v<
        product<
            types<int, double>,
            types<int, double, float>
        >,
        types<
            types<int, int>, types<int, double>, types<int, float>,
            types<double, int>, types<double, double>, types<double, float>
        >
>);

use_definitions takes a class template and a Cartesian product, and applies the template to each of the elements of the product, which yields a list of classes. It throws away the classes that are derived from not_defined. For each remaining definition, use_definitions extracts the method to which the definition applies - by convention, the first argument of the definition template - and adds it to the method.

Finally, we must not forget to register the vector classes themselves.

use_classes<
    abstract_vector, concrete_vector<int>, concrete_vector<double>
> YOMM2_GENSYM;

use_definitions<
    definition,
    product<
        types<addition, subtraction, comparison>,
        types<int, double>,
        types<int, double>
    >
> YOMM2_GENSYM;

Everything is now in place. After calling update as usual, we can exercise the vector methods:

BOOST_AUTO_TEST_CASE(test_vectors) {
    update();

    {
        vector a(new concrete_vector<int>{1, 2});
        vector b(new concrete_vector<double>{3, 4});
        vector actual = a + b;
        BOOST_TEST(actual.type() == typeid(concrete_vector<double>));
        vector expected(new concrete_vector<double>{4, 6});
        bool correct = actual == expected;
        BOOST_TEST(correct);
    }

    {
        vector a(new concrete_vector<int>{1, 2});
        vector b(new concrete_vector<int>{3, 4});
        vector actual = a - b;
        BOOST_TEST(actual.type() == typeid(concrete_vector<int>));
        vector expected(new concrete_vector<double>{-2, -2});
        bool correct = actual == expected;
        BOOST_TEST(correct);
    }
}

Writing a user-friendly instantiation function

Now we have everything we need to write a single template that allows users to initialize the library for the types they need. For that we use std::tuple to lump teh registration object together:

template<typename... NumericTypes>
using use_vector_library = std::tuple<
    use_classes<
        abstract_vector, concrete_vector<NumericTypes>...
    >,
    use_definitions<
        definition,
        product<
            types<addition, subtraction, comparison>,
            types<NumericTypes...>,
            types<NumericTypes...>
        >
    >
>;

All the user of the library needs to do now is to create an instance (object) of use_vector_library, instantiated with the required types:

use_vector_library<int, double> init_vectors;

This example fulfills goals (1) and (2) assigned to library designers. In the following, more complex example, we will see how to fulfill goal (3) as well.

A matrix library

Matrices fall into categories depending on their layouts: square, symmetric, diagonal, etc. Operations involving matrices greatly depend on the matrix type(s). Transposing a symmetric matrix is a null operation. There are optimized algorithms for inverting diagonal matrices. It is not necessary to store all the zeroes in a diagonal matrix, storing the diagonal suffices. And, for a zero or identity matrix we need not store any numbers at all, apart from the matrix’s dimensions.

The following class diagram illustrates a possible design for a matrix template library that can be used with any numeric type T:

classDiagram
any~T~ <|-- ordinary~T~
any~T~ <|-- any_square~T~
any_square~T~ <|-- square~T~
any_square~T~ <|-- any_symmetric~T~
any_symmetric~T~ <|-- symmetric~T~
any_symmetric~T~ <|-- any_diagonal~T~
any_diagonal~T~ <|-- diagonal~T~
any_diagonal~T~ <|-- identity~T~
any_square~T~ <|-- any_triangular~T~
any_triangular~T~ <|-- upper_triangular~T~
any_triangular~T~ <|-- lower_triangular~T~
any~T~
ordinary~T~: size_t rows
ordinary~T~: size_t columns
ordinary~T~: vector~T~ entries
square~T~: size_t order
square~T~: vector~T~ entries
any_triangular~T~: size_t order
any_triangular~T~: vector~T~ entries
symmetric~T~: size_t order
symmetric~T~: vector~T~ entries
diagonal~T~: vector~T~ entries
identity~T~: size_t order

For brevity, we will consider only ordinary (i.e. non-square), square and symmetric matrices. Also, we will omit all maths and storage management, to focus on dealing with types.

Here is the implementation of these classes, and the intermediary abstract classes:

// ---------------------------------------------------------------------------
// matrix root class, parameterized by underlying numeric type

template<typename T>
struct ordinary;

template<typename T>
struct any : ref_count {
    virtual void concrete() = 0;
    using element_type = T;
    using abstract_type = any<T>;
    using concrete_type = ordinary<T>;
};

// ---------------------------------------------------------------------------
// matrix subtypes

template<typename T>
struct ordinary : any<T> {
    void concrete() override {
    }
};

template<typename T>
struct square;

template<typename T>
struct any_square : any<T> {
    using abstract_type = any_square<T>;
    using concrete_type = square<T>;
};

template<typename T>
struct square : any_square<T> {
    void concrete() override {
    }
};

template<typename T>
struct symmetric;

template<typename T>
struct any_symmetric : any_square<T> {
    static constexpr int specificity = 2;
    using abstract_type = any_symmetric<T>;
    using concrete_type = symmetric<T>;
};

template<typename T>
struct symmetric : any_symmetric<T> {
    void concrete() override {
    }
};

The concrete virtual function is just a way of making the any_ classes abstract. Each class contain two aliases: abstract_type, the nearest abstract class; and concrete_type, the nearest concrete class. Their purpose will become clear later in the tutorial.

Operations on Matrices

We will implement three operations: transposition, addition, and scaling (i.e. the multiplication of a matrix by a scalar). They provide a good selection of cases, on which other operations can be modeled.

We want matrix operations to return a result of the most specialized type. Transposing a matrix should return a matrix of the same category, and the matrix itself if it is symmetric. Adding two symmetric matrices should yield a symmetric matrix, not just a square matrix. Adding a symmetric matrix and a square matrix should return a square matrix; etc. Thus the operations must be specialized, depending on the type of the arguments.

At this point, as library designers, we are facing a difficult choice: specialize the operations at compile time (using templates), or at runtime (using virtual functions or open methods). Both ways have pros and cons. The compile-time approach results in (slightly) faster operations, and, more importantly, better interfaces. For example, square matrices have a determinant, which ordinary matrices do not have. Also, type errors can be diagnosed by the compiler. However, compile-time specialization is possible only if the types of the matrices is known in advance. It may not always be the case.

Thus, as library designers, we face a dilemma: which users to serve best?

Here is the good news: combining open methods and templates makes it possible for the user to choose between full static typing and full type erasure, and any level of erasure between the two extremes.

Static Operations

Let’s start by implementing the three operations in a fully typed manner:

// transposition

template<typename T>
auto operator~(const handle<ordinary<T>>& m) {
    return handle(new ordinary<T>(...));
}

template<typename T>
auto operator~(const handle<square<T>>& m) {
    return handle(new square<T>(...));
}

template<typename T>
auto operator~(const handle<symmetric<T>>& m) {
    return m;
}

// addition

template<typename>
struct is_matrix_aux : std::false_type {};

template<template<typename> typename M, typename T>
struct is_matrix_aux<M<T>> : std::is_base_of<any<T>, M<T>> {};

template<typename T>
constexpr bool is_matrix = is_matrix_aux<T>::value;

template<typename T, typename U>
auto operator+(const handle<ordinary<T>>& a, const handle<ordinary<U>>& b) {
    return handle(new ordinary<std::common_type_t<T, U>>(...));
}

template<typename T, typename U, std::enable_if_t<is_matrix<U>, int>...>
auto operator+(const handle<ordinary<T>>& a, const handle<U>& b) {
    return handle(
        new ordinary<std::common_type_t<T, typename U::element_type>>(...));
}

template<typename T, typename U, std::enable_if_t<is_matrix<U>, int>...>
auto operator+(const handle<U>& a, const handle<ordinary<T>>& b) {
    return handle(
        new ordinary<std::common_type_t<T, typename U::element_type>>(...));
}

template<typename T, typename U>
auto operator+(const handle<square<T>>& a, const handle<square<U>>& b) {
    return handle(new square<std::common_type_t<T, U>>(...));
}

template<typename T, typename U>
auto operator+(const handle<symmetric<T>>& a, const handle<symmetric<U>>& b) {
    return handle(new symmetric<std::common_type_t<T, U>>(...));
}

template<typename T, typename U>
auto operator+(const handle<symmetric<T>>& a, const handle<square<U>>& b) {
    return handle(new square<std::common_type_t<T, U>>(...));
}

template<typename T, typename U>
auto operator+(const handle<square<T>>& a, const handle<symmetric<U>>& b) {
    return handle(new square<std::common_type_t<T, U>>(...));
}

// scaling

template<typename T, typename U>
auto operator*(T a, const handle<ordinary<U>>& b) {
    return handle(new ordinary<std::common_type_t<T, U>>(...));
}

template<typename T, typename U>
auto operator*(T a, const handle<square<U>>& b) {
    return handle(new square<std::common_type_t<T, U>>(...));
}

template<typename T, typename U>
auto operator*(T a, const handle<symmetric<U>>& b) {
    return handle(new symmetric<std::common_type_t<T, U>>(...));
}

Let’s exercise the operators with a test case:

BOOST_AUTO_TEST_CASE(test_static_operations) {
    {
        handle<ordinary<double>> o(new ordinary<double>);
        handle<ordinary<double>> ot = ~o;

        handle<square<double>> sq(new square<double>);
        handle<square<double>> tsq = ~sq;

        handle<symmetric<double>> sy(new symmetric<double>);
        handle<symmetric<double>> tsy = ~sy;
        BOOST_TEST(sy.get() == tsy.get());
    }

    {
        handle<ordinary<int>> a(new ordinary<int>);
        handle<ordinary<double>> b(new ordinary<double>);
        handle<ordinary<double>> p = a + b;
    }

    {
        handle<square<double>> a, b;
        handle<square<double>> p = a + b;
    }

    {
        handle<symmetric<double>> a(new symmetric<double>);
        handle<square<double>> b(new square<double>);
        handle<symmetric<double>> aa = a + a;
        handle<square<double>> ab = a + b;
        // handle<symmetric<double>> x = a + b; // wrong: square is-not-a symmetric
        // handle<square<double>> y = a + a;    // wrong: symmetric is-not-a square
        handle<any_square<double>> y =
            a + a; // OK: symmetric is-a any-square - but useless
    }

    {
        handle<ordinary<int>> m(new ordinary<int>);
        handle<ordinary<double>> p = 2. * m;
    }

    {
        handle<square<int>> m(new square<int>);
        handle<square<double>> p = 2. * m;
    }

    {
        handle<symmetric<int>> m(new symmetric<int>);
        handle<symmetric<double>> p = 2. * m;
    }
}

Note that a symmetric matrix cannot be converted to a square matrix. This is normal, because they don’t store the elements in the same way: symmetric stores only half of the elements. A symmetric can be converted to abstract class any_square, but, at this point, this is useless, because the operators are not available for the abstract classes.

Polymorphic Transposition

Let’s change this. First let’s implement transposition for abstract matrix classes.

struct YOMM2_SYMBOL(transpose);

template<typename Matrix>
using transpose =
    method<YOMM2_SYMBOL(transpose), handle<Matrix>(virtual_<Matrix&>)>;

template<
    template<typename> typename Matrix, typename T,
    std::enable_if_t<
        std::is_base_of_v<any<T>, Matrix<T>>  // 1
            && std::is_abstract_v<Matrix<T>>, // 2
        int>...>
auto operator~(const handle<Matrix<T>>& m) {
    return transpose<Matrix<T>>::fn(*m.get());
}

We must make sure that the new operator~ is generated only if the deduced Matrix template satisfy two conditions:

  1. Matrix is actually one of our matrix templates; without (1), any template would match.
  2. Matrix<T> is abstract; without (2), we would have an ambiguity with our own operator~ defined above for concrete matrix types.

C++20 concepts offer a much cleaner solution for taming template instantiations.

Now we must generate definitions for the transpose<Matrix<T>> methods. In a full implementation of a matrix libray, we would also for other, similar methods, e.g. negate. For that we use a templatized definition container called unary_definition. We want to instantiate it for every combination of method template, abstract base class, concrete implementation of the abstract class, and this for the required numeric types. E.g.:

On the other hand, we do not want to generate definitions for combinations that are useless, e.g.:

We can automate this by forming the Cartesian product M x A x C x T, where:

M = { unary method templates: transpose, negate, etc }
A = { abstract class templates }
C = { concrete class templates }
T = { required numeric types }

…and selecting the combinations that satisfy the condition:

This time, we need to make a Cartesian product that involves templates as well as types. For this, YOMM2 provides two constructs:

We can now implement unary_definition:

template<
    typename Method, typename Abstract, typename Concrete, typename T,
    typename = std::bool_constant<std::is_base_of_v<
        typename Abstract::template fn<T>,
        typename Concrete::template fn<T>>> // see note 1
    >
struct unary_definition : not_defined {}; // see note 2

template<
    template<typename> typename Abstract, template<typename> typename Concrete,
    typename T>
struct unary_definition<
    template_<transpose>, // see note 3
    template_<Abstract>, template_<Concrete>, T,
    std::true_type // see note 4
    > {
    using method = transpose<Abstract<T>>;
    static auto fn(Concrete<T>& m) {
        return ~handle<Concrete<T>>(&m); // see note 5
    }
};
  1. The default template argument evaluates to std::true_type if Concrete<T> derives from Abstract<T>, and std::false_type otherwise.
  2. In the general case, unary_definition derives from not_defined, meaning that it will be filtered out by use_definitions.
  3. Specialize for method template transpose.
  4. Specialize for the cases where the concrete class derives from the abstract class.
  5. fn falls back on operator~ for concrete classes.

Let’s create a few convenient aliases for templates lists:

using abstract_matrix_templates = templates<any, any_square, any_symmetric>;
using concrete_matrix_templates = templates<ordinary, square, symmetric>;

We will also need a list of all the templates - abstract and concrete. For that we can use boost::mp11::mp_append:

using matrix_templates = boost::mp11::mp_append<
    abstract_matrix_templates, concrete_matrix_templates>;

We can now use product to create all the combinations, and pass them to use_definitions. Let’s do this for int and double matrices:

use_definitions<
    unary_definition,
    product<
        templates<transpose>, abstract_matrix_templates,
        concrete_matrix_templates, types<double, int>>>
    YOMM2_GENSYM;

We must not forget to register all the matrix classes, for both numeric types. For that, we can use another YOMM2 helper: apply_product. It takes a templates list, and any number of types list; forms the Cartesian product; and, for each resulting combination, applies the first element - a template - to the other elements. The result is a types list of template instantiations.

We saw that use_classes takes a list of classes. It also works with a list of types lists. We can thus inject the result of apply_product into use_classes:

YOMM2_STATIC(use_classes<apply_product<matrix_templates, types<double, int>>>);
BOOST_AUTO_TEST_CASE(test_dynamic_transpose) {
    update();

    handle<any<double>> o(new ordinary<double>);
    handle<any<double>> ot = ~o;
    BOOST_TEST(ot.type() == typeid(ordinary<double>));

    handle<any<double>> sq(new square<double>);
    handle<any<double>> tsq = ~sq;
    BOOST_TEST(tsq.type() == typeid(square<double>));

    handle<any_square<double>> sy(new symmetric<double>);
    handle<any_square<double>> tsy = ~sy;
    BOOST_TEST(sy.get() == tsy.get());
}

Polymorphic Addition

Binary operations are a little more difficult to implement because the return type depends on the type of the arguments. For example, adding a any_symmetric<double> and a square<double> should yield a any_square<double>.

Fortunately, we can leverage the static operator+ to deduce the return type of additions that involve one or two abstract types. This is what the nested abstract_type and concrete_type aliases are for. We can convert the types of the arguments to their nearest concrete type: for a any_symmetric<double>, it is a symmetric<double>; for a square<double>, it is square<double>. Then we can “pretend” to add two handles to objects of the two types, and extract the result’s type using decltype. And finally, we convert that type to its abstract base class, using abstract_type.

template<typename M1, typename M2>
using binary_result_type =
    typename decltype(std::declval<const handle<typename M1::concrete_type>>() + std::declval<const handle<typename M2::concrete_type>>())::
        body_type::abstract_type;

static_assert(std::is_same_v<
              binary_result_type<any_symmetric<double>, any_square<int>>,
              any_square<double>>);

static_assert(std::is_same_v<
              binary_result_type<any_symmetric<double>, square<int>>,
              any_square<double>>);

// etc

We can now declare the add method:

struct YOMM2_SYMBOL(add);

template<typename M1, typename M2>
using add = method<
    YOMM2_SYMBOL(add),
    handle<binary_result_type<M1, M2>>(virtual_<M1&>, virtual_<M2&>)>;

template<
    typename M1, typename M2,
    typename = std::enable_if_t<
        is_matrix<M1> && is_matrix<M2> &&
            (std::is_abstract_v<M1> || std::is_abstract_v<M2>),
        void>>
auto operator+(const handle<M1>& a, const handle<M2>& b) {
    return add<typename M1::abstract_type, typename M2::abstract_type>::fn(
        *a, *b);
}

Let’s generate definitions for the add<T> methods, and other binary methods (multiply, etc). We use the same approach as for transposition, but this time we have two parameters. We use another container definition, binary_definition; we instantiate it for (e.g.):

We can automate this by forming the Cartesian product M x A1 x C1 x T1 x A2 x C2 x T2, where:

M = { unary method templates: transpose, negate, etc }
A1, A2 = { abstract class templates }
C1, C2 = { concrete class templates }
T1, T2 = { required numeric types }

…and selecting the combinations that satisfy the condition:

Since the condition is a little more complex, and we need it for several methods, let’s factorize it:

template<
    typename A1, typename C1, typename T1, typename A2, typename C2,
    typename T2>

struct enable_binary_definition;

template<
    template<typename> typename A1, template<typename> typename C1, typename T1,
    template<typename> typename A2, template<typename> typename C2, typename T2>
struct enable_binary_definition<
    template_<A1>, template_<C1>, T1, template_<A2>, template_<C2>, T2>
    : std::bool_constant<
          std::is_base_of_v<A1<T1>, C1<T1>> &&
          std::is_base_of_v<A2<T2>, C2<T2>> &&
          std::is_base_of_v<
              binary_result_type<A1<T1>, A2<T2>>,
              binary_result_type<C1<T1>, C2<T2>>>> {};

We can now implement binary_defintion:

template<
    typename Method, typename A1, typename C1, typename T1, typename A2,
    typename C2, typename T2,
    typename = typename enable_binary_definition<A1, C1, T1, A2, C2, T2>::type>
struct binary_definition : not_defined {};

template<
    template<typename> typename A1, template<typename> typename D1, typename T1,
    template<typename> typename A2, template<typename> typename D2, typename T2>
struct binary_definition<
    template_<add>, template_<A1>, template_<D1>, T1, template_<A2>,
    template_<D2>, T2, std::true_type> {
    using method = add<A1<T1>, A2<T2>>;
    static auto fn(D1<T1>& a, D2<T2>& b) {
        return handle<D1<T1>>(&a) + handle<D2<T2>>(&b);
    }
};

…and use it to generate definitions for polymorphic addition of matrices of int and double:

use_definitions<
    binary_definition,
    product<
        templates<add>, abstract_matrix_templates, concrete_matrix_templates,
        types<double, int>, abstract_matrix_templates,
        concrete_matrix_templates, types<double, int>>>
    YOMM2_GENSYM;
BOOST_AUTO_TEST_CASE(test_dynamic_operations) {
    update();

    {
        handle<any<int>> a(new ordinary<int>);
        handle<any<double>> b(new ordinary<double>);
        handle<any<double>> s = a + b;
        BOOST_TEST(s.type() == typeid(ordinary<double>));
    }

    {
        handle<any_square<double>> a(new symmetric<double>);
        handle<any_square<double>> b(new square<double>);
        handle<any_square<double>> s = a + b;
        BOOST_TEST(s.type() == typeid(square<double>));
    }

    {
        handle<symmetric<double>> a(new symmetric<double>);
        handle<any_square<int>> b(new symmetric<int>);
        handle<any_square<double>> s = a + b;
        BOOST_TEST(s.type() == typeid(symmetric<double>));
    }
}

Now we need to create a mechanism that enable users to instantiate just the method definitions they want. This is more complicated than for the vector library, because users need to be able to select from three sets: the underlying numeric types, the matrix types, and the methods. For example:

template<typename>
struct template_of;

template<template<typename> typename M, typename T>
struct template_of<M<T>> {
    using type = template_<M>;
};

static_assert(
    std::is_same_v<
        template_of<square<int>::abstract_type>::type, template_<any_square>>);

As for the methods, different method definitions need different Cartesian products. We can handle this with traits:

template<template<typename...> typename Definition>
struct unary_definition_traits {
    template<
        typename AbstractMatrixTemplates, typename ConcreteMatrixTemplates,
        typename NumericTypes>
    using fn = use_definitions<
        unary_definition,
        product<
            templates<Definition>, AbstractMatrixTemplates,
            ConcreteMatrixTemplates, NumericTypes>>;
};

template<template<typename...> typename Definition>
struct binary_definition_traits {
    template<
        typename AbstractMatrixTemplates, typename ConcreteMatrixTemplates,
        typename NumericTypes>
    using fn = use_definitions<
        binary_definition,
        product<
            templates<Definition>, AbstractMatrixTemplates,
            ConcreteMatrixTemplates, NumericTypes, AbstractMatrixTemplates,
            ConcreteMatrixTemplates, NumericTypes>>;
};

template<template<typename...> typename Definition>
struct definition_traits;

template<>
struct definition_traits<transpose> : unary_definition_traits<transpose> {};

template<>
struct definition_traits<add> : binary_definition_traits<add> {};

We can now write use_polymorphic_matrices and its nested templates:

template<template<typename> typename... Ms>
struct use_polymorphic_matrices {
    using abstract_matrix_templates = types<
        typename template_of<typename Ms<double>::abstract_type>::type...>;
    using concrete_matrix_templates = templates<Ms...>;
    template<typename... Ts>
    struct of {
        using numeric_types = types<Ts...>;
        template<template<typename...> typename... Ops>
        struct with
            : use_classes<
                  apply_product<abstract_matrix_templates, types<Ts...>>,
                  apply_product<concrete_matrix_templates, types<Ts...>>>,
              definition_traits<Ops>::template fn<
                  abstract_matrix_templates, concrete_matrix_templates,
                  numeric_types>... {};
        with<transpose, add> all;
    };

    of<double> all;

    template<template<typename...> typename... Ops>
    struct with : of<double>::template with<Ops...> {};
};

template<>
struct use_polymorphic_matrices<>
    : use_polymorphic_matrices<ordinary, square, symmetric> {};