matrix.h 12.8 KB
Newer Older
Tim Gymnich's avatar
Tim Gymnich committed
1
2
3
4
5
#pragma once

#include "global.h"

#include <vector>
Tim Gymnich's avatar
Tim Gymnich committed
6
#include <type_traits>
Tim Gymnich's avatar
Tim Gymnich committed
7
8
9
10
11

namespace pcpo {

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

Tim Gymnich's avatar
Tim Gymnich committed
17
18
19
20
21
22
23
24
25
26
27
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++) {
Tim Gymnich's avatar
Tim Gymnich committed
28
            std::vector<T> vector(width,value);
Tim Gymnich's avatar
Tim Gymnich committed
29
30
31
32
33
34
35
36
37
            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
38
    Matrix<T>() = default;
Tim Gymnich's avatar
Tim Gymnich committed
39
40
41
42
43
44
45
46
47
    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++) {
Tim Gymnich's avatar
Tim Gymnich committed
48
            std::vector<T> vector(width,0);
Tim Gymnich's avatar
Tim Gymnich committed
49
            vector[i] = 1;
Tim Gymnich's avatar
Tim Gymnich committed
50
            vectors.push_back(vector);
Tim Gymnich's avatar
Tim Gymnich committed
51
52
53
54
55
        }
    };

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

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

Tim Gymnich's avatar
Tim Gymnich committed
72
    Matrix(std::vector<T> const &values, int rows, int columns) {
Tim Gymnich's avatar
Tim Gymnich committed
73
        assert(int(values.size()) == rows * columns);
Tim Gymnich's avatar
Tim Gymnich committed
74
        std::vector<std::vector<T>> result;
Tim Gymnich's avatar
Tim Gymnich committed
75
76
        result.reserve(rows);
        for (int row = 0; row < rows; row++) {
Tim Gymnich's avatar
Tim Gymnich committed
77
            std::vector<T> rowVector;
Tim Gymnich's avatar
Tim Gymnich committed
78
79
80
81
82
83
84
85
86
87
88
            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
89
90
91
92
93
94
    // 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; };
95
96
    /// Returns true when the matrix is empty
    bool empty() const { return getWidth() == 0 && getHeight() == 0; };
Tim Gymnich's avatar
Tim Gymnich committed
97
98
99
100
101
102
103
104

    // 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++) {
105
                result.value(j,i) = value(i,j);
Tim Gymnich's avatar
Tim Gymnich committed
106
107
108
109
110
111
            }
        }
        return result;
    };

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

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

152
    /// Basis of the linear span of the column vectors
Tim Gymnich's avatar
Tim Gymnich committed
153
    static Matrix<T> span(Matrix<T> const& matrix, bool transposed = false) {
Tim Gymnich's avatar
Tim Gymnich committed
154
        std::vector<std::vector<T>> rows;
Tim Gymnich's avatar
Tim Gymnich committed
155
156
157
        // if matrix is already transposed don't do it again
        Matrix<T> t = transposed ? matrix : matrix.transpose();
        Matrix<T> te = t.echelonForm();
158
        int rank = te.getRank();
159
        rows.reserve(rank);
160
161
        for (int row = 0; row < rank; row++) {
            rows.push_back(te.row(row));
Tim Gymnich's avatar
Tim Gymnich committed
162
        }
163
        return Matrix(rows).transpose();
Tim Gymnich's avatar
Tim Gymnich committed
164
    }
Tim Gymnich's avatar
Tim Gymnich committed
165

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

169
    /// Converts the matrix to a 1D Vector by stacking the column vectors
170
    std::vector<T> toVector() const {
Tim Gymnich's avatar
Tim Gymnich committed
171
        std::vector<T> result;
172
        result.reserve(getWidth() * getHeight());
173
174
175
176
        for (int column = 0; column < getWidth(); column++) {
            for (int row = 0; row < getHeight(); row++) {
                result.push_back(value(row,column));
            }
177
        }
178
        return result;
179
180
181
182
183
184
    }

    /// 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 {
185
        assert(rows > 0 && columns > 0);
Tim Gymnich's avatar
patch    
Tim Gymnich committed
186
        // FIXME: Performance
187
        Matrix<T> t = transpose();
Tim Gymnich's avatar
patch    
Tim Gymnich committed
188
        return Matrix(t.vectors.front(), rows, columns).transpose();
Tim Gymnich's avatar
Tim Gymnich committed
189
190
    };

Tim Gymnich's avatar
Tim Gymnich committed
191
192
    std::vector<Matrix<T>> reshapeColumns(int height, int width) const {
        std::vector<Matrix<T>> result;
193
        for (int c = 0; c < getWidth(); c++) {
Tim Gymnich's avatar
Tim Gymnich committed
194
            result.push_back(Matrix(column(c), height, width).transpose());
195
        }
Tim Gymnich's avatar
Tim Gymnich committed
196
        return result;
197
    }
198
199
200
201
202

    /// 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
203
        assert(row < getHeight() && column < getWidth());
204
205
206
207
208
209
210
        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
211
        assert(row < getHeight() && column < getWidth());
212
213
214
215
        return vectors[row][column];
    };


Tim Gymnich's avatar
Tim Gymnich committed
216
    /// 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
217
    /// @param i Index of the row to return.
Tim Gymnich's avatar
Tim Gymnich committed
218
    std::vector<T>& row(int i) {
219
220
221
        assert(i < getHeight());
        return vectors[i];
    };
Tim Gymnich's avatar
Tim Gymnich committed
222

Tim Gymnich's avatar
Tim Gymnich committed
223
    std::vector<T> row(int i) const {
224
225
226
227
        assert(i < getHeight());
        return vectors[i];
    };

