Journal [C++14 ] Expressions template pour les nuls

Posté par  (site web personnel) . Licence CC By‑SA.
Étiquettes :
45
31
mai
2016

Sommaire

Expression Templates pour les nuls

Après un contact prolongé avec Joël Falcou, et pas mal de nœuds aux cerveaux pour arriver à émuler le comportement du broadcasting de Numpy avec des expression templates dans Pythran, j'ai eu l'envie soudaine, un peu folle peut-être, de réécrire un moteur d'expressions template en C++14, mais pour faire simple et didactique.

Mais de quoi zy parle

Les expressions templates sont une maintenant assez vieille technique de C++ qui permet par exemple d'éviter de créer des objets intermédiaires lourds quand on calcule sur des tableaux en utilisant des expressions complexes.

Par exemple, pour le type std::vector<double>, si on surcharge naïvement (et de façon fort risquée, pas touche aux conteneurs de la lib standard comme ça) les opérateurs +, * etc avec un code du genre :

std::vector<double> operator+(std::vector<double> const& left, std::vector<double> const& right) {
    assert(left.size() == right.size() && "same size");
    std::vector<double> out(left.size());
    std::transform(left.begin(), left.end(), right.begin(), out.begin(),
                   std::plus<double>());
    return out;
}

et bien l'évaluation d'une expression du genre a + b * c va créer un tableau inutile (le résultat de b * c) et la move semantic ne nous sauve pas totalement. Et on va faire deux boucles, une par opération.

Avec les expressions template, on va plutôt écrire :

add<std::vector<double>, std::vector<double>> operator+(std::vector<double> const& left, std::vector<double> const& right) {
    return {left, right};
}

Avec :

template<typename L, typename R>
struct add {
    L const& left_;
    R const& right_;
    add(L const& left, R const& right) : left-(left), right_(right) {}
    // ...
};

Ce qui nous donnera comme type de retour de l'expression a + b * c :

add<std::vector<double>,
    mul<std::vector<double>,
        std::vector<double>
       >
   >

En ajoutant un opérateur [] aux types add<...> et mul<...> qui fait suivre l'appel aux fils, et une fonction membre size() qui renvoie, par exemple, la plus petite des tailles des deux fils :

template<typename L, typename R>
struct add {
    L const& left_;
    R const& right_;
    add(L const& left, R const& right) : left-(left), right_(right) {}
    auto operator[](size_t i) const { return  left_[i] + right_[i]; }
    auto size() const { return std::min(left_.size(), right_.size(); }
};

on peut alors écrire le code bâtard mais néanmoins relativement facile à comprendre :

auto expr = a + b * c;
auto n = expr.size();
std::vector<double> res(n);
for(size_t i = 0; i < n; ++i)
  res[i] = expr[i];

Et voilà, on a une unique boucle, pas de tableau intermédiaire, c'est la fête et vive les expressions template. Cette technique est très largement utilisée :

Mais elle est assez lourde à mettre en œuvre avec les outils pre-C++11. Je vous propose une solution que je croyais innovante, mais après discussion avec mon mentor sus-nommé, elle n'est qu'ingénieuse, ce qui est somme toute déjà pas mal.

Expressions template à la sauce C++14

Le parcours de l'arbre d'expression de bas en haut à coup d'appels récursifs, ça rappelle fury le design pattern visitor. C'est d'ailleurs une des options proposées par Boost.Proto, le gros pavé pour ceux qui aiment manger des expressions template avec de la crème le matin.

Du coup, proposition didactique (pleins de détails et optimisations possibles sont écartés pour la clarté du propos) : on va associer à chaque type de nœud de notre arbre un tag (une structure vide, un type) et une fonction générique dont le but sera de lancer la mécanique de parcours de l'arbre. Exemple :

struct add_tag {};
template<class A, class B>
auto add(A&& a, B&& b) {
  return [=](auto visitor) { return visitor(add_tag{}, a(visitor), b(visitor));};
  }

Bon ça parait un peu mystérieux comme ça, le gain en clarté n'est pas (clair ?) évident, mais on renvoie juste, quand on fait appel à la fonction add, une fonction générique capable d'appliquer un visitor sur les fils, de bas en haut.

On peut faire la même chose pour les feuilles de l'arbre, ici en distinguant les références vers des objets complexes et les constantes :

struct cst_tag {};
template<class T>
auto cst(T expr) {
 return [=](auto visitor) { return visitor(cst_tag{}, expr); };
}

struct ref_tag {};
template<class T>
auto ref(T& expr) {
 return [&](auto visitor) { return visitor(ref_tag{}, expr); };
}

Ça fait peu ou prou la même chose, mais avec un tag différent.

Maintenant que la mécanique est en place, on peut code de façon élégante nos visiteurs ! Par exemple pour obtenir le i-ème élément:

struct evaluator {
  evaluator(size_t i) : i_(i) {}

