matrix.h 11.8 KB
Newer Older
Tim Gymnich's avatar
Tim Gymnich committed
1
2
3
4
5
6
7
8
9
10
11
12
#pragma once

#include "global.h"

#include <vector>

using namespace std;

namespace pcpo {

/// Matrix. Row and column are indexed beginning at 0
template<typename T> class Matrix {
Tim Gymnich's avatar
Tim Gymnich committed
13
14
15
16
17
protected:
    vector<vector<T>> vectors;
    int width;
    int height;

Tim Gymnich's avatar
Tim Gymnich committed
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
public:
    //MARK: - Initializers
    /// Creates a matrix with dimensions  height x width and initalizes its values to `value`
    /// @param width Width of the matrix
    /// @param height Height of the matrix
    /// @param value Initial value for each element
    Matrix<T>(int height, int width, T value) {
        this->width = width;
        this->height = height;
        this->vectors.reserve(width);
        for (int i = 0; i < height; i++) {
            vector<T> vector(width,value);
            vectors.push_back(vector);
        }
    };

    /// Creates a matrix with dimensions  height x width and initalizes its values to 0
    /// @param width Width of the matrix
    /// @param height Height of the matrix
    Matrix<T>(int height, int width) : Matrix<T>(height, width, 0) {};

Tim Gymnich's avatar
Tim Gymnich committed
39
    Matrix<T>() = default;
Tim Gymnich's avatar
Tim Gymnich committed
40
41
42
43
44
45
46
47
48
49
50
    Matrix(Matrix const& matrix) = default;

    /// Creates an identity matrix with dimension eye x eye
    /// @param eye dimension of the matrix
    Matrix(int eye){
        this->width = eye;
        this->height = eye;
        this->vectors.reserve(width);
        for (int i = 0; i < height; i++) {
            vector<T> vector(width,0);
            vector[i] = 1;
Tim Gymnich's avatar
Tim Gymnich committed
51
            vectors.push_back(vector);
Tim Gymnich's avatar
Tim Gymnich committed
52
53
54
55
56
57
        }
    };

    /// Creates a matrix from a 2D vector
    /// @param vectors 2D vector containing columns with rows
    Matrix(vector<vector<T>> const &vectors) {
58
59
        assert(all_of(vectors.begin(), vectors.end(), [&vectors](vector<T> vec){ return vec.size() == vectors.front().size(); }));
        this->width = vectors.empty() ? 0 : vectors.front().size();
60
        this->height = vectors.size();
Tim Gymnich's avatar
Tim Gymnich committed
61
62
63
64
65
66
67
        this->vectors = vectors;
    };

    /// Creates a vector from a std::vector
    /// @param vector the vector
    Matrix(vector<T> const &vector) {
        std::vector<std::vector<T>> vectors = {vector};
68
69
        this->width = vector.size();
        this->height = vector.empty() ? 0 : 1;
Tim Gymnich's avatar
Tim Gymnich committed
70
71
72
        this->vectors = vectors;
    };

Tim Gymnich's avatar
Tim Gymnich committed
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
    Matrix(vector<T> const &values, int rows, int columns) {
        assert(int(values.size()) == rows * columns);
        vector<vector<T>> result;
        result.reserve(rows);
        for (int row = 0; row < rows; row++) {
            vector<T> rowVector;
            rowVector.reserve(columns);
            for (int column = 0; column < columns; column++) {
                rowVector.push_back(values[row * rows + column]);
            }
            result.push_back(rowVector);
        }
        this->vectors = result;
        this->width = columns;
        this->height = rows;
    };

Tim Gymnich's avatar
Tim Gymnich committed
90
91
92
93
94
95
    // MARK: - Properties

    /// The height of the matrix (number of rows)
    int getHeight() const { return height; };
    /// The width of the matrix (number of columns)
    int getWidth() const { return width; };
96
97
    /// Returns true when the matrix is empty
    bool empty() const { return getWidth() == 0 && getHeight() == 0; };
Tim Gymnich's avatar
Tim Gymnich committed
98
99
100
101
102
103
104
105

    // MARK: - Matrix operations

    /// Transposes the matrix
    Matrix transpose() const {
        Matrix result = Matrix(width, height);
        for (int i = 0; i < height; i++) {
            for (int j = 0; j < width; j++) {
106
                result.value(j,i) = value(i,j);
Tim Gymnich's avatar
Tim Gymnich committed
107
108
109
110
111
112
113
114
115
116
117
118
            }
        }
        return result;
    };

    /// Transforms the matrix to reduced row echelon form
    Matrix echelon() const {
        Matrix result = Matrix(*this);
        int pivot = 0;
        for (int row = 0; row < height; row++) {
            if (pivot >= width) { return result; }
            int i = row;
119
            while (result.value(i,pivot) == 0) {
Tim Gymnich's avatar
Tim Gymnich committed
120
121
122
123
124
125
                if (++i >= height) {
                    i = row;
                    if (++pivot >= width) { return result; }
                }
            }
            result.swap_rows(i, row);
126
            result.divide_row(row, result.value(row,pivot));
Tim Gymnich's avatar
Tim Gymnich committed
127
128
            for (i = 0; i < height; i++) {
                if (i != row) {
129
                    result.add_multiple_row(i, row, -result.value(i,pivot));
Tim Gymnich's avatar
Tim Gymnich committed
130
131
132
133
134
135
                }
            }
        }
        return result;
    };

Tim Gymnich's avatar
Tim Gymnich committed
136
137
138
139
140
141
    /// The rank of the matrix
    int getRank() const {
        Matrix e = echelon();
        int rank = 0;
        for (int row = 0; row < height; row++) {
            for (int column = 0; column < width; column++) {
142
                if ((e.value(row,column) == 0) && (column == width - 1)) {
Tim Gymnich's avatar
Tim Gymnich committed
143
                    return rank;
144
                } else if (e.value(row,column) != 0) {
Tim Gymnich's avatar
Tim Gymnich committed
145
146
147
148
149
150
151
                    break;
                }
            }
            rank++;
        }
        return rank;
    }
Tim Gymnich's avatar
Tim Gymnich committed
152

153
    /// Basis of the linear span of the column vectors
Tim Gymnich's avatar
Tim Gymnich committed
154
    static Matrix<T> span(Matrix<T> const& matrix) {
Tim Gymnich's avatar
Tim Gymnich committed
155
        vector<vector<T>> columns;
Tim Gymnich's avatar
Tim Gymnich committed
156
        int rank = matrix.getRank();
157
        for (int col = 0; col < rank; col++) {
Tim Gymnich's avatar
Tim Gymnich committed
158
            columns.push_back(matrix.column(col));
Tim Gymnich's avatar
Tim Gymnich committed
159
160
161
        }
        return Matrix(columns).transpose();
    }
Tim Gymnich's avatar
Tim Gymnich committed
162

163
    /// Computes the null space for the column vectors
Tim Gymnich's avatar
Tim Gymnich committed
164
    static Matrix<T> null(Matrix<T> const& matrix);
165

166
    /// Converts the matrix to a 1D Vector by stacking the column vectors
167
    std::vector<T> toVector() const {
168
169
170
171
172
        vector<T> result;
        result.reserve(getWidth() * getHeight());
        for (vector<T> vector: vectors) {
            result.insert(result.end(), vector.begin(), vector.end());
        }
173
        return result;
174
175
176
177
178
179
    }

    /// Converts a 1D Vector to a Matrix with given dimensions
    /// @param rows number of rows
    /// @param columns number of columns
    Matrix<T> reshape(int rows, int columns) const {
180
181
182
        assert(rows > 0 && columns > 0);
        Matrix<T> t = transpose();
        return Matrix(t.vectors.front(), rows, columns);
Tim Gymnich's avatar
Tim Gymnich committed
183
184
    };

185
186
    vector<Matrix<T>> reshapeColumns(int height, int width) const {
        vector<Matrix<T>> result;
187
        for (int c = 0; c < getWidth(); c++) {
Tim Gymnich's avatar
Tim Gymnich committed
188
            result.push_back(Matrix(column(c), height, width));
189
        }
Tim Gymnich's avatar
Tim Gymnich committed
190
        return result;
191
    }
192
193
194
195
196

    /// Returns the value at row i and column j
    /// @param row
    /// @param column
    T& value(int row, int column) {
Tim Gymnich's avatar
Tim Gymnich committed
197
        assert(row < getHeight() && column < getWidth());
198
199
200
201
202
203
204
        return vectors[row][column];
    };

    /// Returns the value at row i and column j
    /// @param row
    /// @param column
    T const& value(int row, int column) const {
Tim Gymnich's avatar
Tim Gymnich committed
205
        assert(row < getHeight() && column < getWidth());
206
207
208
209
        return vectors[row][column];
    };


Tim Gymnich's avatar
Tim Gymnich committed
210
    /// Returns a vector with the elements of the row at index i. The returned row can be modified.
Tim Gymnich's avatar
Tim Gymnich committed
211
    /// @param i Index of the row to return.
212
213
214
215
    vector<T>& row(int i) {
        assert(i < getHeight());
        return vectors[i];
    };
Tim Gymnich's avatar
Tim Gymnich committed
216

Tim Gymnich's avatar
Tim Gymnich committed
217
    /// Returns the column at index i. The returned column cannot be modified.
Tim Gymnich's avatar
Tim Gymnich committed
218
219
    /// @param i Index of the column to return
    vector<T> column(int i) const {
220
        assert(i < getWidth());
Tim Gymnich's avatar
Tim Gymnich committed
221
222
223
224
225
226
227
228
        vector<T> row;
        row.reserve(width);
        for (vector<T> const& x : vectors) {
            row.push_back(x[i]);
        }
        return row;
    }

Tim Gymnich's avatar
Tim Gymnich committed
229
230
231
232
233
234
235
    void setColumn(vector<T> const& vector, int column) {
        assert(int(vector.size()) == height);
        for (int row = 0; row < height; row++) {
            value(row,column) = vector[row];
        }
    }

Tim Gymnich's avatar
Tim Gymnich committed
236
237
    // MARK: - Operators

238
    T& operator()(int i, int j) { return value(i,j); };
Tim Gymnich's avatar
Tim Gymnich committed
239
240
241
242
243
244
245

    Matrix<T>& operator *=(Matrix<T> const& rhs) {
        assert(width == rhs.height);
        Matrix result = Matrix(height,rhs.width);
        for (int i = 0; i < height; i++) {
            for (int k = 0; k < width; k++) {
                for (int j = 0; j < rhs.width; j++) {
246
                    result.value(i,j) += value(i,k) * rhs.value(k,j);
Tim Gymnich's avatar
Tim Gymnich committed
247
248
249
250
251
252
253
254
255
256
257
258
                }
            }
        }
        this->vectors = result.vectors;
        this->width = result.width;
        this->height = result.height;
        return *this;
    };

    Matrix<T>& operator *=(T scalar) {
        for (int i = 0; i < height; i++) {
            for (int j = 0; j < width; j++) {
259
                value(i,j) *= scalar;
Tim Gymnich's avatar
Tim Gymnich committed
260
261
262
263
264
265
266
267
268
269
            }
        }
        return *this;
    };


    Matrix<T>& operator +=(Matrix<T> const& rhs) {
        assert(rhs.width == width && rhs.height == height);
        for (int i=0;i<height;i++) {
            for (int j = 0; j < width; j++) {
270
                value(i,j) += rhs.value(i,j);
Tim Gymnich's avatar
Tim Gymnich committed
271
272
273
274
275
276
277
278
            }
        }
        return *this;
    };

    Matrix<T>& operator +=(T scalar) {
        for (int i = 0; i < height; i++) {
            for (int j = 0; j < width; j++) {
279
                value(i,j) += scalar;
Tim Gymnich's avatar
Tim Gymnich committed
280
281
282
283
284
285
286
287
288
            }
        }
        return *this;
    };

    Matrix<T>& operator -=(Matrix<T> const& rhs) {
        assert(rhs.width == width && rhs.height == height);
        for (int i = 0; i < height; i++) {
            for (int j= 0; j < width; j++) {
289
                value(i,j) -= rhs.value(i,j);
Tim Gymnich's avatar
Tim Gymnich committed
290
291
292
293
294
295
296
297
            }
        }
        return *this;
    };

    Matrix<T>& operator -=(int scalar) {
        for (int i = 0; i < height; i++) {
            for (int j = 0; j < width; j++) {
298
                value(i,j) -= scalar;
Tim Gymnich's avatar
Tim Gymnich committed
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
            }
        }
        return *this;
    };

    bool operator==(Matrix<T> const& rhs) const {
        return rhs.vectors == rhs.vectors && width == rhs.width && height == rhs.height;
    };

protected:

    // MARK: - Echelon helpers

    /// Swaps two rows
    /// @param a index of the first row
    /// @param b index of the second row
    void swap_rows(int a, int b) { vectors[a].swap(vectors[b]); };

    /// Divides a row by a constant
    /// @param row  index of the row to divide
    /// @param quotient quotient to divide the row by
    void divide_row(int row, int quotient) {
        for (int column = 0; column < width; column++) {
322
            value(row,column) /= quotient;
Tim Gymnich's avatar
Tim Gymnich committed
323
324
325
326
327
328
329
330
331
        }
    };

    /// Adds a multiple of row b to a row a
    /// @param a Row to add a multiple of b to
    /// @param b Row to be added to row a
    /// @param factor Factor to multiply row b with when adding it to row a
    void add_multiple_row(int a, int b, int factor) {
        for (int column = 0; column < width; column++) {
332
            value(a,column) += value(b,column) * factor;
Tim Gymnich's avatar
Tim Gymnich committed
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
        }
    };

    // MARK: - Utils

    /// Greates common divisor.
    /// @param lhs
    /// @param rhs
    static int ggT(int lhs, int rhs) {
        int h;
        if (lhs == 0) { return abs(rhs); }
        if (rhs == 0) { return abs(lhs); }

        do {
            h = lhs % rhs;
            lhs = rhs;
            rhs = h;
        } while (rhs != 0);

        return abs(lhs);
    };

    /// Least common multiple.
    /// @param lhs
    /// @param rhs
    static int kgV(int lhs, int rhs) { return (lhs * rhs) / ggT(lhs, rhs); };

    
};


template <typename T>
inline Matrix<T> operator*(Matrix<T> lhs, Matrix<T> const& rhs) { return lhs *= rhs; };
template <typename T>
inline Matrix<T> operator*(Matrix<T> lhs, T scalar) { return lhs *= scalar; };
template <typename T>
inline Matrix<T> operator+(Matrix<T> lhs, Matrix<T> const& rhs) { return lhs += rhs; };
template <typename T>
inline Matrix<T> operator+(Matrix<T> lhs, T scalar) { return lhs += scalar; };
template <typename T>
inline Matrix<T> operator-(Matrix<T> lhs, Matrix<T> const& rhs) { return  lhs -= rhs; };
template <typename T>
inline Matrix<T> operator-(Matrix<T> lhs, T scalar) { return lhs -= scalar; };

Tim Gymnich's avatar
Tim Gymnich committed
377
template <typename T>
Tim Gymnich's avatar
Tim Gymnich committed
378
llvm::raw_ostream& operator<<(llvm::raw_ostream& os, Matrix<T> const& matrix) {
Tim Gymnich's avatar
Tim Gymnich committed
379
    for (int row = 0; row < matrix.getHeight(); row++) {
Tim Gymnich's avatar
Tim Gymnich committed
380
        os << "[ ";
Tim Gymnich's avatar
Tim Gymnich committed
381
        for (int column = 0; column < matrix.getWidth(); column++) {
Tim Gymnich's avatar
Tim Gymnich committed
382
            if (column == matrix.getWidth() - 1) {
Tim Gymnich's avatar
Tim Gymnich committed
383
                os << llvm::format("%d", matrix.value(row,column));
Tim Gymnich's avatar
Tim Gymnich committed
384
            } else {
Tim Gymnich's avatar
Tim Gymnich committed
385
                os << llvm::format("%-6d", matrix.value(row,column));
Tim Gymnich's avatar
Tim Gymnich committed
386
            }
Tim Gymnich's avatar
Tim Gymnich committed
387
        }
Tim Gymnich's avatar
Tim Gymnich committed
388
        os << " ]\n";
Tim Gymnich's avatar
Tim Gymnich committed
389
390
391
392
    }
    return os;
};

Tim Gymnich's avatar
Tim Gymnich committed
393
394
}