a19e3c80
David Mayerich
matrix edits for ...
|
1
2
|
#ifndef STIM_MATRIX_H
#define STIM_MATRIX_H
|
57729e5b
David Mayerich
changed directory...
|
3
4
5
6
|
//#include "rts/vector.h"
#include <string.h>
#include <iostream>
|
91589064
David Mayerich
added a global ma...
|
7
|
#include <fstream>
|
cce7daf9
Pavel Govyadinov
added glObj.h to ...
|
8
|
#include <stim/math/vector.h>
|
308a743c
David Mayerich
fixed class compa...
|
9
|
#include <stim/math/vec3.h>
|
1c92dffe
David Mayerich
simplified CUDA i...
|
10
|
//#include <stim/cuda/cudatools/callable.h>
|
57729e5b
David Mayerich
changed directory...
|
11
|
|
8a86bd56
David Mayerich
changed rts names...
|
12
|
namespace stim{
|
57729e5b
David Mayerich
changed directory...
|
13
|
|
7d01bb90
David Mayerich
matlab matrix sav...
|
14
15
16
17
18
19
20
21
22
23
|
enum mat4Format {
mat4_float64,
mat4_float32,
mat4_int32,
mat4_int16,
mat4_uint16,
mat4_uint8,
mat4_float //floating point type, determined automatically
};
|
a11c7efe
David Mayerich
fixed stim::filen...
|
24
|
static size_t mat4Format_size(mat4Format f){
|
7d01bb90
David Mayerich
matlab matrix sav...
|
25
26
27
28
29
30
31
32
33
34
|
switch(f){
case mat4_float64: return 8;
case mat4_float32:
case mat4_int32: return 4;
case mat4_int16:
case mat4_uint16: return 2;
case mat4_uint8: return 1;
default: return 0;
}
}
|
91589064
David Mayerich
added a global ma...
|
35
|
|
a11c7efe
David Mayerich
fixed stim::filen...
|
36
|
static void save_mat4(char* data, std::string filename, std::string varname, size_t sx, size_t sy, mat4Format format){
|
91589064
David Mayerich
added a global ma...
|
37
38
39
40
41
42
43
44
45
|
//save the matrix file here (use the mat4 function above)
//data format: https://maxwell.ict.griffith.edu.au/spl/matlab-page/matfile_format.pdf (page 32)
int MOPT = 0; //initialize the MOPT type value to zero
int m = 0; //little endian
int o = 0; //reserved, always 0
int p = format;
int t = 0;
MOPT = m * 1000 + o * 100 + p * 10 + t; //calculate the type value
|
a11c7efe
David Mayerich
fixed stim::filen...
|
46
47
|
int mrows = (int)sx;
int ncols = (int)sy;
|
91589064
David Mayerich
added a global ma...
|
48
49
|
int imagf = 0; //assume real (for now)
varname.push_back('\0'); //add a null to the string
|
a11c7efe
David Mayerich
fixed stim::filen...
|
50
|
int namlen = (int)varname.size(); //calculate the name size
|
91589064
David Mayerich
added a global ma...
|
51
52
53
54
55
56
57
58
59
60
61
|
size_t bytes = sx * sy * mat4Format_size(format);
std::ofstream outfile(filename, std::ios::binary);
outfile.write((char*)&MOPT, 4);
outfile.write((char*)&mrows, 4);
outfile.write((char*)&ncols, 4);
outfile.write((char*)&imagf, 4);
outfile.write((char*)&namlen, 4);
outfile.write((char*)&varname[0], namlen);
outfile.write((char*)data, bytes); //write the matrix data
outfile.close();
|
91589064
David Mayerich
added a global ma...
|
62
63
|
}
|
6c4afcac
David Mayerich
introduced a gene...
|
64
65
|
template <class T>
class matrix {
|
57729e5b
David Mayerich
changed directory...
|
66
|
//the matrix will be stored in column-major order (compatible with OpenGL)
|
6c4afcac
David Mayerich
introduced a gene...
|
67
68
69
70
|
T* M; //pointer to the matrix data
size_t R; //number of rows
size_t C; //number of colums
|
7d01bb90
David Mayerich
matlab matrix sav...
|
71
72
73
|
size_t bytes() {
return R * C * sizeof(T); //return the number of bytes of matrix data
}
|
65c1aac0
David Mayerich
fixed memory leak...
|
74
|
/*void init(size_t rows, size_t cols){
|
6c4afcac
David Mayerich
introduced a gene...
|
75
76
|
R = rows;
C = cols;
|
65c1aac0
David Mayerich
fixed memory leak...
|
77
78
79
80
81
82
83
84
85
86
87
|
if (R == 0 || C == 0) M = NULL;
else
M = (T*)malloc(R * C * sizeof(T)); //allocate space for the matrix
}*/
T get(const size_t row, const size_t col) const {
if (row >= R || col >= C) {
std::cout << "ERROR: row or column out of range." << std::endl;
exit(1);
}
return M[col * R + row];
|
6c4afcac
David Mayerich
introduced a gene...
|
88
|
}
|
57729e5b
David Mayerich
changed directory...
|
89
|
|
6c4afcac
David Mayerich
introduced a gene...
|
90
|
T& at(size_t row, size_t col){
|
65c1aac0
David Mayerich
fixed memory leak...
|
91
92
93
94
|
if (row >= R || col >= C) {
std::cout << "ERROR: row or column out of range." << std::endl;
exit(1);
}
|
6c4afcac
David Mayerich
introduced a gene...
|
95
|
return M[col * R + row];
|
57729e5b
David Mayerich
changed directory...
|
96
97
|
}
|
6c4afcac
David Mayerich
introduced a gene...
|
98
|
public:
|
65c1aac0
David Mayerich
fixed memory leak...
|
99
100
101
102
103
104
|
matrix() {
R = 0;
C = 0;
M = NULL;
}
|
6c4afcac
David Mayerich
introduced a gene...
|
105
|
matrix(size_t rows, size_t cols) {
|
65c1aac0
David Mayerich
fixed memory leak...
|
106
107
|
R = rows;
C = cols;
|
ef1f6218
David Mayerich
fixed some memory...
|
108
109
110
|
M = NULL;
if (R * C > 0)
M = (T*) malloc(R * C * sizeof(T));
|
b2d10ff1
Pavel Govyadinov
final edits, expa...
|
111
112
|
}
|
ef1f6218
David Mayerich
fixed some memory...
|
113
|
matrix(size_t rows, size_t cols, const T* data) {
|
65c1aac0
David Mayerich
fixed memory leak...
|
114
115
|
R = rows;
C = cols;
|
ef1f6218
David Mayerich
fixed some memory...
|
116
117
118
|
M = NULL;
if (R * C > 0)
M = (T*)malloc(R * C * sizeof(T));
|
dd5aab2f
David Mayerich
fixed errors in s...
|
119
|
memcpy(M, data, R * C * sizeof(T));
|
8e56a0a7
Pavel Govyadinov
Added the propose...
|
120
|
}
|
b2d10ff1
Pavel Govyadinov
final edits, expa...
|
121
|
|
6c4afcac
David Mayerich
introduced a gene...
|
122
|
matrix(const matrix<T>& cpy){
|
ef1f6218
David Mayerich
fixed some memory...
|
123
124
125
|
M = NULL;
if (cpy.R * cpy.C > 0)
M = (T*)malloc(cpy.R * cpy.C * sizeof(T));
|
65c1aac0
David Mayerich
fixed memory leak...
|
126
127
128
129
|
memcpy(M, cpy.M, cpy.R * cpy.C * sizeof(T));
R = cpy.R;
C = cpy.C;
|
6c4afcac
David Mayerich
introduced a gene...
|
130
|
}
|
19a3973b
David Mayerich
updates to matrix...
|
131
|
|
6c4afcac
David Mayerich
introduced a gene...
|
132
|
~matrix() {
|
6c4afcac
David Mayerich
introduced a gene...
|
133
|
if(M) free(M);
|
65c1aac0
David Mayerich
fixed memory leak...
|
134
135
|
M = NULL;
R = C = 0;
|
19a3973b
David Mayerich
updates to matrix...
|
136
137
|
}
|
65c1aac0
David Mayerich
fixed memory leak...
|
138
|
size_t rows() const {
|
6c4afcac
David Mayerich
introduced a gene...
|
139
|
return R;
|
57729e5b
David Mayerich
changed directory...
|
140
141
|
}
|
65c1aac0
David Mayerich
fixed memory leak...
|
142
|
size_t cols() const {
|
6c4afcac
David Mayerich
introduced a gene...
|
143
144
|
return C;
}
|
57729e5b
David Mayerich
changed directory...
|
145
|
|
1c92dffe
David Mayerich
simplified CUDA i...
|
146
|
T& operator()(size_t row, size_t col) {
|
6c4afcac
David Mayerich
introduced a gene...
|
147
148
149
|
return at(row, col);
}
|
65c1aac0
David Mayerich
fixed memory leak...
|
150
151
|
matrix<T>& operator=(const T rhs) {
//init(R, C);
|
6c4afcac
David Mayerich
introduced a gene...
|
152
153
154
155
156
157
158
|
size_t N = R * C;
for(size_t n=0; n<N; n++)
M[n] = rhs;
return *this;
}
|
65c1aac0
David Mayerich
fixed memory leak...
|
159
160
161
162
163
164
165
166
167
168
|
matrix<T>& operator=(const matrix<T>& rhs){
if (this != &rhs) { //if the matrix isn't self-assigned
T* new_matrix = new T[rhs.R * rhs.C]; //allocate new resources
memcpy(new_matrix, rhs.M, rhs.R * rhs.C * sizeof(T)); //copy the matrix
delete[] M; //delete the previous array
M = new_matrix;
R = rhs.R;
C = rhs.C;
}
|
57729e5b
David Mayerich
changed directory...
|
169
170
|
return *this;
}
|
42bb075c
Jiaming Guo
last changes for ...
|
171
|
|
6c4afcac
David Mayerich
introduced a gene...
|
172
|
//element-wise operations
|
65c1aac0
David Mayerich
fixed memory leak...
|
173
|
matrix<T> operator+(const T rhs) const {
|
6c4afcac
David Mayerich
introduced a gene...
|
174
175
176
|
matrix<T> result(R, C); //create a result matrix
size_t N = R * C;
|
42bb075c
Jiaming Guo
last changes for ...
|
177
|
for(int i=0; i<N; i++)
|
6c4afcac
David Mayerich
introduced a gene...
|
178
179
180
|
result.M[i] = M[i] + rhs; //calculate the operation and assign to result
return result;
|
42bb075c
Jiaming Guo
last changes for ...
|
181
|
}
|
57729e5b
David Mayerich
changed directory...
|
182
|
|
65c1aac0
David Mayerich
fixed memory leak...
|
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
|
matrix<T> operator+(const matrix<T> rhs) const {
if (R != rhs.R || C != rhs.C) {
std::cout << "ERROR: addition is only defined for matrices that are the same size." << std::endl;
exit(1);
}
matrix<T> result(R, C); //create a result matrix
size_t N = R * C;
for (int i = 0; i < N; i++)
result.M[i] = M[i] + rhs.M[i]; //calculate the operation and assign to result
return result;
}
matrix<T> operator-(const T rhs) const {
|
6c4afcac
David Mayerich
introduced a gene...
|
198
199
200
|
return operator+(-rhs); //add the negative of rhs
}
|
65c1aac0
David Mayerich
fixed memory leak...
|
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
|
matrix<T> operator-(const matrix<T> rhs) const {
return operator+(-rhs);
}
matrix<T> operator-() const {
matrix<T> result(R, C); //create a result matrix
size_t N = R * C;
for (int i = 0; i < N; i++)
result.M[i] = -M[i]; //calculate the operation and assign to result
return result;
}
matrix<T> operator*(const T rhs) const {
|
6c4afcac
David Mayerich
introduced a gene...
|
216
217
218
219
220
221
222
223
224
|
matrix<T> result(R, C); //create a result matrix
size_t N = R * C;
for(int i=0; i<N; i++)
result.M[i] = M[i] * rhs; //calculate the operation and assign to result
return result;
}
|
65c1aac0
David Mayerich
fixed memory leak...
|
225
|
matrix<T> operator/(const T rhs) const {
|
6c4afcac
David Mayerich
introduced a gene...
|
226
227
228
229
230
|
matrix<T> result(R, C); //create a result matrix
size_t N = R * C;
for(int i=0; i<N; i++)
result.M[i] = M[i] / rhs; //calculate the operation and assign to result
|
487a9b49
David Mayerich
added the ability...
|
231
|
|
6c4afcac
David Mayerich
introduced a gene...
|
232
233
|
return result;
}
|
57729e5b
David Mayerich
changed directory...
|
234
|
|
6c4afcac
David Mayerich
introduced a gene...
|
235
|
//matrix multiplication
|
65c1aac0
David Mayerich
fixed memory leak...
|
236
|
matrix<T> operator*(const matrix<T> rhs) const {
|
6c4afcac
David Mayerich
introduced a gene...
|
237
238
239
240
241
|
if(C != rhs.R){
std::cout<<"ERROR: matrix multiplication is undefined for matrices of size ";
std::cout<<"[ "<<R<<" x "<<C<<" ] and [ "<<rhs.R<<" x "<<rhs.C<<"]"<<std::endl;
exit(1);
}
|
57729e5b
David Mayerich
changed directory...
|
242
|
|
6c4afcac
David Mayerich
introduced a gene...
|
243
244
|
matrix<T> result(R, rhs.C); //create the output matrix
T inner; //stores the running inner product
|
1c92dffe
David Mayerich
simplified CUDA i...
|
245
246
247
|
size_t c, r, i;
for(c = 0; c < rhs.C; c++){
for(r = 0; r < R; r++){
|
6c4afcac
David Mayerich
introduced a gene...
|
248
|
inner = (T)0;
|
1c92dffe
David Mayerich
simplified CUDA i...
|
249
|
for(i = 0; i < C; i++){
|
65c1aac0
David Mayerich
fixed memory leak...
|
250
|
inner += get(r, i) * rhs.get(i, c);
|
6c4afcac
David Mayerich
introduced a gene...
|
251
252
253
254
|
}
result.M[c * R + r] = inner;
}
}
|
57729e5b
David Mayerich
changed directory...
|
255
256
257
|
return result;
}
|
6c4afcac
David Mayerich
introduced a gene...
|
258
259
260
261
|
//returns a pointer to the raw matrix data (in column major format)
T* data(){
return M;
}
|
8e4f8364
David Mayerich
started a new opt...
|
262
|
|
6c4afcac
David Mayerich
introduced a gene...
|
263
|
//return a transposed matrix
|
65c1aac0
David Mayerich
fixed memory leak...
|
264
|
matrix<T> transpose() const {
|
6c4afcac
David Mayerich
introduced a gene...
|
265
266
267
268
269
270
271
|
matrix<T> result(C, R);
size_t c, r;
for(c = 0; c < C; c++){
for(r = 0; r < R; r++){
result.M[r * C + c] = M[c * R + r];
}
}
|
8e4f8364
David Mayerich
started a new opt...
|
272
273
274
|
return result;
}
|
ef1f6218
David Mayerich
fixed some memory...
|
275
276
277
278
279
280
|
// Reshapes the matrix in place
void reshape(size_t rows, size_t cols) {
R = rows;
C = cols;
}
|
65c1aac0
David Mayerich
fixed memory leak...
|
281
282
283
284
285
286
287
288
289
290
291
|
///Calculate and return the determinant of the matrix
T det() const {
if (R != C) {
std::cout << "ERROR: a determinant can only be calculated for a square matrix." << std::endl;
exit(1);
}
if (R == 1) return M[0]; //if the matrix only contains one value, return it
int r, c, ri, cia, cib;
T a = 0;
T b = 0;
|
f7cf19e4
David Mayerich
fixed warnings, t...
|
292
|
for (c = 0; c < (int)C; c++) {
|
65c1aac0
David Mayerich
fixed memory leak...
|
293
294
|
for (r = 0; r < R; r++) {
ri = r;
|
f7cf19e4
David Mayerich
fixed warnings, t...
|
295
296
|
cia = (r + c) % (int)C;
cib = ((int)C - 1 - r) % (int)C;
|
65c1aac0
David Mayerich
fixed memory leak...
|
297
298
299
300
301
302
303
|
a += get(ri, cia);
b += get(ri, cib);
}
}
return a - b;
}
|
d299c720
David Mayerich
added row and col...
|
304
305
306
307
308
309
310
311
|
/// Sum all elements in the matrix
T sum() const {
size_t N = R * C; //calculate the number of elements in the matrix
T s = (T)0; //allocate a register to store the sum
for (size_t n = 0; n < N; n++) s += M[n]; //perform the summation
return s;
}
|
1c92dffe
David Mayerich
simplified CUDA i...
|
312
|
/// Sort rows of the matrix by the specified indices
|
65c1aac0
David Mayerich
fixed memory leak...
|
313
|
matrix<T> sort_rows(size_t* idx) const {
|
1c92dffe
David Mayerich
simplified CUDA i...
|
314
315
316
317
318
319
320
321
322
323
324
|
matrix<T> result(C, R); //create the output matrix
size_t r, c;
for (c = 0; c < C; c++) { //for each column
for (r = 0; r < R; r++) { //for each row element
result.M[c * R + r] = M[c * R + idx[r]]; //copy each element of the row into its new position
}
}
return result;
}
/// Sort columns of the matrix by the specified indices
|
7d01bb90
David Mayerich
matlab matrix sav...
|
325
|
matrix<T> sort_cols(size_t* idx, size_t data_type = mat4_float) const {
|
1c92dffe
David Mayerich
simplified CUDA i...
|
326
327
328
329
330
331
332
333
|
matrix<T> result(C, R);
size_t c;
for (c = 0; c < C; c++) { //for each column
memcpy(&result.M[c * R], &M[idx[c] * R], sizeof(T) * R); //copy the entire column from this matrix to the appropriate location
}
return result;
}
|
d299c720
David Mayerich
added row and col...
|
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
|
/// Return the column specified by index i
matrix<T> col(size_t i) {
matrix<T> c(R, 1); //create a single column matrix
memcpy(c.data(), &data()[R*i], C * sizeof(T)); //copy the column
return c;
}
/// Return the row specified by index i
matrix<T> row(size_t i) {
matrix<T> r(1, C); //create a single row matrix
for (size_t c = 0; c < C; c++)
r(0, c) = at(i, c);
return r;
}
|
65c1aac0
David Mayerich
fixed memory leak...
|
349
|
std::string toStr() const {
|
57729e5b
David Mayerich
changed directory...
|
350
351
|
std::stringstream ss;
|
1c92dffe
David Mayerich
simplified CUDA i...
|
352
|
for(int r = 0; r < R; r++) {
|
b2d10ff1
Pavel Govyadinov
final edits, expa...
|
353
|
ss << "| ";
|
1c92dffe
David Mayerich
simplified CUDA i...
|
354
|
for(int c=0; c<C; c++) {
|
6c4afcac
David Mayerich
introduced a gene...
|
355
|
ss << M[c * R + r] << " ";
|
57729e5b
David Mayerich
changed directory...
|
356
|
}
|
b2d10ff1
Pavel Govyadinov
final edits, expa...
|
357
|
ss << "|" << std::endl;
|
57729e5b
David Mayerich
changed directory...
|
358
|
}
|
57729e5b
David Mayerich
changed directory...
|
359
360
|
return ss.str();
}
|
19a3973b
David Mayerich
updates to matrix...
|
361
|
|
65c1aac0
David Mayerich
fixed memory leak...
|
362
363
|
void csv(std::ostream& out) const {
//std::stringstream csvss;
|
c495d065
David Mayerich
Fixed stim::matri...
|
364
|
for (size_t i = 0; i < R; i++) {
|
65c1aac0
David Mayerich
fixed memory leak...
|
365
366
367
368
|
out << std::fixed << M[i];
for (size_t j = 1; j < C; j++)
out << ", " << std::fixed << M[j * R + i];
out << std::endl;
|
1c92dffe
David Mayerich
simplified CUDA i...
|
369
|
}
|
65c1aac0
David Mayerich
fixed memory leak...
|
370
371
372
373
374
375
376
377
|
//return csvss.str();
}
std::string csv() const {
std::stringstream csvss;
int digits = std::numeric_limits<double>::max_digits10;
csvss.precision(digits);
csv(csvss);
|
1c92dffe
David Mayerich
simplified CUDA i...
|
378
379
380
|
return csvss.str();
}
|
65c1aac0
David Mayerich
fixed memory leak...
|
381
382
|
|
1c92dffe
David Mayerich
simplified CUDA i...
|
383
|
//save the data as a CSV file
|
65c1aac0
David Mayerich
fixed memory leak...
|
384
|
void csv(std::string filename) const {
|
1c92dffe
David Mayerich
simplified CUDA i...
|
385
386
387
388
389
|
ofstream basisfile(filename.c_str());
basisfile << csv();
basisfile.close();
}
|
dd5aab2f
David Mayerich
fixed errors in s...
|
390
|
static matrix<T> I(size_t N) {
|
dd5aab2f
David Mayerich
fixed errors in s...
|
391
392
393
394
395
396
397
398
|
matrix<T> result(N, N); //create the identity matrix
memset(result.M, 0, N * N * sizeof(T)); //set the entire matrix to zero
for (size_t n = 0; n < N; n++) {
result(n, n) = (T)1; //set the diagonal component to 1
}
return result;
}
|
65c1aac0
David Mayerich
fixed memory leak...
|
399
400
401
402
403
404
405
406
407
408
409
410
411
|
//loads a matrix from a stream in CSV format
void csv(std::istream& in) {
size_t c, r;
T v;
for (r = 0; r < R; r++) {
for (c = 0; c < C; c++) {
in >> v;
if (in.peek() == ',') in.seekg(1, std::ios::cur);
at(r, c) = v;;
}
}
}
|
7d01bb90
David Mayerich
matlab matrix sav...
|
412
|
// saves the matrix as a Level-4 MATLAB file
|
91589064
David Mayerich
added a global ma...
|
413
414
415
416
417
418
419
420
421
422
|
void mat4(std::string filename, std::string name = std::string("unknown"), mat4Format format = mat4_float) {
if (format == mat4_float) {
if (sizeof(T) == 4) format = mat4_float32;
else if (sizeof(T) == 8) format = mat4_float64;
else {
std::cout << "stim::matrix ERROR - incorrect format specified" << std::endl;
exit(1);
}
}
stim::save_mat4((char*)M, filename, name, rows(), cols(), format);
|
7d01bb90
David Mayerich
matlab matrix sav...
|
423
|
}
|
57729e5b
David Mayerich
changed directory...
|
424
425
|
};
|
57729e5b
David Mayerich
changed directory...
|
426
427
|
} //end namespace rts
|
57729e5b
David Mayerich
changed directory...
|
428
429
|
#endif
|