1 /* Matrix.cpp 2 * 3 * Copyright (C) 1992-2021 Paul Boersma 4 * 5 * This code is free software; you can redistribute it and/or modify 6 * it under the terms of the GNU General Public License as published by 7 * the Free Software Foundation; either version 2 of the License, or (at 8 * your option) any later version. 9 * 10 * This code is distributed in the hope that it will be useful, but 11 * WITHOUT ANY WARRANTY; without even the implied warranty of 12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. 13 * See the GNU General Public License for more details. 14 * 15 * You should have received a copy of the GNU General Public License 16 * along with this work. If not, see <http://www.gnu.org/licenses/>. 17 */ 18 19 #include "Matrix.h" 20 #include "NUM2.h" 21 #include "Formula.h" 22 #include "Eigen.h" 23 24 #include "oo_DESTROY.h" 25 #include "Matrix_def.h" 26 #include "oo_COPY.h" 27 #include "Matrix_def.h" 28 #include "oo_EQUAL.h" 29 #include "Matrix_def.h" 30 #include "oo_CAN_WRITE_AS_ENCODING.h" 31 #include "Matrix_def.h" 32 #include "oo_WRITE_TEXT.h" 33 #include "Matrix_def.h" 34 #include "oo_WRITE_BINARY.h" 35 #include "Matrix_def.h" 36 #include "oo_READ_BINARY.h" 37 #include "Matrix_def.h" 38 #include "oo_DESCRIPTION.h" 39 #include "Matrix_def.h" 40 41 Thing_implement (Matrix, SampledXY, 2); 42 43 void structMatrix :: v_info () { 44 structDaata :: v_info (); 45 double minimum = 0.0, maximum = 0.0; 46 Matrix_getWindowExtrema (this, 1, our nx, 1, our ny, & minimum, & maximum); 47 MelderInfo_writeLine (U"xmin: ", our xmin); 48 MelderInfo_writeLine (U"xmax: ", our xmax); 49 MelderInfo_writeLine (U"Number of columns: ", our nx); 50 MelderInfo_writeLine (U"dx: ", our dx, U" (-> sampling rate ", 1.0 / our dx, U" )"); 51 MelderInfo_writeLine (U"x1: ", our x1); 52 MelderInfo_writeLine (U"ymin: ", our ymin); 53 MelderInfo_writeLine (U"ymax: ", our ymax); 54 MelderInfo_writeLine (U"Number of rows: ", our ny); 55 MelderInfo_writeLine (U"dy: ", our dy, U" (-> sampling rate ", 1.0 / our dy, U" )"); 56 MelderInfo_writeLine (U"y1: ", our y1); 57 MelderInfo_writeLine (U"Minimum value: ", minimum); 58 MelderInfo_writeLine (U"Maximum value: ", maximum); 59 } 60 61 void structMatrix :: v_readText (MelderReadText text, int formatVersion) { 62 if (formatVersion < 0) { 63 our xmin = texgetr64 (text); 64 our xmax = texgetr64 (text); 65 our ymin = texgetr64 (text); 66 our ymax = texgetr64 (text); 67 our nx = texgeti32 (text); 68 our ny = texgeti32 (text); 69 our dx = texgetr64 (text); 70 our dy = texgetr64 (text); 71 our x1 = texgetr64 (text); 72 our y1 = texgetr64 (text); 73 } else { 74 Matrix_Parent :: v_readText (text, formatVersion); 75 } 76 Melder_require (our xmin <= our xmax, U"xmin should be less than or equal to xmax."); 77 Melder_require (our ymin <= our ymax, U"ymin should be less than or equal to ymax."); 78 Melder_require (our nx >= 1, U"nx should be at least 1."); 79 Melder_require (our ny >= 1, U"ny should be at least 1."); 80 Melder_require (our dx > 0.0, U"dx should be greater than 0.0."); 81 Melder_require (our dy > 0.0, U"dy should be greater than 0.0."); 82 our z = matrix_readText_r64 (our ny, our nx, text, "z"); 83 } 84 85 double structMatrix :: v_getValueAtSample (integer isamp, integer ilevel, int unit) { 86 const double value = our z [ilevel] [isamp]; 87 return ( isdefined (value) ? our v_convertStandardToSpecialUnit (value, ilevel, unit) : undefined ); 88 } 89 90 double structMatrix :: v_getMatrix (integer irow, integer icol) { 91 if (irow < 1 || irow > our ny) 92 return 0.0; 93 if (icol < 1 || icol > our nx) 94 return 0.0; 95 return z [irow] [icol]; 96 } 97 98 double structMatrix :: v_getFunction2 (double x, double y) { 99 const double rrow = (y - our y1) / our dy + 1.0; 100 const double rcol = (x - our x1) / our dx + 1.0; 101 const integer irow = Melder_ifloor (rrow), icol = Melder_ifloor (rcol); 102 const double drow = rrow - irow, dcol = rcol - icol; 103 const double z1 = irow < 1 || irow > our ny || icol < 1 || icol > our nx ? 0.0 : z [irow] [icol]; 104 const double z2 = irow < 0 || irow >= our ny || icol < 1 || icol > our nx ? 0.0 : z [irow + 1] [icol]; 105 const double z3 = irow < 1 || irow > our ny || icol < 0 || icol >= our nx ? 0.0 : z [irow] [icol + 1]; 106 const double z4 = irow < 0 || irow >= our ny || icol < 0 || icol >= our nx ? 0.0 : z [irow + 1] [icol + 1]; 107 return (1.0 - drow) * (1.0 - dcol) * z1 + drow * (1.0 - dcol) * z2 + (1.0 - drow) * dcol * z3 + drow * dcol * z4; 108 } 109 110 void Matrix_init 111 (Matrix me, double xmin, double xmax, integer nx, double dx, double x1, 112 double ymin, double ymax, integer ny, double dy, double y1) 113 { 114 Sampled_init (me, xmin, xmax, nx, dx, x1); 115 my ymin = ymin; 116 my ymax = ymax; 117 my ny = ny; 118 my dy = dy; 119 my y1 = y1; 120 my z = zero_MAT (my ny, my nx); 121 } 122 123 autoMatrix Matrix_create 124 (double xmin, double xmax, integer nx, double dx, double x1, 125 double ymin, double ymax, integer ny, double dy, double y1) 126 { 127 try { 128 autoMatrix me = Thing_new (Matrix); 129 Matrix_init (me.get(), xmin, xmax, nx, dx, x1, ymin, ymax, ny, dy, y1); 130 return me; 131 } catch (MelderError) { 132 Melder_throw (U"Matrix object not created."); 133 } 134 } 135 136 autoMatrix Matrix_createSimple (integer numberOfRows, integer numberOfColumns) { 137 try { 138 autoMatrix me = Thing_new (Matrix); 139 Matrix_init (me.get(), 140 0.5, numberOfColumns + 0.5, numberOfColumns, 1.0, 1.0, 141 0.5, numberOfRows + 0.5, numberOfRows , 1.0, 1.0); 142 return me; 143 } catch (MelderError) { 144 Melder_throw (U"Matrix object not created."); 145 } 146 } 147 148 double Matrix_columnToX (Matrix me, double column) { return my x1 + (column - 1.0) * my dx; } // FIXME inline and use Sampled 149 150 double Matrix_rowToY (Matrix me, double row) { return my y1 + (row - 1.0) * my dy; } 151 152 double Matrix_xToColumn (Matrix me, double x) { return (x - my x1) / my dx + 1.0; } 153 154 integer Matrix_xToLowColumn (Matrix me, double x) { return Melder_ifloor (Matrix_xToColumn (me, x)); } 155 156 integer Matrix_xToHighColumn (Matrix me, double x) { return Melder_iceiling (Matrix_xToColumn (me, x)); } 157 158 integer Matrix_xToNearestColumn (Matrix me, double x) { return Melder_iround (Matrix_xToColumn (me, x)); } 159 160 double Matrix_yToRow (Matrix me, double y) { return (y - my y1) / my dy + 1.0; } 161 162 integer Matrix_yToLowRow (Matrix me, double y) { return Melder_ifloor (Matrix_yToRow (me, y)); } 163 164 integer Matrix_yToHighRow (Matrix me, double y) { return Melder_iceiling (Matrix_yToRow (me, y)); } 165 166 integer Matrix_yToNearestRow (Matrix me, double y) { return Melder_iround (Matrix_yToRow (me, y)); } 167 168 integer Matrix_getWindowSamplesX (Matrix me, double xmin, double xmax, integer *ixmin, integer *ixmax) { 169 *ixmin = 1 + Melder_iceiling ((xmin - my x1) / my dx); 170 *ixmax = 1 + Melder_ifloor ((xmax - my x1) / my dx); 171 if (*ixmin < 1) 172 *ixmin = 1; 173 if (*ixmax > my nx) 174 *ixmax = my nx; 175 if (*ixmin > *ixmax) 176 return 0; 177 return *ixmax - *ixmin + 1; 178 } 179 180 integer Matrix_getWindowSamplesY (Matrix me, double ymin, double ymax, integer *iymin, integer *iymax) { 181 *iymin = 1 + Melder_iceiling ((ymin - my y1) / my dy); 182 *iymax = 1 + Melder_ifloor ((ymax - my y1) / my dy); 183 if (*iymin < 1) 184 *iymin = 1; 185 if (*iymax > my ny) 186 *iymax = my ny; 187 if (*iymin > *iymax) 188 return 0; 189 return *iymax - *iymin + 1; 190 } 191 192 integer Matrix_getWindowExtrema (Matrix me, integer ixmin, integer ixmax, integer iymin, integer iymax, 193 double *minimum, double *maximum) 194 { 195 if (ixmin == 0) // default = all 196 ixmin = 1; 197 if (ixmax == 0) // default = all 198 ixmax = my nx; 199 if (iymin == 0) // default = all 200 iymin = 1; 201 if (iymax == 0) // default = all 202 iymax = my ny; 203 if (ixmin > ixmax || iymin > iymax) 204 return 0; 205 MelderExtremaWithInit extrema; 206 for (integer iy = iymin; iy <= iymax; iy ++) 207 for (integer ix = ixmin; ix <= ixmax; ix ++) 208 extrema.update (my z [iy] [ix]); 209 *minimum = extrema.min; 210 *maximum = extrema.max; 211 return (ixmax - ixmin + 1) * (iymax - iymin + 1); 212 } 213 214 double Matrix_getValueAtXY (Matrix me, double x, double y) { 215 const double row_real = (y - my y1) / my dy + 1.0; 216 const double col_real = (x - my x1) / my dx + 1.0; 217 /* 218 * We imagine a unit square around every (xi, yi) point in the matrix. 219 * For (x, y) values outside the union of these squares, the z value is undefined. 220 */ 221 if (row_real < 0.5 || row_real > my ny + 0.5) 222 return undefined; 223 if (col_real < 0.5 || col_real > my nx + 0.5) 224 return undefined; 225 /* 226 * Determine the four nearest (xi, yi) points. 227 */ 228 integer bottomRow = Melder_ifloor (row_real); // 0 <= bottomRow <= my ny 229 integer topRow = bottomRow + 1; // 1 <= topRow <= my ny + 1 230 integer leftCol = Melder_ifloor (col_real); // 0 <= leftCol <= my nx 231 integer rightCol = leftCol + 1; // 1 <= rightCol <= my nx + 1 232 double drow = row_real - bottomRow; // 0.0 <= drow < 1.0 233 double dcol = col_real - leftCol; // 0.0 <= dcol < 1.0 234 /* 235 * If adjacent points exist 236 * (i.e., both row numbers are between 1 and my ny, 237 * or both column numbers are between 1 and my nx), 238 * we do linear interpolation. 239 * If not, we do constant extrapolation, 240 * which can be simulated by an interpolation between equal z values. 241 */ 242 if (bottomRow < 1) 243 bottomRow = 1; // 1 <= bottomRow <= my ny 244 if (topRow > my ny) 245 topRow = my ny; // 1 <= topRow <= my ny 246 if (leftCol < 1) 247 leftCol = 1; // 1 <= leftCol <= my nx 248 if (rightCol > my nx) 249 rightCol = my nx; // 1 <= rightCol <= my nx 250 return (1.0 - drow) * (1.0 - dcol) * my z [bottomRow] [leftCol] + 251 drow * (1.0 - dcol) * my z [topRow] [leftCol] + 252 (1.0 - drow) * dcol * my z [bottomRow] [rightCol] + 253 drow * dcol * my z [topRow] [rightCol]; 254 } 255 256 double Matrix_getSum (Matrix me) { 257 return NUMsum (my z.all()); 258 } 259 260 double Matrix_getNorm (Matrix me) { 261 return NUMnorm (my z.get(), 2.0); 262 } 263 264 void Matrix_drawRows (Matrix me, Graphics g, double xmin, double xmax, double ymin, double ymax, 265 double minimum, double maximum) 266 { 267 Function_unidirectionalAutowindow (me, & xmin, & xmax); 268 if (ymax <= ymin) { 269 ymin = my ymin; 270 ymax = my ymax; 271 } 272 integer ixmin, ixmax, iymin, iymax; 273 (void) Matrix_getWindowSamplesX (me, xmin, xmax, & ixmin, & ixmax); 274 (void) Matrix_getWindowSamplesY (me, ymin, ymax, & iymin, & iymax); 275 if (maximum <= minimum) 276 (void) Matrix_getWindowExtrema (me, ixmin, ixmax, iymin, iymax, & minimum, & maximum); 277 if (maximum <= minimum) { 278 minimum -= 1.0; 279 maximum += 1.0; 280 } 281 if (xmin >= xmax) 282 return; 283 Graphics_setInner (g); 284 for (integer iy = iymin; iy <= iymax; iy ++) { 285 Graphics_setWindow (g, xmin, xmax, 286 minimum - (iy - iymin) * (maximum - minimum), 287 maximum + (iymax - iy) * (maximum - minimum)); 288 Graphics_function (g, & my z [iy] [0], ixmin, ixmax, 289 Matrix_columnToX (me, ixmin), Matrix_columnToX (me, ixmax)); 290 } 291 Graphics_unsetInner (g); 292 if (iymin < iymax) 293 Graphics_setWindow (g, xmin, xmax, my y1 + (iymin - 1.5) * my dy, my y1 + (iymax - 0.5) * my dy); 294 } 295 296 void Matrix_drawOneContour (Matrix me, Graphics g, double xmin, double xmax, double ymin, double ymax, 297 double height) 298 { 299 bool xreversed = xmin > xmax, yreversed = ymin > ymax; 300 if (xmax == xmin) { xmin = my xmin; xmax = my xmax; } 301 if (ymax == ymin) { ymin = my ymin; ymax = my ymax; } 302 if (xreversed) { double temp = xmin; xmin = xmax; xmax = temp; } 303 if (yreversed) { double temp = ymin; ymin = ymax; ymax = temp; } 304 integer ixmin, ixmax, iymin, iymax; 305 (void) Matrix_getWindowSamplesX (me, xmin, xmax, & ixmin, & ixmax); 306 (void) Matrix_getWindowSamplesY (me, ymin, ymax, & iymin, & iymax); 307 if (xmin == xmax || ymin == ymax) return; 308 Graphics_setInner (g); 309 Graphics_setWindow (g, xreversed ? xmax : xmin, xreversed ? xmin : xmax, yreversed ? ymax : ymin, yreversed ? ymin : ymax); 310 Graphics_contour (g, my z.part (iymin, iymax, ixmin, ixmax), 311 Matrix_columnToX (me, ixmin), Matrix_columnToX (me, ixmax), 312 Matrix_rowToY (me, iymin), Matrix_rowToY (me, iymax), 313 height); 314 Graphics_rectangle (g, xmin, xmax, ymin, ymax); 315 Graphics_unsetInner (g); 316 } 317 318 void Matrix_drawContours (Matrix me, Graphics g, double xmin, double xmax, double ymin, double ymax, 319 double minimum, double maximum) 320 { 321 double border [1 + 8]; 322 if (xmax == xmin) { xmin = my xmin; xmax = my xmax; } 323 if (ymax == ymin) { ymin = my ymin; ymax = my ymax; } 324 integer ixmin, ixmax, iymin, iymax; 325 (void) Matrix_getWindowSamplesX (me, xmin, xmax, & ixmin, & ixmax); 326 (void) Matrix_getWindowSamplesY (me, ymin, ymax, & iymin, & iymax); 327 if (maximum <= minimum) 328 (void) Matrix_getWindowExtrema (me, ixmin, ixmax, iymin, iymax, & minimum, & maximum); 329 if (maximum <= minimum) { minimum -= 1.0; maximum += 1.0; } 330 for (integer iborder = 1; iborder <= 8; iborder ++) 331 border [iborder] = minimum + iborder * (maximum - minimum) / (8 + 1); 332 if (xmin == xmax || ymin == ymax) return; 333 Graphics_setInner (g); 334 Graphics_setWindow (g, xmin, xmax, ymin, ymax); 335 Graphics_altitude (g, my z.part (iymin, iymax, ixmin, ixmax), 336 Matrix_columnToX (me, ixmin), Matrix_columnToX (me, ixmax), 337 Matrix_rowToY (me, iymin), Matrix_rowToY (me, iymax), 338 8, border); 339 Graphics_rectangle (g, xmin, xmax, ymin, ymax); 340 Graphics_unsetInner (g); 341 } 342 343 void Matrix_paintContours (Matrix me, Graphics g, double xmin, double xmax, double ymin, double ymax, 344 double minimum, double maximum) 345 { 346 double border [1 + 30]; 347 Function_unidirectionalAutowindow (me, & xmin, & xmax); 348 if (ymax <= ymin) { ymin = my ymin; ymax = my ymax; } 349 integer ixmin, ixmax, iymin, iymax; 350 (void) Matrix_getWindowSamplesX (me, xmin, xmax, & ixmin, & ixmax); 351 (void) Matrix_getWindowSamplesY (me, ymin, ymax, & iymin, & iymax); 352 if (maximum <= minimum) 353 (void) Matrix_getWindowExtrema (me, ixmin, ixmax, iymin, iymax, & minimum, & maximum); 354 if (maximum <= minimum) { minimum -= 1.0; maximum += 1.0; } 355 for (integer iborder = 1; iborder <= 30; iborder ++) 356 border [iborder] = minimum + iborder * (maximum - minimum) / (30 + 1); 357 if (xmin >= xmax || ymin >= ymax) return; 358 Graphics_setInner (g); 359 Graphics_setWindow (g, xmin, xmax, ymin, ymax); 360 Graphics_grey (g, my z.part (iymin, iymax, ixmin, ixmax), 361 Matrix_columnToX (me, ixmin), Matrix_columnToX (me, ixmax), 362 Matrix_rowToY (me, iymin), Matrix_rowToY (me, iymax), 363 30, border); 364 Graphics_rectangle (g, xmin, xmax, ymin, ymax); 365 Graphics_unsetInner (g); 366 } 367 368 static void cellArrayOrImage (Matrix me, Graphics g, double xmin, double xmax, double ymin, double ymax, 369 double minimum, double maximum, bool interpolate) 370 { 371 Function_unidirectionalAutowindow (me, & xmin, & xmax); 372 if (ymax <= ymin) { 373 ymin = my ymin; 374 ymax = my ymax; 375 } 376 integer ixmin, ixmax, iymin, iymax; 377 (void) Matrix_getWindowSamplesX (me, xmin - 0.49999 * my dx, xmax + 0.49999 * my dx, 378 & ixmin, & ixmax); 379 (void) Matrix_getWindowSamplesY (me, ymin - 0.49999 * my dy, ymax + 0.49999 * my dy, 380 & iymin, & iymax); 381 if (maximum <= minimum) 382 (void) Matrix_getWindowExtrema (me, ixmin, ixmax, iymin, iymax, & minimum, & maximum); 383 if (maximum <= minimum) { 384 minimum -= 1.0; 385 maximum += 1.0; 386 } 387 if (xmin >= xmax || ymin >= ymax) 388 return; 389 Graphics_setInner (g); 390 Graphics_setWindow (g, xmin, xmax, ymin, ymax); 391 if (interpolate) 392 Graphics_image (g, my z.part (iymin, iymax, ixmin, ixmax), 393 Sampled_indexToX (me, ixmin - 0.5), Sampled_indexToX (me, ixmax + 0.5), 394 SampledXY_indexToY (me, iymin - 0.5), SampledXY_indexToY (me, iymax + 0.5), 395 minimum, maximum); 396 else 397 Graphics_cellArray (g, my z.part (iymin, iymax, ixmin, ixmax), 398 Sampled_indexToX (me, ixmin - 0.5), Sampled_indexToX (me, ixmax + 0.5), 399 SampledXY_indexToY (me, iymin - 0.5), SampledXY_indexToY (me, iymax + 0.5), 400 minimum, maximum); 401 Graphics_rectangle (g, xmin, xmax, ymin, ymax); 402 Graphics_unsetInner (g); 403 } 404 405 void Matrix_paintImage (Matrix me, Graphics g, double xmin, double xmax, double ymin, double ymax, 406 double minimum, double maximum) 407 { 408 cellArrayOrImage (me, g, xmin, xmax, ymin, ymax, minimum, maximum, true); 409 } 410 411 void Matrix_paintCells (Matrix me, Graphics g, double xmin, double xmax, double ymin, double ymax, 412 double minimum, double maximum) 413 { 414 cellArrayOrImage (me, g, xmin, xmax, ymin, ymax, minimum, maximum, false); 415 } 416 417 void Matrix_paintSurface (Matrix me, Graphics g, double xmin, double xmax, double ymin, double ymax, 418 double minimum, double maximum, double elevation, double azimuth) 419 { 420 Function_unidirectionalAutowindow (me, & xmin, & xmax); 421 if (ymax <= ymin) { ymin = my ymin; ymax = my ymax; } 422 integer ixmin, ixmax, iymin, iymax; 423 (void) Matrix_getWindowSamplesX (me, xmin, xmax, & ixmin, & ixmax); 424 (void) Matrix_getWindowSamplesY (me, ymin, ymax, & iymin, & iymax); 425 if (maximum <= minimum) 426 (void) Matrix_getWindowExtrema (me, ixmin, ixmax, iymin, iymax, & minimum, & maximum); 427 if (maximum <= minimum) { minimum -= 1.0; maximum += 1.0; } 428 Graphics_setInner (g); 429 Graphics_setWindow (g, -1.0, 1.0, minimum, maximum); 430 Graphics_surface (g, my z.part (iymin, iymax, ixmin, ixmax), 431 Matrix_columnToX (me, ixmin), Matrix_columnToX (me, ixmax), 432 Matrix_rowToY (me, iymin), Matrix_rowToY (me, iymax), 433 minimum, maximum, elevation, azimuth); 434 Graphics_unsetInner (g); 435 } 436 437 void Matrix_playMovie (Matrix me, Graphics g) { 438 Melder_require (my ny >= 2, 439 me, U": cannot play a movie for a Matrix with less than 2 rows."); 440 if (my xmin == my xmax || my ymin == my ymax) 441 return; 442 autoVEC column = raw_VEC (my ny); 443 double minimum = 0.0, maximum = 1.0; 444 Matrix_getWindowExtrema (me, 1, my nx, 1, my ny, & minimum, & maximum); 445 if (minimum == maximum) { 446 minimum -= 0.5; 447 maximum += 0.5; 448 } 449 for (integer icol = 1; icol <= my nx; icol ++) { 450 column.all() <<= my z.column (icol); 451 Graphics_beginMovieFrame (g, & Melder_WHITE); 452 Graphics_setWindow (g, my ymin, my ymax, minimum, maximum); 453 Graphics_function (g, column.asArgumentToFunctionThatExpectsOneBasedArray(), 1, my ny, my ymin, my ymax); 454 Graphics_endMovieFrame (g, 0.03); 455 } 456 } 457 458 autoMatrix Matrix_readAP (MelderFile file) { 459 try { 460 autofile f = Melder_fopen (file, "rb"); 461 int16 header [256]; 462 for (integer i = 0; i < 256; i ++) 463 header [i] = bingeti16LE (f); 464 const integer numberOfFrames = header [34]; 465 const integer numberOfWordsPerFrame = header [35]; 466 const double samplingFrequency = double (header [100]); // converting up (from 16 to 54 bits) 467 //const integer numberOfSamplesPerFrame = header [110]; 468 Melder_casual (U"Sampling frequency ", samplingFrequency); 469 autoMatrix me = Matrix_create (0.0, (double) numberOfFrames, numberOfFrames, 1.0, 0.5, 470 0.0, (double) numberOfWordsPerFrame, numberOfWordsPerFrame, 1.0, 0.5); 471 Melder_casual (U"... Loading ", numberOfFrames, U" frames", 472 U" of ", numberOfWordsPerFrame, U" words each ..."); 473 for (integer i = 1; i <= my nx; i ++) 474 for (integer j = 1; j <= my ny; j ++) 475 my z [j] [i] = bingeti16LE (f); // converting up (from 16 to 54 bits) 476 477 /* 478 Get pitch frequencies. 479 */ 480 for (integer i = 1; i <= my nx; i ++) 481 if (my z [1] [i] != 0.0) 482 my z [1] [i] = - samplingFrequency / my z [1] [i]; 483 484 f.close (file); 485 return me; 486 } catch (MelderError) { 487 Melder_throw (U"Matrix object not read from AP file ", file); 488 } 489 } 490 491 autoMatrix Matrix_appendRows (Matrix me, Matrix thee, ClassInfo klas) { 492 try { 493 autoMatrix him = Thing_newFromClass (klas).static_cast_move<structMatrix>(); 494 Matrix_init (him.get(), 495 std::min (my xmin, thy xmin), 496 std::max (my xmax, thy xmax), 497 std::max (my nx, thy nx), 498 my dx, 499 std::min (my x1, thy x1), 500 my ymin, my ymax + (thy ymax - thy ymin), my ny + thy ny, my dy, my y1 501 ); 502 for (integer irow = 1; irow <= my ny; irow ++) 503 for (integer icol = 1; icol <= my nx; icol ++) 504 his z [irow] [icol] = my z [irow] [icol]; 505 for (integer irow = 1; irow <= thy ny; irow ++) 506 for (integer icol = 1; icol <= thy nx; icol ++) 507 his z [irow + my ny] [icol] = thy z [irow] [icol]; 508 return him; 509 } catch (MelderError) { 510 Melder_throw (me, U" & ", thee, U": rows not appended."); 511 } 512 } 513 514 autoMatrix Matrix_readFromRawTextFile (MelderFile file) { // BUG: not Unicode-compatible (use of fscanf) 515 try { 516 autofile f = Melder_fopen (file, "rb"); 517 518 /* 519 Count columns. 520 */ 521 integer numberOfColumns = 0; 522 for (;;) { 523 int kar = fgetc (f); 524 if (kar == EOF || Melder_isVerticalSpace ((char32) kar)) 525 break; 526 if (Melder_isHorizontalSpace ((char32) kar)) 527 continue; 528 numberOfColumns ++; 529 do { 530 kar = fgetc (f); 531 } while (kar != EOF && ! Melder_isHorizontalOrVerticalSpace ((char32) kar)); 532 if (kar == EOF || Melder_isVerticalSpace ((char32) kar)) 533 break; 534 } 535 if (numberOfColumns == 0) 536 Melder_throw (U"File empty"); 537 538 /* 539 Count elements. 540 */ 541 rewind (f); 542 integer numberOfElements = 0; 543 for (;;) { 544 double element; 545 if (fscanf (f, "%lf", & element) < 1) 546 break; // zero or end-of-file 547 numberOfElements ++; 548 } 549 550 /* 551 Check if all columns are complete. 552 */ 553 if (numberOfElements == 0 || numberOfElements % numberOfColumns != 0) 554 Melder_throw (U"The number of elements (", numberOfElements, U") is not a multiple of the number of columns (", numberOfColumns, U")."); 555 556 /* 557 Create simple matrix. 558 */ 559 integer numberOfRows = numberOfElements / numberOfColumns; 560 autoMatrix me = Matrix_createSimple (numberOfRows, numberOfColumns); 561 562 /* 563 Read elements. 564 */ 565 rewind (f); 566 for (integer irow = 1; irow <= numberOfRows; irow ++) 567 for (integer icol = 1; icol <= numberOfColumns; icol ++) 568 fscanf (f, "%lf", & my z [irow] [icol]); 569 570 f.close (file); 571 return me; 572 } catch (MelderError) { 573 Melder_throw (U"Matrix object not read from raw text file ", file, U"."); 574 } 575 } 576 577 static bool isSymmetric (Matrix me) { 578 for (integer irow = 1; irow <= my ny - 1; irow ++) 579 for (integer icol = irow + 1; icol <= my nx; icol ++) 580 if (my z [irow] [icol] != my z [icol] [irow]) 581 return false; 582 return true; 583 } 584 585 void Matrix_eigen (Matrix me, autoMatrix *out_eigenvectors, autoMatrix *out_eigenvalues) { 586 try { 587 Melder_require (my nx == my ny, 588 U"The number of rows (here ", my ny, U") should be equal to the number of columns (here ", my nx, U")."); 589 Melder_require (isSymmetric (me), 590 U"The matrix should be symmetric."); 591 autoEigen eigen = Thing_new (Eigen); 592 Eigen_initFromSymmetricMatrix (eigen.get(), my z.get()); 593 autoMatrix eigenvectors = Data_copy (me); 594 autoMatrix eigenvalues = Matrix_create (1.0, 1.0, 1, 1.0, 1.0, my ymin, my ymax, my ny, my dy, my y1); 595 for (integer i = 1; i <= my nx; i ++) { 596 eigenvalues -> z [i] [1] = eigen -> eigenvalues [i]; 597 for (integer j = 1; j <= my nx; j ++) 598 eigenvectors -> z [i] [j] = eigen -> eigenvectors [j] [i]; 599 } 600 *out_eigenvectors = eigenvectors.move(); 601 *out_eigenvalues = eigenvalues.move(); 602 } catch (MelderError) { 603 Melder_throw (me, U": eigenstructure not computed."); 604 } 605 } 606 607 autoMatrix Matrix_power (Matrix me, integer power) { 608 try { 609 Melder_require (my nx == my ny, 610 U"The number of rows (here ", my ny, U") should be equal to the number of columns (here ", my nx, U")."); 611 autoMatrix thee = Data_copy (me); 612 autoMatrix him = Data_copy (me); 613 for (integer ipow = 2; ipow <= power; ipow ++) { 614 std::swap (his z, thy z); 615 for (integer irow = 1; irow <= my ny; irow ++) { 616 for (integer icol = 1; icol <= my nx; icol ++) { 617 thy z [irow] [icol] = 0.0; 618 for (integer i = 1; i <= my nx; i ++) 619 thy z [irow] [icol] += his z [irow] [i] * my z [i] [icol]; 620 } 621 } 622 } 623 return thee; 624 } catch (MelderError) { 625 Melder_throw (me, U": power not computed."); 626 } 627 } 628 629 void Matrix_writeToMatrixTextFile (Matrix me, MelderFile file) { 630 try { 631 autofile f = Melder_fopen (file, "w"); 632 fprintf (f, "\"ooTextFile\"\n\"Matrix\"\n%s %s %s %s %s\n%s %s %s %s %s\n", 633 Melder8_double (my xmin), Melder8_double (my xmax), Melder8_integer (my nx), 634 Melder8_double (my dx), Melder8_double (my x1), 635 Melder8_double (my ymin), Melder8_double (my ymax), Melder8_integer (my ny), 636 Melder8_double (my dy), Melder8_double (my y1)); 637 for (integer i = 1; i <= my ny; i ++) { 638 for (integer j = 1; j <= my nx; j ++) { 639 if (j > 1) 640 fprintf (f, " "); 641 fprintf (f, "%s", Melder8_double (my z [i] [j])); 642 } 643 fprintf (f, "\n"); 644 } 645 f.close (file); 646 } catch (MelderError) { 647 Melder_throw (me, U": not written to Matrix text file."); 648 } 649 } 650 651 void Matrix_writeToHeaderlessSpreadsheetFile (Matrix me, MelderFile file) { 652 try { 653 autofile f = Melder_fopen (file, "w"); 654 for (integer i = 1; i <= my ny; i ++) { 655 for (integer j = 1; j <= my nx; j ++) { 656 if (j > 1) 657 fprintf (f, "\t"); 658 fprintf (f, "%s", Melder8_single (my z [i] [j])); 659 } 660 fprintf (f, "\n"); 661 } 662 f.close (file); 663 } catch (MelderError) { 664 Melder_throw (me, U": not saved as tab-separated file ", file); 665 } 666 } 667 668 void Matrix_formula (Matrix me, conststring32 expression, Interpreter interpreter, Matrix target) { 669 try { 670 Formula_compile (interpreter, me, expression, kFormula_EXPRESSION_TYPE_NUMERIC, true); 671 Formula_Result result; 672 if (! target) 673 target = me; 674 for (integer irow = 1; irow <= my ny; irow ++) { 675 for (integer icol = 1; icol <= my nx; icol ++) { 676 Formula_run (irow, icol, & result); 677 target -> z [irow] [icol] = result. numericResult; 678 } 679 } 680 } catch (MelderError) { 681 Melder_throw (me, U": formula not completed."); 682 } 683 } 684 685 void Matrix_formula_part (Matrix me, double xmin, double xmax, double ymin, double ymax, 686 conststring32 expression, Interpreter interpreter, Matrix target) 687 { 688 try { 689 Function_unidirectionalAutowindow (me, & xmin, & xmax); 690 if (ymax <= ymin) { 691 ymin = my ymin; 692 ymax = my ymax; 693 } 694 integer ixmin, ixmax, iymin, iymax; 695 (void) Matrix_getWindowSamplesX (me, xmin, xmax, & ixmin, & ixmax); 696 (void) Matrix_getWindowSamplesY (me, ymin, ymax, & iymin, & iymax); 697 Formula_compile (interpreter, me, expression, kFormula_EXPRESSION_TYPE_NUMERIC, true); 698 Formula_Result result; 699 if (! target) 700 target = me; 701 for (integer irow = iymin; irow <= iymax; irow ++) { 702 for (integer icol = ixmin; icol <= ixmax; icol ++) { 703 Formula_run (irow, icol, & result); 704 target -> z [irow] [icol] = result. numericResult; 705 } 706 } 707 } catch (MelderError) { 708 Melder_throw (me, U": formula not completed."); 709 } 710 } 711 712 void Matrix_scaleAbsoluteExtremum (Matrix me, double scale) { 713 double extremum = 0.0; 714 for (integer i = 1; i <= my ny; i ++) 715 for (integer j = 1; j <= my nx; j ++) 716 if (fabs (my z [i] [j]) > extremum) 717 extremum = fabs (my z [i] [j]); 718 if (extremum != 0.0) { 719 double factor = scale / extremum; 720 for (integer i = 1; i <= my ny; i ++) 721 for (integer j = 1; j <= my nx; j ++) 722 my z [i] [j] *= factor; 723 } 724 } 725 726 autoMatrix TableOfReal_to_Matrix (TableOfReal me) { 727 try { 728 autoMatrix thee = Matrix_createSimple (my numberOfRows, my numberOfColumns); 729 for (integer i = 1; i <= my numberOfRows; i ++) 730 for (integer j = 1; j <= my numberOfColumns; j ++) 731 thy z [i] [j] = my data [i] [j]; 732 return thee; 733 } catch (MelderError) { 734 Melder_throw (me, U": not converted to Matrix."); 735 } 736 } 737 738 autoTableOfReal Matrix_to_TableOfReal (Matrix me) { 739 try { 740 autoTableOfReal thee = TableOfReal_create (my ny, my nx); 741 for (integer i = 1; i <= my ny; i ++) 742 for (integer j = 1; j <= my nx; j ++) 743 thy data [i] [j] = my z [i] [j]; 744 return thee; 745 } catch (MelderError) { 746 Melder_throw (me, U": not converted to TableOfReal."); 747 } 748 } 749 750 autoMatrix Table_to_Matrix (Table me) { 751 try { 752 autoMatrix thee = Matrix_createSimple (my rows.size, my numberOfColumns); 753 for (integer icol = 1; icol <= my numberOfColumns; icol ++) 754 Table_numericize_Assert (me, icol); 755 for (integer irow = 1; irow <= my rows.size; irow ++) { 756 TableRow row = my rows.at [irow]; 757 for (integer icol = 1; icol <= my numberOfColumns; icol ++) 758 thy z [irow] [icol] = row -> cells [icol]. number; 759 } 760 return thee; 761 } catch (MelderError) { 762 Melder_throw (me, U": not converted to Matrix."); 763 } 764 } 765 766 /* End of file Matrix.cpp */ 767