finitediff
SNoperators.h
Go to the documentation of this file.
1 /*
2 Copyright 2017 Laurent Claessens
3 contact : laurent@claessens-donadello.eu
4 
5 This program is free software: you can redistribute it and/or modify
6 it under the terms of the GNU General Public License as published by
7 the Free Software Foundation, either version 3 of the License, or
8 (at your option) any later version.
9 
10 This program is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU General Public License for more details.
14 
15 You should have received a copy of the GNU General Public License
16 along with this program. If not, see <http://www.gnu.org/licenses/>.
17 */
18 
19 
20 /*
21  This file contains operators *,== between different types of matrices.
22 */
23 
24 #ifndef OPERATORS_H__064802_
25 #define OPERATORS_H__064802_
26 
27 #include "SNpermutation.h"
28 #include "SNgaussian.h"
29 #include "SNmultiGaussian.h"
30 #include "SNidentity.h"
31 #include "SNlowerTriangular.h"
32 #include "SNupperTriangular.h"
33 #include "MathUtilities.h"
34 #include "../exceptions/SNexceptions.cpp"
35 
36 
37 // PRODUCTS ------------------------------------------
38 
39 
52 template <class U,class V,unsigned int s,unsigned int t>
54 {
55  tooGenericWarning("Warning : using a very generic product 'SNgeneric * SNgeneric'. Can't you be more specific ?");
57  SNmatrix<U,s> ans;
58  for (unsigned int i=0;i<s;i++)
59  {
60  for (unsigned int j=0;j<s;j++)
61  {
62  ans.at(i,j)=matrixProductComponent(A,B,i,j);
63  }
64  }
65  return ans; //relies on RVO.
66 }
67 
68 // SNgaussian * SNgeneric
69 
87 template <class U,class V,unsigned int s,unsigned int t>
88 SNmatrix<U,s> operator*
89 (const SNgaussian<U,s>& A, const SNgeneric<V,t>& B)
90 {
91 
93 
94  SNmatrix<U,s> ans;
95 
96  const unsigned int c=A.getColumn();
97 
98  copyFirstLines(ans,B,c);
99 
100  // for the other lines, one has only two non vanishing elements
101  // in the gaussian.
102  for (m_num line=c+1;line<s;++line)
103  {
104  for (m_num col=0;col<s;++col)
105  {
106  ans.at(line,col)=B.get(line,col)+A.get(line,c)*B.get(c,col);
107  }
108  }
109  return ans;
110 }
111 
112 
113 // SNgaussian * SNgaussian
114 
115 template <class U,class V,unsigned int s,unsigned int t>
116 SNmultiGaussian<U,s> operator*
117 (const SNgaussian<U,s>& A, const SNgaussian<V,t>& B)
118 {
120  unsigned int tp_size=s;
121 
123  if (A.getColumn()==B.getColumn())
124  {
125  ans.setLastColumn(A.getColumn());
126  m_num col=A.getColumn();
127  for (m_num line = col+1 ; line< tp_size ;++line )
128  {
129  ans.at(line,col)=A.get(line,col)+B.get(line,col);
130  }
131  }
132  else if (A.getColumn()>B.getColumn())
133  {
134  ans.setLastColumn(A.getColumn());
135 
136  // The `_at` function in in `SNmultigauss` automatically
137  // returns '0' when the column number is large than
138  // // `getLastColumn`. In other words, these are
139  // "special values" and attempting to access them with `_at`
140  // throws a `SNchangeNotAllowedException`.
141  for (m_num col=0;col <= ans.getLastColumn() ;++col)
142  {
143  if (col==A.getColumn())
144  {
145  for (m_num line=col+1;line<s;++line)
146  {
147  ans.at(line,col)=A.get(line,col);
148  }
149  }
150  else if (col==B.getColumn())
151  {
152  for (m_num line=B.getColumn()+1;line<A.getColumn()+1;++line)
153  {
154  ans.at(line,col)=B.get(line,col);
155  }
156  for (m_num line=A.getColumn()+1;line<s;++line)
157  {
158  ans.at(line,col)=A.get(line,A.getColumn())*B.get(A.getColumn(),col)+B.get(line,B.getColumn());
159  }
160  }
161  else
162  {
163  for (m_num line=col+1;line<s;++line)
164  {
165  ans.at(line,col)=0;
166  }
167  }
168  }
169  }
170  else
171  {
172  ans.setLastColumn(B.getColumn());
173 
174  for (m_num c=0;c <= ans.getLastColumn();c++)
175  {
176  for (m_num l=c+1;l<s;l++)
177  {
178  ans.at(l,c)=A.get(l,c)+B.get(l,c);
179  }
180  }
181  }
182  return ans;
183 }
184 
185 // SNmultiGaussian * SNgaussian
186 
187 template <class U,class V,unsigned int s,unsigned int t>
188 SNmultiGaussian<U,s> operator*
190 {
192  SNmultiGaussian<U,s> ans(M);
193 
194  if (M.getLastColumn() >= G.getColumn())
195  {
196  ans.setLastColumn(M.getLastColumn());
197 
198  m_num col=G.getColumn();
199  for (m_num line=col+1;line<s;++line)
200  {
201  ans.at(line,col)=matrixProductComponent(M,G,line,col);
202  }
203  }
204  else
205  {
206  ans.setLastColumn(G.getColumn());
207 
208  for (m_num l=G.getColumn()+1;l<s;l++)
209  {
210  ans.at(l,G.getColumn())+=G.get(l,G.getColumn());
211  }
212  }
213  return ans;
214 }
215 
216 // SNmultiGaussian * SNgeneric
217 
241 template <class U,class V,unsigned int s,unsigned int t>
242 SNmatrix<U,s> operator*
244 {
245 
247  const unsigned int tp_size=M.getSize(); // for homogeneous notations.
248  const m_num last_col=M.getLastColumn();
249 
250  SNmatrix<U,s> ans;
251 
252  // copy the first line
253 
254  for (m_num col=0;col<tp_size;++col)
255  {
256  ans.at(0,col)=E.get(0,col);
257  }
258 
259  // when the line number is smaller than 'last_col'
260  for (m_num line=1;line <= last_col;++line)
261  {
262  for (m_num col=0;col<tp_size;++col)
263  {
264  U acc=0;
265  for (m_num k=0; k < line;++k)
266  {
267  acc+=M.get(line,k)*E.get(k,col);
268  }
269  ans.at(line,col)=acc+E.get(line,col);
270  }
271  }
272 
273  // for the last lines
274  for (m_num line=last_col+1;line<tp_size;++line)
275  {
276  for (m_num col=0;col<tp_size;++col)
277  {
278  U acc=0;
279  for (m_num k=0;k <= last_col;++k)
280  {
281  acc+=M.get(line,k)*E.get(k,col);
282  }
283  ans.at(line,col)=acc+E.get(line,col);
284  }
285  }
286  return ans;
287 }
288 
289 
290 // SNmultiGaussian * SNmultigaussian
291 
292 template <class U,class V,unsigned int s,unsigned int t>
293 SNmultiGaussian<U,s> operator*
295 {
298  ans.setLastColumn(std::max(A.getLastColumn(),B.getLastColumn()));
299 
300  for (m_num col=0;col<ans.getLastColumn();++col)
301  {
302  for (m_num line=col+1;line<s;++line)
303  {
304  U acc=0;
305  // TODO : non optimal because the first and last products are 1*something and something*1.
306  for (m_num k=col;k <= line;++k)
307  {
308  acc+=(A.get(line,k)*B.get(k,col));
309  }
310  ans.at(line,col)=acc;
311  }
312  }
313  return ans;
314 }
315 
316 // SNgaussian * SNlowerTriangular
317 
318 template <class U,class V,unsigned int s,unsigned int t>
319 SNlowerTriangular<U,s> operator*
321 {
322 
323  // gaussian * lower trig -> lower trig
324  //
325  // Copy the first 'c' lines (only under the diagonal)
326  // Then compute the product for the other lines (only under the diagonal)
327 
329  unsigned int size=A.getSize();
330  unsigned int c=A.getColumn();
332 
333  for (unsigned int i=0;i<c+1;i++)
334  {
335  for (unsigned int j=0;j<i+1;j++)
336  {
337  ans.at(i,j)=B.get(i,j);
338  }
339  }
340  for (unsigned int i=c+1;i<size;i++)
341  {
342  for (unsigned int j=0;j<i+1;j++)
343  {
344  ans.at(i,j)=matrixProductComponent(A,B,i,j);
345  }
346  }
347  return ans;
348 }
349 
350 // SNgaussian * SNmultiGaussian
351 
376 template <class U,class V,unsigned int s,unsigned int t>
377 SNmultiGaussian<U,s> operator*
379 {
380 
382 
383  const unsigned int tp_size=G.getSize(); // for homogeneity
384  const m_num col=G.getColumn();
385  const m_num last_col=M.getLastColumn();
386 
388  ans.setLastColumn( std::max(col,last_col) );
389 
390  ans.setFirstLines(M,col);
391 
392  for (m_num i=col+1;i<tp_size;++i) // loop over the next lines
393  {
394  for (m_num j=0;j<i;++j)
395  {
396  ans.at(i,j)=M.get(i,j)+G.get(i,col)*M.get(col,j);
397  }
398  }
399  return ans;
400 }
401 
402 // Mpermutation * Mpermutation
403 
408 template <unsigned int tp_size>
409 Mpermutation<tp_size> operator*
411 {
412  Mpermutation<tp_size> new_perm;
413  for (unsigned int i=0;i<tp_size;++i)
414  {
415  new_perm.at(i)=p1.image( p2.image(i) );
416  }
417  return new_perm;
418 }
419 
429 template <unsigned int tp_size>
431 {
433  for (unsigned int k=0;k<tp_size;++k)
434  {
435  tmp.at(k)=A.image( B.image(k) );
436  }
437  return tmp;
438 }
439 
440 // SUM ---------------------------------------
441 
450 template <class U,unsigned int s,class V,unsigned int t>
452 {
453 
454  // TODO : since this sum is not commutative, maybe one has to implement it
455  // as member function.
456 
457  SNmatrix<U,s> new_matrix(A);
458  for (unsigned int k=0;k<s*s;k++)
459  {
460  new_matrix.data.at(k)+=B.data.at(k);
461  }
462  return new_matrix;
463 }
464 
465 // DIFFERENCE ---------------------------------------
466 
467 template <class U,unsigned int s,class V,unsigned int t>
469 {
471  SNmatrix<U,s> new_matrix(A);
472  for (m_num k=0;k<s;k++)
473  {
474  new_matrix.at(k,k)-=1;
475  }
476  return new_matrix;
477 }
478 
479 
480 // EQUALITIES ---------------------------------------
481 
482 template <class U,unsigned int s,class V,unsigned int t>
483 bool operator==(const SNgeneric<U,s>& A,const SNgeneric<V,t>& B)
484 {
485  return componentWiseeEquality(A,B);
486 }
487 
488 template <class U,unsigned int s,class V,unsigned int t>
489 bool operator==(const SNmatrix<U,s>& A,const SNmatrix<V,t>& B)
490 {
492  return A.data==B.data;
493 }
494 
495 template <class U,unsigned int s,class V,unsigned int t>
497 {
498  return B==A;
499 }
500 
501 
505 template <unsigned int tp_size>
507 {
508  return A.data==B.data;
509 }
510 
514 template <unsigned int tp_size>
516 {
517  for (unsigned int k=0;k<tp_size;++k)
518  {
519  if (A(k)!=B(k) )
520  {
521  return false;
522  }
523  }
524  return true;
525 }
526 
527 #endif
virtual T & at(const m_num, const m_num) final
Definition: SNgeneric.h:218
Definition: SNupperTriangular.h:34
virtual unsigned int image(const unsigned int k) const =0
Definition: SNgaussian.h:61
Represent a square numerical matrix.
Definition: SNline.h:29
void checkSizeCompatibility(const SNgeneric< U, s > &A, const SNgeneric< V, t > &B)
Definition: MathUtilities.h:35
This is the base class for the other matrices types.
Definition: MathUtilities.h:32
m_num getLastColumn() const
Definition: SNmultiGaussian.h:209
void setLastColumn(const m_num &lc)
Set the number of the last non trivial column.
Definition: SNmultiGaussian.h:215
Definition: m_num.h:34
Definition: SNidentity.h:28
A permutation is a bijection of a finite subset of N.
Definition: MgenericPermutation.h:48
Definition: SNmultiGaussian.h:42
std::array< unsigned int, tp_size > data
Definition: Mpermutation.h:52
virtual T get(const m_num, const m_num) const final
Definition: SNgeneric.h:211
SNmatrix< U, s > operator+(const SNmatrix< U, s > &A, const SNmatrix< V, t > &B)
Definition: SNoperators.h:451
m_num getColumn() const
Definition: SNgaussian.h:144
void setFirstLines(const SNgeneric< U, s > &other, const m_num max_l)
copies the first max_l lines from other to this.
Definition: SNmultiGaussian.h:228
unsigned int & at(const unsigned int k)
return by reference the image of &#39;k&#39; by the permutation
Definition: Mpermutation.h:78
unsigned int image(const unsigned int k) const override
Definition: Mpermutation.h:129
SNmatrix< U, s > operator-(const SNgeneric< U, s > &A, const SNidentity< V, t > &B)
Definition: SNoperators.h:468
Represents a lower triangular matrix (the diagonal can be non zero).
Definition: SNlowerTriangular.h:37
void tooGenericWarning(const std::string &message)
Display a small warning.
Definition: Utilities.cpp:34
bool operator==(const SNgeneric< U, s > &A, const SNgeneric< V, t > &B)
Definition: SNoperators.h:483
void copyFirstLines(SNgeneric< U, t > &ans, const SNgeneric< T, tp_size > &A, const m_num max_l)
copy the first lines of a matrix in an other one.
Definition: MathUtilities.h:114
virtual unsigned int getSize() const final
Definition: SNgeneric.h:203
bool componentWiseeEquality(const T &A, const U &B)
Definition: MathUtilities.h:45
std::array< T, tp_size *tp_size > data
Definition: SNmatrix.h:95
SNmatrix< U, s > operator*(const SNgeneric< U, s > &A, const SNgeneric< V, t > &B)
SNgeneric * SNgeneric.
Definition: SNoperators.h:53
This class represents a permutation (not a matrix).
Definition: MelementaryPermutation.h:26
T matrixProductComponent(const SNgeneric< T, s > &A, const SNgeneric< U, t > &B, unsigned int i, unsigned int j)
compute the component (i,j) of the matrix product A*B
Definition: MathUtilities.h:78