  template<class T>
  auto operator()(lazy::cst_tag, T c) { return c; }

  template<class T>
  auto operator()(lazy::ref_tag, T& r) { return r[i_]; }

  template<class A, class B>
  auto operator()(lazy::mul_tag, A a, B b) { return a * b; }

  template<class A, class B>
  auto operator()(lazy::add_tag, A a, B b) { return a + b; }

  private:
  size_t i_;
};

On utilise le tag pour spécialiser l'appel et déterminer le nœud à visiter. L'appel se fait simplement par :

auto expr = add(cst(12), mul(ref(a), ref(b)));
auto expr_0 = expr(evaluator(0));

On peut même spécifier un comportement par défaut, en jouant sur l'ordre de résolution des surcharges :

struct size {

  template<class T>
  auto operator()(lazy::cst_tag, T c) { return std::numeric_limits<size_t>::max(); }

  template<class T>
  auto operator()(lazy::ref_tag, T& r) { return r.size(); }

  template<class T, class A, class B>
  auto operator()(T, A a, B b) { return std::min(a,  b); }

};

À la différence de la première approche, plus besoin de modifier les types proxy add, mul etc. On peut mettre toute la logique à un endroit, et l'étendre de manière non intrusive. Au final, si on sait vouloir stocker notre expression finale dans un std::vector<...>, on peut même écrire de façon assez élégante :

template<class E>
auto eval(E const & expr) {
  size_t n = expr(size());
  std::vector<decltype(expr(evaluator(0)))> res(n);
  for(size_t i = 0; i < n; ++i)
    res[i] = expr(evaluator(i));
  return res;
}

qui fait le bonheur des petits et des grands :

auto expr = add(cst(12), mul(ref(a), ref(b)));
auto evaluated = eval(expr);

Performance

Un petit coup de clang++ -std=c++14 -Ofast -march=native sur un code utilisant l'expression précédente a fait apparaitre cette séquence d'assembleur:

400ab0:       c4 c1 7a 6f 0c 1e       vmovdqu (%r14,%rbx,1),%xmm1
400ab6:       c4 c1 7a 6f 54 1e 10    vmovdqu 0x10(%r14,%rbx,1),%xmm2
400abd:       c4 c1 7a 6f 5c 1e 20    vmovdqu 0x20(%r14,%rbx,1),%xmm3
400ac4:       c4 c1 7a 6f 64 1e 30    vmovdqu 0x30(%r14,%rbx,1),%xmm4
400acb:       c4 c2 71 40 0c 1f       vpmulld (%r15,%rbx,1),%xmm1,%xmm1
400ad1:       c4 c2 69 40 54 1f 10    vpmulld 0x10(%r15,%rbx,1),%xmm2,%xmm2
400ad8:       c4 c2 61 40 5c 1f 20    vpmulld 0x20(%r15,%rbx,1),%xmm3,%xmm3
400adf:       c4 c2 59 40 64 1f 30    vpmulld 0x30(%r15,%rbx,1),%xmm4,%xmm4
400ae6:       c5 f1 fe c8             vpaddd %xmm0,%xmm1,%xmm1
400aea:       c5 e9 fe d0             vpaddd %xmm0,%xmm2,%xmm2
400aee:       c5 e1 fe d8             vpaddd %xmm0,%xmm3,%xmm3
400af2:       c5 d9 fe e0             vpaddd %xmm0,%xmm4,%xmm4
400af6:       c4 c1 7a 7f 4c 1d 00    vmovdqu %xmm1,0x0(%r13,%rbx,1)
400afd:       c4 c1 7a 7f 54 1d 10    vmovdqu %xmm2,0x10(%r13,%rbx,1)
400b04:       c4 c1 7a 7f 5c 1d 20    vmovdqu %xmm3,0x20(%r13,%rbx,1)
400b0b:       c4 c1 7a 7f 64 1d 30    vmovdqu %xmm4,0x30(%r13,%rbx,1)
400b12:       48 83 c3 40             add    $0x40,%rbx
400b16:       48 81 fb 00 03 00 00    cmp    $0x300,%rbx
400b1d:       75 91                   jne    400ab0 <main+0x70>

Ce qui est franchement cool : le code a été vectorisé, déroulé et il est… vachement propre. On retrouve bien là le principe de costless abstraction si cher au C++. J'♥.

Répéter, c'est apprendre

La petite lib que j'ai écrite pour tester ces idées n'ira pas sur le grand ninternet, c'était juste pour le lulz. Et puis après coup, Il m'a pointé vers https://isocpp.org/blog/2016/05/cppcon-2015-expression-templates-past-present-future-joel-falcou qui contient bien plus de contenu que ce maigre article, qui ne sera finalement qu'une petite introduction et qui aura eu l'avantage de me forcer de coucher sur le clavier (avant d'aller moi-même me coucher) mes idées. Et n'est-ce pas ça, le but d'un journal ?

