1 /* TableOfReal_extensions.cpp 2 * 3 * Copyright (C) 1993-2021 David Weenink, 2017 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. See the GNU 13 * 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 /* 20 djmw 20000202 17 typos in F1-3,L1-3 table corrected 21 djmw 20030707 Changed TableOfReal_drawVectors interface. 22 djmw 20031030 Added TableOfReal_appendColumns 23 djmw 20031216 Interface change in TableOfReal_choleskyDecomposition 24 djmw 20040108 Corrected a memory leak in TableOfReal_drawBoxPlots 25 djmw 20040211 Modified TableOfReal_copyLabels behaviour: copy NULL-labels too. 26 djmw 20040511 Removed TableOfReal_extractRowsByRowLabel,TableOfReal_selectColumnsWhereRow SVD_to_TableOfReal(SVD me,integer from,integer to)27 djmw 20040617 Removed column selection bug in TableOfReal_drawRowsAsHistogram 28 djmw 20041105 Added TableOfReal_createFromVanNieropData_25females 29 djmw 20041115 TableOfReal_drawScatterPlot: plotting a NULL-label crashed Praat. 30 djmw 20041213 Added TableOfReal_createFromWeeninkData. 31 djmw 20050221 TableOfReal_meansByRowLabels, extra reduce parameter. 32 djmw 20050222 TableOfReal_drawVectors didn't draw rowlabels. 33 djmw 20050512 TableOfReal TableOfReal_meansByRowLabels crashed if first label in sorted was NULL. 34 djmw 20051116 TableOfReal_drawScatterPlot draw reverse permited by choosing xmin > xmax and/or ymin>ymax 35 djmw 20060301 TableOfReal_meansByRowLabels extra medianize 36 djmw 20060626 Extra NULL argument for ExecRE. 37 djmw 20070822 wchar 38 djmw 20070902 Better error messages (object type and name feedback) 39 djmw 20070614 updated to version 1.30 of regular expressions. 40 djmw 20071202 Melder_warning<n> 41 djmw 20080122 float -> double 42 djmw 20081119 +TableOfReal_areAllCellsDefined 43 djmw 20090506 +setInner for _drawScatterPlotMatrix 44 djmw 20091009 +TableOfReal_drawColumnAsDistribution 45 djmw 20100222 Corrected a bug in TableOfReal_copyOneRowWithLabel which caused label corruption if 46 from and to table were equal and rows were equal too. 47 djmw 20111110 Use autostringvector 48 djmw 20111123 Always use Melder_wcscmp 49 */ 50 51 #include "Correlation.h" 52 #include "Covariance.h" 53 #include "Graphics_extensions.h" 54 #include "Matrix_extensions.h" 55 #include "NUMlapack.h" 56 #include "NUM2.h" 57 #include "SVD.h" 58 #include "Table_extensions.h" 59 #include "TableOfReal_extensions.h" 60 #include "TableOfReal_and_Permutation.h" 61 #include "regularExp.h" 62 #include "Formula.h" 63 64 #define EMPTY_STRING(s) (! (s) || s [0] == '\0') 65 66 #define Graphics_ARROW 1 67 #define Graphics_TWOWAYARROW 2 68 #define Graphics_LINE 3 69 70 static autoTableOfReal TableOfReal_TableOfReal_columnCorrelations (TableOfReal me, TableOfReal thee, bool center, bool normalize); 71 static autoTableOfReal TableOfReal_TableOfReal_rowCorrelations (TableOfReal me, TableOfReal thee, bool center, bool normalize); 72 73 integer TableOfReal_getColumnIndexAtMaximumInRow (TableOfReal me, integer rowNumber) { 74 return ( rowNumber > 0 && rowNumber <= my numberOfRows ? 75 NUMmaxPos (my data.row (rowNumber)) : 0 ); 76 } TablesOfReal_to_GSVD(TableOfReal me,TableOfReal thee)77 78 conststring32 TableOfReal_getColumnLabelAtMaximumInRow (TableOfReal me, integer rowNumber) { 79 integer columnNumber = TableOfReal_getColumnIndexAtMaximumInRow (me, rowNumber); 80 return my v_getColStr (columnNumber); 81 } 82 83 void TableOfReal_copyOneRowWithLabel (TableOfReal me, TableOfReal thee, integer myrow, integer thyrow) { 84 try { 85 if (me == thee && myrow == thyrow) 86 return; 87 Melder_require (myrow > 0 && myrow <= my numberOfRows && thyrow > 0 && thyrow <= thy numberOfRows && my numberOfColumns == thy numberOfColumns, 88 U"The dimensions do not fit."); 89 90 thy rowLabels [thyrow] = Melder_dup (my rowLabels [myrow].get()); 91 thy data.row (thyrow) <<= my data.row (myrow); 92 } catch (MelderError) { 93 Melder_throw (me, U": row ", myrow, U" not copied to ", thee); 94 } 95 } 96 97 bool TableOfReal_hasRowLabels (TableOfReal me) { 98 if (! my rowLabels) 99 return false; 100 for (integer i = 1; i <= my numberOfRows; i ++) 101 if (EMPTY_STRING (my rowLabels [i])) 102 return false; 103 return true; 104 } 105 106 bool TableOfReal_hasColumnLabels (TableOfReal me) { 107 if (! my columnLabels) 108 return false; 109 for (integer i = 1; i <= my numberOfColumns; i ++) 110 if (EMPTY_STRING (my columnLabels [i])) 111 return false; 112 return true; 113 } 114 115 autoTableOfReal TableOfReal_createIrisDataset () { 116 const double iris [150] [4] = { 117 {5.1, 3.5, 1.4, 0.2}, {4.9, 3.0, 1.4, 0.2}, {4.7, 3.2, 1.3, 0.2}, {4.6, 3.1, 1.5, 0.2}, {5.0, 3.6, 1.4, 0.2}, 118 {5.4, 3.9, 1.7, 0.4}, {4.6, 3.4, 1.4, 0.3}, {5.0, 3.4, 1.5, 0.2}, {4.4, 2.9, 1.4, 0.2}, {4.9, 3.1, 1.5, 0.1}, 119 {5.4, 3.7, 1.5, 0.2}, {4.8, 3.4, 1.6, 0.2}, {4.8, 3.0, 1.4, 0.1}, {4.3, 3.0, 1.1, 0.1}, {5.8, 4.0, 1.2, 0.2}, 120 {5.7, 4.4, 1.5, 0.4}, {5.4, 3.9, 1.3, 0.4}, {5.1, 3.5, 1.4, 0.3}, {5.7, 3.8, 1.7, 0.3}, {5.1, 3.8, 1.5, 0.3}, 121 {5.4, 3.4, 1.7, 0.2}, {5.1, 3.7, 1.5, 0.4}, {4.6, 3.6, 1.0, 0.2}, {5.1, 3.3, 1.7, 0.5}, {4.8, 3.4, 1.9, 0.2}, 122 {5.0, 3.0, 1.6, 0.2}, {5.0, 3.4, 1.6, 0.4}, {5.2, 3.5, 1.5, 0.2}, {5.2, 3.4, 1.4, 0.2}, {4.7, 3.2, 1.6, 0.2}, 123 {4.8, 3.1, 1.6, 0.2}, {5.4, 3.4, 1.5, 0.4}, {5.2, 4.1, 1.5, 0.1}, {5.5, 4.2, 1.4, 0.2}, {4.9, 3.1, 1.5, 0.2}, 124 {5.0, 3.2, 1.2, 0.2}, {5.5, 3.5, 1.3, 0.2}, {4.9, 3.6, 1.4, 0.1}, {4.4, 3.0, 1.3, 0.2}, {5.1, 3.4, 1.5, 0.2}, 125 {5.0, 3.5, 1.3, 0.3}, {4.5, 2.3, 1.3, 0.3}, {4.4, 3.2, 1.3, 0.2}, {5.0, 3.5, 1.6, 0.6}, {5.1, 3.8, 1.9, 0.4}, 126 {4.8, 3.0, 1.4, 0.3}, {5.1, 3.8, 1.6, 0.2}, {4.6, 3.2, 1.4, 0.2}, {5.3, 3.7, 1.5, 0.2}, {5.0, 3.3, 1.4, 0.2}, 127 {7.0, 3.2, 4.7, 1.4}, {6.4, 3.2, 4.5, 1.5}, {6.9, 3.1, 4.9, 1.5}, {5.5, 2.3, 4.0, 1.3}, {6.5, 2.8, 4.6, 1.5}, 128 {5.7, 2.8, 4.5, 1.3}, {6.3, 3.3, 4.7, 1.6}, {4.9, 2.4, 3.3, 1.0}, {6.6, 2.9, 4.6, 1.3}, {5.2, 2.7, 3.9, 1.4}, 129 {5.0, 2.0, 3.5, 1.0}, {5.9, 3.0, 4.2, 1.5}, {6.0, 2.2, 4.0, 1.0}, {6.1, 2.9, 4.7, 1.4}, {5.6, 2.9, 3.6, 1.3}, 130 {6.7, 3.1, 4.4, 1.4}, {5.6, 3.0, 4.5, 1.5}, {5.8, 2.7, 4.1, 1.0}, {6.2, 2.2, 4.5, 1.5}, {5.6, 2.5, 3.9, 1.1}, 131 {5.9, 3.2, 4.8, 1.8}, {6.1, 2.8, 4.0, 1.3}, {6.3, 2.5, 4.9, 1.5}, {6.1, 2.8, 4.7, 1.2}, {6.4, 2.9, 4.3, 1.3}, 132 {6.6, 3.0, 4.4, 1.4}, {6.8, 2.8, 4.8, 1.4}, {6.7, 3.0, 5.0, 1.7}, {6.0, 2.9, 4.5, 1.5}, {5.7, 2.6, 3.5, 1.0}, 133 {5.5, 2.4, 3.8, 1.1}, {5.5, 2.4, 3.7, 1.0}, {5.8, 2.7, 3.9, 1.2}, {6.0, 2.7, 5.1, 1.6}, {5.4, 3.0, 4.5, 1.5}, 134 {6.0, 3.4, 4.5, 1.6}, {6.7, 3.1, 4.7, 1.5}, {6.3, 2.3, 4.4, 1.3}, {5.6, 3.0, 4.1, 1.3}, {5.5, 2.5, 4.0, 1.3}, 135 {5.5, 2.6, 4.4, 1.2}, {6.1, 3.0, 4.6, 1.4}, {5.8, 2.6, 4.0, 1.2}, {5.0, 2.3, 3.3, 1.0}, {5.6, 2.7, 4.2, 1.3}, 136 {5.7, 3.0, 4.2, 1.2}, {5.7, 2.9, 4.2, 1.3}, {6.2, 2.9, 4.3, 1.3}, {5.1, 2.5, 3.0, 1.1}, {5.7, 2.8, 4.1, 1.3}, 137 {6.3, 3.3, 6.0, 2.5}, {5.8, 2.7, 5.1, 1.9}, {7.1, 3.0, 5.9, 2.1}, {6.3, 2.9, 5.6, 1.8}, {6.5, 3.0, 5.8, 2.2}, 138 {7.6, 3.0, 6.6, 2.1}, {4.9, 2.5, 4.5, 1.7}, {7.3, 2.9, 6.3, 1.8}, {6.7, 2.5, 5.8, 1.8}, {7.2, 3.6, 6.1, 2.5}, 139 {6.5, 3.2, 5.1, 2.0}, {6.4, 2.7, 5.3, 1.9}, {6.8, 3.0, 5.5, 2.1}, {5.7, 2.5, 5.0, 2.0}, {5.8, 2.8, 5.1, 2.4}, 140 {6.4, 3.2, 5.3, 2.3}, {6.5, 3.0, 5.5, 1.8}, {7.7, 3.8, 6.7, 2.2}, {7.7, 2.6, 6.9, 2.3}, {6.0, 2.2, 5.0, 1.5}, 141 {6.9, 3.2, 5.7, 2.3}, {5.6, 2.8, 4.9, 2.0}, {7.7, 2.8, 6.7, 2.0}, {6.3, 2.7, 4.9, 1.8}, {6.7, 3.3, 5.7, 2.1}, 142 {7.2, 3.2, 6.0, 1.8}, {6.2, 2.8, 4.8, 1.8}, {6.1, 3.0, 4.9, 1.8}, {6.4, 2.8, 5.6, 2.1}, {7.2, 3.0, 5.8, 1.6}, 143 {7.4, 2.8, 6.1, 1.9}, {7.9, 3.8, 6.4, 2.0}, {6.4, 2.8, 5.6, 2.2}, {6.3, 2.8, 5.1, 1.5}, {6.1, 2.6, 5.6, 1.4}, 144 {7.7, 3.0, 6.1, 2.3}, {6.3, 3.4, 5.6, 2.4}, {6.4, 3.1, 5.5, 1.8}, {6.0, 3.0, 4.8, 1.8}, {6.9, 3.1, 5.4, 2.1}, 145 {6.7, 3.1, 5.6, 2.4}, {6.9, 3.1, 5.1, 2.3}, {5.8, 2.7, 5.1, 1.9}, {6.8, 3.2, 5.9, 2.3}, {6.7, 3.3, 5.7, 2.5}, 146 {6.7, 3.0, 5.2, 2.3}, {6.3, 2.5, 5.0, 1.9}, {6.5, 3.0, 5.2, 2.0}, {6.2, 3.4, 5.4, 2.3}, {5.9, 3.0, 5.1, 1.8} 147 }; 148 149 try { 150 autoTableOfReal me = TableOfReal_create (150, 4); 151 152 TableOfReal_setColumnLabel (me.get(), 1, U"sl"); 153 TableOfReal_setColumnLabel (me.get(), 2, U"sw"); 154 TableOfReal_setColumnLabel (me.get(), 3, U"pl"); 155 TableOfReal_setColumnLabel (me.get(), 4, U"pw"); 156 for (integer i = 1; i <= 150; i ++) { 157 const integer kind = (i - 1) / 50 + 1; 158 char32 const *label = kind == 1 ? U"1" : kind == 2 ? U"2" : U"3"; 159 for (integer j = 1; j <= 4; j ++) 160 my data [i] [j] = iris [i - 1] [j - 1]; 161 TableOfReal_setRowLabel (me.get(), i, label); 162 } 163 Thing_setName (me.get(), U"iris"); 164 return me; 165 } catch (MelderError) { 166 Melder_throw (U"TableOfReal from iris data not created."); 167 } 168 } 169 170 autoStrings TableOfReal_extractRowLabels (TableOfReal me) { 171 try { 172 autoStrings thee = Thing_new (Strings); 173 174 if (my numberOfRows > 0) { 175 thy strings = autoSTRVEC (my numberOfRows); 176 thy numberOfStrings = my numberOfRows; 177 for (integer i = 1; i <= my numberOfRows; i ++) { 178 const conststring32 label = my rowLabels [i] ? my rowLabels [i].get() : U"?"; 179 thy strings [i] = Melder_dup (label); 180 } 181 } 182 return thee; 183 } catch (MelderError) { 184 Melder_throw (me, U": row labels not extracted."); 185 } 186 } 187 188 autoStrings TableOfReal_extractColumnLabels (TableOfReal me) { 189 try { 190 autoStrings thee = Thing_new (Strings); 191 192 if (my numberOfColumns > 0) { 193 thy strings = autoSTRVEC (my numberOfColumns); 194 thy numberOfStrings = my numberOfColumns; 195 196 for (integer i = 1; i <= my numberOfColumns; i ++) { 197 const conststring32 label = my columnLabels [i] ? my columnLabels [i].get() : U"?"; 198 thy strings [i] = Melder_dup (label); 199 } 200 } 201 return thee; 202 } catch (MelderError) { 203 Melder_throw (me, U": column labels not extracted."); 204 } 205 } 206 207 autoTableOfReal TableOfReal_transpose (TableOfReal me) { 208 try { 209 autoTableOfReal thee = TableOfReal_create (my numberOfColumns, my numberOfRows); 210 thy data.all() <<= my data.transpose(); 211 thy columnLabels.all() <<= my rowLabels.all(); 212 thy rowLabels.all() <<= my columnLabels.all(); 213 return thee; 214 } catch (MelderError) { 215 Melder_throw (me, U": not transposed."); 216 } 217 } 218 219 void TableOfReal_to_PatternList_and_Categories (TableOfReal me, integer fromrow, integer torow, integer fromcol, integer tocol, 220 autoPatternList *out_p, autoCategories *out_c) 221 { 222 try { 223 integer ncols = my numberOfColumns, nrows = my numberOfRows; 224 fromrow = fromrow == 0 ? 1 : fromrow; 225 torow = torow == 0 ? nrows : torow; 226 Melder_require (fromrow > 0 && fromrow <= torow && torow <= nrows, 227 U"Invalid row selection."); 228 fromcol = fromcol == 0 ? 1 : fromcol; 229 tocol = tocol == 0 ? ncols : tocol; 230 231 Melder_require (fromcol > 0 && fromcol <= tocol && tocol <= ncols, 232 U"Invalid column selection."); 233 234 nrows = torow - fromrow + 1; 235 ncols = tocol - fromcol + 1; 236 autoPatternList ap = PatternList_create (nrows, ncols); 237 autoCategories ac = Categories_create (); 238 239 integer row = 1; 240 for (integer i = fromrow; i <= torow; i ++, row ++) { 241 const char32 *s = my rowLabels [i] ? my rowLabels [i].get() : U"?"; 242 autoSimpleString item = SimpleString_create (s); 243 ac -> addItem_move (item.move()); 244 integer col = 1; 245 for (integer j = fromcol; j <= tocol; j ++, col ++) 246 ap -> z [row] [col] = my data [i] [j]; 247 } 248 if (out_p) 249 *out_p = ap.move(); 250 if (out_c) 251 *out_c = ac.move(); 252 } catch (MelderError) { 253 Melder_throw (U"PatternList and Categories not created from TableOfReal."); 254 } 255 } 256 257 void TableOfReal_getColumnExtrema (TableOfReal me, integer col, double *out_min, double *out_max) { 258 Melder_require (col > 0 && col <= my numberOfColumns, 259 U"Invalid column number."); 260 if (out_min) 261 *out_min = NUMmin (my data.column (col)); 262 if (out_max) 263 *out_max = NUMmax (my data.column (col)); 264 } 265 266 void TableOfReal_drawRowsAsHistogram (TableOfReal me, Graphics g, constINTVECVU const& rowNumbers, integer colb, integer cole, double ymin, 267 double ymax, double xoffsetFraction, double interbarFraction, double interbarsFraction, constVECVU const& greys, bool garnish) 268 { 269 if (colb == 0) // undefined? 270 colb = 1; // sensible default: all 271 if (cole == 0) // undefined? 272 cole = my numberOfColumns; // sensible default: all 273 274 Melder_require (NUMisSorted4 (1_integer, colb, cole, my numberOfColumns), 275 U"Invalid column numbers"); 276 277 for (integer i = 1; i <= rowNumbers.size; i ++) { 278 const integer rowNumber = rowNumbers [i]; 279 Melder_require (rowNumber > 0 && rowNumber <= my numberOfRows, 280 U"Invalid row (", rowNumber, U")."); 281 if (ymin >= ymax) { 282 double min, max; 283 NUMextrema (my data.row (rowNumber).part (colb, cole), & min, & max); 284 if (i > 1) { 285 if (min < ymin) 286 ymin = min; 287 if (max > ymax) 288 ymax = max; 289 } else { 290 ymin = min; 291 ymax = max; 292 } 293 } 294 } 295 296 Graphics_setWindow (g, 0.0, 1.0, ymin, ymax); 297 Graphics_setInner (g); 298 299 const integer ncols = cole - colb + 1, nrows = rowNumbers.size; 300 const double bar_width = 1.0 / (ncols * nrows + 2.0 * xoffsetFraction + (ncols - 1) * interbarsFraction + ncols * (nrows - 1) * interbarFraction); 301 const double dx = (interbarsFraction + nrows + (nrows - 1) * interbarFraction) * bar_width; 302 303 for (integer i = 1; i <= nrows; i ++) { 304 const integer rowNumber = rowNumbers [i]; 305 const double xb = xoffsetFraction * bar_width + (i - 1) * (1.0 + interbarFraction) * bar_width; 306 307 const double grey = greys [1 + (i - 1) % greys.size]; // cycle through the colours 308 double x1 = xb; 309 for (integer j = colb; j <= cole; j ++) { 310 const double x2 = x1 + bar_width; 311 double y2 = my data [rowNumber] [j]; 312 if (y2 > ymin) { 313 if (y2 > ymax) 314 y2 = ymax; 315 Graphics_setGrey (g, grey); 316 Graphics_fillRectangle (g, x1, x2, ymin, y2); 317 Graphics_setGrey (g, 0.0); // black 318 Graphics_rectangle (g, x1, x2, ymin, y2); 319 } 320 x1 += dx; 321 } 322 } 323 324 Graphics_unsetInner (g); 325 326 if (garnish) { 327 double xb = (xoffsetFraction + 0.5 * (nrows + (nrows - 1) * interbarFraction)) * bar_width; 328 for (integer j = colb; j <= cole; j ++) { 329 if (my columnLabels [j]) 330 Graphics_markBottom (g, xb, false, false, false, my columnLabels [j].get()); 331 xb += dx; 332 } 333 Graphics_drawInnerBox (g); 334 Graphics_marksLeft (g, 2, true, true, false); 335 } 336 } 337 338 void TableOfReal_drawBiplot (TableOfReal me, Graphics g, double xmin, double xmax, double ymin, double ymax, double sv_splitfactor, int labelsize, bool garnish) { 339 const integer nr = my numberOfRows, nc = my numberOfColumns, nPoints = nr + nc; 340 const double fontsize = Graphics_inqFontSize (g); 341 342 autoSVD svd = SVD_create (nr, nc); 343 344 svd -> u.all() <<= my data.all(); 345 centreEachColumn_MAT_inout (svd -> u.get()); 346 347 SVD_compute (svd.get()); 348 const integer numberOfZeroed = SVD_zeroSmallSingularValues (svd.get(), 0.0); 349 350 const integer nmin = std::min (nr, nc) - numberOfZeroed; 351 Melder_require (nmin > 1, 352 U"There should be at least two (independent) columns in the table."); 353 354 autoVEC x = raw_VEC (nPoints); 355 autoVEC y = raw_VEC (nPoints); 356 357 const double lambda1 = pow (svd -> d [1], sv_splitfactor); 358 const double lambda2 = pow (svd -> d [2], sv_splitfactor); 359 x.part (1, nr) <<= svd -> u.column (1) * lambda1; 360 y.part (1, nr) <<= svd -> u.column (2) * lambda2; 361 362 const double lambda3 = svd -> d [1] / lambda1; 363 const double lambda4 = svd -> d [2] / lambda2; 364 x.part (nr + 1, nPoints) <<= svd -> v.column (1) * lambda3; 365 y.part (nr + 1, nPoints) <<= svd -> v.column (2) * lambda4; 366 367 if (xmax <= xmin) 368 NUMextrema (x.get(), & xmin, & xmax); 369 370 if (xmax <= xmin) { 371 xmax += 1.0; 372 xmin -= 1.0; 373 } 374 if (ymax <= ymin) 375 NUMextrema (y.get(), & ymin, & ymax); 376 if (ymax <= ymin) { 377 ymax += 1.0; 378 ymin -= 1.0; 379 } 380 381 Graphics_setWindow (g, xmin, xmax, ymin, ymax); 382 Graphics_setInner (g); 383 if (labelsize > 0) 384 Graphics_setFontSize (g, labelsize); 385 Graphics_setTextAlignment (g, kGraphics_horizontalAlignment::CENTRE, Graphics_HALF); 386 387 for (integer i = 1; i <= nPoints; i ++) { 388 char32 const *label; 389 if (i <= nr) { 390 label = my rowLabels [i].get(); 391 if (! label) 392 label = U"?__r_"; 393 } else { 394 label = my columnLabels [i - nr].get(); 395 if (! label) 396 label = U"?__c_"; 397 } 398 Graphics_text (g, x [i], y [i], label); 399 } 400 401 Graphics_unsetInner (g); 402 403 if (garnish) { 404 Graphics_drawInnerBox (g); 405 Graphics_marksLeft (g, 2, true, true, false); 406 Graphics_marksBottom (g, 2, true, true, false); 407 } 408 409 if (labelsize > 0) 410 Graphics_setFontSize (g, fontsize); 411 } 412 413 void TableOfReal_drawBoxPlots (TableOfReal me, Graphics g, integer rowmin, integer rowmax, integer colmin, integer colmax, double ymin, double ymax, bool garnish) { 414 if (rowmax < rowmin || rowmax < 1) { 415 rowmin = 1; 416 rowmax = my numberOfRows; 417 } 418 if (rowmin < 1) 419 rowmin = 1; 420 if (rowmax > my numberOfRows) 421 rowmax = my numberOfRows; 422 423 const integer numberOfRows = rowmax - rowmin + 1; 424 if (colmax < colmin || colmax < 1) { 425 colmin = 1; 426 colmax = my numberOfColumns; 427 } 428 if (colmin < 1) 429 colmin = 1; 430 if (colmax > my numberOfColumns) 431 colmax = my numberOfColumns; 432 433 if (ymax <= ymin) { 434 MelderRealRange yrange = NUMextrema (my data.part (rowmin, rowmax, colmin, colmax)); 435 ymin = yrange.min; 436 ymax = yrange.max; 437 } 438 439 Graphics_setWindow (g, colmin - 0.5, colmax + 0.5, ymin, ymax); 440 Graphics_setInner (g); 441 442 autoVEC data = raw_VEC (numberOfRows); 443 for (integer j = colmin; j <= colmax; j ++) { 444 const double x = j, r = 0.05, w = 0.2; 445 integer ndata = 0; 446 447 for (integer i = 1; i <= numberOfRows; i ++) { 448 if (isdefined (my data [rowmin + i - 1] [j])) 449 data [ ++ ndata] = my data [rowmin + i - 1] [j]; 450 } 451 Graphics_boxAndWhiskerPlot (g, data.get(), x, r, w, ymin, ymax); 452 } 453 Graphics_unsetInner (g); 454 if (garnish) { 455 Graphics_drawInnerBox (g); 456 for (integer j = colmin; j <= colmax; j ++) { 457 if (my columnLabels && my columnLabels [j] && my columnLabels [j] [0]) 458 Graphics_markBottom (g, j, false, true, false, my columnLabels [j].get()); 459 } 460 Graphics_marksLeft (g, 2, true, true, false); 461 } 462 } 463 464 void TableOfReal_copyLabels (TableOfReal me, TableOfReal thee, int rowOrigin, int columnOrigin) { 465 if (rowOrigin == 1) { 466 Melder_require (my numberOfRows == thy numberOfRows, 467 U"Both tables must have the same number of rows."); 468 thy rowLabels.all() <<= my rowLabels.all(); 469 } else if (rowOrigin == -1) { 470 Melder_require (my numberOfColumns == thy numberOfRows, 471 U"Both tables must have the same number of columns."); 472 thy rowLabels.all() <<= my columnLabels.all(); 473 } 474 if (columnOrigin == 1) { 475 Melder_require (my numberOfColumns == thy numberOfColumns, 476 U"Both tables must have the same number of columns."); 477 thy columnLabels.all() <<= my columnLabels.all(); 478 } else if (columnOrigin == -1) { 479 Melder_require (my numberOfRows == thy numberOfColumns, 480 U"Both tables must have the same number of rows."); 481 thy columnLabels.all() <<= my rowLabels.all(); 482 } 483 } 484 485 void TableOfReal_setLabelsFromCollectionItemNames (TableOfReal me, Collection thee, bool setRowLabels, bool setColumnLabels) { 486 try { 487 if (setRowLabels) { 488 Melder_assert (my numberOfRows == thy size); 489 for (integer i = 1; i <= my numberOfRows; i ++) { 490 const conststring32 name = Thing_getName (thy at [i]); 491 TableOfReal_setRowLabel (me, i, name); 492 } 493 } 494 if (setColumnLabels) { 495 Melder_assert (my numberOfColumns == thy size); 496 for (integer i = 1; i <= my numberOfColumns; i ++) { 497 const conststring32 name = Thing_getName (thy at [i]); 498 TableOfReal_setColumnLabel (me, i, name); 499 } 500 } 501 } catch (MelderError) { 502 Melder_throw (me, U": labels not changed."); 503 } 504 } 505 506 void TableOfReal_centreColumns (TableOfReal me) { 507 centreEachColumn_MAT_inout (my data.get()); 508 } 509 510 void TableOfReal_Categories_setRowLabels (TableOfReal me, Categories thee) { 511 try { 512 Melder_require (my numberOfRows == thy size, 513 U"The number of items in both objects should be equal."); 514 /* 515 Create without change. 516 */ 517 autoCategories categories_copy = Data_copy (thee); 518 /* 519 Change without error. 520 */ 521 for (integer i = 1; i <= my numberOfRows; i ++) 522 my rowLabels [i] = categories_copy->at [i] -> string.move(); 523 } catch (MelderError) { 524 Melder_throw (me, U": row labels not set from categories."); 525 } 526 } 527 528 void TableOfReal_centreColumns_byRowLabel (TableOfReal me) { 529 conststring32 label = my rowLabels [1].get(); 530 integer index = 1; 531 for (integer i = 2; i <= my numberOfRows; i ++) { 532 const conststring32 li = my rowLabels [i].get(); 533 if (! Melder_equ (li, label)) { 534 centreEachColumn_MAT_inout (my data.horizontalBand (index, i - 1)); 535 label = li; 536 index = i; 537 } 538 } 539 centreEachColumn_MAT_inout (my data.horizontalBand (index, my numberOfRows)); 540 } 541 542 double TableOfReal_getRowSum (TableOfReal me, integer rowNumber) { 543 Melder_require (rowNumber > 0 && rowNumber <= my numberOfRows, 544 U"Row number not in valid range."); 545 return NUMsum (my data.row (rowNumber)); 546 } 547 548 double TableOfReal_getColumnSum (TableOfReal me, integer columnNumber) { 549 Melder_require (columnNumber > 0 && columnNumber <= my numberOfColumns, 550 U"Column number not in valid range."); 551 return NUMsum (my data.column (columnNumber)); 552 } 553 554 double TableOfReal_getRowSumByLabel (TableOfReal me, conststring32 rowLabel) { 555 const integer rowNumber = TableOfReal_rowLabelToIndex (me, rowLabel); 556 Melder_require (rowNumber > 0, 557 U"There is no \"", rowLabel, U"\" row label."); 558 return TableOfReal_getRowSum (me, rowNumber); 559 } 560 561 double TableOfReal_getColumnSumByLabel (TableOfReal me, conststring32 columnLabel) { 562 const integer columnNumber = TableOfReal_columnLabelToIndex (me, columnLabel); 563 Melder_require (columnNumber > 0, 564 U"There is no \"", columnLabel, U"\" column label."); 565 return TableOfReal_getColumnSum (me, columnNumber); 566 } 567 568 double TableOfReal_getGrandSum (TableOfReal me) { 569 return NUMsum (my data.all()); 570 } 571 572 void TableOfReal_centreRows (TableOfReal me) { 573 centreEachRow_MAT_inout (my data.get()); 574 } 575 576 void TableOfReal_doubleCentre (TableOfReal me) { 577 doubleCentre_MAT_inout (my data.get()); 578 } 579 580 void TableOfReal_normalizeColumns (TableOfReal me, double norm) { 581 MATnormalizeRows_inplace (my data.transpose(), 2.0, norm); 582 } 583 584 void TableOfReal_normalizeRows (TableOfReal me, double norm) { 585 MATnormalizeRows_inplace (my data.get(), 2.0, norm); 586 } 587 588 void TableOfReal_standardizeColumns (TableOfReal me) { 589 if (my numberOfRows <= 1) { 590 my data.all() <<= 0.0; 591 return; 592 } 593 for (integer icol = 1; icol <= my numberOfColumns; icol ++) { 594 MelderGaussianStats stats = NUMmeanStdev (my data.column (icol)); 595 for (integer irow = 1; irow <= my numberOfRows; irow ++) 596 my data [irow] [icol] = (my data [irow] [icol] - stats.mean) / stats.stdev; 597 } 598 } 599 600 void TableOfReal_standardizeRows (TableOfReal me) { 601 if (my numberOfColumns <= 1) { 602 my data.all() <<= 0.0; 603 return; 604 } 605 for (integer irow = 1; irow <= my numberOfRows; irow ++) { 606 MelderGaussianStats stats = NUMmeanStdev (my data.row (irow)); 607 for (integer icol = 1; icol <= my numberOfColumns; icol ++) 608 my data [irow] [icol] = (my data [irow] [icol] - stats.mean) / stats.stdev; 609 } 610 } 611 612 void TableOfReal_normalizeTable (TableOfReal me, double norm) { 613 MATnormalize_inplace (my data.get(), 2.0, norm); 614 } 615 616 double TableOfReal_getTableNorm (TableOfReal me) { 617 return NUMnorm (my data.get(), 2.0); 618 } 619 620 bool TableOfReal_isNonNegative (TableOfReal me) { 621 return NUMisNonNegative (my data.all()); 622 } 623 624 void TableOfReal_drawScatterPlotMatrix (TableOfReal me, Graphics g, integer colb, integer cole, double fractionWhite) { 625 626 if (colb == 0 && cole == 0) { 627 colb = 1; 628 cole = my numberOfColumns; 629 } else if (cole < colb || colb < 1 || cole > my numberOfColumns) 630 return; 631 632 const integer numberOfColumns = cole - colb + 1; 633 if (numberOfColumns == 1) 634 return; 635 autoVEC colmin = raw_VEC (numberOfColumns); 636 autoVEC colmax = raw_VEC (numberOfColumns); 637 638 for (integer j = 1; j <= numberOfColumns; j ++) { 639 colmin [j] = NUMmin (my data.column (colb + j - 1)); 640 colmax [j] = NUMmax (my data.column (colb + j - 1)); 641 } 642 643 for (integer j = 1; j <= numberOfColumns; j ++) { 644 double extra = fractionWhite * fabs (colmax [j] - colmin [j]); 645 if (extra == 0.0) 646 extra = 0.5; 647 colmin [j] -= extra; 648 colmax [j] += extra; 649 } 650 651 Graphics_setWindow (g, 0.0, numberOfColumns, 0.0, numberOfColumns); 652 Graphics_setInner (g); 653 Graphics_line (g, 0.0, numberOfColumns, numberOfColumns, numberOfColumns); 654 Graphics_line (g, 0.0, 0.0, 0.0, numberOfColumns); 655 Graphics_setTextAlignment (g, kGraphics_horizontalAlignment::CENTRE, Graphics_HALF); 656 657 for (integer icol = 1; icol <= numberOfColumns; icol ++) { 658 const integer ycol = colb + icol - 1; 659 Graphics_line (g, 0.0, numberOfColumns - icol, numberOfColumns, numberOfColumns - icol); 660 Graphics_line (g, icol, numberOfColumns, icol, 0.0); 661 for (integer irow = 1; irow <= numberOfColumns; irow ++) { 662 const integer xcol = colb + irow - 1; 663 if (icol == irow) { 664 conststring32 mark = my columnLabels [xcol].get(); 665 char32 label [40]; 666 if (! mark) { 667 Melder_sprint (label, 40, U"Column ", xcol); 668 mark = label; 669 } 670 Graphics_text (g, irow - 0.5, numberOfColumns - icol + 0.5, mark); 671 } else { 672 for (integer k = 1; k <= my numberOfRows; k ++) { 673 const double x = irow - 1 + (my data [k] [xcol] - colmin [irow]) / (colmax [irow] - colmin [irow]); 674 const double y = numberOfColumns - icol + (my data [k] [ycol] - colmin [icol]) / (colmax [icol] - colmin [icol]); 675 conststring32 mark = EMPTY_STRING (my rowLabels [k]) ? U"+" : my rowLabels [k].get(); 676 Graphics_text (g, x, y, mark); 677 } 678 } 679 } 680 } 681 Graphics_unsetInner (g); 682 } 683 684 void TableOfReal_drawAsScalableSquares (TableOfReal me, Graphics g, integer rowmin, integer rowmax, integer colmin, integer colmax, kGraphicsMatrixOrigin origin, double cellSizeFactor, kGraphicsMatrixCellDrawingOrder fillOrder, bool garnish) { 685 try { 686 NUMfixIndicesInRange (1, my numberOfRows, & rowmin, & rowmax); 687 NUMfixIndicesInRange (1, my numberOfColumns, & colmin, & colmax); 688 autoMatrix thee = TableOfReal_to_Matrix (me); 689 Graphics_setWindow (g, colmin - 0.5, colmax + 0.5, rowmin - 0.5, rowmax + 0.5); 690 Graphics_setInner (g); 691 Matrix_drawAsSquares_inside (thee.get(), g, colmin - 0.5, colmax + 0.5, rowmin - 0.5, rowmax + 0.5, origin, cellSizeFactor, fillOrder); 692 Graphics_unsetInner (g); 693 if (garnish) { 694 Graphics_drawInnerBox (g); 695 Graphics_marksBottomEvery (g, 1.0, 1.0, false, true, false); 696 Graphics_marksLeftEvery (g, 1.0, 1.0, false, true, false); 697 } 698 } catch (MelderError) { 699 Melder_clearError (); // drawing errors shall be ignored 700 } 701 } 702 703 void TableOfReal_drawScatterPlot (TableOfReal me, Graphics g, integer icx, integer icy, integer rowb, integer rowe, double xmin, double xmax, double ymin, double ymax, 704 integer labelSize, bool useRowLabels, conststring32 label, bool garnish) 705 { 706 const double fontSize = Graphics_inqFontSize (g); 707 708 Melder_require (icx >= 1 && icx <= my numberOfColumns, 709 U"The horizontal column number should be in the range from 1 to ", my numberOfColumns, U"."); 710 Melder_require (icy >= 1 && icy <= my numberOfColumns, 711 U"The vertical column number should be in the range from 1 to ", my numberOfColumns, U"."); 712 if (rowb < 1) 713 rowb = 1; 714 if (rowe > my numberOfRows) 715 rowe = my numberOfRows; 716 if (rowe <= rowb) { 717 rowb = 1; 718 rowe = my numberOfRows; 719 } 720 if (xmax == xmin) { 721 MelderRealRange xrange = NUMextrema (my data.part (rowb, rowe, icx, icx)); 722 const double tmp = ( xrange.max == xrange.min ? 0.5 : 0.0 ); 723 xmin = xrange.min - tmp; 724 xmax = xrange.max + tmp; 725 } 726 if (ymax == ymin) { 727 MelderRealRange yrange = NUMextrema (my data.part (rowb, rowe, icy, icy)); 728 const double tmp = ( yrange.max == yrange.min ? 0.5 : 0.0 ); 729 ymin = yrange.min - tmp; 730 ymax = yrange.max + tmp; 731 } 732 Graphics_setWindow (g, xmin, xmax, ymin, ymax); 733 Graphics_setInner (g); 734 Graphics_setTextAlignment (g, kGraphics_horizontalAlignment::CENTRE, Graphics_HALF); 735 Graphics_setFontSize (g, labelSize); 736 737 for (integer i = rowb; i <= rowe; i ++) { 738 const double x = my data [i] [icx], y = my data [i] [icy]; 739 if (((xmin < xmax && x >= xmin && x <= xmax) || (xmin > xmax && x <= xmin && x >= xmax)) && 740 ((ymin < ymax && y >= ymin && y <= ymax) || (ymin > ymax && y <= ymin && y >= ymax))) { 741 conststring32 mark = useRowLabels ? my rowLabels [i].get() : label; 742 if (mark) 743 Graphics_text (g, x, y, mark); 744 } 745 } 746 747 Graphics_setFontSize (g, fontSize); 748 Graphics_unsetInner (g); 749 750 if (garnish) { 751 Graphics_drawInnerBox (g); 752 if (ymin < ymax) { 753 if (my columnLabels [icx]) 754 Graphics_textBottom (g, true, my columnLabels [icx].get()); 755 Graphics_marksBottom (g, 2, true, true, false); 756 } else { 757 if (my columnLabels [icx]) 758 Graphics_textTop (g, true, my columnLabels [icx].get()); 759 Graphics_marksTop (g, 2, true, true, false); 760 } 761 if (xmin < xmax) { 762 if (my columnLabels [icy]) 763 Graphics_textLeft (g, true, my columnLabels [icy].get()); 764 Graphics_marksLeft (g, 2, true, true, false); 765 } else { 766 if (my columnLabels [icy]) 767 Graphics_textRight (g, true, my columnLabels [icy].get()); 768 Graphics_marksRight (g, 2, true, true, false); 769 } 770 } 771 } 772 773 774 #pragma mark - class TableOfRealList 775 776 autoTableOfReal TableOfRealList_sum (TableOfRealList me) { 777 try { 778 if (my size <= 0) 779 return autoTableOfReal(); 780 781 autoTableOfReal thee = Data_copy (my at [1]); 782 783 for (integer i = 2; i <= my size; i ++) { 784 TableOfReal him = my at [i]; 785 Melder_require (thy numberOfRows == his numberOfRows && thy numberOfColumns == his numberOfColumns 786 && NUMequal (thy rowLabels.get(), his rowLabels.get()) 787 && NUMequal (thy columnLabels.get(), his columnLabels.get()), 788 U"Dimensions or labels differ for table ", i, U"."); 789 thy data.get() += his data.get(); 790 } 791 return thee; 792 } catch (MelderError) { 793 Melder_throw (me, U": sum not created."); 794 } 795 } 796 797 bool TableOfRealList_haveIdenticalDimensions (TableOfRealList me) { 798 if (my size < 2) 799 return true; 800 const TableOfReal t1 = my at [1]; 801 for (integer i = 2; i <= my size; i ++) { 802 const TableOfReal t = my at [i]; 803 if (t -> numberOfColumns != t1 -> numberOfColumns || t -> numberOfRows != t1 -> numberOfRows) 804 return false; 805 } 806 return true; 807 } 808 809 double TableOfReal_getColumnQuantile (TableOfReal me, integer columnNumber, double quantile) { 810 try { 811 if (columnNumber < 1 || columnNumber > my numberOfColumns) 812 return undefined; 813 autoVEC values = column_VEC (my data.get(), columnNumber); 814 sort_VEC_inout (values.get()); 815 return NUMquantile (values.get(), quantile); 816 } catch (MelderError) { 817 return undefined; 818 } 819 } 820 821 autoTableOfReal TableOfReal_create_sandwell1987 () { 822 try { 823 /* 824 The following numbers are the 21 (approximate) data points in Fig, 2 in Sandwell (1987). 825 They were measured by djmw from a printed picture blown up by 800%. 826 The following numbers are in cm measured from the left (x) and the bottom (y) of the figure. 827 Vertical scale: 8.25 cm in the picture is 12 units, y [1] is at y = 0.0. 828 Horizontal scale: 17.75 cm is 10 units, x [1] is at x = 0.0. 829 */ 830 const double x [22] = { 0.0, 0.9, 2.15, 3.5, 4.75, 5.3, 6.15, 7.15, 7.95, 8.85, 9.95, 10.15, 10.3, 11.5, 12.4, 13.3, 14.2, 15.15, 16.0, 16.85, 17.25, 18.15 }; 831 const double y [22] = { 0.0, 4.2, 3.5, 4.2, 5.65, 10.1, 8.5, 7.8, 7.1, 6.4, 5.65, 0.6, 5.65, 4.2, 5.65, 7.1, 6.75, 6.35, 4.2, 2.05, 4.95, 4.25 }; 832 const integer numberOfSamples = 21; 833 autoTableOfReal thee = TableOfReal_create (numberOfSamples, 2); 834 for (integer isample = 1; isample <= numberOfSamples; isample ++) { 835 thy data [isample] [1] = (x [isample] - x [1]) * 10.0 / 17.75; 836 thy data [isample] [2] = (y [isample] - y [1]) * 12.0 / 8.25; 837 } 838 return thee; 839 } catch (MelderError) { 840 Melder_throw (U"Sandwell (1987) table not created."); 841 } 842 } 843 844 static autoTableOfReal TableOfReal_createPolsVanNieropData (int choice, bool include_levels) { 845 try { 846 autoTable table = Table_create_polsVanNierop1973 (); 847 /* 848 Default: Pols 50 males, first part of the table. 849 */ 850 const integer ncols = include_levels ? 6 : 3; 851 integer ib = 1,nrows = 50 * 12; 852 if (choice == 2) { /* Van Nierop 25 females */ 853 ib = nrows + 1; 854 nrows = 25 * 12; 855 } 856 857 autoTableOfReal thee = TableOfReal_create (nrows, ncols); 858 859 for (integer i = 1; i <= nrows; i ++) { 860 const TableRow row = table -> rows.at [ib + i - 1]; 861 TableOfReal_setRowLabel (thee.get(), i, row -> cells [4]. string.get()); 862 for (integer j = 1; j <= 3; j ++) { 863 thy data [i] [j] = Melder_atof (row -> cells [4 + j]. string.get()); 864 if (include_levels) 865 thy data [i] [3 + j] = Melder_atof (row -> cells [7 + j]. string.get()); 866 } 867 } 868 for (integer j = 1; j <= 3; j ++) { 869 conststring32 label = table -> columnHeaders [4 + j]. label.get(); 870 TableOfReal_setColumnLabel (thee.get(), j, label); 871 if (include_levels) { 872 label = table -> columnHeaders [7 + j]. label.get(); 873 TableOfReal_setColumnLabel (thee.get(), 3 + j, label); 874 } 875 } 876 return thee; 877 } catch (MelderError) { 878 Melder_throw (U"TableOfReal from Pols & Van Nierop data not created."); 879 } 880 } 881 882 autoTableOfReal TableOfReal_create_pols1973 (bool include_levels) { 883 return TableOfReal_createPolsVanNieropData (1, include_levels); 884 } 885 886 autoTableOfReal TableOfReal_create_vanNierop1973 (bool include_levels) { 887 return TableOfReal_createPolsVanNieropData (2, include_levels); 888 } 889 890 autoTableOfReal TableOfReal_create_weenink1983 (int option) { 891 try { 892 const integer nvowels = 12, ncols = 3, nrows = 10 * nvowels; 893 894 autoTable table = Table_create_weenink1983 (); 895 896 integer ib = ( option == 1 ? 1 : option == 2 ? 11 : 21 ); /* m f c*/ 897 ib = (ib - 1) * nvowels + 1; 898 899 autoTableOfReal thee = TableOfReal_create (nrows, ncols); 900 for (integer i = 1; i <= nrows; i ++) { 901 const TableRow row = table -> rows.at [ib + i - 1]; 902 TableOfReal_setRowLabel (thee.get(), i, row -> cells [5]. string.get()); 903 for (integer j = 1; j <= 3; j ++) 904 thy data [i] [j] = Melder_atof (row -> cells [6 + j]. string.get()); /* Skip F0 */ 905 } 906 for (integer j = 1; j <= 3; j ++) { 907 const conststring32 label = table -> columnHeaders [6 + j]. label.get(); 908 TableOfReal_setColumnLabel (thee.get(), j, label); 909 } 910 return thee; 911 } catch (MelderError) { 912 Melder_throw (U"TableOfReal from Weenink data not created."); 913 } 914 } 915 916 autoTableOfReal TableOfReal_randomizeRows (TableOfReal me) { 917 try { 918 autoPermutation p = Permutation_create (my numberOfRows); 919 Permutation_permuteRandomly_inplace (p.get(), 0, 0); 920 autoTableOfReal thee = TableOfReal_Permutation_permuteRows (me, p.get()); 921 return thee; 922 } catch (MelderError) { 923 Melder_throw (me, U": randomized rows not created"); 924 } 925 } 926 927 autoTableOfReal TableOfReal_bootstrap (TableOfReal me) { 928 try { 929 autoTableOfReal thee = TableOfReal_create (my numberOfRows, my numberOfColumns); 930 for (integer icol = 1; icol <= my numberOfColumns; icol ++) 931 TableOfReal_setColumnLabel (thee.get(), icol, my columnLabels [icol].get()); 932 /* 933 Select randomly from table with replacement. Because of replacement 934 you do not get back the original data set. A random fraction, 935 typically 1/e (37%), are replaced by duplicates. 936 */ 937 for (integer thyRow = 1; thyRow <= thy numberOfRows; thyRow ++) { 938 const integer myRow = NUMrandomInteger (1, my numberOfRows); 939 thy data.row (thyRow) <<= my data.row (myRow); 940 TableOfReal_setRowLabel (thee.get(), thyRow, my rowLabels [myRow].get()); 941 } 942 return thee; 943 } catch (MelderError) { 944 Melder_throw (me, U": bootstrapped data not created."); 945 } 946 } 947 948 void TableOfReal_changeRowLabels (TableOfReal me, 949 conststring32 search, conststring32 replace, integer maximumNumberOfReplaces, 950 integer *nmatches, integer *nstringmatches, bool use_regexp) { 951 try { 952 autoSTRVEC rowLabels = string32vector_searchAndReplace (my rowLabels.get(), 953 search, replace, maximumNumberOfReplaces, nmatches, nstringmatches, use_regexp); 954 my rowLabels = std::move (rowLabels); 955 } catch (MelderError) { 956 Melder_throw (me, U": row labels not changed."); 957 } 958 } 959 960 void TableOfReal_changeColumnLabels (TableOfReal me, 961 conststring32 search, conststring32 replace, integer maximumNumberOfReplaces, 962 integer *nmatches, integer *nstringmatches, bool use_regexp) { 963 try { 964 autoSTRVEC columnLabels = string32vector_searchAndReplace (my columnLabels.get(), 965 search, replace, maximumNumberOfReplaces, nmatches, nstringmatches, use_regexp); 966 my columnLabels = std::move (columnLabels); 967 } catch (MelderError) { 968 Melder_throw (me, U": column labels not changed."); 969 } 970 } 971 972 integer TableOfReal_getNumberOfLabelMatches (TableOfReal me, conststring32 search, bool columnLabels, bool use_regexp) { 973 integer nmatches = 0, numberOfLabels = my numberOfRows; 974 char32 **labels = my rowLabels.peek2(); 975 regexp *compiled_regexp = nullptr; 976 977 if (! search || search [0] == U'\0') 978 return 0; 979 if (columnLabels) { 980 numberOfLabels = my numberOfColumns; 981 labels = my columnLabels.peek2(); 982 } 983 try { 984 if (use_regexp) 985 compiled_regexp = CompileRE_throwable (search, 0); 986 for (integer i = 1; i <= numberOfLabels; i ++) { 987 if (! labels [i]) 988 continue; 989 if (use_regexp) { 990 if (ExecRE (compiled_regexp, nullptr, labels [i], nullptr, false, U'\0', U'\0', nullptr, nullptr)) 991 nmatches ++; 992 } else if (str32equ (labels [i], search)) 993 nmatches ++; 994 } 995 if (use_regexp) 996 free (compiled_regexp); 997 } catch (MelderError) { 998 if (use_regexp) 999 free (compiled_regexp); 1000 Melder_throw (me, U"Can't get number of matches."); 1001 } 1002 return nmatches; 1003 } 1004 1005 void TableOfReal_drawVectors (TableOfReal me, Graphics g, integer colx1, integer coly1, integer colx2, integer coly2, double xmin, double xmax, double ymin, double ymax, int vectype, int labelsize, bool garnish) { 1006 const integer nx = my numberOfColumns, ny = my numberOfRows; 1007 const double fontsize = Graphics_inqFontSize (g); 1008 1009 Melder_require (colx1 > 0 && colx1 <= nx && coly1 > 0 && coly1 <= nx, 1010 U"The index in the \"From\" column(s) should be in range [1, ", nx, U"]."); 1011 Melder_require (colx2 > 0 && colx2 <= nx && coly2 > 0 && coly2 <= nx, 1012 U"The index in the \"To\" column(s) should be in range [1, ", nx, U"]."); 1013 1014 if (xmin >= xmax) { 1015 MelderRealRange x1 = NUMextrema (my data.column (colx1)); 1016 MelderRealRange x2 = NUMextrema (my data.column (colx2)); 1017 xmin = std::min (x1.min, x2.min); 1018 xmax = std::max (x1.max, x2.max); 1019 } 1020 if (ymin >= ymax) { 1021 MelderRealRange y1 = NUMextrema (my data.column (coly1)); 1022 MelderRealRange y2 = NUMextrema (my data.column (coly2)); 1023 ymin = std::min (y1.min, y2.min); 1024 ymax = std::max (y1.max, y2.max); 1025 } 1026 if (xmin == xmax) { 1027 if (ymin == ymax) return; 1028 xmin -= 0.5; 1029 xmax += 0.5; 1030 } 1031 if (ymin == ymax) { 1032 ymin -= 0.5; 1033 ymax += 0.5; 1034 } 1035 1036 Graphics_setWindow (g, xmin, xmax, ymin, ymax); 1037 Graphics_setInner (g); 1038 Graphics_setTextAlignment (g, kGraphics_horizontalAlignment::CENTRE, Graphics_HALF); 1039 1040 if (labelsize > 0) 1041 Graphics_setFontSize (g, labelsize); 1042 1043 for (integer i = 1; i <= ny; i ++) { 1044 const double x1 = my data [i] [colx1]; 1045 const double y1 = my data [i] [coly1]; 1046 const double x2 = my data [i] [colx2]; 1047 const double y2 = my data [i] [coly2]; 1048 const conststring32 mark = EMPTY_STRING (my rowLabels [i]) ? U"" : my rowLabels [i].get(); 1049 if (vectype == Graphics_LINE) 1050 Graphics_line (g, x1, y1, x2, y2); 1051 else if (vectype == Graphics_TWOWAYARROW) { 1052 Graphics_arrow (g, x1, y1, x2, y2); 1053 Graphics_arrow (g, x2, y2, x1, y1); 1054 } else /*if (vectype == Graphics_ARROW) */ 1055 Graphics_arrow (g, x1, y1, x2, y2); 1056 if (labelsize > 0) 1057 Graphics_text (g, x1, y1, mark); 1058 } 1059 if (labelsize > 0) 1060 Graphics_setFontSize (g, fontsize); 1061 1062 Graphics_unsetInner (g); 1063 if (garnish) { 1064 Graphics_drawInnerBox (g); 1065 Graphics_marksLeft (g, 2, true, true, false); 1066 Graphics_marksBottom (g, 2, true, true, false); 1067 } 1068 } 1069 1070 void TableOfReal_drawColumnAsDistribution (TableOfReal me, Graphics g, integer column, double minimum, double maximum, integer nBins, double freqMin, double freqMax, bool cumulative, bool garnish) { 1071 if (column < 1 || column > my numberOfColumns) 1072 return; 1073 autoMatrix thee = TableOfReal_to_Matrix (me); 1074 Matrix_drawDistribution (thee.get(), g, column - 0.5, column + 0.5, 0.0, 0.0, minimum, maximum, nBins, freqMin, freqMax, cumulative, garnish); 1075 if (garnish && my columnLabels [column]) 1076 Graphics_textBottom (g, true, my columnLabels [column].get()); 1077 } 1078 1079 autoTableOfReal TableOfReal_sortRowsByIndex (TableOfReal me, constINTVEC index, bool reverse) { 1080 try { 1081 Melder_require (my rowLabels, 1082 U"No labels to sort"); 1083 MelderIntegerRange range = NUMextrema (index); 1084 Melder_require (range.first > 0 && range.first <= my numberOfRows && range.last > 0 && range.last <= my numberOfRows, 1085 U"One or more indices out of range [1, ", my numberOfRows, U"]."); 1086 autoTableOfReal thee = TableOfReal_create (my numberOfRows, my numberOfColumns); 1087 for (integer i = 1; i <= my numberOfRows; i ++) { 1088 const integer myindex = reverse ? i : index [i]; 1089 const integer thyindex = reverse ? index [i] : i; 1090 thy rowLabels [thyindex] = Melder_dup (my rowLabels [myindex].get()); 1091 thy data.row (thyindex) <<= my data.row (myindex); 1092 } 1093 thy columnLabels.all() <<= my columnLabels.all(); 1094 return thee; 1095 } catch (MelderError) { 1096 Melder_throw (me, U": not sorted by row index."); 1097 } 1098 } 1099 1100 autoINTVEC TableOfReal_getSortedIndexFromRowLabels (TableOfReal me) { 1101 try { 1102 return newINTVECindex (my rowLabels.get()); 1103 } catch (MelderError) { 1104 Melder_throw (me, U": no sorted index created."); 1105 } 1106 } 1107 1108 autoTableOfReal TableOfReal_sortOnlyByRowLabels (TableOfReal me) { 1109 try { 1110 autoPermutation index = TableOfReal_to_Permutation_sortRowLabels (me); 1111 autoTableOfReal thee = TableOfReal_Permutation_permuteRows (me, index.get()); 1112 return thee; 1113 } catch (MelderError) { 1114 Melder_throw (me, U": not sorted by row labels."); 1115 } 1116 } 1117 1118 static void NUMaverageBlock_byColumns_inplace (MAT a, integer rb, integer re, integer cb, integer ce, bool medians) { 1119 Melder_assert (rb > 0 && rb <= a.nrow); 1120 Melder_assert (rb <= re && re <= a.nrow); 1121 Melder_assert (cb > 0 && cb <= a.ncol); 1122 Melder_assert (cb <= ce && ce <= a.ncol); 1123 const integer n = re - rb + 1; 1124 if (n < 2) 1125 return; 1126 autoVEC tmp = raw_VEC (n); 1127 for (integer j = cb; j <= ce; j ++) { 1128 integer k = 1; 1129 for (integer i = rb; i <= re; i ++, k ++) 1130 tmp [k] = a [i] [j]; 1131 double average; 1132 if (medians) { 1133 sort_VEC_inout (tmp.get()); 1134 average = NUMquantile (tmp.get(), 0.5); 1135 } else { 1136 average = NUMmean (tmp.get()); 1137 } 1138 for (integer i = rb; i <= re; i ++) 1139 a [i] [j] = average; 1140 } 1141 } 1142 1143 autoTableOfReal TableOfReal_meansByRowLabels (TableOfReal me, bool expand, bool useMedians) { 1144 try { 1145 autoTableOfReal thee; 1146 autoINTVEC index = TableOfReal_getSortedIndexFromRowLabels (me); 1147 autoTableOfReal sorted = TableOfReal_sortRowsByIndex (me, index.get(), false); 1148 1149 integer indexi = 1, indexr = 0; 1150 conststring32 label = sorted -> rowLabels [1].get(); 1151 for (integer i = 2; i <= my numberOfRows; i ++) { 1152 const conststring32 li = sorted -> rowLabels [i].get(); 1153 if (Melder_cmp (li, label) != 0) { 1154 NUMaverageBlock_byColumns_inplace (sorted -> data.get(), indexi, i - 1, 1, my numberOfColumns, useMedians); 1155 if (! expand) { 1156 indexr ++; 1157 TableOfReal_copyOneRowWithLabel (sorted.get(), sorted.get(), indexi, indexr); 1158 } 1159 label = li; 1160 indexi = i; 1161 } 1162 } 1163 1164 NUMaverageBlock_byColumns_inplace (sorted -> data.get(), indexi, my numberOfRows, 1, my numberOfColumns, useMedians); 1165 1166 if (expand) { 1167 /* 1168 Now revert the sorting 1169 */ 1170 thee = TableOfReal_sortRowsByIndex (sorted.get(), index.get(), true); 1171 } else { 1172 indexr ++; 1173 TableOfReal_copyOneRowWithLabel (sorted.get(), sorted.get(), indexi, indexr); 1174 thee = TableOfReal_create (indexr, my numberOfColumns); 1175 for (integer i = 1; i <= indexr; i ++) 1176 TableOfReal_copyOneRowWithLabel (sorted.get(), thee.get(), i, i); 1177 thy columnLabels.all() <<= sorted -> columnLabels.all(); 1178 } 1179 return thee; 1180 } catch (MelderError) { 1181 Melder_throw (me, U": means by row labels not created."); 1182 } 1183 } 1184 1185 autoTableOfReal TableOfReal_rankColumns (TableOfReal me, integer fromColumn, integer toColumn) { 1186 try { 1187 fixAndCheckColumnRange (& fromColumn, & toColumn, my data.get(), 1); 1188 autoTableOfReal thee = Data_copy (me); 1189 MATrankColumns (thy data.get(), fromColumn, toColumn); 1190 return thee; 1191 } catch (MelderError) { 1192 Melder_throw (me, U": column ranks not created."); 1193 } 1194 } 1195 1196 /* 1197 s [lo] = precursor<number> 1198 s [lo+1] = precursor<number+1> 1199 ... 1200 s [hi] = precursor<number+hi-lo> 1201 */ 1202 void TableOfReal_setSequentialColumnLabels (TableOfReal me, integer from, integer to, conststring32 precursor, integer number, integer increment) { 1203 from = ( from == 0 ? 1 : from ); 1204 to = ( to == 0 ? my numberOfColumns : to ); 1205 Melder_require (from > 0 && from <= to && to <= my numberOfColumns, 1206 U"Wrong column indices."); 1207 for (integer i = from; i <= to; i ++, number += increment) 1208 my columnLabels [i] = Melder_dup (Melder_cat (precursor, number)); 1209 } 1210 1211 void TableOfReal_setSequentialRowLabels (TableOfReal me, integer from, integer to, conststring32 precursor, integer number, integer increment) { 1212 from = ( from == 0 ? 1 : from ); 1213 to = ( to == 0 ? my numberOfRows : to ); 1214 Melder_require (from > 0 && from <= to && to <= my numberOfRows, 1215 U"Wrong row indices."); 1216 for (integer i = from; i <= to; i ++, number += increment) 1217 my rowLabels [i] = Melder_dup (Melder_cat (precursor, number)); 1218 } 1219 1220 /* For the inheritors */ 1221 autoTableOfReal TableOfReal_to_TableOfReal (TableOfReal me) { 1222 try { 1223 autoTableOfReal thee = TableOfReal_create (my numberOfRows, my numberOfColumns); 1224 thy data.all() <<= my data.all(); 1225 TableOfReal_copyLabels (me, thee.get(), 1, 1); 1226 return thee; 1227 } catch (MelderError) { 1228 Melder_throw (me, U": not copied."); 1229 } 1230 } 1231 1232 autoTableOfReal TableOfReal_choleskyDecomposition (TableOfReal me, bool upper, bool inverse) { 1233 try { 1234 integer n = my numberOfColumns, lda = my numberOfRows, info; 1235 1236 Melder_require (n == lda, 1237 U"The table should be a square symmetric table."); 1238 1239 autoTableOfReal thee = Data_copy (me); 1240 1241 if (upper) { 1242 for (integer i = 2; i <= n; i ++) 1243 for (integer j = 1; j < i; j ++) 1244 thy data [i] [j] = 0.0; 1245 } else { 1246 for (integer i = 1; i < n; i ++) 1247 for (integer j = i + 1; j <= n; j ++) 1248 thy data [i] [j] = 0.0; 1249 } 1250 const char *uplo = ( upper ? "L" : "U" ); 1251 NUMlapack_dpotf2_ (uplo, n, & thy data [1] [1], lda, & info); 1252 Melder_require (info == 0, 1253 U"dpotf2 fails"); 1254 1255 if (inverse) { 1256 NUMlapack_dtrtri_ (uplo, "N", n, & thy data [1] [1], lda, &info); 1257 Melder_require (info == 0, 1258 U"dtrtri fails"); 1259 } 1260 return thee; 1261 } catch (MelderError) { 1262 Melder_throw (me, U": Cholesky decomposition not performed."); 1263 } 1264 } 1265 1266 autoTableOfReal TableOfReal_appendColumns (TableOfReal me, TableOfReal thee) { 1267 try { 1268 const integer ncols = my numberOfColumns + thy numberOfColumns; 1269 Melder_require (my numberOfRows == thy numberOfRows, 1270 U"The numbers of rows should be equal."); 1271 /* 1272 Stricter label checking??? 1273 append only if 1274 (my rowLabels [i] == thy rowlabels [i], i=1..my numberOfRows) or 1275 (my rowLabels [i] == 'empty', i=1..my numberOfRows) or 1276 (thy rowLabels [i] == 'empty', i=1..my numberOfRows); 1277 'empty': nullptr or \w* 1278 */ 1279 autoTableOfReal him = TableOfReal_create (my numberOfRows, ncols); 1280 his rowLabels.all() <<= my rowLabels.all(); 1281 his columnLabels.part (1, my numberOfColumns) <<= my columnLabels.all(); 1282 his columnLabels.part (my numberOfColumns + 1, ncols) <<= thy columnLabels.all(); 1283 his data.verticalBand (1, my numberOfColumns) <<= my data.all(); 1284 his data.verticalBand (my numberOfColumns + 1, ncols) <<= thy data.all(); 1285 1286 integer labeldiffs = 0; 1287 for (integer i = 1; i <= my numberOfRows; i ++) 1288 if (! Melder_equ (my rowLabels [i].get(), thy rowLabels [i].get())) 1289 labeldiffs ++; 1290 if (labeldiffs > 0) 1291 Melder_warning (labeldiffs, U" row labels differed."); 1292 return him; 1293 } catch (MelderError) { 1294 Melder_throw (U"TableOfReal with appended columns not created."); 1295 } 1296 } 1297 1298 autoTableOfReal TableOfRealList_appendColumnsMany (TableOfRealList me) { 1299 try { 1300 Melder_require (my size > 0, 1301 U"No tables selected."); 1302 const TableOfReal first = my at [1]; 1303 const integer nrows = first -> numberOfRows; 1304 integer ncols = first -> numberOfColumns; 1305 for (integer itab = 2; itab <= my size; itab ++) { 1306 const TableOfReal next = my at [itab]; 1307 ncols += next -> numberOfColumns; 1308 Melder_require (next -> numberOfRows == nrows, 1309 U"Numbers of rows in item ", itab, U" differs from previous."); 1310 } 1311 autoTableOfReal him = TableOfReal_create (nrows, ncols); 1312 /* Unsafe: new attributes not initialized. */ 1313 his rowLabels.all() <<= first -> rowLabels.all(); 1314 1315 integer hisColumnIndex = 0; 1316 for (integer itab = 1; itab <= my size; itab ++) { 1317 const TableOfReal table = my at [itab]; 1318 for (integer icol = 1; icol <= table -> numberOfColumns; icol ++) { 1319 hisColumnIndex ++; 1320 TableOfReal_setColumnLabel (him.get(), hisColumnIndex, table -> columnLabels [icol].get()); 1321 his data.column (hisColumnIndex) <<= table -> data.column (icol); 1322 } 1323 } 1324 Melder_assert (hisColumnIndex == his numberOfColumns); 1325 return him; 1326 } catch (MelderError) { 1327 Melder_throw (U"TableOfReal with appended columns not created."); 1328 } 1329 } 1330 1331 double TableOfReal_normalityTest_BHEP (TableOfReal me, double *h, double *out_tnb, double *out_lnmu, double *out_lnvar, bool *out_singularCovariance) { 1332 try { 1333 /* 1334 Henze & Wagner (1997), A new approach to the BHEP tests for multivariate normality, 1335 Journal of Multivariate Analysis 62, 1-23. 1336 */ 1337 autoCovariance thee = TableOfReal_to_Covariance (me); 1338 autoVEC weights; 1339 const double probability = Covariance_TableOfReal_normalityTest_BHEP (thee.get(), me, weights.get(), h, out_tnb, out_lnmu, out_lnvar, out_singularCovariance); 1340 return probability; 1341 } catch (MelderError) { 1342 Melder_throw (me, U": cannot determine normality."); 1343 } 1344 } 1345 1346 autoTableOfReal TableOfReal_TableOfReal_crossCorrelations (TableOfReal me, TableOfReal thee, bool by_columns, bool center, bool normalize) { 1347 return by_columns ? TableOfReal_TableOfReal_columnCorrelations (me, thee, center, normalize) : 1348 TableOfReal_TableOfReal_rowCorrelations (me, thee, center, normalize); 1349 } 1350 1351 autoTableOfReal TableOfReal_TableOfReal_rowCorrelations (TableOfReal me, TableOfReal thee, bool centre, bool normalize) { 1352 try { 1353 Melder_require (my numberOfColumns == thy numberOfColumns, 1354 U"Both tables should have the same number of columns."); 1355 autoTableOfReal him = TableOfReal_create (my numberOfRows, thy numberOfRows); 1356 autoMAT my_data = copy_MAT (my data.get()); 1357 autoMAT thy_data = copy_MAT (thy data.get()); 1358 if (centre) { 1359 centreEachRow_MAT_inout (my_data.get()); 1360 centreEachRow_MAT_inout (thy_data.get()); 1361 } 1362 if (normalize) { 1363 MATnormalizeRows_inplace (my_data.get(), 2.0, 1.0); 1364 MATnormalizeRows_inplace (thy_data.get(), 2.0, 1.0); 1365 } 1366 his rowLabels.all() <<= my rowLabels.all(); 1367 his columnLabels.all() <<= thy rowLabels.all(); 1368 mul_MAT_out (his data.get(), my_data.get(), thy_data.transpose()); 1369 return him; 1370 } catch (MelderError) { 1371 Melder_throw (U"TableOfReal with row correlations not created."); 1372 } 1373 } 1374 1375 autoTableOfReal TableOfReal_TableOfReal_columnCorrelations (TableOfReal me, TableOfReal thee, bool center, bool normalize) { 1376 try { 1377 Melder_require (my numberOfRows == thy numberOfRows, 1378 U"Both tables should have the same number of rows."); 1379 autoTableOfReal him = TableOfReal_create (my numberOfColumns, thy numberOfColumns); 1380 autoMAT my_data = copy_MAT (my data.get()); 1381 autoMAT thy_data = copy_MAT (thy data.get()); 1382 if (center) { 1383 centreEachColumn_MAT_inout (my_data.get()); 1384 centreEachColumn_MAT_inout (thy_data.get()); 1385 } 1386 if (normalize) { 1387 MATnormalizeColumns_inplace (my_data.get(), 2.0, 1.0); 1388 MATnormalizeColumns_inplace (thy_data.get(), 2.0, 1.0); 1389 } 1390 his rowLabels.all() <<= my columnLabels.all(); 1391 his columnLabels.all() <<= thy columnLabels.all(); 1392 mul_MAT_out (his data.get(), my_data.transpose(), thy_data.get()); 1393 return him; 1394 } catch (MelderError) { 1395 Melder_throw (U"TableOfReal with column correlations not created."); 1396 } 1397 } 1398 1399 autoMatrix TableOfReal_to_Matrix_interpolateOnRectangularGrid (TableOfReal me, double xmin, double xmax, double nx, double ymin, double ymax, integer ny, int /* method */) { 1400 try { 1401 Melder_require (my numberOfColumns > 2 && my numberOfColumns > 2, 1402 U"There should be at least three colums and three rows."); 1403 autoVEC weights = newVECbiharmonic2DSplineInterpolation_getWeights (my data.column (1), my data.column (2), my data.column (3)); 1404 const double dx = (xmax - xmin) / nx, dy = (ymax - ymin) / ny; 1405 autoMatrix thee = Matrix_create (xmin, xmax, nx, dx, xmin + 0.5 * dx, 1406 ymin, ymax, ny, dy, ymin + 0.5 * dy); 1407 for (integer irow = 1; irow <= ny; irow ++) { 1408 const double yp = thy y1 + (irow - 1) * dy; 1409 for (integer icol = 1; icol <= nx; icol ++) { 1410 const double xp = thy x1 + (icol - 1) * dx; 1411 thy z [irow] [icol] = NUMbiharmonic2DSplineInterpolation (my data.column (1), my data.column (2), weights.get(), xp, yp); 1412 } 1413 } 1414 return thee; 1415 } catch (MelderError) { 1416 Melder_throw (me, U": interpolation not finished."); 1417 } 1418 } 1419 1420 static void TableOfReal_permuteRowElements_inplace (TableOfReal me, integer fromRow, integer toRow) { 1421 if (fromRow == 0 && toRow == 0) { 1422 fromRow = 1; 1423 toRow = my numberOfRows; 1424 } 1425 Melder_require (fromRow <= toRow && fromRow >= 1 && toRow <= my numberOfRows, 1426 U"The row range should be larger than 0 and smaller than ", my numberOfRows + 1); 1427 autoPermutation permutation = Permutation_create (my numberOfColumns); 1428 autoVEC rowData = raw_VEC (my numberOfColumns); 1429 for (integer irow = fromRow; irow <= toRow; irow ++) { 1430 Permutation_permuteRandomly_inplace (permutation.get(), 0, 0); 1431 rowData.get() <<= my data.row (irow); 1432 for (integer icol = 1; icol <= my numberOfColumns; icol ++) 1433 my data [irow][icol] = rowData [Permutation_getValueAtIndex (permutation.get(), icol)]; 1434 } 1435 } 1436 1437 static void TableOfReal_permuteColumnElements_inplace (TableOfReal me, integer fromColumn, integer toColumn) { 1438 if (fromColumn == 0 && toColumn == 0) { 1439 fromColumn = 1; 1440 toColumn = my numberOfColumns; 1441 } 1442 Melder_require (fromColumn <= toColumn && fromColumn >= 1 && toColumn <= my numberOfColumns, 1443 U"The column range should be larger than 0 and smaller than ", my numberOfColumns + 1); 1444 autoPermutation permutation = Permutation_create (my numberOfRows); 1445 autoVEC columnData = raw_VEC (my numberOfRows); 1446 for (integer icol = fromColumn; icol <= toColumn; icol ++) { 1447 Permutation_permuteRandomly_inplace (permutation.get(), 0, 0); 1448 columnData.get() <<= my data.column (icol); 1449 for (integer irow = 1; irow <= my numberOfRows; irow ++) 1450 my data [irow] [icol] = columnData [Permutation_getValueAtIndex (permutation.get(), irow)]; 1451 } 1452 } 1453 1454 static void TableOfReal_shuffleCombinedRows (TableOfReal xShuffled, TableOfReal yShuffled, TableOfReal x, TableOfReal y) { 1455 Melder_require (x -> numberOfColumns == y -> numberOfColumns && 1456 xShuffled -> numberOfColumns == yShuffled -> numberOfColumns && 1457 x -> numberOfColumns == yShuffled -> numberOfColumns, 1458 U"All TableOfReals should have the same number of columns."); 1459 Melder_require (x -> numberOfRows == xShuffled -> numberOfRows && y -> numberOfRows == yShuffled -> numberOfRows, 1460 U"The number of rows of the first and the third TableOfReal should be equal and the number of rows of the " 1461 "second and forth TableOfReal should also be equal."); 1462 const integer numberOfRowsCombined = x -> numberOfRows + y -> numberOfRows; 1463 autoPermutation permutation = Permutation_create (numberOfRowsCombined); 1464 Permutation_permuteRandomly_inplace (permutation.get(), 0, 0); 1465 for (integer irow = 1; irow <= numberOfRowsCombined; irow ++) { 1466 const integer rowIndex = Permutation_getValueAtIndex (permutation.get(), irow); 1467 constVEC fromRow = ( rowIndex > x -> numberOfRows ? y -> data.row (rowIndex - x -> numberOfRows) : 1468 x -> data.row (rowIndex) ); 1469 VEC toRow = ( irow > xShuffled -> numberOfRows ? yShuffled -> data.row (irow - xShuffled -> numberOfRows) : 1470 xShuffled -> data.row (irow) ); 1471 toRow <<= fromRow; 1472 } 1473 } 1474 1475 static autoVEC vectorizeLowerMinusDiagonal_VEC (constMATVU m) { 1476 Melder_assert (m.nrow == m.ncol); 1477 integer index = 0; 1478 const integer size = m.nrow * (m.nrow - 1) / 2; 1479 autoVEC vech = raw_VEC (size); 1480 for (integer irow = 2; irow <= m.nrow; irow ++) 1481 for (integer icol = 1; icol < irow; icol ++) 1482 vech [++ index] = m [irow] [icol]; 1483 Melder_assert (index == size); 1484 return vech; 1485 } 1486 1487 static autoVEC vectorizeLowerPlusDiagonal_VEC (constMATVU m) { 1488 Melder_assert (m.nrow == m.ncol); 1489 integer index = 0; 1490 const integer size = m.nrow * (m.nrow + 1) / 2; 1491 autoVEC vech = raw_VEC (size); 1492 for (integer irow = 1; irow <= m.nrow; irow ++) 1493 for (integer icol = 1; icol <= irow; icol ++) 1494 vech [++ index] = m [irow] [icol]; 1495 Melder_assert (index == size); 1496 return vech; 1497 } 1498 1499 static double TableOfReal_computeTestStatistic_WuEtAl1215 (TableOfReal me, bool useCorrelation) { 1500 autoCovariance thee = TableOfReal_to_Covariance (me); 1501 double testStatistic; 1502 if (useCorrelation) { 1503 autoCorrelation him = SSCP_to_Correlation (thee.get()); 1504 autoVEC vech = vectorizeLowerMinusDiagonal_VEC (his data.get()); 1505 testStatistic = 1.0 - sqrt (my numberOfColumns) / NUMnorm (vech.get(), 2.0); 1506 } else { 1507 autoVEC vech = vectorizeLowerPlusDiagonal_VEC (thy data.get()); 1508 testStatistic = 1.0 - NUMtrace (thy data.get()) / (sqrt (my numberOfColumns) * NUMnorm (vech.get(), 2.0)); 1509 } 1510 return testStatistic; 1511 } 1512 1513 static double TableOfReal_computeTestStatistic_WuEtAl16 (TableOfReal me, bool useCorrelation) { 1514 autoCovariance thee = TableOfReal_to_Covariance (me); 1515 autoCorrelation him; 1516 constMATVU data = thy data.get(); 1517 if (useCorrelation) { 1518 him = SSCP_to_Correlation (thee.get()); 1519 data = his data.get(); 1520 } 1521 autoVEC vech = vectorizeLowerMinusDiagonal_VEC (data); 1522 const double testStatistic = 1.0 - NUMsum (vech.get()) / (NUMnorm (vech.get(), 2.0) * sqrt (vech.size)); 1523 return testStatistic; 1524 } 1525 1526 static double TableOfReal_computeTestStatistic_WuEtAl17 (TableOfReal me, TableOfReal thee, bool useCorrelation) { 1527 autoCovariance mycov = TableOfReal_to_Covariance (me); 1528 autoCovariance thycov = TableOfReal_to_Covariance (thee); 1529 autoVEC myvech,thyvech; 1530 if (useCorrelation) { 1531 autoCorrelation mycor = SSCP_to_Correlation (mycov.get()); 1532 myvech = vectorizeLowerMinusDiagonal_VEC (mycor -> data.get()); 1533 autoCorrelation thycor = SSCP_to_Correlation (thycov.get()); 1534 thyvech = vectorizeLowerMinusDiagonal_VEC (thycor -> data.get()); 1535 } else { 1536 myvech = vectorizeLowerPlusDiagonal_VEC (mycov -> data.get()); 1537 thyvech = vectorizeLowerPlusDiagonal_VEC (thycov -> data.get()); 1538 } 1539 const double testStatistic = 1.0 - NUMinner (myvech.get(), thyvech.get()) / (NUMnorm (myvech.get(), 2.0) * NUMnorm (thyvech.get(), 2.0)); 1540 return testStatistic; 1541 } 1542 1543 double TableOfReal_testSphericityOfCovariance (TableOfReal me, integer numberOfPermutations, bool useCorrelation) { 1544 try { 1545 autoTableOfReal thee = Data_copy (me); 1546 double testStatistic0 = TableOfReal_computeTestStatistic_WuEtAl1215 (me, useCorrelation); 1547 integer countLargerOrEqual = 0; 1548 for (integer iperm = 1; iperm <= numberOfPermutations; iperm ++) { 1549 TableOfReal_permuteRowElements_inplace (thee.get(), 0, 0); 1550 TableOfReal_permuteColumnElements_inplace (thee.get(), 0, 0); 1551 const double testStatistic = TableOfReal_computeTestStatistic_WuEtAl1215 (thee.get(), useCorrelation); 1552 if (testStatistic >= testStatistic0) 1553 countLargerOrEqual ++; 1554 } 1555 return (1.0 + countLargerOrEqual) / (1.0 + numberOfPermutations); 1556 } catch (MelderError) { 1557 Melder_throw (me, U": could not determine probability for sphericity."); 1558 } 1559 } 1560 1561 double TableOfReal_testCovarianceEqualsIdentityMatrix (TableOfReal me, integer numberOfPermutations, bool useCorrelation) { 1562 try { 1563 autoTableOfReal mycopy = Data_copy (me); 1564 double testStatistic0 = TableOfReal_computeTestStatistic_WuEtAl1215 (me, useCorrelation); 1565 integer countLargerOrEqual = 0; 1566 for (integer iperm = 1; iperm <= numberOfPermutations; iperm ++) { 1567 TableOfReal_permuteColumnElements_inplace (mycopy.get(), 0, 0); 1568 const double testStatistic = TableOfReal_computeTestStatistic_WuEtAl1215 (mycopy.get(), useCorrelation); 1569 if (testStatistic >= testStatistic0) 1570 countLargerOrEqual ++; 1571 } 1572 return (1.0 + countLargerOrEqual) / (1.0 + numberOfPermutations); 1573 } catch (MelderError) { 1574 Melder_throw (me, U": cannot determine probability for identity."); 1575 } 1576 } 1577 1578 double TableOfReal_testCovarianceCompoundSymmetry (TableOfReal me, integer numberOfPermutations, bool useCorrelation) { 1579 try { 1580 autoTableOfReal mycopy = Data_copy (me); 1581 double testStatistic0 = TableOfReal_computeTestStatistic_WuEtAl16 (me, useCorrelation); 1582 integer countLargerOrEqual = 0; 1583 for (integer iperm = 1; iperm <= numberOfPermutations; iperm ++) { 1584 TableOfReal_permuteRowElements_inplace (mycopy.get(), 0, 0); 1585 const double testStatistic = TableOfReal_computeTestStatistic_WuEtAl16 (mycopy.get(), useCorrelation); 1586 if (testStatistic >= testStatistic0) 1587 countLargerOrEqual ++; 1588 } 1589 return (1.0 + countLargerOrEqual) / (1.0 + numberOfPermutations); 1590 } catch (MelderError) { 1591 Melder_throw (me, U": cannot determine compound symmetry."); 1592 } 1593 } 1594 1595 double TableOfReal_testEqualityOfCovariances (TableOfReal me, TableOfReal thee, integer numberOfPermutations, bool useCorrelation) { 1596 try { 1597 Melder_require (my numberOfColumns == thy numberOfColumns, 1598 U"The number of columns of both TableOfReal should be equal."); 1599 autoTableOfReal myCopy = Data_copy (me); 1600 autoTableOfReal thyCopy = Data_copy (thee); 1601 double testStatistic0 = TableOfReal_computeTestStatistic_WuEtAl17 (me, thee, useCorrelation); 1602 integer countLargerOrEqual = 0; 1603 for (integer iperm = 1; iperm <= numberOfPermutations; iperm ++) { 1604 double testStatistic; 1605 if (iperm % 2 == 1) { 1606 TableOfReal_shuffleCombinedRows (myCopy.get(), thyCopy.get(), me, thee); 1607 testStatistic = TableOfReal_computeTestStatistic_WuEtAl17 (myCopy.get(), thyCopy.get(), useCorrelation); 1608 } else { 1609 TableOfReal_shuffleCombinedRows (me, thee, myCopy.get(), thyCopy.get()); 1610 testStatistic = TableOfReal_computeTestStatistic_WuEtAl17 (me, thee, useCorrelation); 1611 } 1612 if (testStatistic >= testStatistic0) 1613 countLargerOrEqual ++; 1614 } 1615 return (1.0 + countLargerOrEqual) / (1.0 + numberOfPermutations); 1616 } catch (MelderError) { 1617 Melder_throw (me, U" & ", thee, U": could not determine probability for equality."); 1618 } 1619 } 1620 1621 #undef EMPTY_STRING 1622 1623 /* End of file TableOfReal_extensions.c */ 1624