VERB_code_2.2  2
 All Classes Namespaces Files Functions Variables Typedefs Macros Pages
Matrix.h
Go to the documentation of this file.
1 /** Matrix 1D, 2D and 3D and operations with them
2  *
3  * \file Matrix.h
4  *
5  * File has 1D-class, 2D-class and 3D-class of matrixes and variouse functions to work with them.
6  *
7  * \author Developed under supervision of the PI Yuri Shprits
8  *
9  */
10 
11 #include <assert.h>
12 #include <string>
13 #include <string.h>
14 #include <fstream>
15 #include <memory.h>
16 #include <math.h>
17 #include <map>
18 
19 #ifndef matrix_array_MATRIX_H
20 #define matrix_array_MATRIX_H
21 
22 #include "./Interpolation/linear.h"
23 #include "./Interpolation/spline.h"
25 #include "./Interpolation/polint.h"
26 #include "./Interpolation/ratint.h"
27 #include "./Interpolation/akima.h"
28 
29 
30 
31 #include "../Exceptions/error.h"
32 
33 using namespace std;
34 
35 /**
36  * Matrix 1D class
37  *
38  * Matrixes and operations.
39  */
40 template <typename T>
41 class Matrix1D {
42 public:
43  T *matrix_array; ///< Array to keep the values
44 
45  bool initialized; ///< Flag, equal true if initialized
46  int size_x; ///< size x
47  string name; ///< name of the Matrix
48 
49  // constructors and destructors
50  Matrix1D() { initialized = false; };
51  Matrix1D( int size_x );
52  Matrix1D( int x_size , string name);
53  Matrix1D( const Matrix1D<T> &M );
54  ~Matrix1D();
55 
56  void AllocateMemory( int size_x );
57 
58  // Operators
59  inline T& operator[](int i); ///< Return the i-th value of matrix
60  inline T& operator[](int i) const;
61  inline T& operator()(int x); ///< Return the x-th value of matrix
62  inline T& Value (int x) { return operator()(x); } ///< Return the (x,y) value of matrix
63  inline Matrix1D<T>& MatrixArray () { return *this; } ///< Return pointer to the instance of the class.
64  inline T* MatrixArrayPointer () { return matrix_array; } ///< Return pointer to the instance of the class.
65 
66  // unary
67  inline const Matrix1D& operator+() const { return *this; }
68  inline const Matrix1D operator-() const { return ((*this)*(-1)); }
69 
70  // The following operators modify the matrix they applied to
71  inline Matrix1D& operator= (const Matrix1D<T> &M);
72  inline Matrix1D& operator= (const T val);
73 
74  // \todo Some of the matrix operators still need to be implimented
75  // I didn't have time yet to write these functions - these are matrix opearations
76  //inline Matrix1D& operator*= (const Matrix1D<T> &M); // reserved for something good
77  //inline Matrix1D& operator/= (const Matrix1D<T> &M); // reserved for something good
78  //inline Matrix1D& operator+= (const Matrix1D<T> &M);
79  //inline Matrix1D& operator-= (const Matrix1D<T> &M);
80  //inline Matrix1D& operator*= (const T Val);
81  //inline Matrix1D& operator/= (const T Val);
82  //inline Matrix1D& operator+= (const T Val); ///< Add the Val to each matrix element, stores result in the matrix it's applied to
83  //inline Matrix1D& operator-= (const T Val); ///< Substract the Val from each matrix element, stores result in the matrix it's applied to
84 
85  //inline Matrix1D& times_equal (const Matrix1D<T> &M); ///< Arraywise multiplication (A.*B), stores result in the matrix it's applied to
86  //inline Matrix1D& divide_equal (const Matrix1D<T> &M); ///< Arraywise division (A./B), stores result in the matrix it's applied to
87 
88  // The following operators save the result to a new matrix
89  //inline Matrix1D operator* (const Matrix1D<T> &M) const; // reserved for something good
90  //inline Matrix1D operator/ (const Matrix1D<T> &M) const; // reserved for something good
91  //inline Matrix1D operator+ (const Matrix1D<T> &M) const;
92  //inline Matrix1D operator- (const Matrix1D<T> &M) const;
93  inline Matrix1D operator* (const T Val) const;
94  inline Matrix1D operator/ (const T Val) const;
95  //inline Matrix1D operator+ (const T Val) const; ///< Add the Val to each matrix element, stores result in a new matrix
96  //inline Matrix1D operator- (const T Val) const; ///< Substract the Val from each matrix element, stores result in a new matrix
97 
98  inline Matrix1D times (const Matrix1D<T> &M) const; ///< Arraywise multiplication (A.*B), stores result in a new matrix
99  inline Matrix1D divide (const Matrix1D<T> &M) const; ///< Arraywise division (A./B), stores result in a new matrix
100 
101  inline T dot (const Matrix1D<T> &M) const; ///< Dot product
102  inline T norm () const; ///< Norm
103 
104 
105  // interpolations
106  void Akima_interpolation(Matrix1D<T> &old_function, Matrix1D<T> &old_grid, Matrix1D<T> &new_grid, int extrapolation_type = 0, double lb = 0, double ub = 0);
107  void Akima_interpolation2(Matrix1D<T> &old_function, Matrix1D<T> &old_grid, Matrix1D<T> &new_grid, int extrapolation_type = 0, double lb = 0, double ub = 0);
108  void Interpolate(Matrix1D<T> &old_function, Matrix1D<T> &old_grid, Matrix1D<T> &new_grid);
109  void Spline(Matrix1D<T> &old_function, Matrix1D<T> &old_grid, Matrix1D<T> &new_grid, double lb, double ub, double lin_spline_coef = 0, double max_second_der = 0);
110  void Polint(Matrix1D<T> &old_function, Matrix1D<T> &old_grid, Matrix1D<T> &new_grid);
111  void Ratint(Matrix1D<T> &old_function, Matrix1D<T> &old_grid, Matrix1D<T> &new_grid);
112  void Polilinear(Matrix1D<T> &old_function, Matrix1D<T> &old_grid, Matrix1D<T> &new_grid, double lb, double ub);
113 
114  void Spline2(Matrix1D<T> &old_function, Matrix1D<T> &old_grid, Matrix1D<T> &new_grid, double lb, double ub, double lin_spline_coef = 0, double max_second_der = 0);
115 
116 
117  //T max();
118  //T maxabs();
119 
120  // writeToFileting
121  void writeToFile(string filename);
122  void writeToFile(string filename, Matrix1D<T> &grid_x);
123  void readFromFile(string filename);
124  void readFromFile(string filename, Matrix1D<T> &grid_x);
125 };
126 
127 /**
128  * Matrix 2D class
129  *
130  * Matrixes and operations.
131  */
132 template <typename T> class Matrix2D {
133 private:
134  /// Matrix array (array of links to other arrays). Final likns pointed to the memory addresses of the plane array. Matrix[x][y] can be used.
135  /// Also, all rows saved in the memory one after anouther as a big array. So Matrix[x+x_size*y] can be also used.
137 public:
138  bool initialized; ///< Flag, equal true if initialized
139  int size_x, size_y; ///< size x, size_y
140  string name; ///< name of the Matrix
141 
142  // Constructors and destructors
143  Matrix2D() { initialized = false; };
144  Matrix2D( const Matrix2D<T> &M );
145  Matrix2D( int size_x, int size_y );
146  ~Matrix2D();
147 
148  void AllocateMemory(int size_x, int size_y);
149 
150  // Operators
151 
152  inline T* operator[](int i) { return matrix_array[i]; } ///< Return the i-th pointer to 1d-array. Next [j] can be applied, so we have regular [i][j].
153  inline T* operator[](int i) const { return matrix_array[i]; }
154  inline T& operator()(int x, int y) { return matrix_array[0][x*size_y + y]; } ///< Return the (x,y)-th value of matrix
155  inline T& Value (int x, int y) { return operator()(x, y); } ///< Return the (x,y) value of matrix
156  inline Matrix2D<T>& MatrixArray () { return *this; } ///< Return pointer to the instance of the class.
157 
158  // unary
159  inline const Matrix2D& operator+() const { return *this; }
160  inline const Matrix2D operator-() const { return ((*this)*(-1)); }
161 
162  // The following operators modify the matrix they applied to
163  inline Matrix2D& operator= (const Matrix2D<T> &M);
164  inline Matrix2D& operator= (const T val);
165 
166  //inline Matrix2D& operator*= (const Matrix2D<T> &M); // reserved for something good
167  //inline Matrix2D& operator/= (const Matrix2D<T> &M); // reserved for something good
168  //inline Matrix2D& operator+= (const Matrix2D<T> &M);
169  //inline Matrix2D& operator-= (const Matrix2D<T> &M);
170  //inline Matrix2D& operator*= (const T Val);
171  //inline Matrix2D& operator/= (const T Val);
172  //inline Matrix2D& operator+= (const T Val); ///< Add the Val to each matrix element, stores result in the matrix it's applied to
173  //inline Matrix2D& operator-= (const T Val); ///< Substract the Val from each matrix element, stores result in the matrix it's applied to
174 
175  //inline Matrix2D& times_equal (const Matrix2D<T> &M); ///< Arraywise multiplication (A.*B), stores result in the matrix it's applied to
176  //inline Matrix2D& divide_equal (const Matrix2D<T> &M); ///< Arraywise division (A./B), stores result in the matrix it's applied to
177 
178  // The following operators save the result to a new matrix
179  //inline Matrix2D operator* (const Matrix2D<T> &M) const; // reserved for something good
180  //inline Matrix2D operator/ (const Matrix2D<T> &M) const; // reserved for something good
181  //inline Matrix2D operator+ (const Matrix2D<T> &M) const;
182  //inline Matrix2D operator- (const Matrix2D<T> &M) const;
183  inline Matrix2D operator* (const T Val) const;
184  inline Matrix2D operator/ (const T Val) const;
185  //inline Matrix2D operator+ (const T Val) const; ///< Add the Val to each matrix element, stores result in a new matrix
186  //inline Matrix2D operator- (const T Val) const; ///< Substract the Val from each matrix element, stores result in a new matrix
187 
188  inline Matrix2D times (const Matrix2D<T> &M) const; ///< Arraywise multiplication (A.*B), stores result in a new matrix
189  inline Matrix2D divide (const Matrix2D<T> &M) const; ///< Arraywise division (A./B), stores result in a new matrix
190 
191 
192  // It returns maximum between values from class psd2DSlice and argument (VC::zero_f in that case)
193  Matrix2D max_of(T val);
194 
195  void Interpolate(Matrix2D<T> &old_function, Matrix2D<T> &old_grid_x, Matrix2D<T> &old_grid_y, Matrix2D<T> &new_grid_x, Matrix2D<T> &new_grid_y);
196 
197  // Return corresponding index of 1d array
198  inline int index1d(int x, int y) const;
199 
200  // writeToFileting
201  void writeToFile(string filename);
202  void writeToFile(string filename, Matrix2D<T> &grid_x, Matrix2D<T> &grid_y);
203  void readFromFile(string filename);
204  void readFromFile(string filename, Matrix2D<T> &grid_x, Matrix2D<T> &grid_y);
205 };
206 
207 /**
208  * Matrix 3D class
209  *
210  * Matrixes and operations.
211  */
212 template <typename T> class Matrix3D {
213 private:
214  /// Plane array of values. All rows saved in the memory one after anouther as a big array.
216  /// Matrix array (array of links to other arrays). Final likns pointed to the memory addresses of the plane array. Matrix[x][y] can be used.
218 public:
219  bool initialized; ///< Flag, equal true if initialized
220  int size_x, size_y, size_z; ///< size x, size y, size z
221  string name; ///< name of the Matrix
222 
223  // constructors and destructors
224  /// Default constructor. Do nothing.
225  Matrix3D() { initialized = false; };
226  Matrix3D( const Matrix3D<T> &M );
227  Matrix3D( int size_x, int size_y, int size_z );
228  ~Matrix3D();
229 
230  void AllocateMemory(int size_x, int size_y, int size_z);
231 
232  // Operators
233  inline T** operator[] (int i); ///< Return the i-th pointer to 2d-array. Next [j][k] can be applied, so we have regular [i][j][k].
234  inline T** operator[] (int i) const { return matrix_array[i]; }
235  inline T& operator() (int x, int y, int z); ///< Return the (x,y,z) value of matrix
236  inline T& Value (int x, int y, int z) { return operator()(x, y, z); } ///< Return the (x,y,z) value of matrix
237  inline Matrix3D<T>& MatrixArray () { return *this; } ///< Return pointer to the instance of the class.
238 
239  // The following operators modify the matrix they applied to
240  inline Matrix3D& operator= (const Matrix3D<T> &M);
241 // inline Matrix3D& operator= (const Matrix2D<T> &M); // confusing function
242  inline Matrix3D& operator= (const T Val);
243 
244  // unary
245  inline const Matrix3D& operator+() const { return *this; }
246  inline const Matrix3D operator-() const { return ((*this)*(-1)); }
247 
248  //inline Matrix3D& operator*= (const Matrix3D<T> &M); // reserved for something good
249  //inline Matrix3D& operator/= (const Matrix3D<T> &M); // reserved for something good
250  inline Matrix3D& operator+= (const Matrix3D<T> &M);
251  inline Matrix3D& operator-= (const Matrix3D<T> &M);
252  inline Matrix3D& operator*= (const T Val);
253  inline Matrix3D& operator/= (const T Val);
254  inline Matrix3D& operator+= (const T Val); ///< Add the Val to each matrix element, stores result in the matrix it's applied to
255  inline Matrix3D& operator-= (const T Val); ///< Substract the Val from each matrix element, stores result in the matrix it's applied to
256 
257  inline Matrix3D& times_equal (const Matrix3D<T> &M); ///< Arraywise multiplication (A.*B), stores result in the matrix it's applied to
258  inline Matrix3D& divide_equal (const Matrix3D<T> &M); ///< Arraywise division (A./B), stores result in the matrix it's applied to
259 
260  // The following operators save the result to a new matrix
261  //inline Matrix3D operator* (const Matrix3D<T> &M) const; // reserved for something good
262  //inline Matrix3D operator/ (const Matrix3D<T> &M) const; // reserved for something good
263  inline Matrix3D operator+ (const Matrix3D<T> &M) const;
264  inline Matrix3D operator- (const Matrix3D<T> &M) const;
265  inline Matrix3D operator* (const T Val) const;
266  inline Matrix3D operator/ (const T Val) const;
267  //inline Matrix3D operator+ (const T Val) const; ///< Add the Val to each matrix element, stores result in a new matrix
268  //inline Matrix3D operator- (const T Val) const; ///< Substract the Val from each matrix element, stores result in a new matrix
269 
270  inline Matrix3D times (const Matrix3D<T> &M) const; ///< Arraywise multiplication (A.*B), stores result in a new matrix
271  inline Matrix3D divide (const Matrix3D<T> &M) const; ///< Arraywise division (A./B), stores result in a new matrix
272 
273  // Saving (loading) of a matrix into (from) file
274  void writeToFile(string filename); ///< Save matrix to a file
275  void writeToFile(string filename,
276  Matrix3D<T> &grid_x, Matrix3D<T> &grid_y, Matrix3D<T> &grid_z); ///< Save matrix to a file, including grid
277  void readFromFile(string filename); ///< Load matrix from a file
278  void readFromFile(string filename,
279  Matrix3D<T> &grid_x, Matrix3D<T> &grid_y, Matrix3D<T> &grid_z); ///< Load matrix to a file
280 
281  // Some other stuff
282  string change_ind; ///< Variables useful for tracking of changes (time of change can be stored here)
283 
284  inline int index1d(int x, int y, int z); ///< Returns index of the element (x,y,z) in 1d array
285 
286  T max();
287  T maxabs();
288  Matrix3D<T> abs();
289 
290  // slices - get 2D slice from 3D array
291  Matrix2D<T> xSlice(int p_x) const;
292  Matrix2D<T> ySlice(int p_y) const;
293  Matrix2D<T> zSlice(int p_z) const;
294 
295 };
296 
297 
298 /** Diagonal matrix.
299  * This method of storage for matrices is convenient for diagonal (spread) matrices.
300  * Stored as map (diagonal number, 1d diagonal array)
301  * The USED diagonals of the matrix are stored in 1d arrays.
302  */
303 typedef map <int , Matrix1D<double> > DiagMatrix;
304 
305 /** Model matrix (or related matrices)
306  * It is based on Diagonal matrix and have methods for conversion from 3D or 2D PSD (and related) arrays into 1d array of unknown elements
307  */
309 public:
310 
312  int x_size, y_size, total_size;
313  // flag, if needs to be recalculated
314  string change_ind; ///< Variables useful for changes tracking (store here time when changed)
315 
316  // Constructors
317  CalculationMatrix() { this->initialized = false; }
318  // !!! CalculationMatrix(int x_size, int y_size = 0, int z_size = 0, int n_of_diags = 1);
319  CalculationMatrix(int x_size, int y_size = 1, int z_size = 1, int n_of_diags = 1);
320 
321  // Initialization
322  // !!! void Initialize(int x_size, int y_size = 0, int z_size = 0, int n_of_diags = 1);
323  void Initialize(int x_size, int y_size = 1, int z_size = 1, int n_of_diags = 1);
324 
325  // Returns 1d index for multiple dimension array (2D or 3D)
326  int index1d(int x, int y = 0, int z = 0);
327 
328  // Save to a file
329  void writeToFile(string filename);
330 
331  // Operators
332  Matrix1D<double> operator* (Matrix1D<double> &V) const;
333 
334 };
335 
336 
337 #endif