Tim Gymnich's avatar
Tim Gymnich committed
228
    /// Returns the column at index i. The returned column cannot be modified.
Tim Gymnich's avatar
Tim Gymnich committed
229
    /// @param i Index of the column to return
Tim Gymnich's avatar
Tim Gymnich committed
230
    std::vector<T> column(int i) const {
231
        assert(i < getWidth());
Tim Gymnich's avatar
Tim Gymnich committed
232
        std::vector<T> row;
Tim Gymnich's avatar
Tim Gymnich committed
233
        row.reserve(width);
Tim Gymnich's avatar
Tim Gymnich committed
234
        for (std::vector<T> const& x : vectors) {
Tim Gymnich's avatar
Tim Gymnich committed
235
236
237
238
239
            row.push_back(x[i]);
        }
        return row;
    }

Tim Gymnich's avatar
Tim Gymnich committed
240
    void setColumn(std::vector<T> const& vector, int column) {
Tim Gymnich's avatar
Tim Gymnich committed
241
242
243
244
245
246
        assert(int(vector.size()) == height);
        for (int row = 0; row < height; row++) {
            value(row,column) = vector[row];
        }
    }

Tim Gymnich's avatar
Tim Gymnich committed
247
248
    // MARK: - Operators

249
    T& operator()(int i, int j) { return value(i,j); };
Tim Gymnich's avatar
Tim Gymnich committed
250
251
252
253
254
255
256

    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++) {
257
                    result.value(i,j) += value(i,k) * rhs.value(k,j);
Tim Gymnich's avatar
Tim Gymnich committed
258
259
260
261
262
263
264
265
266
267
268
269
                }
            }
        }
        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++) {
270
                value(i,j) *= scalar;
Tim Gymnich's avatar
Tim Gymnich committed
271
272
273
274
275
276
277
278
279
280
            }
        }
        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++) {
281
                value(i,j) += rhs.value(i,j);
Tim Gymnich's avatar
Tim Gymnich committed
282
283
284
285
286
287
288
289
            }
        }
        return *this;
    };

    Matrix<T>& operator +=(T scalar) {
        for (int i = 0; i < height; i++) {
            for (int j = 0; j < width; j++) {
290
                value(i,j) += scalar;
Tim Gymnich's avatar
Tim Gymnich committed
291
292
293
294
295
296
297
298
299
            }
        }
        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++) {
300
                value(i,j) -= rhs.value(i,j);
Tim Gymnich's avatar
Tim Gymnich committed
301
302
303
304
305
306
307
308
            }
        }
        return *this;
    };

    Matrix<T>& operator -=(int scalar) {
        for (int i = 0; i < height; i++) {
            for (int j = 0; j < width; j++) {
309
                value(i,j) -= scalar;
Tim Gymnich's avatar
Tim Gymnich committed
310
311
312
313
314
315
316
317
318
            }
        }
        return *this;
    };

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

319
320
321
322
    void print() const {
        dbgs(4) << *this;
    }

Tim Gymnich's avatar
Tim Gymnich committed
323
324
325
326
327
328
329
330
331
332
333
334
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
335
    void divide_row(int row, T quotient) {
Tim Gymnich's avatar
Tim Gymnich committed
336
        for (int column = 0; column < width; column++) {
337
            value(row,column) /= quotient;
Tim Gymnich's avatar
Tim Gymnich committed
338
339
340
341
342
343
344
        }
    };

    /// 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
345
    void add_multiple_row(int a, int b, T factor) {
Tim Gymnich's avatar
Tim Gymnich committed
346
        for (int column = 0; column < width; column++) {
347
            value(a,column) += value(b,column) * factor;
Tim Gymnich's avatar
Tim Gymnich committed
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
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
        }
    };

    // 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
392
template <typename T>
Tim Gymnich's avatar
Tim Gymnich committed
393
llvm::raw_ostream& operator<<(llvm::raw_ostream& os, Matrix<T> const& matrix) {
Tim Gymnich's avatar
Tim Gymnich committed
394
    for (int row = 0; row < matrix.getHeight(); row++) {
Tim Gymnich's avatar
Tim Gymnich committed
395
        os << "[ ";
Tim Gymnich's avatar
Tim Gymnich committed
396
        for (int column = 0; column < matrix.getWidth(); column++) {
Tim Gymnich's avatar
Tim Gymnich committed
397
398
            if constexpr (std::is_floating_point_v<T>) {
                if (column == matrix.getWidth() - 1) {
Tim Gymnich's avatar
Tim Gymnich committed
399
                    os << llvm::format("%g", matrix.value(row,column));
Tim Gymnich's avatar
Tim Gymnich committed
400
                } else {
Tim Gymnich's avatar
Tim Gymnich committed
401
                    os << llvm::format("%-6g", matrix.value(row,column));
Tim Gymnich's avatar
Tim Gymnich committed
402
                }
Tim Gymnich's avatar
Tim Gymnich committed
403
            } else {
Tim Gymnich's avatar
Tim Gymnich committed
404
405
406
407
408
                if (column == matrix.getWidth() - 1) {
                    os << llvm::format("%d", matrix.value(row,column));
                } else {
                    os << llvm::format("%-6d", matrix.value(row,column));
                }
Tim Gymnich's avatar
Tim Gymnich committed
409
            }
Tim Gymnich's avatar
Tim Gymnich committed
410
        }
Tim Gymnich's avatar
Tim Gymnich committed
411
        os << " ]\n";
Tim Gymnich's avatar
Tim Gymnich committed
412
413
414
415
    }
    return os;
};

Tim Gymnich's avatar
Tim Gymnich committed
416
417
}