a Code for the Combination of Indirect and Direct Constraints on High Energy Physics Models Logo
Expanded.h
Go to the documentation of this file.
1 /*
2  * To change this license header, choose License Headers in Project Properties.
3  * To change this template file, choose Tools | Templates
4  * and open the template in the editor.
5  */
6 
7 /*
8  * File: Expanded.h
9  * Author: enrico
10  *
11  * Created on 16 novembre 2017, 12.02
12  */
13 
14 #ifndef EXPANDED_H
15 #define EXPANDED_H
16 
17 #include <cmath>
18 #include <vector>
19 #include <iostream>
20 #include <stdexcept>
21 
22 #include "gslpp_complex.h"
23 #include "gslpp_vector_double.h"
24 #include "gslpp_vector_complex.h"
25 #include "gslpp_matrix_complex.h"
26 #include "gslpp_matrix_double.h"
27 
28 using namespace gslpp;
29 
30 #define isMat(X,M) (std::is_same<X,matrix<M> >::value)
31 #define isVec(X,V) (std::is_same<X,vector<V> >::value)
32 #define isSc(X,S) (std::is_same<X,S>::value)
33 #define isS(X) (isSc(X,double) || isSc(X,complex))
34 #define isV(X) (isVec(X,double) || isVec(X,complex))
35 #define isM(X) (isMat(X,double) || isMat(X,complex))
36 #define isComp(X) (isSc(X,complex) || isVec(X,complex) || isMat(X,complex))
37 #define sameType(X,Y) ( (isS(X) && isS(Y)) || (isV(X) && isV(Y)) || (isM(X) && isM(Y)))
38 
55 template <class T> class Expanded {
56 public:
60  Expanded<T>();
61 
68  Expanded<T>(std::vector<T> & dinp, int minord2_i = 0);
69 
78  Expanded<T>(std::vector<std::vector<T> > & dinp, int minord2_i = 0, int minord1_i = 0);
79 
80  template <class Q>
81  typename std::enable_if<
82  isSc(T, double) ||
83  (isSc(T, complex) && isComp(Q)) ||
84  (isMat(T, double) && !isS(Q)) ||
85  (isMat(T, complex) && (isVec(Q, complex) || isMat(Q, complex))), Expanded<Q> >::type
86  operator*(const Expanded<Q>& z) const {
87  int min1 = minord1 + z.getMin1();
88  int min2 = minord2 + z.getMin2();
89  int up1 = std::min(n1, z.getN1());
90  std::vector<std::vector<Q> > res;
91  for (int i1 = 0; i1 < up1; i1++) {
92  int up2 = std::min(n2.at(i1), z.getN2().at(i1));
93  res.push_back(std::vector<Q>(up2, z.Zero()));
94  for (int i2 = 0; i2 < up2; i2++)
95  for (int j1 = 0; j1 <= i1; j1++)
96  for (int j2 = 0; j2 <= i2; j2++)
97  res[i1][i2] += data[j1][j2] * z.getOrd(i2 - j2 + z.getMin2(), i1 - j1 + z.getMin1());
98  }
99  return Expanded<Q>(res, min2, min1);
100  }
101 
102  template <class Q>
103  typename std::enable_if<
104  (!isSc(T, double) && isSc(Q, double)) ||
105  (!isS(T) && !isMat(T, double) && isMat(Q, double)) ||
106  ((isMat(T, complex) || isVec(T, complex)) && isSc(Q, complex)) ||
107  (isVec(T, complex) && isMat(Q, complex)), Expanded<T> >::type
108  operator*(const Expanded<Q>& z) const {
109  int min1 = minord1 + z.getMin1();
110  int min2 = minord2 + z.getMin2();
111  int up1 = std::min(n1, z.getN1());
112  std::vector<std::vector<T> > res;
113  for (int i1 = 0; i1 < up1; i1++) {
114  int up2 = std::min(n2.at(i1), z.getN2().at(i1));
115  res.push_back(std::vector<T>(up2, Zero()));
116  for (int i2 = 0; i2 < up2; i2++)
117  for (int j1 = 0; j1 <= i1; j1++)
118  for (int j2 = 0; j2 <= i2; j2++)
119  res[i1][i2] += data[j1][j2] * z.getOrd(i2 - j2 + z.getMin2(), i1 - j1 + z.getMin1());
120  }
121  return Expanded<T>(res, min2, min1);
122  }
123 
124  template <class Q>
125  typename std::enable_if<
126  (isVec(T, double) && isSc(Q, complex)) ||
127  (isVec(T, double) && isMat(Q, complex)), Expanded<vector<complex> > >::type
128  operator*(const Expanded<Q>& z) const {
129  int min1 = minord1 + z.getMin1();
130  int min2 = minord2 + z.getMin2();
131  int up1 = std::min(n1, z.getN1());
132  std::vector<std::vector<vector<complex> > > res;
133  for (int i1 = 0; i1 < up1; i1++) {
134  int up2 = std::min(n2.at(i1), z.getN2().at(i1));
135  res.push_back(std::vector<vector<complex> >(up2, Zero()));
136  for (int i2 = 0; i2 < up2; i2++)
137  for (int j1 = 0; j1 <= i1; j1++)
138  for (int j2 = 0; j2 <= i2; j2++)
139  res[i1][i2] += data[j1][j2] * z.getOrd(i2 - j2 + z.getMin2(), i1 - j1 + z.getMin1());
140  }
141  return Expanded<vector<complex> >(res, min2, min1);
142  }
143 
144  template <class Q>
145  typename std::enable_if<
146  (isSc(T, complex) && isVec(Q, double)) ||
147  (isMat(T, complex) && isVec(Q, double)), Expanded<vector<complex> > >::type
148  operator*(const Expanded<Q>& z) const {
149  int min1 = minord1 + z.getMin1();
150  int min2 = minord2 + z.getMin2();
151  int up1 = std::min(n1, z.getN1());
152  std::vector<std::vector<vector<complex> > > res;
153  for (int i1 = 0; i1 < up1; i1++) {
154  int up2 = std::min(n2.at(i1), z.getN2().at(i1));
155  res.push_back(std::vector<vector<complex> >(up2, z.Zero()));
156  for (int i2 = 0; i2 < up2; i2++)
157  for (int j1 = 0; j1 <= i1; j1++)
158  for (int j2 = 0; j2 <= i2; j2++)
159  res[i1][i2] += data[j1][j2] * z.getOrd(i2 - j2 + z.getMin2(), i1 - j1 + z.getMin1());
160  }
161  return Expanded<vector<complex> >(res, min2, min1);
162  }
163 
164  template <class Q>
165  typename std::enable_if<
166  (isVec(T, double) && isVec(Q, complex)) ||
167  (isVec(T, complex) && isVec(Q, double)) ||
168  (isVec(T, complex) && isVec(Q, complex)), Expanded<complex> >::type
169  operator*(const Expanded<Q>& z) const {
170  int min1 = minord1 + z.getMin1();
171  int min2 = minord2 + z.getMin2();
172  int up1 = std::min(n1, z.getN1());
173  std::vector<std::vector<complex> > res;
174  for (int i1 = 0; i1 < up1; i1++) {
175  int up2 = std::min(n2.at(i1), z.getN2().at(i1));
176  res.push_back(std::vector<complex>(up2, 0.));
177  for (int i2 = 0; i2 < up2; i2++)
178  for (int j1 = 0; j1 <= i1; j1++)
179  for (int j2 = 0; j2 <= i2; j2++)
180  res[i1][i2] += data[j1][j2] * z.getOrd(i2 - j2 + z.getMin2(), i1 - j1 + z.getMin1());
181  }
182  return Expanded<complex>(res, min2, min1);
183  }
184 
185  template <class Q>
186  typename std::enable_if<
187  (isMat(T, double) && isSc(Q, complex)), Expanded<matrix<complex> > >::type
188  operator*(const Expanded<Q>& z) const {
189  int min1 = minord1 + z.getMin1();
190  int min2 = minord2 + z.getMin2();
191  int up1 = std::min(n1, z.getN1());
192  std::vector<std::vector<matrix<complex> > > res;
193  for (int i1 = 0; i1 < up1; i1++) {
194  int up2 = std::min(n2.at(i1), z.getN2().at(i1));
195  res.push_back(std::vector<matrix<complex> >(up2, Zero()));
196  for (int i2 = 0; i2 < up2; i2++)
197  for (int j1 = 0; j1 <= i1; j1++)
198  for (int j2 = 0; j2 <= i2; j2++)
199  res[i1][i2] += data[j1][j2] * z.getOrd(i2 - j2 + z.getMin2(), i1 - j1 + z.getMin1());
200  }
201  return Expanded<matrix<complex> >(res, min2, min1);
202  }
203 
204  template <class Q>
205  typename std::enable_if<
206  (isSc(T, complex) && isMat(Q, double)), Expanded<matrix<complex> > >::type
207  operator*(const Expanded<Q>& z) const {
208  int min1 = minord1 + z.getMin1();
209  int min2 = minord2 + z.getMin2();
210  int up1 = std::min(n1, z.getN1());
211  std::vector<std::vector<matrix<complex> > > res;
212  for (int i1 = 0; i1 < up1; i1++) {
213  int up2 = std::min(n2.at(i1), z.getN2().at(i1));
214  res.push_back(std::vector<matrix<complex> >(up2, z.Zero()));
215  for (int i2 = 0; i2 < up2; i2++)
216  for (int j1 = 0; j1 <= i1; j1++)
217  for (int j2 = 0; j2 <= i2; j2++)
218  res[i1][i2] += data[j1][j2] * z.getOrd(i2 - j2 + z.getMin2(), i1 - j1 + z.getMin1());
219  }
220  return Expanded<matrix<complex> >(res, min2, min1);
221  }
222 
223  template <class Q>
224  typename std::enable_if<
225  isVec(T, double) && isVec(Q, double), Expanded<double> >::type
226  operator*(const Expanded<Q>& z) const {
227  int min1 = minord1 + z.getMin1();
228  int min2 = minord2 + z.getMin2();
229  int up1 = std::min(n1, z.getN1());
230  std::vector<std::vector<double> > res;
231  for (int i1 = 0; i1 < up1; i1++) {
232  int up2 = std::min(n2.at(i1), z.getN2().at(i1));
233  res.push_back(std::vector<double>(up2, 0.));
234  for (int i2 = 0; i2 < up2; i2++)
235  for (int j1 = 0; j1 <= i1; j1++)
236  for (int j2 = 0; j2 <= i2; j2++)
237  res[i1][i2] += data[j1][j2] * z.getOrd(i2 - j2 + z.getMin2(), i1 - j1 + z.getMin1());
238  }
239  return Expanded<double>(res, min2, min1);
240  }
241 
242  // Operator +
243 
244  template <class Q>
245  typename std::enable_if<
246  sameType(T, Q) && isComp(T), Expanded<T> >::type
247  operator+(const Expanded<Q>& z) const {
248  int min1 = std::min(minord1, z.getMin1());
249  int min2 = std::min(minord2, z.getMin2());
250  int up1 = std::min(n1 + minord1, z.getN1() + z.getMin1()) - min1;
251  std::vector<std::vector<T> > res;
252  for (int i1 = 0; i1 < up1; i1++) {
253  int up2 = std::min(n2.at(i1) + minord2, z.getN2().at(i1) + z.getMin2()) - min2;
254  res.push_back(std::vector<T>(up2, z.Zero()));
255  for (int i2 = 0; i2 < up2; i2++) {
256  if (i1 + min1 >= minord1 && i2 + min2 >= minord2) res[i1][i2] += data[i1 + min1 - minord1][i2 + min2 - minord2];
257  if (i1 + min1 >= z.getMin1() && i2 + min2 >= z.getMin2()) res[i1][i2] += z.getOrd(i2 + min2, i1 + min1);
258  }
259  }
260  return Expanded<T>(res, min2, min1);
261  }
262 
263  template <class Q>
264  typename std::enable_if<
265  sameType(T, Q) && !isComp(T), Expanded<Q> >::type
266  operator+(const Expanded<Q>& z) const {
267  int min1 = std::min(minord1, z.getMin1());
268  int min2 = std::min(minord2, z.getMin2());
269  int up1 = std::min(n1 + minord1, z.getN1() + z.getMin1()) - min1;
270  std::vector<std::vector<Q> > res;
271  for (int i1 = 0; i1 < up1; i1++) {
272  int up2 = std::min(n2.at(i1) + minord2, z.getN2().at(i1) + z.getMin2()) - min2;
273  res.push_back(std::vector<Q>(up2, z.Zero()));
274  for (int i2 = 0; i2 < up2; i2++) {
275  if (i1 + min1 >= minord1 && i2 + min2 >= minord2) res[i1][i2] += data[i1 + min1 - minord1][i2 + min2 - minord2];
276  if (i1 + min1 >= z.getMin1() && i2 + min2 >= z.getMin2()) res[i1][i2] += z.getOrd(i2 + min2, i1 + min1);
277  }
278  }
279  return Expanded<Q>(res, min2, min1);
280  }
281 
282  // Operator -
283 
284  template <class Q>
285  typename std::enable_if<
286  sameType(T, Q) && isComp(T), Expanded<T> >::type
287  operator-(const Expanded<Q>& z) const {
288  int min1 = std::min(minord1, z.getMin1());
289  int min2 = std::min(minord2, z.getMin2());
290  int up1 = std::min(n1 + minord1, z.getN1() + z.getMin1()) - min1;
291  std::vector<std::vector<T> > res;
292  for (int i1 = 0; i1 < up1; i1++) {
293  int up2 = std::min(n2.at(i1) + minord2, z.getN2().at(i1) + z.getMin2()) - min2;
294  res.push_back(std::vector<T>(up2, z.Zero()));
295  for (int i2 = 0; i2 < up2; i2++) {
296  if (i1 + min1 >= minord1 && i2 + min2 >= minord2) res[i1][i2] += data[i1 + min1 - minord1][i2 + min2 - minord2];
297  if (i1 + min1 >= z.getMin1() && i2 + min2 >= z.getMin2()) res[i1][i2] -= z.getOrd(i2 + min2, i1 + min1);
298  }
299  }
300  return Expanded<T>(res, min2, min1);
301  }
302 
303  template <class Q>
304  typename std::enable_if<
305  sameType(T, Q) && !isComp(T), Expanded<Q> >::type
306  operator-(const Expanded<Q>& z) const {
307  int min1 = std::min(minord1, z.getMin1());
308  int min2 = std::min(minord2, z.getMin2());
309  int up1 = std::min(n1 + minord1, z.getN1() + z.getMin1()) - min1;
310  std::vector<std::vector<Q> > res;
311  for (int i1 = 0; i1 < up1; i1++) {
312  int up2 = std::min(n2.at(i1) + minord2, z.getN2().at(i1) + z.getMin2()) - min2;
313  res.push_back(std::vector<Q>(up2, z.Zero()));
314  for (int i2 = 0; i2 < up2; i2++) {
315  if (i1 + min1 >= minord1 && i2 + min2 >= minord2) res[i1][i2] += data[i1 + min1 - minord1][i2 + min2 - minord2];
316  if (i1 + min1 >= z.getMin1() && i2 + min2 >= z.getMin2()) res[i1][i2] -= z.getOrd(i2 + min2, i1 + min1);
317  }
318  }
319  return Expanded<Q>(res, min2, min1);
320  }
321 
322  Expanded<T> inverse() const {
323  throw std::runtime_error("Expanded::inverse: method not implemented");
324  }
325 
326  Expanded<T> operator-() const {
327  std::vector<std::vector<T> > res;
328  for (int i1 = 0; i1 < n1; i1++) {
329  res.push_back(std::vector<T>(n2.at(i1), Zero()));
330  for (int i2 = 0; i2 < n2.at(i1); i2++)
331  res[i1][i2] = -data[i1][i2];
332  }
333  return Expanded<T>(res, minord2, minord1);
334  }
335 
336  template <class Q>
337  typename std::enable_if< isSc(Q,double) , Expanded<T> >::type
338  operator/(const Expanded<Q>& z) const {
339  return (*this)*z.inverse();
340  }
341 
342  template <class Q>
343  typename std::enable_if< isSc(Q,complex) && isS(T) , Expanded<complex> >::type
344  operator/(const Expanded<Q>& z) const {
345  return (*this)*z.inverse();
346  }
347 
348  template <class Q>
349  typename std::enable_if< isSc(Q,complex) && isV(T) , Expanded<vector<complex> > >::type
350  operator/(const Expanded<Q>& z) const {
351  return (*this)*z.inverse();
352  }
353 
354  template <class Q>
355  typename std::enable_if< isSc(Q,complex) && isM(T) , Expanded<matrix<complex> > >::type
356  operator/(const Expanded<Q>& z) const {
357  return (*this)*z.inverse();
358  }
359 
360  bool operator==(const Expanded<T>& z) const {
361 
362  if (minord1 != z.getMin1() || minord2 != z.getMin2() || n1 != z.getN1()) return false;
363 
364  for (int i1 = 0; i1 < n1; i1++)
365  {
366  if(n2.at(i1) != z.getN2().at(i1)) return false;
367  for (int i2 = 0; i2 < n2.at(i1); i2++)
368  if (z.getOrd(i2 + z.getMin2(), i1 + z.getMin1()) != data[i1][i2]) return false;
369  }
370 
371  return true;
372  }
373 
374  bool operator!=(const Expanded<T>& z) const {
375  return !(z == *this);
376  }
377 
378  /***************** End Expanded-Expanded *****************/
379 
380  /***************** Begin Expanded * UnExpanded *****************/
381 
382  // Return T.
383 
384  template <class Q>
385  typename std::enable_if<
386  (std::is_same<T, Q>::value && !isV(T)) ||
387  isSc(Q, double) ||
388  (!isS(T) && isMat(Q, double)) ||
389  ((isMat(T, complex) || isVec(T, complex)) && isS(Q)) ||
390  (isVec(T, complex) && isMat(Q, complex)) ||
391  (isV(T) && isSc(Q, double)) ||
392  (isMat(T, double) && isSc(Q, double)), Expanded<T> >::type
393  operator*(const Q& z) const {
394  std::vector<std::vector<T> > res;
395  for (int i1 = 0; i1 < n1; i1++) {
396  res.push_back(std::vector<T>(n2.at(i1), Zero()));
397  for (int i2 = 0; i2 < n2.at(i1); i2++)
398  res[i1][i2] = data[i1][i2] * z;
399  }
400  return Expanded<T>(res, minord2, minord1);
401  }
402 
403  // Return double
404 
405  template <class Q>
406  typename std::enable_if<
407  isVec(T, double) && isVec(Q, double), Expanded<double> >::type
408  operator*(const Q& z) const {
409  std::vector<std::vector<double> > res;
410  for (int i1 = 0; i1 < n1; i1++) {
411  res.push_back(std::vector<double>(n2.at(i1), 0.));
412  for (int i2 = 0; i2 < n2.at(i1); i2++)
413  res[i1][i2] = data[i1][i2] * z;
414  }
415  return Expanded<double>(res, minord2, minord1);
416  }
417 
418  // Return complex
419 
420  template <class Q>
421  typename std::enable_if<
422  (isSc(T, double) && isSc(Q, complex)) ||
423  (isV(T) && isV(Q) && !(isVec(T, double) && isVec(Q, double))), Expanded<complex> >::type
424  operator*(const Q& z) const {
425  std::vector<std::vector<complex> > res;
426  for (int i1 = 0; i1 < n1; i1++) {
427  res.push_back(std::vector<complex>(n2.at(i1), 0.));
428  for (int i2 = 0; i2 < n2.at(i1); i2++)
429  res[i1][i2] = data[i1][i2] * z;
430  }
431  return Expanded<complex>(res, minord2, minord1);
432  }
433  // Return vector<double>
434 
435  template <class Q>
436  typename std::enable_if<
437  isVec(Q, double) && (isSc(T, double) || isMat(T, double)), Expanded<vector<double> > >::type
438  operator*(const Q& z) const {
439  std::vector<std::vector<vector<double> > > res;
440  for (int i1 = 0; i1 < n1; i1++) {
441  res.push_back(std::vector<vector<double> >(n2.at(i1), vector<double>(z.size(), 0.)));
442  for (int i2 = 0; i2 < n2.at(i1); i2++)
443  res[i1][i2] = data[i1][i2] * z;
444  }
445  return Expanded<vector<double> >(res, minord2, minord1);
446  }
447 
448 
449  // Return vector<complex>
450 
451  template <class Q>
452  typename std::enable_if<
453  (isVec(Q, complex) && (isS(T) || isM(T))) ||
454  (isVec(Q, double) && (isSc(T, complex) || isMat(T, complex))), Expanded<vector<complex> > >::type
455  operator*(const Q& z) const {
456  std::vector<std::vector<vector<complex> > > res;
457  for (int i1 = 0; i1 < n1; i1++) {
458  res.push_back(std::vector<vector<complex> >(n2.at(i1), vector<complex>(z.size(), 0.)));
459  for (int i2 = 0; i2 < n2.at(i1); i2++)
460  res[i1][i2] = data[i1][i2] * z;
461  }
462  return Expanded<vector<complex> >(res, minord2, minord1);
463  }
464 
465  template <class Q>
466  typename std::enable_if<
467  isVec(T, double) && (isSc(Q, complex) || isMat(Q, complex)), Expanded<vector<complex> > >::type
468  operator*(const Q& z) const {
469  std::vector<std::vector<vector<complex> > > res;
470  for (int i1 = 0; i1 < n1; i1++) {
471  res.push_back(std::vector<vector<complex> >(n2.at(i1), vector<complex>(data[0][0].size(), 0.)));
472  for (int i2 = 0; i2 < n2.at(i1); i2++)
473  res[i1][i2] = data[i1][i2] * z;
474  }
475  return Expanded<vector<complex> >(res, minord2, minord1);
476  }
477 
478  // Return matrix<double>
479 
480  template <class Q>
481  typename std::enable_if<
482  (isSc(T, double) && isMat(Q, double)), Expanded<matrix<double> > >::type
483  operator*(const Q& z) const {
484  std::vector<std::vector<matrix<double> > > res;
485  for (int i1 = 0; i1 < n1; i1++) {
486  res.push_back(std::vector<matrix<double> >(n2.at(i1), matrix<double>(z.size_i(), z.size_j(), 0.)));
487  for (int i2 = 0; i2 < n2.at(i1); i2++)
488  res[i1][i2] = data[i1][i2] * z;
489  }
490  return Expanded<matrix<double> >(res, minord2, minord1);
491  }
492 
493  // Return matrix<complex>
494 
495  template <class Q>
496  typename std::enable_if<
497  (isS(T) && isMat(Q, complex)) ||
498  (isSc(T, complex) && isMat(Q, double)) ||
499  (isMat(T, double) && isMat(Q, complex)), Expanded<matrix<complex> > >::type
500  operator*(const Q& z) const {
501  std::vector<std::vector<matrix<complex> > > res;
502  for (int i1 = 0; i1 < n1; i1++) {
503  res.push_back(std::vector<matrix<complex> >(n2.at(i1), matrix<complex>(z.size_i(), z.size_j(), 0.)));
504  for (int i2 = 0; i2 < n2.at(i1); i2++)
505  res[i1][i2] = data[i1][i2] * z;
506  }
507  return Expanded<matrix<complex> >(res, minord2, minord1);
508  }
509 
510  template <class Q>
511  typename std::enable_if<
512  isMat(T, double) && isSc(Q, complex), Expanded<matrix<complex> > >::type
513  operator*(const Q& z) const {
514  std::vector<std::vector<matrix<complex> > > res;
515  for (int i1 = 0; i1 < n1; i1++) {
516  res.push_back(std::vector<matrix<complex> >(n2.at(i1), matrix<complex>(data[0][0].size_i(), data[0][0].size_j(), 0.)));
517  for (int i2 = 0; i2 < n2.at(i1); i2++)
518  res[i1][i2] = data[i1][i2] * z;
519  }
520  return Expanded<matrix<complex> >(res, minord2, minord1);
521  }
522 
523  /***************** End Expanded * UnExpanded *****************/
524 
525  /***************** Begin Expanded / UnExpanded-Scalar *****************/
526 
527 
528  Expanded<T> operator/(const double& z) const {
529  std::vector<std::vector<T> > res;
530  for (int i1 = 0; i1 < n1; i1++) {
531  res.push_back(std::vector<T>(n2.at(i1), Zero()));
532  for (int i2 = 0; i2 < n2.at(i1); i2++)
533  res[i1][i2] = data[i1][i2] / z;
534  }
535  return Expanded<T>(res, minord2, minord1);
536  }
537 
538  template <class Q>
539  typename std::enable_if<isComp(T) && isSc(Q, complex), Expanded<T> >::type
540  operator/(const Q& z) const {
541  std::vector<std::vector<T> > res;
542  for (int i1 = 0; i1 < n1; i1++) {
543  res.push_back(std::vector<T>(n2.at(i1), Zero()));
544  for (int i2 = 0; i2 < n2.at(i1); i2++)
545  res[i1][i2] = data[i1][i2] / z;
546  }
547  return Expanded<T>(res, minord2, minord1);
548  }
549 
550  template <class Q>
551  typename std::enable_if<isSc(T, double) && isSc(Q, complex), Expanded<complex > >::type
552  operator/(const Q& z) const {
553  std::vector<std::vector<complex> > res;
554  for (int i1 = 0; i1 < n1; i1++) {
555  res.push_back(std::vector<complex>(n2.at(i1), 0.));
556  for (int i2 = 0; i2 < n2.at(i1); i2++)
557  res[i1][i2] = data[i1][i2] / z;
558  }
559  return Expanded<complex>(res, minord2, minord1);
560  }
561 
562  template <class Q>
563  typename std::enable_if<isVec(T, double) && isSc(Q, complex), Expanded<vector<complex> > >::type
564  operator/(const Q& z) const {
565  std::vector<std::vector<vector<complex> > > res;
566  for (int i1 = 0; i1 < n1; i1++) {
567  res.push_back(std::vector<vector<complex> >(n2.at(i1), vector<complex>(data[0][0].size(), 0.)));
568  for (int i2 = 0; i2 < n2.at(i1); i2++)
569  res[i1][i2] = data[i1][i2] / z;
570  }
571  return Expanded<vector<complex> >(res, minord2, minord1);
572  }
573 
574  template <class Q>
575  typename std::enable_if<isMat(T, double) && isSc(Q, complex), Expanded<matrix<complex> > >::type
576  operator/(const Q& z) const {
577  std::vector<std::vector<matrix<complex> > > res;
578  for (int i1 = 0; i1 < n1; i1++) {
579  res.push_back(std::vector<matrix<complex> >(n2.at(i1), matrix<complex>(data[0][0].size_i(), data[0][0].size_j(), 0.)));
580  for (int i2 = 0; i2 < n2.at(i1); i2++)
581  res[i1][i2] = data[i1][i2] / z;
582  }
583  return Expanded<matrix<complex> >(res, minord2, minord1);
584  }
585 
586  /***************** End Expanded / UnExpanded-Scalar *****************/
587 
588  friend Expanded<T> operator*(const double& ue, const Expanded<T>& ex) {
589  return (ex * ue);
590  }
591 
596  Expanded<double> real() const;
597 
602  Expanded<double> abs2() const;
603 
608  Expanded<T> conjugate() const;
609 
615  Expanded<T> transpose() const;
616 
624  Expanded<T> truncate(std::vector<int> n2max, int n1max) {
625  if (n1max < minord1 || n1max >= minord1 + n1)
626  throw std::runtime_error("Expanded::truncate(): order of the outer expansion not present in Expanded");
627  std::vector<std::vector<T> > res;
628  for (int i1 = 0; i1 <= n1max - minord1; i1++) {
629  std::vector<T> tmp;
630  if (n2max.at(i1) >= minord2 + n2.at(i1))
631  throw std::runtime_error("Expanded::truncate(): order of the inner expansion not present in Expanded");
632  for (int j1 = 0; j1 <= n2max.at(i1) - minord2; j1++)
633  tmp.push_back(data[i1][j1]);
634  res.push_back(tmp);
635  }
636  return Expanded<T>(res, minord2, minord1);
637  }
643  Expanded<T> truncate(int n2max) {
644  if(n1 >= 1 || minord1 > 0)
645  throw std::runtime_error("Expanded::truncate(): wrong method, please use truncate(std::vector<int>, int)");
646  return truncate(std::vector<int>(1, n2max), 0);
647  }
648 
657  T Series(std::vector<int> ord2, int ord1, double als2 = 1., double als1 = 1.) {
658  T res;
659  if (ord1 >= minord1)
660  {
661  for (int i = 0; i <= std::min(n1, ord1 - minord1); i++)
662  {
663  if (ord2.at(i) >= minord2)
664  for (int j = 0; j <= std::min(ord2.at(i) - minord2, n2.at(i)); j++)
665  res += data[i][j] * pow(als1, (double) (minord1 + i)) * pow(als2, (double) (minord2 + j));
666  }
667  }
668  return res;
669  }
670 
677  T Series(int ord2, double als2 = 1.) {
678  return Series(std::vector<int> (1, ord2), 0, als2);
679  }
680 
685  T Zero() const {
686  return T();
687  }
688 
694  const T& getOrd(int j) const {
695  return getOrd(j, minord1);
696  }
697 
704  const T& getOrd(int j, int i) const {
705  checkOrd(j, i, "Expanded::getOrd(): order non present in Expanded");
706  return data[i - minord1][j - minord2];
707  }
708 
713  void setOrd(int j, T value) {
714  setOrd(j, minord1, value);
715  }
716 
722  void setOrd(int j, int i, T value) {
723  checkOrd(j, i, "Expanded::setOrd(): order non present in Expanded");
724  data[i - minord1][j - minord2] = value;
725  }
726 
727  void checkOrd(int j, int i, std::string s) const {
728  if (i < minord1 || i >= minord1 + n1 || j < minord2 || j >= minord2 + n2.at(i - minord1))
729  throw std::runtime_error(s);
730  }
731 
740  template <class Q> typename std::enable_if<isMat(T, Q), void>::type
741  setMatrixElement(int j, int i, int h, int k, Q x) {
742  checkOrd(j, i, "Expanded::setMatrixElement(): order non present in Expanded");
743  data[i - minord1][j - minord2].assign(h, k, x);
744  }
752  template <class Q> typename std::enable_if<isMat(T, Q), void>::type
753  setMatrixElement(int j, int h, int k, Q x) {
754  setMatrixElement(j, minord1, h, k, x);
755  }
756 
764  template <class Q> typename std::enable_if<isVec(T, Q), void>::type
765  setVectorElement(int j, int i, int h, Q x) {
766  checkOrd(j, i, "Expanded::setVectorElement(): order non present in Expanded");
767  //THERE IS NO ASSIGN METHOD FOR VECTORS
768  //data[i - minord1][j - minord2].assign(h, x);
769  data[i - minord1][j - minord2](h) = x;
770  }
771 
778  template <class Q> typename std::enable_if<isVec(T, Q), void>::type
779  setVectorElement(int j, int h, Q x) {
780  setVectorElement(j, minord1, h, x);
781  }
782 
787  const int getMin1() const {
788  return minord1;
789  }
790 
795  const int getMin2() const {
796  return minord2;
797  }
798 
803  const int getN1() const {
804  return n1;
805  }
806 
812  const std::vector<int>& getN2() const {
813  return n2;
814  }
815 
816  friend std::ostream& operator<<(std::ostream& output, const Expanded<T> z) {
817  for (int i = 0; i < z.getN1(); i++)
818  for (int j = 0; j < z.getN2().at(i); j++) {
819  output << "Order (" << z.getMin2() + j << "," << z.getMin1() + i << "): " << z.getOrd(j + z.getMin2(), i + z.getMin1());
820  if (i != (z.getN1() - 1) || j != (z.getN2().at(i) - 1))
821  output << std::endl;
822  }
823  return output;
824  }
825 
826 
827 private:
832  T invCoeff(int i, int j) const {
833  throw std::runtime_error("Expanded::invCoeff: method not implemented");
834  }
835 
836  int minord1, minord2, n1;
837  std::vector<int> n2;
838  std::vector<std::vector<T> > data;
839 };
840 
841 
842 // Template specialization
843 
844 // scalar * expanded
845 
846 template <class T>
847 typename std::enable_if<
848 isComp(T), Expanded<T> >::type
849 operator*(const complex& ue, const Expanded<T>& ex) {
850  return (ex * ue);
851 }
852 
853 // vector * expanded scalar
854 
855 template <class T>
856 typename std::enable_if<
857 isS(T), Expanded<vector<complex> > >::type
858 operator*(const vector<complex>& ue, const Expanded<T>& ex) {
859  return (ex * ue);
860 }
861 
862 // vector * vector
863 
864 template <class T, class Q>
865 typename std::enable_if<
866 isV(T) && isV(Q) && !(isVec(T, double) && isVec(Q, double)), Expanded<complex> >::type
867 operator*(const T& ue, const Expanded<Q>& ex) {
868  return (ex * ue);
869 }
870 
871 // vector * matrix
872 
873 template <class T, class Q>
874 typename std::enable_if<
875 isV(T) && isM(Q) && !(isVec(T, double) && isMat(Q, double)), Expanded<vector<complex> > >::type
876 operator*(const T& ue, const Expanded<Q>& ex) {
877  return (ex.transpose() * ue);
878 }
879 
880 // matrix * expanded vector
881 
882 template <class T, class Q>
883 typename std::enable_if<
884 isM(T) && isV(Q) && !(isMat(T, double) && isVec(Q, double)), Expanded<vector<complex> > >::type
885 operator*(const T& ue, const Expanded<Q>& ex) {
886  return (ex * ue.transpose());
887 }
888 
889 // matrix * expanded matrix
890 
891 template <class T, class Q>
892 typename std::enable_if<
893 isM(T) && isM(Q) && !(isMat(T, double) && isMat(Q, double)), Expanded<matrix<complex> > >::type
894 operator*(const T& ue, const Expanded<Q>& ex) {
895  return ((ex.transpose() * ue.transpose()).transpose());
896 }
897 // matrix * expanded scalar
898 
899 template <class T>
900 typename std::enable_if<
901 isS(T), Expanded<matrix<complex> > >::type
902 operator*(const matrix<complex>& ue, const Expanded<T>& ex) {
903  return (ex * ue);
904 }
905 
906 template <> Expanded<complex> Expanded<complex>::conjugate() const;
907 
908 template <> Expanded<matrix<double> > Expanded<matrix<double> >::transpose() const;
909 
910 template <> Expanded<matrix<complex> > Expanded<matrix<complex> >::transpose() const;
911 
912 template <> Expanded<double> Expanded<complex>::real() const;
913 
914 template <> Expanded<double> Expanded<complex>::abs2() const;
915 
916 template <> matrix<double> Expanded<matrix<double> >::Zero() const;
917 
918 template <> matrix<complex> Expanded<matrix<complex> >::Zero() const;
919 
920 template <> vector<double> Expanded<vector<double> >::Zero() const;
921 
922 template<> vector<complex> Expanded<vector<complex> >::Zero() const;
923 
924 /***************** Begin UnExpanded * Expanded *****************/
925 
926 Expanded<complex> operator*(const complex& ue, const Expanded<double>& ex);
927 
928 Expanded<vector<complex> > operator*(const complex& ue, const Expanded<vector<double> >& ex);
929 
930 Expanded<matrix<complex> > operator*(const complex& ue, const Expanded<matrix<double> >& ex);
931 
932 // X * scalar -> commute
933 
934 Expanded<vector<double> > operator*(const vector<double>& ue, const Expanded<double>& ex);
935 
936 Expanded<vector<complex> > operator*(const vector<double>& ue, const Expanded<complex>& ex);
937 
938 Expanded<matrix<double> > operator*(const matrix<double>& ue, const Expanded<double>& ex);
939 
940 Expanded<matrix<complex> > operator*(const matrix<double>& ue, const Expanded<complex>& ex);
941 
942 // both vector -> commute
943 
944 Expanded<double> operator*(const vector<double>& ue, const Expanded<vector<double> >& ex);
945 
946 // vector * matrix
947 
948 Expanded<vector<double> > operator*(const vector<double>& ue, const Expanded<matrix<double> >& ex);
949 
950 // matrix * vector
951 
952 Expanded<vector<double> > operator*(const matrix<double>& ue, const Expanded<vector<double> >& ex);
953 
954 // matrix * matrix
955 
956 Expanded<matrix<double> > operator*(const matrix<double>& ue, const Expanded<matrix<double> >& ex);
957 
958 /***************** End UnExpanded * Expanded *****************/
959 
960 #endif /* EXPANDED_H */
Expanded::real
Expanded< double > real() const
Method to return the real part of a complex Expanded.
Expanded::getOrd
const T & getOrd(int j) const
Get an element of a single expansion.
Definition: Expanded.h:694
gslpp::vector< complex >
A class for constructing and defining operations on complex vectors.
Definition: gslpp_vector_complex.h:33
Expanded::setMatrixElement
std::enable_if< isMat(T, Q), void >::type setMatrixElement(int j, int i, int h, int k, Q x)
Set an element of a matrix in a double Expanded<matrix<*> >
Definition: Expanded.h:741
Expanded::operator/
Expanded< T > operator/(const double &z) const
Definition: Expanded.h:528
Expanded::n2
std::vector< int > n2
Definition: Expanded.h:837
gslpp::matrix< double >
A class for constructing and defining operations on real matrices.
Definition: gslpp_matrix_double.h:48
Expanded::getMin2
const int getMin2() const
Get the minimum order of the inner expansion.
Definition: Expanded.h:795
Expanded::setOrd
void setOrd(int j, int i, T value)
Set an element of a double expansion.
Definition: Expanded.h:722
Expanded::operator/
std::enable_if< isSc(T, double) &&isSc(Q, complex), Expanded< complex > >::type operator/(const Q &z) const
Definition: Expanded.h:552
Expanded::operator+
std::enable_if< sameType(T, Q) &&!isComp(T), Expanded< Q > >::type operator+(const Expanded< Q > &z) const
Definition: Expanded.h:266
Expanded::operator!=
bool operator!=(const Expanded< T > &z) const
Definition: Expanded.h:374
Expanded::operator/
std::enable_if< isComp(T) &&isSc(Q, complex), Expanded< T > >::type operator/(const Q &z) const
Definition: Expanded.h:540
gslpp_matrix_complex.h
Expanded::conjugate
Expanded< T > conjugate() const
Method to return the conjugate value of a Expanded<complex>.
Expanded::operator-
Expanded< T > operator-() const
Definition: Expanded.h:326
Expanded::Series
T Series(int ord2, double als2=1.)
Method to sum an expansion up to a given order.
Definition: Expanded.h:677
Expanded::abs2
Expanded< double > abs2() const
Method to return the squared absolute value of a Expanded<complex>.
gslpp::complex
A class for defining operations on and functions of complex numbers.
Definition: gslpp_complex.h:35
Expanded::operator*
std::enable_if< isVec(T, double) &&isVec(Q, double), Expanded< double > >::type operator*(const Q &z) const
Definition: Expanded.h:408
Expanded::operator*
std::enable_if<(isVec(Q, complex) &&(isS(T)||isM(T)))||(isVec(Q, double) &&(isSc(T, complex)||isMat(T, complex))), Expanded< vector< complex > > >::type operator*(const Q &z) const
Definition: Expanded.h:455
Expanded::getN1
const int getN1() const
Get the number of terms in the outer expansion.
Definition: Expanded.h:803
Expanded::transpose
Expanded< T > transpose() const
Method to transpose a Expanded<matrix<double> > or a Expanded<matrix<complex> >.
Expanded::operator==
bool operator==(const Expanded< T > &z) const
Definition: Expanded.h:360
Expanded::operator*
std::enable_if<(isSc(T, double) &&isSc(Q, complex))||(isV(T) &&isV(Q) &&!(isVec(T, double) &&isVec(Q, double))), Expanded< complex > >::type operator*(const Q &z) const
Definition: Expanded.h:424
Expanded::setMatrixElement
std::enable_if< isMat(T, Q), void >::type setMatrixElement(int j, int h, int k, Q x)
Set an element of a matrix in a single Expanded<matrix<*> >
Definition: Expanded.h:753
gslpp_matrix_double.h
Expanded::operator/
std::enable_if< isSc(Q, double), Expanded< T > >::type operator/(const Expanded< Q > &z) const
Definition: Expanded.h:338
Expanded::operator/
std::enable_if< isSc(Q, complex) &&isS(T), Expanded< complex > >::type operator/(const Expanded< Q > &z) const
Definition: Expanded.h:344
Expanded::Series
T Series(std::vector< int > ord2, int ord1, double als2=1., double als1=1.)
Method to sum a double expansion up to a given order.
Definition: Expanded.h:657
gslpp::pow
complex pow(const complex &z1, const complex &z2)
Definition: gslpp_complex.cpp:395
Expanded
A template class for Taylor double expansion of several objects.
Definition: Expanded.h:55
Expanded::operator-
std::enable_if< sameType(T, Q) &&!isComp(T), Expanded< Q > >::type operator-(const Expanded< Q > &z) const
Definition: Expanded.h:306
Expanded::setVectorElement
std::enable_if< isVec(T, Q), void >::type setVectorElement(int j, int h, Q x)
Set an element of a matrix in a single Expanded<vector<*> >
Definition: Expanded.h:779
Expanded::getMin1
const int getMin1() const
Get the minimum order of the outer expansion.
Definition: Expanded.h:787
Expanded::truncate
Expanded< T > truncate(std::vector< int > n2max, int n1max)
Method to truncate a double expansion.
Definition: Expanded.h:624
gslpp_vector_complex.h
Expanded::operator*
friend Expanded< T > operator*(const double &ue, const Expanded< T > &ex)
Definition: Expanded.h:588
Expanded::setOrd
void setOrd(int j, T value)
Set an element of a single expansion.
Definition: Expanded.h:713
Expanded::operator<<
friend std::ostream & operator<<(std::ostream &output, const Expanded< T > z)
Definition: Expanded.h:816
gslpp::vector< double >
A class for constructing and defining operations on real vectors.
Definition: gslpp_vector_double.h:33
Expanded::operator*
std::enable_if< isSc(T, double)||(isSc(T, complex) &&isComp(Q))||(isMat(T, double) &&!isS(Q))||(isMat(T, complex) &&(isVec(Q, complex)||isMat(Q, complex))), Expanded< Q > >::type operator*(const Expanded< Q > &z) const
Definition: Expanded.h:86
Expanded::operator-
std::enable_if< sameType(T, Q) &&isComp(T), Expanded< T > >::type operator-(const Expanded< Q > &z) const
Definition: Expanded.h:287
Expanded::n1
int n1
Definition: Expanded.h:836
Expanded::operator/
std::enable_if< isSc(Q, complex) &&isM(T), Expanded< matrix< complex > > >::type operator/(const Expanded< Q > &z) const
Definition: Expanded.h:356
Expanded::operator+
std::enable_if< sameType(T, Q) &&isComp(T), Expanded< T > >::type operator+(const Expanded< Q > &z) const
Definition: Expanded.h:247
Expanded::inverse
Expanded< T > inverse() const
Definition: Expanded.h:322
Expanded::getOrd
const T & getOrd(int j, int i) const
Get an element of a double expansion.
Definition: Expanded.h:704
Expanded::getN2
const std::vector< int > & getN2() const
Get the number of terms in the inner expansion.
Definition: Expanded.h:812
Expanded::data
std::vector< std::vector< T > > data
Definition: Expanded.h:838
Expanded::operator/
std::enable_if< isSc(Q, complex) &&isV(T), Expanded< vector< complex > > >::type operator/(const Expanded< Q > &z) const
Definition: Expanded.h:350
gslpp_complex.h
gslpp::operator*
complex operator*(const double &x1, const complex &z2)
Definition: gslpp_complex.cpp:314
Expanded::operator/
std::enable_if< isMat(T, double) &&isSc(Q, complex), Expanded< matrix< complex > > >::type operator/(const Q &z) const
Definition: Expanded.h:576
Expanded::setVectorElement
std::enable_if< isVec(T, Q), void >::type setVectorElement(int j, int i, int h, Q x)
Set an element of a matrix in a double Expanded<vector<*> >
Definition: Expanded.h:765
Expanded::checkOrd
void checkOrd(int j, int i, std::string s) const
Definition: Expanded.h:727
Expanded::operator*
std::enable_if< isMat(T, double) &&isSc(Q, complex), Expanded< matrix< complex > > >::type operator*(const Q &z) const
Definition: Expanded.h:513
Expanded::invCoeff
T invCoeff(int i, int j) const
Return the (i-th, j-th) coefficient of the expansion of 1/Expanded.
Definition: Expanded.h:832
Expanded::truncate
Expanded< T > truncate(int n2max)
Method to truncate a single expansion.
Definition: Expanded.h:643
gslpp
Complex number, vector and matrix manipulation using GSL.
Definition: gslpp_complex.cpp:16
gslpp_vector_double.h
Expanded::Zero
T Zero() const
Return an empty instance of class T.
Definition: Expanded.h:685
gslpp::matrix< complex >
A class for constructing and defining operations on complex matrices.
Definition: gslpp_matrix_complex.h:45
Expanded::operator/
std::enable_if< isVec(T, double) &&isSc(Q, complex), Expanded< vector< complex > > >::type operator/(const Q &z) const
Definition: Expanded.h:564
Expanded::operator*
std::enable_if< isVec(T, double) &&isVec(Q, double), Expanded< double > >::type operator*(const Expanded< Q > &z) const
Definition: Expanded.h:226