1 /* Transition.cpp 2 * 3 * Copyright (C) 1997-2012,2015-2020 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 "Transition.h" 20 #include "NUM2.h" 21 #include "Eigen.h" 22 23 #include "oo_DESTROY.h" 24 #include "Transition_def.h" 25 #include "oo_COPY.h" 26 #include "Transition_def.h" 27 #include "oo_EQUAL.h" 28 #include "Transition_def.h" 29 #include "oo_CAN_WRITE_AS_ENCODING.h" 30 #include "Transition_def.h" 31 #include "oo_WRITE_BINARY.h" 32 #include "Transition_def.h" 33 #include "oo_READ_TEXT.h" 34 #include "Transition_def.h" 35 #include "oo_READ_BINARY.h" 36 #include "Transition_def.h" 37 #include "oo_DESCRIPTION.h" 38 #include "Transition_def.h" 39 40 Thing_implement (Transition, Daata, 0); 41 42 void structTransition :: v_info () { 43 structDaata :: v_info (); 44 MelderInfo_writeLine (U"Number of states: ", numberOfStates); 45 } 46 47 void structTransition :: v_writeText (MelderFile file) { 48 texputi32 (file, numberOfStates, U"numberOfStates"); 49 MelderFile_write (file, U"\nstateLabels []: "); 50 if (numberOfStates < 1) MelderFile_write (file, U"(empty)"); 51 MelderFile_write (file, U"\n"); 52 for (integer i = 1; i <= numberOfStates; i ++) { 53 MelderFile_write (file, U"\""); 54 if (stateLabels [i]) MelderFile_write (file, stateLabels [i].get()); 55 MelderFile_write (file, U"\"\t"); 56 } 57 for (integer i = 1; i <= numberOfStates; i ++) { 58 MelderFile_write (file, U"\nstate [", i, U"]:"); 59 for (integer j = 1; j <= numberOfStates; j ++) { 60 MelderFile_write (file, U"\t", data [i] [j]); 61 } 62 } 63 } 64 65 void Transition_init (Transition me, integer numberOfStates) { 66 if (numberOfStates < 1) 67 Melder_throw (U"Cannot create empty matrix."); 68 my numberOfStates = numberOfStates; 69 my stateLabels = autoSTRVEC (numberOfStates); 70 my data = zero_MAT (my numberOfStates, my numberOfStates); 71 } 72 73 autoTransition Transition_create (integer numberOfStates) { 74 try { 75 autoTransition me = Thing_new (Transition); 76 Transition_init (me.get(), numberOfStates); 77 return me; 78 } catch (MelderError) { 79 Melder_throw (U"Transition not created."); 80 } 81 } 82 83 static void NUMrationalize (double x, integer *numerator, integer *denominator) { 84 double epsilon = 1e-6; 85 *numerator = 1; 86 for (*denominator = 1; *denominator <= 100000; (*denominator) ++) { 87 double numerator_d = x * *denominator, rounded = round (numerator_d); 88 if (fabs (rounded - numerator_d) < epsilon) { 89 *numerator = (integer) rounded; 90 return; 91 } 92 } 93 *denominator = 0; // failure 94 } 95 96 static void print4 (char *buffer, double value, int iformat, int width, int precision) { 97 char formatString [40]; 98 if (iformat == 4) { 99 integer numerator, denominator; 100 NUMrationalize (value, & numerator, & denominator); 101 if (numerator == 0) 102 snprintf (buffer, 40, "0"); 103 else if (denominator > 1) 104 snprintf (buffer, 40, "%s/%s", Melder8_integer (numerator), Melder8_integer (denominator)); 105 else 106 snprintf (buffer, 40, "%.7g", value); 107 } else { 108 snprintf (formatString, 40, "%%%d.%d%c", width, precision, iformat == 1 ? 'f' : iformat == 2 ? 'e' : 'g'); 109 snprintf (buffer, 40, formatString, value); 110 } 111 } 112 113 void Transition_drawAsNumbers (Transition me, Graphics g, int iformat, int precision) { 114 double maxTextWidth = 0.0, maxTextHeight = 0.0; 115 Graphics_setInner (g); 116 Graphics_setWindow (g, 0.5, my numberOfStates + 0.5, 0.0, 1.0); 117 const double leftMargin = Graphics_dxMMtoWC (g, 1.0); 118 const double lineSpacing = Graphics_dyMMtoWC (g, 1.5 * Graphics_inqFontSize (g) * 25.4 / 72.0); 119 Graphics_setTextAlignment (g, Graphics_CENTRE, Graphics_BOTTOM); 120 for (integer col = 1; col <= my numberOfStates; col ++) { 121 if (my stateLabels && my stateLabels [col] && my stateLabels [col] [0]) { 122 Graphics_text (g, col, 1, my stateLabels [col].get()); 123 if (maxTextHeight == 0.0) maxTextHeight = lineSpacing; 124 } 125 } 126 for (integer row = 1; row <= my numberOfStates; row ++) { 127 double y = 1 - lineSpacing * (row - 1 + 0.7); 128 Graphics_setTextAlignment (g, Graphics_RIGHT, Graphics_HALF); 129 if (my stateLabels && my stateLabels [row]) { 130 double textWidth = Graphics_textWidth (g, my stateLabels [row].get()); 131 if (textWidth > maxTextWidth) maxTextWidth = textWidth; 132 Graphics_text (g, 0.5 - leftMargin, y, my stateLabels [row].get()); 133 } 134 Graphics_setTextAlignment (g, Graphics_CENTRE, Graphics_HALF); 135 for (integer col = 1; col <= my numberOfStates; col ++) { 136 char text [40]; 137 print4 (text, my data [row] [col], iformat, 0, precision); 138 Graphics_text (g, col, y, Melder_peek8to32 (text)); 139 } 140 } 141 if (maxTextWidth != 0.0) 142 Graphics_line (g, 0.5 - maxTextWidth - leftMargin, 1, my numberOfStates + 0.5, 1); 143 if (maxTextHeight != 0.0) 144 Graphics_line (g, 0.5, 1 + maxTextHeight, 0.5, 1 - lineSpacing * (my numberOfStates + 0.2)); 145 Graphics_unsetInner (g); 146 } 147 148 static void Transition_transpose (Transition me) { 149 for (integer i = 1; i < my numberOfStates; i ++) { 150 for (integer j = i + 1; j <= my numberOfStates; j ++) { 151 double temp = my data [i] [j]; 152 my data [i] [j] = my data [j] [i]; 153 my data [j] [i] = temp; 154 } 155 } 156 } 157 158 void Transition_eigen (Transition me, autoMatrix *out_eigenvectors, autoMatrix *out_eigenvalues) { 159 bool transposed = false; 160 try { 161 autoEigen eigen = Thing_new (Eigen); 162 Transition_transpose (me); 163 Eigen_initFromSymmetricMatrix (eigen.get(), my data.get()); 164 Transition_transpose (me); 165 transposed = true; 166 autoMatrix eigenvectors = Matrix_createSimple (my numberOfStates, my numberOfStates); 167 autoMatrix eigenvalues = Matrix_createSimple (my numberOfStates, 1); 168 for (integer i = 1; i <= my numberOfStates; i ++) { 169 eigenvalues -> z [i] [1] = eigen -> eigenvalues [i]; 170 for (integer j = 1; j <= my numberOfStates; j ++) 171 eigenvectors -> z [i] [j] = eigen -> eigenvectors [j] [i]; 172 } 173 *out_eigenvectors = eigenvectors.move(); 174 *out_eigenvalues = eigenvalues.move(); 175 } catch (MelderError) { 176 if (transposed) 177 Transition_transpose (me); 178 Melder_throw (me, U": eigenvectors not computed."); 179 } 180 } 181 182 autoTransition Transition_power (Transition me, integer power) { 183 try { 184 autoTransition thee = Data_copy (me); 185 autoTransition him = Data_copy (me); 186 for (integer ipow = 2; ipow <= power; ipow ++) { 187 std::swap (his data.cells, thy data.cells); // OPTIMIZE 188 for (integer irow = 1; irow <= my numberOfStates; irow ++) { 189 for (integer icol = 1; icol <= my numberOfStates; icol ++) { 190 thy data [irow] [icol] = 0.0; 191 for (integer i = 1; i <= my numberOfStates; i ++) { 192 thy data [irow] [icol] += his data [irow] [i] * my data [i] [icol]; 193 } 194 } 195 } 196 } 197 return thee; 198 } catch (MelderError) { 199 Melder_throw (me, U": power not computed."); 200 } 201 } 202 203 autoMatrix Transition_to_Matrix (Transition me) { 204 try { 205 autoMatrix thee = Matrix_createSimple (my numberOfStates, my numberOfStates); 206 for (integer i = 1; i <= my numberOfStates; i ++) 207 for (integer j = 1; j <= my numberOfStates; j ++) 208 thy z [i] [j] = my data [i] [j]; 209 return thee; 210 } catch (MelderError) { 211 Melder_throw (me, U": not converted to Matrix."); 212 } 213 } 214 215 autoTransition Matrix_to_Transition (Matrix me) { 216 try { 217 if (my nx != my ny) 218 Melder_throw (U"Matrix should be square."); 219 autoTransition thee = Transition_create (my nx); 220 for (integer i = 1; i <= my nx; i ++) 221 for (integer j = 1; j <= my nx; j ++) 222 thy data [i] [j] = my z [i] [j]; 223 return thee; 224 } catch (MelderError) { 225 Melder_throw (me, U": not converted to Transition."); 226 } 227 } 228 229 /* End of file Transition.cpp */ 230