  • # Commentaire supprimé

    Posté par  . Évalué à 9.

    Ce commentaire a été supprimé par l’équipe de modération.

    • [^] # Re: Avantage ?

      Posté par  (site web personnel) . Évalué à 7.

      Dans ton exemple c'est évident, mais imagine maintenant un pipeline avec une centaines d'opérations de haut niveau elles même définies par une dizaine d'opérations de bas niveau, tout de suite c'est moins facile à inliner à la main, regarde cet exemple :

      auto preprocess = [](const auto &im){ return malib::gammaCorrect(malib::tonemap(malib::greyscale(im))); };

auto image = malib::max(malib::absdiff(preprocess(im0), preprocess(im1)));

      Un autre example intéressant, imagine que tu ai un pipeline complex et que tu ne veuilles qu'une sous partie de l'image :

      
auto image = malib::regionOfInterest(rectangle, pipelineComplex(im0, im1, im2, im3, im4);

      Tu ne veux pas calculer toute l'image alors que le regionOfInterest limite le calcul à une petite sous partie. Cette technique peut te permettre de compiler un code qui n'effectue le traitement que sur sa sous partie SANS copier une seule image. Note que tu peux aussi spécialiser les templates pour générer du code vectorisé efficace ou des accès mémoire plus sympathiques. (Remarque, ce problème ne se simplifie pas en distribuant regionOfInterest sur les images d'entrée car ton pipelineComplex peut effectuer un traitement global)

    • [^] # Re: Avantage ?

      Posté par  (site web personnel) . Évalué à 3. Dernière modification le 31 mai 2016 à 13:46.

      Je pense que c'est aussi un problème d'API, ou plutôt de compromis API/perf.

      Prenons le cas de GMP, il est chouette. Tu as trois entiers multi-précision a b et c. Tu veux calculer rapidement a + b * c. Tu peux demander à l'utilisateur d'écrire (et ce le cas en C d'ailleurs), mpz_mul(tmp, b, c); mpz_add(out, a, tmp). Le problème c'est que (en plus d'être verbeux) tu introduis un temporaire, et que peut-être qu'il existe une façon de faire un fused multiply add plus efficacement qu'en le cassant en une multiplication et une addition, comme pour la boucle que tu cites où la fusion de boucle a permis, entre autres, de ne pas parcourir deux fois les tableaux. Et tu ne peux pas te permettre d'écrire toutes les fonction add_muladd_add_mul` etc pour chaque optimisation possible.

      Les expression templates en général permettent de résoudre les deux soucis : au niveau API, l'utilisateur écris juste a + b * c et au niveau perf la machinerie (souvent bien plus complexe que ce qui est présentée ici, et bien mieux faite aussi) se charge de l'évaluation efficace de l'expression.

      Et tu as raison, l'évaluation paresseuse fait parti des effets de bord sympa ;-)

    • [^] # Re: Avantage ?

      Posté par  (site web personnel) . Évalué à 2.

      Il doit y en avoir. Peut être est-ce la simplicité de l'exemple qui ne me permet pas de le voir. Eventuellement du lazy loading (l'expression n'est évaluée que lorsqu'elle est nécessaire) ?

      Imagine 5 secondes que a, b et c sont trois matrix de trés grande tailles et que tu enchaînes des opérations beaucoup plus complexe qu'une simple multiplication et soustraction.

      L'ordre dans lequel tu appliques tes opérations (suivant la taille de chacune de tes matrix) peut avoir un impact énorme sur le temps de calcul nécessaire.

      De même pour la nécessité à copier la mémoire ou non, sur des matrix trés grande, ou lorsque tu as besoin de transposition.

      Ce sont des optimisations que presque l’intégralité des bibliothèques désignés pour faire de l'algèbre linéaire mettre en oeuvre en C++: Eigen, Armadillo, blaze, etc, etc.

    • [^] # Re: Avantage ?

      Posté par  (site web personnel) . Évalué à 2.

      L'avantage est que le plus souvent, le vice est poussé au maximum, et la boucle peut-être supprimé totalement :

      Vector a, c;
Matrix b;
Vector d = a + b*c;

      Si tu veux voir la différence, compare l'API d'Eigen (http://eigen.tuxfamily.org/dox/group__QuickRefPage.html) avec celle de BLAS (http://www.netlib.org/blas/blasqr.pdf). Les expressions templates permet d'adapter les calculs effectués au mieux.

    • [^] # Commentaire supprimé

      Posté par  . Évalué à 3.

      Ce commentaire a été supprimé par l’équipe de modération.

  • # C++, où quand le code asm généré est plus lisible

    Posté par  . Évalué à 10.

    Peut-être que je n'ai pas assez fait de C++ pour voir à quel point c'est génial, mais j'ai l'impression que le code assembleur que tu présentes est plus simple que le code C++. L'expressivité d'un langage ne lui permet-il pas justement de simplifier sa compréhension et d'exprimer de manière plus synthétique et plus sémantique un problème donné ? Avec le C++ j'ai l'impression qu'on répond à un problème technique par une complexité au niveau du langage, justifiée par la résolution d'un problème technique.

    Le résultat est très verbeux, peut-on faire plus simple ? Comment se gère ce problème dans les autres langages ?

    • [^] # Re: C++, où quand le code asm généré est plus lisible

      Posté par  (site web personnel) . Évalué à 2.

      D'autres langages vont proposer soit de résoudre le problème au runtime en créant l'arbre d'opération à l’exécution et en générant du code avec un compilateur intégré (JIT), c'est ce que font pas mal d'approches au dessus de numpy, numba, … Tu peux aussi générer du code OpenCL à la volée. Ces approches ont l'avantage de pouvoir en plus avoir une vision sur la taille de la structure, ses caractéristiques et la machine utilisée (Quantité de ram, taille des caches, …), et ainsi pouvoir générer un code plus adapté. En contrepartie c'est encore plus complexe et cela a un coût lors de l’exécution ;)

      D'autres langages vont proposer un outil de méta-programmation lors de la compilation qui est "plus lisible" que les templates C++.

      D'autres langages ont le support de ce type d'optimisation en natif dans le compilateur (regarde Haskell, les fusion et les rewrite rules), mais à ma connaissance cela ne donne pas ce niveau de satisfaction.

      Au final on est dans un domaine de la performance à tout prix et c'est dur, cela fait du code crade et seulement maintenable par des gars qui ont 3 PhD, dont un consacré à la maintenance du code en question, mais bon, c'est ça qui est marrant non ? (Il faut voir que l'exemple donné dans ce journal est "trivial" et "très lisible"… si si… ;(

    • [^] # Re: C++, où quand le code asm généré est plus lisible

      Posté par  (site web personnel) . Évalué à 1.

      "Le résultat est très verbeux, peut-on faire plus simple ? Comment se gère ce problème dans les autres langages ?"

      C'est un problème de meta-programmation. Tu as une forme de code qui tu veux transformer dans une autre forme qui se compile bien. Cette transformation se fait souvent au runtime (tensorflow).

      Je ne connais pas d'autres langage qui gère correctement la metaprogrammation qui fonctionne à la compilation. Il me semble que Lisp le fait, et Rust avec ses "macro".

      "La première sécurité est la liberté"

      • [^] # Re: C++, où quand le code asm généré est plus lisible

        Posté par  . Évalué à 3.

        Le langage de Jonathan Blow pourra sans doute gérer ça de manière encore plus lisible. La génération de code compile-time se fait avec le même langage.

        C'est très long, mais c'est trèèèès intéressant:

        https://www.youtube.com/watch?v=2IBr0XZOPsk

        Dans la vidéo ci-dessus il crée différentes classes Matrix (4x4) à l'aide d'une fonction template qui prend en paramètre un flag qui va définir comment les différentes structures de matrices vont se multiplier.
        Dans un cas concret on peut utiliser différentes matrices 4x4 et savoir à l'avance quels nombres ne vont jamais se multiplier (car toujours à zéro).
        En résumé il crée différentes classes de matrice ayant des fonctions "mult(matrice1, matrice2)" différentes avec le même code.

        Je me suis toujours demandé pourquoi on ne pouvait pas générer de code avec son code (sans avoir a utiliser un outils externe). Et bien, Jonathan Blow aussi et a décider de ne rien attendre du C/C++ pour ça (ce qui semble assez raisonnable vu l'exemple de ce journal).

      • [^] # Re: C++, où quand le code asm généré est plus lisible

        Posté par  . Évalué à 1.

        Je ne connais pas d'autres langage qui gère correctement la metaprogrammation qui fonctionne à la compilation. Il me semble que Lisp le fait, et Rust avec ses "macro".

        D le fait aussi. C’est même sa principale force par rapport à C++ : il le fait avec une syntaxe beaucoup plus facile à appréhender.

        Mes commentaires sont en wtfpl. Une licence sur les commentaires, sérieux ? o_0

    • [^] # Re: C++, où quand le code asm généré est plus lisible

      Posté par  (site web personnel) . Évalué à 2. Dernière modification le 01 juin 2016 à 10:38.

      Peut-être que je n'ai pas assez fait de C++ pour voir à quel point c'est génial, mais j'ai l'impression que le code assembleur que tu présentes est plus simple que le code C++.

      Parce que le code assembleur que tu vois en bas n'est pas représentatif du code templated que tu vois plus haut, ni même comparable.

      Le code C++ templated que tu vois plus haut est valide pour une infinité de combinaison de mul / add et pour une infinité d'expressions différente.

      Alors que le code assembler que tu vois plus bas est valide uniquement pour la simple expression "auto expr = add(cst(12), mul(ref(a), ref(b)));" avec des paramètres donnés.

      Si tu gardes ça à l'esprit, alors non le code C++ n'est pas verbeux, il est même extrêmement concis.
      Ce que tu as écris en C++ se rapproche au final d'un un "micro-compiler" d'expressions mathématiques. Alors que ce que tu as en assembler est juste une simple opération mathématique sur une range.

      C'est là toute la magie de la méta-programmation en C++, faire tout ce que tu peux faire à compile time à compile time pour générer un runtime aussi léger que possible ( en théorie du moins ).
      ```

  • # La vraie nouvelle astuce de C++14

    Posté par  . Évalué à 3.

    La vraie nouvelle astuce pour les expressions templates et C++14, c'est effectivement d'utiliser les lambdas polymorphes avec capture, comme tu le fais ici.
    Par contre moi je ferais pas la récursion dedans, c'est moins générique.

    J'ai personnellement commencé à expérimenter avec ce principe pour une nouvelle bibliothèque pour la différentiation automatique d'ordre supérieur.
    Par contre malheureusement retourner une lambda ne suffit pas; on veut pouvoir surcharger les opérateurs, donc il faut enrober la lambda dans quelque chose (moi en fait j'ai fait l'inverse, j'ai implémenté un tuple avec une lambda).

  • # Compatible avec auto ?

    Posté par  (site web personnel) . Évalué à 0.

    Si je comprends bien, avec une lambda vous solutionnez le problème des références pendouillantes que l'on a avec les libs usuelles comme Eigen sur

    auto temp = a * x + b;
vect r = temp * 42 + y;

    J'avais tenté de résoudre ça avec une approche bien plus tordue.

  • # Petites précisions

    Posté par  . Évalué à 1.

    Superbe dépêche… Pour bien le comprendre, j'ai une petite question (mon C++ est tout rouillé et mes design patterns balbutiants): dans l'exemple

    struct cst_tag {};
template<class T>
auto cst(T expr) {
 return [=](auto visitor) { return visitor(cst_tag{}, expr); };
}

    qu'est-ce que fait cst_tag{}? Si je me reporte au livre du GoF, ça sert à indiquer le type du nœud donc c'est lié aux différentes surcharges de operator () dans evaluator… Mais comme je ne sais pas ce que fait lazy, j'ai du mal à tout comprendre…

Suivre le flux des commentaires

Note : les commentaires appartiennent à celles et ceux qui les ont postés. Nous n’en sommes pas responsables.