The Higher Education and Research forge

Home My Page Projects Code Snippets Project Openings Garamon
Summary Activity SCM

SCM Repository

1 // Copyright (c) 2018 by University Paris-Est Marne-la-Vallee
2 // Mvec.hpp
3 // This file is part of the Garamon for project_namespace.
4 // Authors: Stephane Breuils and Vincent Nozick
5 // Contact: vincent.nozick@u-pem.fr
6 //
7 // Licence MIT
8 // A a copy of the MIT License is given along with this program
10 /// \file Mvec.hpp
11 /// \author Stephane Breuils, Vincent Nozick
12 /// \brief Class to define a multivector and its basic operators in the Geometric algebra of project_namespace.
15 // Anti-doublon
16 #ifndef project_inclusion_guard
17 #define project_inclusion_guard
18 #pragma once
20 // External Includes
21 #include <Eigen/Core>
22 #include <list>
23 #include <iostream>
24 #include <cmath>
25 #include <limits>
27 // Internal Includes
28 #include "project_namespace/Utility.hpp"
29 #include "project_namespace/Constants.hpp"
31 #include "project_namespace/Outer.hpp"
32 #include "project_namespace/Inner.hpp"
33 #include "project_namespace/Geometric.hpp"
35 #include "project_namespace/OuterExplicit.hpp"
36 #include "project_namespace/InnerExplicit.hpp"
37 #include "project_namespace/GeometricExplicit.hpp"
39 /*!
40  * @namespace project_namespace
41  */
42 namespace project_namespace{
45     /// \class Kvec
46     /// \brief class defining a single grade component of a multivector.
47     template<class T>
48     struct Kvec{
50         Eigen::Matrix<T, Eigen::Dynamic, 1> vec;  /*!< dynamic vector of Eigen Library */
52         unsigned int grade; /*!< grade k of the k-vector */
54         /// \brief operator == to test the equality between 2 k-vectors.
55         bool operator == (const Kvec<T>& other) const {
56             return vec == other.vec;
57         }
58     };
62     /// \class Mvec
63     /// \brief class defining multivectors.
64     template<typename T = double>
65     class Mvec {
67     protected:
68         std::list<Kvec<T>> mvData;  /*!< set of k-vectors, mapped by grade */
69         unsigned int gradeBitmap;   /*!< ith bit to 1 if grade i is contained in the multivector */
71     public:
73         /// \brief Default constructor, generate an empty multivector equivalent to the scalar 0.
74         Mvec();
76         /// \brief Copy constructor
77         /// \param mv - the multivector to be copied
78         Mvec(const Mvec& mv);
80         /// \brief Copy constructor of Mvec
81         /// \param mv - the multivector which has to be copied
82         Mvec(Mvec&& mv); // move constructor
85         template <typename U>
86         friend class Mvec;
88         /// \brief Copy constructor of Mvec from different types
89         /// \param mv - the multivector with another template type
90         template<typename U>
91         Mvec<T>(const Mvec<U> &mv);
93         /// \brief Constructor of Mvec from a scalar
94         /// \param val - scalar value
95         template<typename S>
96         Mvec(const S val);
98         /// Destructor
99         ~Mvec();
101         /// \brief Overload the assignment operator. No need any copy when the argument is an R-value, just need to move.
102         /// \param mv - Mvec used as a Rvalue
103         /// \return assign other to this
104         Mvec& operator=(Mvec&& mv);
106         /// \brief Overload the assignment operator
107         /// \param mv - Mvec
108         /// \return assign mv to this object
109         Mvec& operator=(const Mvec& mv);
111         /// \brief defines the addition between two Mvec
112         /// \param mv2 - second operand of type Mvec
113         /// \return this + mv2
114         Mvec operator+(const Mvec &mv2) const;
116         /// \brief defines the addition between a Mvec and a scalar
117         /// \param value - second operand (scalar)
118         /// \return this + scalar
119         template<typename S>
120         Mvec operator+(const S &value) const;
122         /// \brief defines the addition between a scalar and a Mvec
123         /// \param value a scalar
124         /// \param mv - second operand of type Mvec
125         /// \return scalar + mv
126         template<typename U, typename S>
127         friend Mvec<U> operator+(const S &value, const Mvec<U> &mv);
129         /// \brief Overload the += operator, corresponds to this += mv
130         /// \param mv - Mvec to be added to this object
131         /// \return this += mv
132         Mvec& operator+=(const Mvec& mv);
134         /// \brief defines the opposite of a multivector
135         /// \param mv: operand of type Mvec
136         /// \return -mv
137         template<typename U>
138         friend Mvec<U> operator-(const Mvec<U> &mv); // unary operator -mv1
140         /// \brief defines the difference between two Mvec
141         /// \param mv2 - second operand of type Mvec
142         /// \return this - mv2
143         Mvec<T> operator-(const Mvec<T> &mv2) const;
145         /// \brief defines the difference between a Mvec and a scalar
146         /// \param value - second operand (scalar)
147         /// \return this - scalar
148         template<typename S>
149         Mvec operator-(const S &value) const;
151         /// \brief defines the difference between a scalar and a Mvec
152         /// \param value a scalar
153         /// \param mv - second operand of type Mvec
154         /// \return scalar - mvX
155         template<typename U, typename S>
156         friend Mvec<U> operator-(const S &value, const Mvec<U> &mv);
158         /// \brief Overload the -= operator, corresponds to this -= mv
159         /// \param mv - Mvec to be added to this object
160         /// \return this -= mv
161         Mvec& operator-=(const Mvec& mv);
163         /// \brief defines the outer product between two multivectors
164         /// \param mv2 - a multivector
165         /// \return this^mv2
166         Mvec operator^(const Mvec &mv2) const;
168         /// \brief defines the outer product a multivector and a scalar
169         /// \param value - a scalar
170         /// \return this^value
171         template<typename S>
172         Mvec operator^(const S &value) const;
174         /// \brief defines the outer product between a scalar and a multivector
175         /// \param value - a scalar
176         /// \param mv - a multivector
177         /// \return value ^ mv
178         template<typename U, typename S>
179         friend Mvec<U> operator^(const S &value, const Mvec<U> &mv);
181         /// \brief Overload the outer product with operator =, corresponds to this ^= mv
182         /// \param mv - Mvec to be wedged to this object
183         /// \return this ^= mv
184         Mvec& operator^=(const Mvec& mv);
186         /// \brief defines the inner product between two multivectors
187         /// \param mv2 - a multivector
188         /// \return this.mv2
189         Mvec operator|(const Mvec &mv2) const;
191         /// \brief defines the inner product between a multivector and a scalar
192         /// \param value - a scalar
193         /// \return this.value = 0 (empty multivector)
194         template<typename S>
195         Mvec operator|(const S &value) const;
197         /// \brief defines the inner product between a scalar and a multivector
198         /// \param value - a scalar
199         /// \param mv - a multivector
200         /// \return mv.value = 0 (empty multivector)
201         template<typename U, typename S>
202         friend Mvec<U> operator|(const S &value, const Mvec<U> &mv);
204         /// \brief Overload the inner product with operator =, corresponds to this |= mv
205         /// \param mv - Mvec to be computed in the inner product with this object
206         /// \return this |= mv
207         Mvec& operator|=(const Mvec& mv);
209         /// \brief defines the right contraction between two multivectors
210         /// \param mv2 - a multivector
211         /// \return the right contraction : $this \\lfloor mv2$
212         Mvec operator>(const Mvec &mv2) const;
214         /// \brief defines the right contraction between a multivector and a scalar
215         /// \param value - a scalar
216         /// \return $this \\lfloor value$
217         template<typename S>
218         Mvec operator>(const S &value) const;
220         /// \brief defines the right contraction between a scalar and a multivector
221         /// \param value - a scalar
222         /// \param mv - a multivector
223         /// \return $value \\lfloor this = 0$ (empty multivector)
224         template<typename U, typename S>
225         friend Mvec<U> operator>(const S &value, const Mvec<U> &mv);
227         /// \brief defines the left contraction between two multivectors
228         /// \param mv2 - a multivector
229         /// \return the left contraction $this \\rfloor mv2$
230         Mvec operator<(const Mvec &mv2) const;
232         /// \brief defines the left contraction between a multivector and a scalar
233         /// \param value - a scalar
234         /// \return $value \\rfloor this$ = 0 (empty multivector)
235         template<typename S>
236         Mvec operator<(const S &value) const;
238         /// \brief defines the left contraction between a scalar and a multivector
239         /// \param value - a scalar
240         /// \param mv - a multivector
241         /// \return $this \\rfloor value$
242         template<typename U, typename S>
243         friend Mvec<U> operator<(const S &value, const Mvec<U> &mv);
245         /// \brief defines the geometric product between two multivectors
246         /// \param mv2 - a multivector
247         /// \return mv1*mv2
248         Mvec operator*(const Mvec &mv2) const;
250     project_singular_metric_comment_begin
251         /// \brief defines the outer product between two multivectors, where the second multivector is dualized during the product computation.
252         /// \param mv2 - a multivector that will be dualized during the wedge with the calling multivector.
253         /// \return a multivector.
254         Mvec<T> outerPrimalDual(const Mvec<T> &mv2) const;
256         /// \brief defines the outer product between two multivectors, where the first multivector is dualized during the product computation.
257         /// \param mv2 - a primal form of a multivector; the object multivector will be dualized during the wedge with the calling multivector.
258         /// \return a multivector.
259         Mvec<T> outerDualPrimal(const Mvec<T> &mv2) const;
261         /// \brief defines the outer product between two multivectors, where both multivectors are dualized during the product computation.
262         /// \param mv2 - a primal form of a multivector; the object multivector will be dualized during the wedge with the calling multivector.
263         /// \return a multivector.
264         Mvec<T> outerDualDual(const Mvec<T> &mv2) const;
266     project_singular_metric_comment_end
267         
269         /// \brief defines the scalar product between two multivectors (sum of the inner products between same grade pairs from the 2 multivectors)
270         /// \param mv2 - a multivector
271         /// \return a scalar
272         Mvec<T> scalarProduct(const Mvec<T> &mv2) const;
274         /// \brief defines the Hestenes product between two multivectors (Inner product - scalar product)
275         /// \param mv2 - a multivector
276         /// \return a multivector.
277         Mvec<T> hestenesProduct(const Mvec<T> &mv2) const;
279         /// \brief defines the geometric product between a multivector and a scalar
280         /// \param value - a scalar
281         /// \return mv2*value
282         template<typename S>
283         Mvec operator*(const S &value) const;
285         /// \brief defines the geometric product between a scalar and a multivector
286         /// \param value - a scalar
287         /// \param mv - a multivector
288         /// \return value*mv2
289         template<typename U, typename S>
290         friend Mvec<U> operator*(const S &value, const Mvec<U> &mv);
292         /// \brief Overload the geometric product with operator =, corresponds to this *= mv
293         /// \param mv - Mvec to be multiplied to this object
294         /// \return this *= mv
295         Mvec& operator*=(const Mvec& mv);
297         /// \brief defines the geometric product with a multivector and the inverse of a second multivector
298         /// \param mv2 - a multivector
299         /// \return this / mv2
300         Mvec operator/(const Mvec &mv2) const;
302         /// \brief defines the scalar inverse between a multivector and a scalar
303         /// \param value - a scalar
304         /// \return this / value
305         template<typename S>
306         Mvec operator/(const S &value) const;
308         /// \brief defines the inverse product between a scalar and a multivector
309         /// \param value - a scalar
310         /// \param mv - a multivector
311         /// \return value / mv
312         template<typename U, typename S>
313         friend Mvec<U> operator/(const S &value, const Mvec<U> &mv);
315         /// \brief Overload the inverse with operator =, corresponds to this /= mv
316         /// \param mv - Mvec to be inversed to this object
317         /// \return this /= mv
318         Mvec& operator/=(const Mvec& mv);
320         /// \brief the reverse of a multivector, i.e. if mv = a1^a2^...^an, then reverse(mv) = an^...^a2^a1
321         /// \param mv - a multivector
322         /// \return reverse of mv
323         template<typename U>
324         friend Mvec<U> operator~(const Mvec<U> &mv);
326     project_singular_metric_comment_begin
327         /// \brief the dual of a k-vector is defined as $A_k^* = A_k \\lcont I_n^{-1}$, for a multivector, we just dualize all its components.
328         /// \param mv - a multivector
329         /// \return the dual of mv
330         template<typename U>
331         friend Mvec<U> operator!(const Mvec<U> &mv);
332     project_singular_metric_comment_end
333         /// \brief boolean operator that tests the equality between two Mvec
334         /// \param mv2 - second operand of type Mvec
335         /// \return whether two Mvec have the same coefficients
336         inline bool operator==(const Mvec& mv2){
337             if(gradeBitmap != mv2.gradeBitmap)
338                 return false;
340             return mvData == mv2.mvData;   //// listcaca : ca marche que si les listes sont ordonnees
341         }
343         /// \brief compute the inverse of a multivector
344         /// \return - the inverse of the current multivector
345         Mvec<T> inv() const;
348             /// \brief operator to test whether two Mvec have not the same coefficients
349         /// \param mv2 - second operand of type Mvec
350         /// \return boolean that specify the non-equality between two Mvec
351         inline bool operator!=(const Mvec& mv2){
352             return !(mvData == mv2.mvData);
353         }
355         /// \brief Display all the non-null basis blades of this objects
356         /// \param stream - destination stream
357         /// \param mvec - the multivector to be outputed
358         /// \return a stream that contains the list of the non-zero element of the multivector
359         template<typename U>
360         friend std::ostream& operator<< (std::ostream& stream, const Mvec<U> &mvec);
362         /// \brief overload the casting operator, using this function, it is now possible to compute : float a = float(mv);
363         /// \return the scalar part of the multivector
364         operator T () {
365             // if no scalar part, return
367             if( (gradeBitmap & 1) == 0 )
368                 return 0;
370             // assuming now that the scalar part exists, return it
371             return findGrade(0)->vec[0];
372         }
374         /// \brief Overload the [] operator to assign a basis blade to a multivector. As an example, float a = mv[E12] = 42.
375         /// \param idx - the basis vector index related to the query.
376         /// \return the coefficient of the multivector corresponding to the "idx" component.
377         /// \todo handle the cases when the number of indices is higher than the dimension, and when the index is too high
378         T& operator[](const int idx ){
380             const unsigned int grade = xorIndexToGrade[idx];
381             const unsigned int idxHomogeneous = xorIndexToHomogeneousIndex[idx];
383             auto it = mvData.begin();
384             while(it != mvData.end()){
386                 // if grade not reach yet, continue
387                 if(it->grade < grade)
388                     ++it;
389                 else{
391                     // if grade found
392                     if(it->grade == grade)
393                         return it->vec[idxHomogeneous];
395                     // if grade exceed, create it and inster it before the current element
396                     Kvec<T> kvec = {Eigen::Matrix<T, Eigen::Dynamic, 1>::Zero(binomialArray[grade]),grade};
397                     auto it2 = mvData.insert(it,kvec);
398                     gradeBitmap |= 1 << (grade);
399                     return it2->vec[idxHomogeneous];
400                 }
401             }
403             // if the searched element should be added at the end, add it
404             Kvec<T> kvec = {Eigen::Matrix<T, Eigen::Dynamic, 1>::Zero(binomialArray[grade]),grade};
406             auto it2 = mvData.insert(mvData.end(),kvec);
407             gradeBitmap |= 1 << (grade);
408             return it2->vec[idxHomogeneous];
409         }
411         /// \brief Overload the [] operator to copy a basis blade of this multivector. As an example, float a = mv[E12].
412         /// \param idx - the basis vector index related to the query.
413         /// \return the coefficient of the multivector corresponding to the "idx" component.
414         const T& operator[](const int idx) const{
416             const unsigned int grade = xorIndexToGrade[idx];
417             const unsigned int idxHomogeneous = xorIndexToHomogeneousIndex[idx];
419             auto it = mvData.begin();
420             while(it != mvData.end()){
422                 // if grade not reach yet, continue
423                 if(it->grade < grade)
424                     ++it;
425                 else{
427                     // if grade found
428                     if(it->grade == grade)
429                         return it->vec[idxHomogeneous];
431                     // if grade exceed, return a reference on ... humm ...
432                     T val = T(0);
433                     return val;
434                 }
435             }
437             // if the searched element should be added at the end, humm ...
438             T val = T(0);   /// \todo find an elegant way to correct this
439             return val;
440         }
442 /*
443         /// \cond DEV
444         // Variadic operators
445         /// \brief Overload the () operator to assign a basis blade to a multivector. As an example, consider we have a Mvec a and we want to set its e23 component to 4.7, using operator(), it will consists in writing a[2,3]=4.8.
446         /// \tparam Arguments - denotes a variadic list of index
447         /// \param listIndices - the list of indices
448         /// \return the right index in the homogeneous vectorXd
449         /// \todo future work + handle the cases when the number of indices is higher than the dimension, and when the index is too high
450         template<typename... List>
451         T& operator()(List ...listIndices){
452             // operation on these arguments
453             // First determine the grade, simply use the variadic function sizeof...()
454             // This function is correct in terms of indexing
455             const int gd = sizeof...(listIndices); //
457             // Second, from these parameters, compute the index in the corresponding VectorXd
458             const int idx = (Binomial<algebraDimension,gd>::value-1) - computeIdxFromList(algebraDimension,gd,listIndices...);//Binomial<algebraDimension,binomCoef>::value;//binomCoef - sumB(first,listIndices...);//computeOrderingIndex<>(listIndices...);//(Binomial<algebraDimension,gd>::value-1);//-computeOrderingIndex<int,List>(4,2,listIndices...);//computeOrderingIndex<algebraDimension,gd>(listIndices...);//idxVariadic<algebraDimension,gd,first,listIndices...>::value; //computeOrderingIndex<algebraDimension,gd,listIndices...>();//idxVariadic<algebraDimension,gd,listIndices...>::value; // (Binomial<algebraDimension,gd>::value-1)- idxVariadic<algebraDimension,gd,listIndices...>::value
460             if(mvData.count(gd)>0)
461                 return (mvData.at(gd))(idx);
463             mvData[gd] = Eigen::Matrix<T, Eigen::Dynamic, 1>::Zero(Binomial<algebraDimension,gd>::value);
465             return mvData[gd](idx);
466         }
467         /// \endcond
468 */
470         /// \cond DEV
471         /// \brief returns a multivector with one component to 1
472         /// \param grade : grade of the component to enable
473         /// \param index : index of the parameter in the k-vector (k = grade)
474         inline Mvec componentToOne(const unsigned int grade, const int index){
475             Kvec<T> kvec = {.vec=Eigen::Matrix<T, Eigen::Dynamic, 1>::Zero(binomialArray[grade]),.grade=grade};
476             kvec.vec[index] = T(1);
477             Mvec mv1;
478             mv1.mvData.push_back(kvec);
479             mv1.gradeBitmap = 1 << (grade);
481             return mv1;
482         }
483         /// \endcond // do not comment this functions
485         /// \cond DEV
486         /// \brief create a VectorXd if it has not yet been created
487         /// \param grade - grade of the considered kvector
488         /// \return nothing
489         inline typename std::list<Kvec<T>>::iterator createVectorXdIfDoesNotExist(const unsigned int grade){
491             auto it = mvData.begin();
492             while(it != mvData.end()){
494                 // if grade not reach yet, continue
495                 if(it->grade < grade)
496                     ++it;
497                 else{
499                     // if grade found
500                     if(it->grade == grade)
501                         return it;
503                     // if grade exceed, create it and inster it before the current element
504                     Kvec<T> kvec = {.vec=Eigen::Matrix<T, Eigen::Dynamic, 1>::Zero(binomialArray[grade]),.grade=grade};
505                     auto it2 = mvData.insert(it,kvec);
506                     gradeBitmap |= 1 << (grade);
507                     return it2;
508                 }
509             }
511             // if the searched element should be added at the end, add it
512             Kvec<T> kvec = {.vec=Eigen::Matrix<T, Eigen::Dynamic, 1>::Zero(binomialArray[grade]),.grade=grade};
513             auto it2 = mvData.insert(mvData.end(),kvec);
514             gradeBitmap |= 1 << (grade);
515             return it2;
516         }
519         /// \cond DEV
520         /// \brief modify the element of the multivector whose grade is "grade"
521         /// \param grade : the grade to access
522         /// \param indexVectorXd : index of the k-vector (with k = "grade")
523         /// \return the element of the Mv whose grade is grade and index in the VectorXd is indexVectorXd
524         inline T& at(const int grade, const int indexVectorXd){
525             auto it = createVectorXdIfDoesNotExist(grade);
526             return it->vec.coeffRef(indexVectorXd);  
527         }
528         /// \endcond // do not comment this functions
530         /// \cond DEV
531         /// \brief return the element of the multivector whose grade is "grade"
532         /// \param grade : the grade to access
533         /// \param indexVectorXd : index of the k-vector (with k = "grade")
534         /// \return the element of the Mv whose grade is grade and index in the VectorXd is indexVectorXd
535         inline T at(const int grade, const int indexVectorXd) const{
537             if((gradeBitmap & (1<<(grade))) == 0 )
538                 return T(0);
540             auto it = findGrade(grade);
541             return it->vec.coeff(indexVectorXd);
542         }
543         /// \endcond // do not comment this functions
545         /// \brief the L2-norm of the mv is sqrt( abs( mv.mv ) )
546         /// \return the L2-norm of the multivector (as a double)
547         T inline norm() const {
548             return sqrt( fabs( (*this).scalarProduct( this->reverse() ) ));
549         };
551         /// \brief the L2-norm over 2 of the mv is mv.mv
552         /// \return the L2-norm of the multivector (as a double)
553         T inline quadraticNorm() const {
554             return ( this->reverse() | (*this) );
555         };
557     project_singular_metric_comment_begin
558         /// \brief compute the dual of a multivector (i.e mv* = reverse(mv) * Iinv)
559         /// \return - the dual of the multivector
560         Mvec<T> dual() const;
561     project_singular_metric_comment_end
562         /// \brief compute the reverse of a multivector
563         /// \return - the reverse of the multivector
564         Mvec<T> reverse() const;
566         /// \brief search in the multivector for a Kvec of grade "grade"
567         /// \return return a const iterator on the Kvec if exist, else return mvData.end()
568         inline typename std::list<Kvec<T>>::const_iterator findGrade(const unsigned int & gradeToFind) const {
569             for(auto it = mvData.begin();it != mvData.end();++it){
570                 if(it->grade==gradeToFind){
571                     return it;
572                 }
573             }
574             return mvData.end();
575         }
577         /// \brief search in the multivector for a Kvec of grade "grade"
578         /// \return return an iterator on the Kvec if exist, else return mvData.end()
579         inline typename std::list<Kvec<T>>::iterator findGrade(const unsigned int & gradeToFind) {
580             for(auto it = mvData.begin();it != mvData.end();++it){
581                 if(it->grade==gradeToFind){
582                     return it;
583                 }
584             }
585             return mvData.end();
586         }
588         /// \brief return the (highest) grade of the multivector
589         /// \return the highest grade of the multivector
590         inline int grade() const {
591             return (mvData.begin() == mvData.end()) ? 0 : mvData.rbegin()->grade;
592         }
594         /// \brief return the all non-zero grades of the multivector (several grades for non-homogeneous multivectors)
595         /// \return a list of the grades encountered in the multivector
596         std::vector<unsigned int> grades() const;
598         /// \brief returns a multivector that contains all the components of this multivector whose grade is i
599         /// \return an empty Mvec if the requested element is not part of the multivector, or the multivector that contains only this element if present in the current multivector.
600         Mvec grade(const int i) const;
602         /// \brief tell whether the multivector has grade component
603         /// \param grade - grade of the considered kvector
604         /// \return true if multivector has grade component, false else
605         inline bool isGrade(const unsigned int grade) const{
606             return ( (gradeBitmap & (1<<(grade+1))) != 0 );
607         }
609         /// \brief partially or completely erase the content of a multivector
610         /// \param grade : if < 0, erase all the multivector, else just erase the part of grade "grade".
611         void clear(const int grade = -1);
613         /// \brief check is a mutivector is empty, i.e. corresponds to 0.
614         /// \return True if the multivector is empty, else False.
615         inline bool isEmpty() const {
616             return mvData.empty();
617         }
619         /// \brief A multivector is homogeneous if all its components have the same grade.
620         /// \return True if the multivector is homogeneous, else False.
621         inline bool isHomogeneous() const {
622             return mvData.size() < 2; // only one element, or zero element on the list
623         }
625         /// \brief inplace simplify the multivector such that all the values with a magnitude lower than a epsilon in the Mv are set to 0.
626         /// \param epsilon - threshold, with default value the epsilon of the float/double/long double type from numeric_limits
627         void roundZero(const T epsilon = std::numeric_limits<T>::epsilon());
629         /// \brief Specify if two multivectors have the same grade
630         /// \param mv1 - first multivector
631         /// \param mv2 - second multivector
632         /// \return true if the two multivectors have the same grade, else return false
633         inline bool sameGrade(const Mvec<T>& mv) const{
634             return grade() == mv.grade();
635         }
637         /// \brief Display the multivector data (per grade value)
638         void display() const;
640         /// \cond DEV
641         /// \brief a function to output all the element of a k-vector in the multivector indexing order.
642         /// \tparam T - type of the multivector
643         /// \param stream - stream that will contain the result
644         /// \param mvec - multivector to be processed
645         /// \param gradeMV - the considered grade
646         /// \param moreThanOne - true if it the first element to display (should we put a '+' before)
647         template<typename U>
648         friend void traverseKVector(std::ostream &stream, const Eigen::Matrix<U, Eigen::Dynamic, 1>, unsigned int gradeMV, bool& moreThanOne);
649         /// \endcond // do not comment this functions
651 /*
652         /// \cond DEV
653         /// \brief functions that enables to  to be able to extract a component of a multivector as a multivector. As an example, mv1.e(1,3) will create a multivector whose 1,3 component will be the component of mv1.
654         /// \tparam Arguments - denotes a variadic list of index
655         /// \param listIndices - the list of indices
656         /// \return a homogeneous multivector filled with zero except the listIndices index
657         /// \todo future work + handle the cases when the number of indices is higher than the dimension, and when the index is too high
658         template<typename... List>
659         Mvec e(List ...listIndices){
660             // operation on these arguments
661             // First determine the grade, simply use the variadic function sizeof...()
662             const int gd = sizeof...(listIndices); //
664             // Second, from these parameters, compute the index in the corresponding VectorXd
665             const int idx = (Binomial<algebraDimension,gd>::value-1) - computeIdxFromList(algebraDimension,gd,listIndices...);
667             Mvec res;
668             res.mvData[gd] = Eigen::Matrix<T, Eigen::Dynamic, 1>::Zero(Binomial<algebraDimension,gd>::value);
669             res.mvData[gd](idx) = mvData[gd](idx);
671             return res;
672         }
673         /// \endcond // do not comment this functions
674 */
677         /// \brief functions that enables to extract one component of this multivector and initialize another mv with this component
678         /// \tparam Arguments - denotes a variadic list of index
679         /// \param grade : the considered grade
680         /// \param sizeOfKVector: C(dimension, grade)
681         /// \param indexInKvector: the considerd index
682         /// \return a new mv with the right component at the right place
683         Mvec extractOneComponent(const int grade, const int sizeOfKVector, const int indexInKvector) const {
684             Mvec mv;
685             auto itThisMV = this->findGrade(grade);
686             if(itThisMV==this->mvData.end()){
687                 return mv;
688             }
689             mv = mv.componentToOne(grade,indexInKvector);
690             mv.mvData.begin()->vec.coeffRef(indexInKvector) = itThisMV->vec.coeff(indexInKvector);
691             return mv;
692         }
695         /// \brief set of functions that return multivectors containing only the value '1' for the specified component.
696 project_multivector_one_component
698     };  // end of the class definition
701     /* ------------------------------------------------------------------------------------------------ */
704     template<typename T>
705     Mvec<T>::Mvec():gradeBitmap(0)
706     {}
709     template<typename T>
710     Mvec<T>::Mvec(const Mvec& mv) : mvData(mv.mvData), gradeBitmap(mv.gradeBitmap)
711     {
712         //std::cout << "copy constructor " << std::endl;
713     }
716     template<typename T>
717     Mvec<T>::Mvec(Mvec<T>&& multivector) : mvData(std::move(multivector.mvData)), gradeBitmap(multivector.gradeBitmap)
718     {
719         // the move constructor
720         //std::cout << "move constructor" << std::endl;
721     }
724     template<typename T>
725     template<typename U>
726     Mvec<T>::Mvec(const Mvec<U> &mv) : gradeBitmap(mv.gradeBitmap)
727     {
728         for(auto it=mv.mvData.begin(); it != mv.mvData.end(); ++it){
729             Kvec<T> kvec;
730             kvec.grade = it->grade;
731             kvec.vec = Eigen::Matrix<T, Eigen::Dynamic, 1>::Zero(binomialArray[it->grade]);
732             for(unsigned int i=0; i<it->vec.size(); ++i)
733                 kvec.vec.coeffRef(i) = T(it->vec.coeff(i));
734             mvData.push_back(kvec);
735         }
736     }
739     template<typename T>
740     template<typename U>
741     Mvec<T>::Mvec(const U val) {
742         if(val != U(0)) {
743             gradeBitmap = 1;
744             mvData.push_back({.vec=Eigen::Matrix<T, Eigen::Dynamic, 1>(1),.grade=0});
745             mvData.begin()->vec.coeffRef(0) = val;
746         }
747     }
750     template<typename T>
751     Mvec<T>::~Mvec()
752     {}
755     template<typename T>
756     Mvec<T>& Mvec<T>::operator=(const Mvec& mv){
757         if(&mv == this) return *this;
758         gradeBitmap = mv.gradeBitmap;
759         mvData = mv.mvData;
760         return *this;
761     }
764     template<typename T>
765     Mvec<T>& Mvec<T>::operator=(Mvec&& mv){
766         mvData = std::move(mv.mvData);
767         gradeBitmap = mv.gradeBitmap;
768         return *this;
769     }
772     template<typename T>
773     Mvec<T> Mvec<T>::operator+(const Mvec<T> &mv2) const {
774         Mvec<T> mv3(*this);
775         for(auto & itMv : mv2.mvData) {
776             auto it = mv3.createVectorXdIfDoesNotExist(itMv.grade);
777             it->vec += itMv.vec;
778         }
779         return mv3;
780     }
783     template<typename U, typename S>
784     Mvec<U> operator+(const S &value, const Mvec<U> &mv){
785         return mv + value;
786     }
789     template<typename T>
790     template<typename S>
791     Mvec<T> Mvec<T>::operator+(const S &value) const {
792         Mvec<T> mv(*this);
793         if(value != T(0)) {
794             auto it = mv.createVectorXdIfDoesNotExist(0);
795             it->vec.coeffRef(0) += value;
796         }
797         return mv;
798     }
801     template<typename T>
802     Mvec<T>& Mvec<T>::operator+=(const Mvec& mv){
803         *this = *this + mv;
804         return *this;
805     }
808     template<typename T>
809     Mvec<T> operator-(const Mvec<T> &mv) { // unary -
810         Mvec<T> mv2(mv);
811         for(auto & itMv : mv2.mvData)
812             itMv.vec = -itMv.vec;
813         return mv2;
814     }
817     template<typename T>
818     Mvec<T> Mvec<T>::operator-(const Mvec<T> &mv2) const {
819         Mvec<T> mv3(*this);
820         for(auto & itMv : mv2.mvData) {
821             auto it = mv3.createVectorXdIfDoesNotExist(itMv.grade);
822             it->vec -= itMv.vec;
823         }
824         return mv3;
825     }
828     template<typename T>
829     Mvec<T>& Mvec<T>::operator-=(const Mvec& mv){
830         *this = *this - mv;
831         return *this;
832     }
835     template<typename U, typename S>
836     Mvec<U> operator-(const S &value, const Mvec<U> &mv){
837         return (-mv) + value;
838     }
841     template<typename T>
842     template<typename S>
843     Mvec<T> Mvec<T>::operator-(const S &value) const {
844         Mvec<T> mv(*this);
845         if(value != T(0)) {
846             auto it = mv.createVectorXdIfDoesNotExist(0);
847             it->vec.coeffRef(0) -= value;
848         }
849         return mv;
850     }
853     template<typename T>
854     Mvec<T> Mvec<T>::operator^(const Mvec<T> &mv2) const {
855 #if 0 // only recursive version
856         // Loop over non-empty grade of mv1 and mv2
857         // This version (with recursive call) is only faster than the standard recursive call if
858         // the multivector mv1 and mv2 are homogeneous or near from homogeneous.
859         Mvec<T> mv3;
860         for(const auto & itMv1 : this->mvData)    // all per-grade component of mv1
861             for(const auto & itMv2 : mv2.mvData)  // all per-grade component of mv2
862             {
863                 unsigned int grade_mv3 = itMv1.grade + itMv2.grade;
864                 if(grade_mv3 <= algebraDimension) {
865                     auto itMv3 = mv3.createVectorXdIfDoesNotExist(grade_mv3);
866                     outerProductHomogeneous(itMv1.vec, itMv2.vec, itMv3->vec,
867                                             itMv1.grade, itMv2.grade, grade_mv3);
868                 }
869             }
870         return mv3;
871 #else // use the adaptative pointer function array
872         // Loop over non-empty grade of mv1 and mv2
873         // call the right explicit unrolled outer function using the functions pointer called outerFunctionsContainer
874         Mvec<T> mv3;
875         for(const auto & itMv1 : this->mvData)
876             for(const auto & itMv2 : mv2.mvData){
877                 if(itMv1.grade + itMv2.grade <= (int) algebraDimension ){
878                     auto itMv3 = mv3.createVectorXdIfDoesNotExist(itMv1.grade + itMv2.grade);
879                     outerFunctionsContainer<T>[itMv1.grade][itMv2.grade](itMv1.vec, itMv2.vec, itMv3->vec);
880                 }
881             }
882         return mv3;
883 #endif
884     }
886     template<typename U, typename S>
887     Mvec<U> operator^(const S &value, const Mvec<U> &mv) {
888         return  mv^value;
889     }
892     template<typename T>
893     template<typename S>
894     Mvec<T> Mvec<T>::operator^(const S &value) const {
895         Mvec<T> mv2(*this);
896         for(auto & itMv : mv2.mvData)
897             itMv.vec *= T(value);
898         return mv2;
899     }
902     template<typename T>
903     Mvec<T> &Mvec<T>::operator^=(const Mvec &mv) {
904         *this = *this ^ mv;
905         return *this;
906     }
909     template<typename T>
910     Mvec<T> Mvec<T>::operator|(const Mvec<T> &mv2) const{
911         // Loop over non-empty grade of mv1 and mv2
912         // call the right explicit unrolled inner function using the functions pointer called innerFunctionsContainer
913         Mvec<T> mv3;
914         for(const auto & itMv1 : this->mvData)
915             for(const auto & itMv2 : mv2.mvData){
917                 // perform the inner product
918                 int absGradeMv3 = std::abs((int)(itMv1.grade - itMv2.grade));
919                 auto itMv3 = mv3.createVectorXdIfDoesNotExist(absGradeMv3);
920                 innerFunctionsContainer<T>[itMv1.grade][itMv2.grade](itMv1.vec, itMv2.vec, itMv3->vec);
922                 // check if the result is non-zero
923                 if(!((itMv3->vec.array() != 0.0).any())){
924                     mv3.mvData.erase(itMv3);
925                     mv3.gradeBitmap &= ~(1<<absGradeMv3);
926                 }
927             }
928         return mv3;
929     }
932     template<typename U, typename S>
933     Mvec<U> operator|(const S &value, const Mvec<U> &mv) {
934         return mv | value;
935     }
938     template<typename T>
939     template<typename S>
940     Mvec<T> Mvec<T>::operator|(const S &value) const {
941         return (*this) > value ; // inner between a mv and a scalar gives 0 (empty multivector)
942     }
945     template<typename T>
946     Mvec<T> &Mvec<T>::operator|=(const Mvec &mv) {
947         *this = *this | mv;
948         return *this;
949     }
952     template<typename T>
953     Mvec<T> Mvec<T>::operator>(const Mvec<T> &mv2) const{
954         // Loop over non-empty grade of mv1 and mv2
955         // call the right explicit unrolled inner function using the functions pointer called innerFunctionsContainer
956         Mvec<T> mv3;
957         for(const auto & itMv1 : this->mvData)
958             for(const auto & itMv2 : mv2.mvData){
960                 // right contraction constraint: gradeMv1 >= gradeMv2
961                 if(itMv1.grade < itMv2.grade)
962                     continue;
964                 // perform the inner product
965                 int absGradeMv3 = std::abs((int)(itMv1.grade - itMv2.grade));
966                 auto itMv3 = mv3.createVectorXdIfDoesNotExist(absGradeMv3);
967                 innerFunctionsContainer<T>[itMv1.grade][itMv2.grade](itMv1.vec, itMv2.vec, itMv3->vec);
969                 // check if the result is non-zero
970                 if(!((itMv3->vec.array() != 0.0).any())){
971                     mv3.mvData.erase(itMv3);
972                     mv3.gradeBitmap &= ~(1<<absGradeMv3);
973                 }
974             }
976         return mv3;
977     }
980     template<typename U, typename S>
981     Mvec<U> operator>(const S &value, const Mvec<U> &mv) {
982         if( (mv.gradeBitmap & 1) == 0) return Mvec<U>();
983         else return Mvec<U>( U(mv.mvData.begin()->vec.coeff(0) * value ) );
984     }
987     template<typename T>
988     template<typename S>
989     Mvec<T> Mvec<T>::operator>(const S &value) const {
990         // return mv x value
991         Mvec<T> mv(*this);
992         for(auto & itMv : mv.mvData)
993             itMv.vec *= T(value);
994         return mv;
995     }
998     template<typename T>
999     Mvec<T> Mvec<T>::operator<(const Mvec<T> &mv2) const{
1001         // Loop over non-empty grade of mv1 and mv2
1002         // call the right explicit unrolled inner function using the functions pointer called innerFunctionsContainer
1003         Mvec<T> mv3;
1004         for(const auto & itMv1 : this->mvData)
1005             for(const auto & itMv2 : mv2.mvData){
1007                 // left contraction constraint: gradeMv1 <= gradeMv2
1008                 if(itMv1.grade > itMv2.grade)
1009                     continue;
1011                 // perform the inner product
1012                 int absGradeMv3 = std::abs((int)(itMv1.grade - itMv2.grade));
1013                 auto itMv3 = mv3.createVectorXdIfDoesNotExist(absGradeMv3);
1014                 innerFunctionsContainer<T>[itMv1.grade][itMv2.grade](itMv1.vec, itMv2.vec, itMv3->vec);
1016                 // check if the result is non-zero
1017                 if(!((itMv3->vec.array() != 0.0).any())){
1018                     mv3.mvData.erase(itMv3);
1019                     mv3.gradeBitmap &= ~(1<<absGradeMv3);
1020                 }
1021             }
1022         return mv3;
1023     }
1026     template<typename U, typename S>
1027     Mvec<U> operator<(const S &value, const Mvec<U> &mv) {
1028         // return mv x value
1029         Mvec<U> mv3(mv);
1030         for(auto & itMv : mv3.mvData)
1031             itMv.vec *= U(value);
1032         return mv3;
1033     }
1036     template<typename T>
1037     template<typename S>
1038     Mvec<T> Mvec<T>::operator<(const S &value) const {
1039         if( (gradeBitmap & 1) == 0) return Mvec<T>();
1040         else return Mvec<T>( T(mvData.begin()->vec.coeff(0) * value ) );
1041     }
1044     template<typename T>
1045     Mvec<T> Mvec<T>::scalarProduct(const Mvec<T> &mv2) const{
1046         // Loop over non-empty grade of mv1 and mv2
1047         // call the right explicit unrolled inner function using the functions pointer called innerFunctionsContainer
1048         Mvec<T> mv3;
1049         for(const auto & itMv1 : this->mvData)
1050             for(const auto & itMv2 : mv2.mvData){
1052                 if(itMv1.grade != itMv2.grade)
1053                     continue;
1055                 // perform the inner product
1056                 int absGradeMv3 = 0;
1057                 auto itMv3 = mv3.createVectorXdIfDoesNotExist(absGradeMv3);
1058                 innerFunctionsContainer<T>[itMv1.grade][itMv2.grade](itMv1.vec, itMv2.vec, itMv3->vec);
1060                 // check if the result is non-zero
1061                 if(!((itMv3->vec.array() != 0.0).any())){
1062                     mv3.mvData.erase(itMv3);
1063                     mv3.gradeBitmap &= ~(1<<absGradeMv3);
1064                 }
1065             }
1066         return mv3;
1067     }
1070 project_singular_metric_comment_begin
1071     template<typename T>
1072     Mvec<T> Mvec<T>::outerPrimalDual(const Mvec<T> &mv2) const{
1073         Mvec<T> mv3;
1075         for(const auto & itMv1 : this->mvData)    // all per-grade component of mv1
1076             for(const auto & itMv2 : mv2.mvData)  // all per-grade component of mv2
1077             {
1078                 unsigned int grade_mv3 = itMv1.grade + (algebraDimension-itMv2.grade);
1079                 if(grade_mv3 <= algebraDimension) {
1080                     // handle the scalar product as well as the left contraction
1081                     auto itMv3 = mv3.createVectorXdIfDoesNotExist(grade_mv3);
1082                     outerProductPrimalDual(itMv1.vec, itMv2.vec, itMv3->vec,
1083                                            itMv1.grade, itMv2.grade, (unsigned)(algebraDimension-grade_mv3));
1085                     // check if the result is non-zero
1086                     if(!((itMv3->vec.array() != 0.0).any())){
1087                         mv3.mvData.erase(itMv3);
1088                         mv3.gradeBitmap &= ~(1<<grade_mv3);
1089                     }
1090                 }
1091             }
1093         return mv3;
1094     }
1096     template<typename T>
1097     Mvec<T> Mvec<T>::outerDualPrimal(const Mvec<T> &mv2) const{
1098         Mvec<T> mv3;
1100         for(const auto & itMv1 : this->mvData)    // all per-grade component of mv1
1101             for(const auto & itMv2 : mv2.mvData)  // all per-grade component of mv2
1102             {
1103                 unsigned int grade_mv3 = itMv1.grade + (algebraDimension-itMv2.grade);
1104                 if(grade_mv3 <= algebraDimension) {
1105                     // handle the scalar product as well as the left contraction
1106                     auto itMv3 = mv3.createVectorXdIfDoesNotExist(grade_mv3);
1107                     outerProductDualPrimal(itMv1.vec, itMv2.vec, itMv3->vec,
1108                                            itMv1.grade, itMv2.grade, (unsigned)(algebraDimension-grade_mv3));
1110                     // check if the result is non-zero
1111                     if(!((itMv3->vec.array() != 0.0).any())){
1112                         mv3.mvData.erase(itMv3);
1113                         mv3.gradeBitmap &= ~(1<<grade_mv3);
1114                     }
1115                 }
1116             }
1119         return mv3;
1121     }
1123     template<typename T>
1124     Mvec<T> Mvec<T>::outerDualDual(const Mvec<T> &mv2) const{
1125         Mvec<T> mv3;
1127         for(const auto & itMv1 : this->mvData)    // all per-grade component of mv1
1128             for(const auto & itMv2 : mv2.mvData)  // all per-grade component of mv2
1129             {
1130                 unsigned int grade_mv3 = itMv1.grade + (algebraDimension-itMv2.grade);
1131                 if(grade_mv3 <= algebraDimension) {
1132                     // handle the scalar product as well as the left contraction
1133                     auto itMv3 = mv3.createVectorXdIfDoesNotExist(grade_mv3);
1134                     outerProductDualDual(itMv1.vec, itMv2.vec, itMv3->vec,
1135                                            itMv1.grade, itMv2.grade, (unsigned)(algebraDimension-grade_mv3));
1137                     // check if the result is non-zero
1138                     if(!((itMv3->vec.array() != 0.0).any())){
1139                         mv3.mvData.erase(itMv3);
1140                         mv3.gradeBitmap &= ~(1<<grade_mv3);
1141                     }
1142                 }
1143             }
1146         return mv3;
1148     }
1150 project_singular_metric_comment_end
1152     template<typename T>
1153     Mvec<T> Mvec<T>::hestenesProduct(const Mvec<T> &mv2) const{
1154         // Loop over non-empty grade of mv1 and mv2
1155         // call the right explicit unrolled inner function using the functions pointer called innerFunctionsContainer
1156         Mvec<T> mv3;
1157         for(const auto & itMv1 : this->mvData)
1158             for(const auto & itMv2 : mv2.mvData){
1160                 if(itMv1.grade*itMv2.grade == 0)
1161                     continue;
1163                 // perform the inner product
1164                 int absGradeMv3 = std::abs((int)(itMv1.grade - itMv2.grade));
1165                 auto itMv3 = mv3.createVectorXdIfDoesNotExist(absGradeMv3);
1166                 innerFunctionsContainer<T>[itMv1.grade][itMv2.grade](itMv1.vec, itMv2.vec, itMv3->vec);
1168                 // check if the result is non-zero
1169                 if(!((itMv3->vec.array() != 0.0).any())){
1170                     mv3.mvData.erase(itMv3);
1171                     mv3.gradeBitmap &= ~(1<<absGradeMv3);
1172                 }
1173             }
1174         return mv3;
1175     }
1178     template<typename T>
1179     Mvec<T> Mvec<T>::operator*(const Mvec<T> &mv2) const {
1180         // Loop over non-empty grade of mv1 and mv2
1181         // call the right explicit unrolled product function using the functions pointer
1182         Mvec<T> mv3;
1183         for(const auto & itMv1 : this->mvData)
1184             for(const auto & itMv2 : mv2.mvData){
1185                 project_select_recursive_geometric_product_template
1186                 // outer product block
1187                 unsigned int gradeOuter = itMv1.grade + itMv2.grade;
1188                 if(gradeOuter <=  algebraDimension ){
1189                     auto itMv3 = mv3.createVectorXdIfDoesNotExist(gradeOuter);
1190                     outerFunctionsContainer<T>[itMv1.grade][itMv2.grade](itMv1.vec, itMv2.vec, itMv3->vec);
1191                     if(!((itMv3->vec.array() != 0.0).any())){
1192                         mv3.mvData.erase(itMv3);
1193                         mv3.gradeBitmap &= ~(1<<gradeOuter);
1194                     }
1195                 }
1197                 // inner product block
1198                 unsigned int gradeInner = (unsigned int)std::abs((int)(itMv1.grade-itMv2.grade));
1199                 // when the grade of one of the kvectors is zero, the inner product is the same as the outer product
1200                 if(gradeInner != gradeOuter) {
1201                     auto itMv3 = mv3.createVectorXdIfDoesNotExist(gradeInner);
1202                     innerFunctionsContainer<T>[itMv1.grade][itMv2.grade](itMv1.vec, itMv2.vec, itMv3->vec);
1203                     // check if the result is non-zero
1204                     if(!((itMv3->vec.array() != 0.0).any())){
1205                         mv3.mvData.erase(itMv3);
1206                         mv3.gradeBitmap &= ~(1<<gradeInner);
1207                     }
1209                     // geometric product part
1210                     int gradeMax = std::min(((2*algebraDimension)-gradeOuter)+1,gradeOuter);
1211                     for (int gradeResult = gradeInner+2; gradeResult < gradeMax; gradeResult+=2) {
1212                         auto itMv3 = mv3.createVectorXdIfDoesNotExist(gradeResult);
1213                         geometricFunctionsContainer<T>[gradeResult][itMv1.grade][itMv2.grade](itMv1.vec, itMv2.vec, itMv3->vec);
1214                         // check if the result is non-zero
1215                         if(!((itMv3->vec.array() != 0.0).any())){
1216                             mv3.mvData.erase(itMv3);
1217                             mv3.gradeBitmap &= ~(1<<gradeResult);
1218                         }
1219                     }
1220                 }
1221             }
1222         return mv3;
1223     }
1226     template<typename U, typename S>
1227     Mvec<U> operator*(const S &value, const Mvec<U> &mv){
1228         return mv * value;
1229     }
1232     template<typename T>
1233     template<typename S>
1234     Mvec<T> Mvec<T>::operator*(const S &value) const{
1235         Mvec<T> mv(*this);
1236         for(auto & itMv : mv.mvData)
1237             itMv.vec *= T(value);
1238         return mv;
1239     }
1242     template<typename T>
1243     Mvec<T> &Mvec<T>::operator*=(const Mvec &mv) {
1244         *this = *this * mv;
1245         return *this;
1246     }
1249     template<typename T>
1250     Mvec<T> Mvec<T>::operator/(const Mvec<T> &mv2) const {
1251         return *this * mv2.inv();
1252     }
1255     template<typename U, typename S>
1256     Mvec<U> operator/(const S &value, const Mvec<U> &mv){
1257         return value * mv.inv();
1258     }
1261     template<typename T>
1262     template<typename S>
1263     Mvec<T> Mvec<T>::operator/(const S &value) const {
1264         Mvec<T> mv(*this);
1265         for(auto & itMv : mv.mvData)
1266             itMv.vec /= value;
1267         return mv;
1268     }
1271     template<typename T>
1272     Mvec<T> &Mvec<T>::operator/=(const Mvec &mv) {
1273         *this = *this / mv;
1274         return *this;
1275     }
1278     template<typename T>
1279     Mvec<T> operator~(const Mvec<T> &mv){
1280         return mv.reverse();
1281     }
1284     template<typename T>
1285     Mvec<T> Mvec<T>::inv() const {
1286         T n = this->quadraticNorm();
1287         if(n<std::numeric_limits<T>::epsilon() && n>-std::numeric_limits<T>::epsilon())
1288             return Mvec<T>(); // return 0, this is was gaviewer does.
1289         return this->reverse() / n;
1290     };
1293     project_singular_metric_comment_begin
1294     template<typename U>
1295     Mvec<U> operator!(const Mvec<U> &mv) {
1296         return mv.dual();
1297     }
1298 project_singular_metric_comment_end
1300     /// \brief defines the left contraction between two multivectors
1301     /// \param mv1 - a multivector
1302     /// \param mv2 - second operand, corresponds to a multivector
1303     /// \return mv1 left contraction mv2
1304     template<typename T>
1305     Mvec<T> leftContraction(const Mvec<T> &mv1, const Mvec<T> &mv2){
1306         return mv1 < mv2;
1307     }
1310     /// \brief defines the right contraction between two multivectors
1311     /// \param mv1 - a multivector
1312     /// \param mv2 - second operand, corresponds to a multivector
1313     /// \return mv1 right contraction mv2
1314     template<typename T>
1315     Mvec<T> rightContraction(const Mvec<T> &mv1, const Mvec<T> &mv2){
1316         return mv1 > mv2;
1317     }
1320     /// \brief returns a multivector that only contains the coefficient associated to the pseudoscalar.
1321     /// \return an empty Mvec if the requested element is not part of the multivector, or the multivector that contains only this element if present in the current multivector.
1322     template<typename T>
1323     Mvec<T> I() {
1324         Mvec<T> mvec;
1325         mvec[project_pseudo_scalar] = 1.0;
1326         return mvec;
1327     }
1329 project_singular_metric_comment_begin
1330     /// \brief return the inverse of the pseudo scalar.
1331     /// \return returns a multivector corresponding to the inverse of the pseudo scalar.
1332     template<typename T>
1333     Mvec<T> Iinv() {
1334         Mvec<T> mvec;
1335         mvec[project_pseudo_scalar] = pseudoScalarInverse; // we do not return just a scalar of type T because the grade should be dimension and not 0.
1336         return mvec;
1337     }
1338 project_singular_metric_comment_end
1340 project_singular_metric_comment_begin
1341     // compute the dual of a multivector (i.e mv* = mv.reverse() * Iinv)
1342     template<typename T>
1343     Mvec<T> Mvec<T>::dual() const {
1344         Mvec<T> mvResult;
1345         // for each k-vectors of the multivector
1346         for(auto itMv=mvData.rbegin(); itMv!=mvData.rend(); ++itMv){
1348             // create the dual k-vector
1349             Kvec<T> kvec ={itMv->vec,algebraDimension-(itMv->grade)};
1351             // some elements need to be permuted
1352             for(unsigned int i=0;i<binomialArray[itMv->grade];++i)
1353                 kvec.vec.coeffRef(dualPermutations[itMv->grade][i]) = itMv->vec.coeff(i);
1355             // the inner product may involve some constant multiplucation for the dual elements
1356             kvec.vec = kvec.vec.cwiseProduct(dualCoefficients[itMv->grade].template cast<T>());
1358             // add the k-vector to the resulting multivector
1359             mvResult.mvData.push_back(kvec);
1360             mvResult.gradeBitmap |= (1 << kvec.grade);
1361         }
1362         return mvResult;
1363     };
1364 project_singular_metric_comment_end
1366     // \brief compute the reverse of a multivector
1367     // \return - the reverse of the multivector
1368     template<typename T>
1369     Mvec<T> Mvec<T>::reverse() const {
1370         Mvec<T> mv(*this);
1371         for(auto & itMv : mv.mvData)
1372             if(signReversePerGrade[itMv.grade] == -1)
1373                 itMv.vec *= -1;
1374         return mv;
1375     }
1378     template<typename T>
1379     Mvec<T> Mvec<T>::grade(const int i) const{
1381         auto it = findGrade(i);
1382         Mvec<T> mv;
1384         // if not found, return the empty multivector
1385         if(it == mvData.end())
1386             return mv;
1388         // else return the grade 'i' data
1389         mv.mvData.push_back(*it);
1390         mv.gradeBitmap = 1 << (i);
1392         return mv;
1393     }
1396     template<typename T>
1397     std::vector<unsigned int> Mvec<T>::grades() const {
1398         std::vector<unsigned int> mvGrades;
1399         for(unsigned int i=0; i< algebraDimension+1; ++i){
1401             if( (gradeBitmap & (1<<i)))
1402                 mvGrades.push_back(i);
1403         }
1405         return mvGrades;
1406     }
1409     template<typename T>
1410     void Mvec<T>::roundZero(const T epsilon) {
1411         // loop over each k-vector of the multivector
1412         auto itMv=mvData.begin();
1413         while(itMv != mvData.end()){
1414             // loop over each element of the k-vector
1415             for(unsigned int i=0; i<(unsigned int)itMv->vec.size(); ++i)
1416                 if(fabs(itMv->vec.coeff(i)) <= epsilon)
1417                     itMv->vec.coeffRef(i) = 0.0;
1418             // if the k-vector is full of 0, remove it
1419             if(!((itMv->vec.array() != 0.0).any())){
1420                 gradeBitmap = gradeBitmap - (1 << itMv->grade);
1421                 mvData.erase(itMv++);
1422             }
1423             else ++itMv;
1424         }
1425     }
1428     template<typename T>
1429     void Mvec<T>::clear(const int grade) {
1431         // full erase
1432         if(grade < 0){
1433             mvData.clear();
1434             gradeBitmap = 0;
1435             return;
1436         }
1438         // partial erase
1439         auto iter = findGrade(grade);
1440         if(iter != mvData.end()) {
1441             gradeBitmap = gradeBitmap - (1 << iter->grade);
1442             iter = mvData.erase(iter);
1443         }
1444     }
1447     template<typename T>
1448     void Mvec<T>::display() const {
1449         if(mvData.size() == 0)
1450             std::cout << " null grade , null value " <<std::endl;
1452         for(auto itmv=mvData.begin(); itmv!=mvData.end(); ++itmv){
1453             std::cout << "  grade   : " << itmv->grade  << std::endl;
1454             std::cout << "  kvector : " << itmv->vec.transpose() << std::endl;
1455         }
1456         std::cout << std::endl;
1457     }
1460     /// \cond DEV
1461     template<typename U>
1462     void traverseKVector(std::ostream &stream, const Eigen::Matrix<U, Eigen::Dynamic, 1> kvector, unsigned int gradeMV, bool& moreThanOne ){
1464         // version with XOR indices
1465         // for the current grade, generate all the XOR indices
1466         std::vector<bool> booleanXorIndex(algebraDimension);
1467         std::fill(booleanXorIndex.begin(), booleanXorIndex.end(), false);
1468         // build the first xorIndex of grade 'grade' (with 'grade' values to true).
1469         std::fill(booleanXorIndex.begin(), booleanXorIndex.begin() + gradeMV, true);
1470         unsigned int positionInKVector = 0;
1472         do {
1473             // convert the vector of bool to the string containing all the basis vectors
1474             std::string basisBlades={};
1475             for(unsigned int i=0; i<algebraDimension; ++i) {
1476                 if(booleanXorIndex[i]) {
1477                     basisBlades += basisVectors[i];
1478                 }
1479             }
1481             if(kvector.coeff(positionInKVector)!= 0) {
1482                 if(!(moreThanOne)){
1483                     stream<< kvector.coeff(positionInKVector) << "*e"+ basisBlades;
1484                     moreThanOne = true;
1485                 }else{
1486                     if(kvector.coeff(positionInKVector)>0)
1487                         stream<< " + " << kvector.coeff(positionInKVector) << "*e" + basisBlades;
1488                     else
1489                         stream<< " - " << -kvector.coeff(positionInKVector) << "*e" + basisBlades;
1490                 }
1491             }
1493             positionInKVector++;
1495         } while(std::prev_permutation(booleanXorIndex.begin(), booleanXorIndex.end())); // compute next permutation of the true values of booleanXorIndex
1497     }
1498     /// \endcond
1501     template<typename U>
1502     std::ostream &operator<<(std::ostream &stream, const Mvec<U> &mvec) {
1503         if(mvec.mvData.size()==0){
1504             stream << " 0 ";
1505             return stream;
1506         }
1508         bool moreThanOne = false;
1509         auto mvIterator = mvec.mvData.begin();
1511         // if the iterator corresponds to the scalar
1512         if(mvIterator->grade == 0){
1513             stream << mvIterator->vec.coeff(0);
1514             moreThanOne = true;
1515             ++mvIterator;
1516         }
1518         // for all other k-vectors of mvec
1519         for(; mvIterator!=mvec.mvData.end(); ++mvIterator){
1521             // call the function that covers a single k-vector
1522             traverseKVector(stream,mvIterator->vec,mvIterator->grade, moreThanOne);
1523         }
1525         if(!moreThanOne)
1526             stream << " 0 ";
1528         return stream;
1529     }
1531     // static unit multivector functions (returns a multivector with only one component)
1532 project_static_multivector_one_component
1534     //    template<typename U>
1535     //    void recursiveTraversalMultivector(std::ostream &stream, const Mvec<U> &mvec, unsigned int currentGrade, int currentIndex, std::vector<int> listBasisBlades, unsigned int lastIndex, unsigned int gradeMV, bool& moreThanOne);
1538     void temporaryFunction1();
1540 }     /// End of Namespace
1542 #endif // project_inclusion_guard