1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126112711281129113011311132113311341135113611371138113911401141114211431144114511461147114811491150115111521153115411551156115711581159116011611162116311641165116611671168116911701171117211731174117511761177117811791180118111821183118411851186118711881189119011911192119311941195119611971198119912001201120212031204120512061207120812091210121112121213121412151216121712181219122012211222122312241225122612271228122912301231123212331234123512361237123812391240124112421243124412451246124712481249125012511252125312541255125612571258125912601261126212631264126512661267126812691270127112721273127412751276127712781279128012811282128312841285128612871288128912901291129212931294129512961297129812991300130113021303130413051306130713081309131013111312131313141315131613171318131913201321132213231324132513261327132813291330133113321333133413351336133713381339134013411342134313441345134613471348134913501351135213531354135513561357135813591360136113621363136413651366136713681369137013711372137313741375137613771378137913801381138213831384138513861387138813891390139113921393139413951396139713981399140014011402140314041405140614071408140914101411141214131414141514161417141814191420142114221423142414251426142714281429143014311432143314341435143614371438143914401441144214431444144514461447144814491450145114521453145414551456145714581459146014611462146314641465146614671468146914701471147214731474147514761477147814791480148114821483148414851486148714881489149014911492149314941495149614971498149915001501150215031504150515061507150815091510151115121513151415151516151715181519152015211522152315241525152615271528152915301531153215331534153515361537153815391540154115421543154415451546154715481549155015511552155315541555155615571558155915601561156215631564156515661567156815691570157115721573157415751576157715781579158015811582158315841585158615871588158915901591159215931594159515961597159815991600160116021603160416051606160716081609161016111612161316141615161616171618161916201621162216231624162516261627162816291630163116321633163416351636163716381639 |
- /* HMM.cpp
- *
- * Copyright (C) 2010-2017 David Weenink, 2015,2017 Paul Boersma
- *
- * This code is free software; you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation; either version 2 of the License, or (at
- * your option) any later version.
- *
- * This code is distributed in the hope that it will be useful, but
- * WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
- * General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this work. If not, see <http://www.gnu.org/licenses/>.
- */
- /*
- djmw 20110304 Thing_new
- */
- #include "Distributions_and_Strings.h"
- #include "HMM.h"
- #include "Index.h"
- #include "NUM2.h"
- #include "Strings_extensions.h"
- #include "oo_DESTROY.h"
- #include "HMM_def.h"
- #include "oo_COPY.h"
- #include "HMM_def.h"
- #include "oo_EQUAL.h"
- #include "HMM_def.h"
- #include "oo_CAN_WRITE_AS_ENCODING.h"
- #include "HMM_def.h"
- #include "oo_WRITE_TEXT.h"
- #include "HMM_def.h"
- #include "oo_WRITE_BINARY.h"
- #include "HMM_def.h"
- #include "oo_READ_TEXT.h"
- #include "HMM_def.h"
- #include "oo_READ_BINARY.h"
- #include "HMM_def.h"
- #include "oo_DESCRIPTION.h"
- #include "HMM_def.h"
- Thing_implement (HMM, Daata, 0);
- Thing_implement (HMMState, Daata, 0);
- Thing_implement (HMMStateList, Ordered, 0);
- Thing_implement (HMMObservation, Daata, 0);
- Thing_implement (HMMObservationList, Ordered, 0);
- Thing_implement (HMMBaumWelch, Daata, 0);
- Thing_implement (HMMViterbi, Daata, 0);
- Thing_implement (HMMObservationSequence, Table, 0);
- Thing_implement (HMMObservationSequenceBag, Collection, 0);
- Thing_implement (HMMStateSequence, Strings, 0);
- /*
- Whenever a routine returns ln(p), the result for p=0 is -INFINITY.
- On IEEE floating point hardware this number behaves reasonably.
- This means that when the variable q equals INFINITY, q + a -> INFINITY,
- where a is a finite number.
- */
- // helpers
- int NUMget_line_intersection_with_circle (double xc, double yc, double r, double a, double b, double *x1, double *y1, double *x2, double *y2);
- autoHMMObservation HMMObservation_create (conststring32 label, integer numberOfComponents, integer dimension, integer storage);
- integer HMM_HMMObservationSequence_getLongestSequence (HMM me, HMMObservationSequence thee, integer symbolNumber);
- integer StringsIndex_getLongestSequence (StringsIndex me, integer index, integer *pos);
- integer Strings_getLongestSequence (Strings me, char32 *string, integer *pos);
- autoHMMState HMMState_create (conststring32 label);
- autoHMMBaumWelch HMMBaumWelch_create (integer nstates, integer nsymbols, integer capacity);
- void HMMBaumWelch_getGamma (HMMBaumWelch me);
- autoHMMBaumWelch HMM_forward (HMM me, integer *obs, integer nt);
- void HMMBaumWelch_reInit (HMMBaumWelch me);
- void HMM_HMMBaumWelch_getXi (HMM me, HMMBaumWelch thee, integer *obs);
- void HMM_HMMBaumWelch_reestimate (HMM me, HMMBaumWelch thee);
- void HMM_HMMBaumWelch_addEstimate (HMM me, HMMBaumWelch thee, integer *obs);
- void HMM_HMMBaumWelch_forward (HMM me, HMMBaumWelch thee, integer *obs);
- void HMM_HMMBaumWelch_backward (HMM me, HMMBaumWelch thee, integer *obs);
- void HMM_HMMViterbi_decode (HMM me, HMMViterbi thee, integer *obs);
- double HMM_getProbabilityOfObservations (HMM me, integer *obs, integer numberOfTimes);
- autoTableOfReal StringsIndex_to_TableOfReal_transitions (StringsIndex me, int probabilities);
- autoStringsIndex HMM_HMMStateSequence_to_StringsIndex (HMM me, HMMStateSequence thee);
- autoHMMViterbi HMMViterbi_create (integer nstates, integer ntimes);
- autoHMMViterbi HMM_to_HMMViterbi (HMM me, integer *obs, integer ntimes);
- // evaluate the numbers given to probabilities
- static autoVEC NUMwstring_to_probs (conststring32 s, integer nwanted) {
- autoVEC numbers = VEC_createFromString (s);
- if (numbers.size != nwanted)
- Melder_throw (U"You supplied ", numbers.size, U", while ", nwanted, U" numbers needed.");
- longdouble sum = 0.0;
- for (integer i = 1; i <= numbers.size; i ++) {
- if (numbers [i] < 0.0)
- Melder_throw (U"Numbers have to be positive.");
- sum += numbers [i];
- }
- if (sum <= 0.0)
- Melder_throw (U"All probabilities cannot be zero.");
- for (integer i = 1; i <= numbers.size; i ++)
- numbers [i] /= sum;
- return numbers;
- }
- int NUMget_line_intersection_with_circle (double xc, double yc, double r, double a, double b, double *x1, double *y1, double *x2, double *y2) {
- double ca = a * a + 1.0, bmyc = (b - yc);
- double cb = 2.0 * (a * bmyc - xc);
- double cc = bmyc * bmyc + xc * xc - r * r;
- integer nroots = NUMsolveQuadraticEquation (ca, cb, cc, x1, x2);
- if (nroots == 1) {
- *y1 = a * *x1 + b;
- *x2 = *x1;
- *y2 = *y1;
- } else if (nroots == 2) {
- if (*x1 > *x2) {
- double tmp = *x1;
- *x1 = *x2;
- *x2 = tmp;
- }
- *y1 = *x1 * a + b;
- *y2 = *x2 * a + b;
- }
- return nroots;
- }
- // D(l_1,l_2)=1/n( log p(O_2|l_1) - log p(O_2|l_2)
- static double HMM_HMM_getCrossEntropy_asym (HMM me, HMM thee, integer observationLength) {
- autoHMMObservationSequence os = HMM_to_HMMObservationSequence (thee, 0, observationLength);
- double ce = HMM_HMMObservationSequence_getCrossEntropy (me, os.get());
- if (isundef (ce)) {
- return ce;
- }
- double ce2 = HMM_HMMObservationSequence_getCrossEntropy (thee, os.get());
- if (isundef (ce2)) {
- return ce2;
- }
- return ce - ce2;
- }
- /**************** HMMObservation ******************************/
- static void HMMObservation_init (HMMObservation me, conststring32 label, integer numberOfComponents, integer dimension, integer storage) {
- my label = Melder_dup (label);
- my gm = GaussianMixture_create (numberOfComponents, dimension, storage);
- }
- autoHMMObservation HMMObservation_create (conststring32 label, integer numberOfComponents, integer dimension, integer storage) {
- try {
- autoHMMObservation me = Thing_new (HMMObservation);
- HMMObservation_init (me.get(), label, numberOfComponents, dimension, storage);
- return me;
- } catch (MelderError) {
- Melder_throw (U"HMMObservation not created.");
- }
- }
- integer Strings_getLongestSequence (Strings me, char32 *string, integer *pos) {
- integer length = 0, longest = 0, lpos = 0;
- for (integer i = 1; i <= my numberOfStrings; i ++) {
- if (Melder_equ (my strings [i].get(), string)) {
- if (length == 0) {
- lpos = i;
- }
- length ++;
- } else {
- if (length > 0) {
- if (length > longest) {
- longest = length; *pos = lpos;
- }
- length = 0;
- }
- }
- }
- return length;
- }
- integer StringsIndex_getLongestSequence (StringsIndex me, integer index, integer *pos) {
- integer length = 0, longest = 0, lpos = 0;
- for (integer i = 1; i <= my numberOfItems; i ++) {
- if (my classIndex [i] == index) {
- if (length == 0)
- lpos = i;
- length ++;
- } else {
- if (length > 0) {
- if (length > longest) {
- longest = length;
- *pos = lpos;
- }
- length = 0;
- }
- }
- }
- return length;
- }
- /**************** HMMState ******************************/
- static void HMMState_init (HMMState me, conststring32 label) {
- my label = Melder_dup (label);
- }
- autoHMMState HMMState_create (conststring32 label) {
- try {
- autoHMMState me = Thing_new (HMMState);
- HMMState_init (me.get(), label);
- return me;
- } catch (MelderError) {
- Melder_throw (U"HMMState not created.");
- }
- }
- void HMMState_setLabel (HMMState me, char32 *label) {
- my label = Melder_dup (label);
- }
- /**************** HMMBaumWelch ******************************/
- void structHMMBaumWelch :: v_destroy () noexcept {
- for (integer it = 1; it <= capacity; it ++) {
- NUMmatrix_free (xi [it], 1, 1);
- }
- NUMvector_free (xi, 1);
- NUMvector_free (scale, 1);
- NUMmatrix_free (beta, 1, 1);
- NUMmatrix_free (alpha, 1, 1);
- NUMmatrix_free (gamma, 1, 1);
- NUMmatrix_free (aij_num, 0, 1);
- NUMmatrix_free (aij_denom, 0, 1);
- NUMmatrix_free (bik_num, 1, 1);
- NUMmatrix_free (bik_denom, 1, 1);
- }
- autoHMMBaumWelch HMMBaumWelch_create (integer nstates, integer nsymbols, integer capacity) {
- try {
- autoHMMBaumWelch me = Thing_new (HMMBaumWelch);
- my numberOfTimes = my capacity = capacity;
- my numberOfStates = nstates;
- my numberOfSymbols = nsymbols;
- my alpha = NUMmatrix<double> (1, nstates, 1, capacity);
- my beta = NUMmatrix<double> (1, nstates, 1, capacity);
- my scale = NUMvector<double> (1, capacity);
- my xi = NUMvector<double **> (1, capacity);
- my aij_num = NUMmatrix<double> (0, nstates, 1, nstates + 1);
- my aij_denom = NUMmatrix<double> (0, nstates, 1, nstates + 1);
- my bik_num = NUMmatrix<double> (1, nstates, 1, nsymbols);
- my bik_denom = NUMmatrix<double> (1, nstates, 1, nsymbols);
- my gamma = NUMmatrix<double> (1, nstates, 1, capacity);
- for (integer it = 1; it <= capacity; it ++)
- my xi [it] = NUMmatrix<double> (1, nstates, 1, nstates);
- return me;
- } catch (MelderError) {
- Melder_throw (U"HMMBaumWelch not created.");
- }
- }
- void HMMBaumWelch_getGamma (HMMBaumWelch me) {
- for (integer it = 1; it <= my numberOfTimes; it ++) {
- double sum = 0.0;
- for (integer is = 1; is <= my numberOfStates; is ++) {
- my gamma [is] [it] = my alpha [is] [it] * my beta [is] [it];
- sum += my gamma [is] [it];
- }
- for (integer is = 1; is <= my numberOfStates; is ++)
- my gamma [is] [it] /= sum;
- }
- }
- /**************** HMMViterbi ******************************/
- autoHMMViterbi HMMViterbi_create (integer nstates, integer ntimes) {
- try {
- autoHMMViterbi me = Thing_new (HMMViterbi);
- my numberOfTimes = ntimes;
- my numberOfStates = nstates;
- my viterbi = NUMmatrix <double> (1, nstates, 1, ntimes);
- my bp = NUMmatrix <integer> (1, nstates, 1 , ntimes);
- my path = NUMvector <integer> (1, ntimes);
- return me;
- } catch (MelderError) {
- Melder_throw (U"HMMViterbi not created.");
- }
- }
- /******************* HMMObservationSequence & HMMStateSequence ***/
- autoHMMObservationSequence HMMObservationSequence_create (integer numberOfItems, integer dataLength) {
- try {
- autoHMMObservationSequence me = Thing_new (HMMObservationSequence);
- Table_initWithoutColumnNames (me.get(), numberOfItems, dataLength + 1);
- return me;
- } catch (MelderError) {
- Melder_throw (U"HMMObservationSequence not created.");
- }
- }
- integer HMMObservationSequence_getNumberOfObservations (HMMObservationSequence me) {
- return my rows.size;
- }
- void HMMObservationSequence_removeObservation (HMMObservationSequence me, integer index) {
- Table_removeRow (me, index);
- }
- autoStrings HMMObservationSequence_to_Strings (HMMObservationSequence me) {
- try {
- integer numberOfStrings = my rows.size;
- autoStrings thee = Thing_new (Strings);
- thy strings = autostring32vector (numberOfStrings);
- for (integer i = 1; i <= numberOfStrings; i ++)
- thy strings [i] = Melder_dup (Table_getStringValue_Assert (me, i, 1));
- thy numberOfStrings = numberOfStrings;
- return thee;
- } catch (MelderError) {
- Melder_throw (me, U": no Strings created.");
- }
- }
- autoHMMObservationSequence Strings_to_HMMObservationSequence (Strings me) {
- try {
- autoHMMObservationSequence thee = HMMObservationSequence_create (my numberOfStrings, 0);
- for (integer i = 1; i <= my numberOfStrings; i ++)
- Table_setStringValue (thee.get(), i, 1, my strings [i].get());
- return thee;
- } catch (MelderError) {
- Melder_throw (me, U": no HMMObservationSequence created.");
- }
- }
- autoStringsIndex HMMObservationSequence_to_StringsIndex (HMMObservationSequence me) {
- try {
- autoStrings s = HMMObservationSequence_to_Strings (me);
- autoStringsIndex thee = Strings_to_StringsIndex (s.get());
- return thee;
- } catch (MelderError) {
- Melder_throw (me, U": no StringsIndex created.");
- }
- }
- integer HMM_HMMObservationSequence_getLongestSequence (HMM me, HMMObservationSequence thee, integer symbolNumber) {
- autoStringsIndex si = HMM_HMMObservationSequence_to_StringsIndex (me, thee);
- // TODO
- (void) symbolNumber;
- return 1;
- }
- integer HMMObservationSequenceBag_getLongestSequence (HMMObservationSequenceBag me) {
- integer longest = 0;
- for (integer i = 1; i <= my size; i ++) {
- HMMObservationSequence thee = my at [i];
- if (thy rows.size > longest)
- longest = thy rows.size;
- }
- return longest;
- }
- autoHMMStateSequence HMMStateSequence_create (integer numberOfItems) {
- try {
- autoHMMStateSequence me = Thing_new (HMMStateSequence);
- my strings = autostring32vector (numberOfItems);
- return me;
- } catch (MelderError) {
- Melder_throw (U"HMMStateSequence not created.");
- }
- }
- autoStrings HMMStateSequence_to_Strings (HMMStateSequence me) {
- try {
- autoStrings thee = Thing_new (Strings);
- my structStrings :: v_copy (thee.get());
- return thee;
- } catch (MelderError) {
- Melder_throw (me, U": no Strings created.");
- }
- }
- /**************** HMM ******************************/
- void structHMM :: v_info () {
- structDaata :: v_info ();
- MelderInfo_writeLine (U"Number of states: ", numberOfStates);
- for (integer i = 1; i <= numberOfStates; i ++) {
- HMMState hmms = our states->at [i];
- MelderInfo_writeLine (U" ", hmms -> label.get());
- }
- MelderInfo_writeLine (U"Number of symbols: ", numberOfObservationSymbols);
- for (integer i = 1; i <= numberOfObservationSymbols; i ++) {
- HMMObservation hmms = our observationSymbols->at [i];
- MelderInfo_writeLine (U" ", hmms -> label.get());
- }
- }
- static void HMM_init (HMM me, integer numberOfStates, integer numberOfObservationSymbols, int leftToRight) {
- my numberOfStates = numberOfStates;
- my numberOfObservationSymbols = numberOfObservationSymbols;
- my componentStorage = 1;
- my leftToRight = leftToRight;
- my states = HMMStateList_create ();
- my observationSymbols = HMMObservationList_create ();
- my transitionProbs = NUMmatrix<double> (0, numberOfStates, 1, numberOfStates + 1);
- my emissionProbs = NUMmatrix<double> (1, numberOfStates, 1, numberOfObservationSymbols);
- }
- autoHMM HMM_create (int leftToRight, integer numberOfStates, integer numberOfObservationSymbols) {
- try {
- autoHMM me = Thing_new (HMM);
- HMM_init (me.get(), numberOfStates, numberOfObservationSymbols, leftToRight);
- HMM_setDefaultStates (me.get());
- HMM_setDefaultObservations (me.get());
- HMM_setDefaultTransitionProbs (me.get());
- HMM_setDefaultStartProbs (me.get());
- HMM_setDefaultEmissionProbs (me.get());
- return me;
- } catch (MelderError) {
- Melder_throw (U"HMM not created.");
- }
- }
- void HMM_setDefaultStates (HMM me) {
- for (integer i = 1; i <= my numberOfStates; i ++) {
- autoHMMState hmms = HMMState_create (Melder_cat (U"S", i));
- HMM_addState_move (me, hmms.move());
- }
- }
- autoHMM HMM_createFullContinuousModel (int leftToRight, integer numberOfStates, integer numberOfObservationSymbols, integer numberOfFeatureStreams, integer *dimensionOfStream, integer *numberOfGaussiansforStream) {
- (void) leftToRight;
- (void) numberOfStates;
- (void) numberOfObservationSymbols;
- (void) numberOfFeatureStreams;
- (void) dimensionOfStream;
- (void) numberOfGaussiansforStream;
- return autoHMM();
- }
- autoHMM HMM_createContinuousModel (int leftToRight, integer numberOfStates, integer numberOfObservationSymbols, integer numberOfMixtureComponentsPerSymbol, integer componentDimension, integer componentStorage) {
- try {
- autoHMM me = Thing_new (HMM);
- HMM_init (me.get(), numberOfStates, numberOfObservationSymbols, leftToRight);
- my numberOfMixtureComponents = numberOfMixtureComponentsPerSymbol;
- my componentDimension = componentDimension;
- my componentStorage = componentStorage;
- for (integer i = 1; i <= numberOfStates; i ++) {
- autoHMMState state = HMMState_create (Melder_cat (U"S", i));
- HMM_addState_move (me.get(), state.move());
- }
- for (integer j = 1; j <= numberOfObservationSymbols; j ++) {
- autoHMMObservation obs = HMMObservation_create (Melder_cat (U"s", j), numberOfMixtureComponentsPerSymbol, componentDimension, componentStorage);
- HMM_addObservation_move (me.get(), obs.move());
- }
- HMM_setDefaultTransitionProbs (me.get());
- HMM_setDefaultStartProbs (me.get());
- HMM_setDefaultEmissionProbs (me.get());
- HMM_setDefaultMixingProbabilities (me.get());
- return me;
- } catch (MelderError) {
- Melder_throw (U"Continuous model HMM not created.");
- }
- }
- // for a simple non-hidden model leave either states empty or symbols empty !!!
- autoHMM HMM_createSimple (int leftToRight, conststring32 states_string, conststring32 symbols_string) {
- try {
- autostring32vector states = STRVECtokenize (states_string);
- autostring32vector symbols = STRVECtokenize (symbols_string);
- autoHMM me = Thing_new (HMM);
- Melder_require (states.size > 0 || symbols.size > 0,
- U"The states and symbols should not both be empty.");
-
- if (symbols.size <= 0) {
- symbols = STRVECclone (states.get());
- my notHidden = 1;
- }
- if (states.size <= 0) {
- states = STRVECclone (symbols.get());
- my notHidden = 1;
- }
- HMM_init (me.get(), states.size, symbols.size, leftToRight);
- for (integer istate = 1; istate <= states.size; istate ++) {
- autoHMMState state = HMMState_create (states [istate].get());
- HMM_addState_move (me.get(), state.move());
- }
- for (integer isymbol = 1; isymbol <= symbols.size; isymbol ++) {
- autoHMMObservation symbol = HMMObservation_create (symbols [isymbol].get(), 0, 0, 0);
- HMM_addObservation_move (me.get(), symbol.move());
- }
- HMM_setDefaultTransitionProbs (me.get());
- HMM_setDefaultStartProbs (me.get());
- HMM_setDefaultEmissionProbs (me.get());
- return me;
- } catch (MelderError) {
- Melder_throw (U"Simple HMM not created.");
- }
- }
- void HMM_setDefaultObservations (HMM me) {
- conststring32 def = my notHidden ? U"S" : U"s";
- for (integer i = 1; i <= my numberOfObservationSymbols; i ++) {
- autoHMMObservation hmms = HMMObservation_create (Melder_cat (def, i), 0, 0, 0);
- HMM_addObservation_move (me, hmms.move());
- }
- }
- void HMM_setDefaultTransitionProbs (HMM me) {
- for (integer i = 1; i <= my numberOfStates; i ++) {
- double p = ( my leftToRight ? 1.0 / (my numberOfStates - i + 1.0) : 1.0 / my numberOfStates );
- for (integer j = 1; j <= my numberOfStates; j ++)
- my transitionProbs [i] [j] = ( my leftToRight && j < i ? 0.0 : p );
- }
- // leftToRight must have end state!
- if (my leftToRight) my transitionProbs [my numberOfStates] [my numberOfStates] =
- my transitionProbs [my numberOfStates] [my numberOfStates + 1] = 0.5;
- }
- void HMM_setDefaultStartProbs (HMM me) {
- double p = 1.0 / my numberOfStates;
- for (integer j = 1; j <= my numberOfStates; j ++)
- my transitionProbs [0] [j] = p;
- }
- void HMM_setDefaultEmissionProbs (HMM me) {
- double p = 1.0 / my numberOfObservationSymbols;
- for (integer i = 1; i <= my numberOfStates; i ++) {
- for (integer j = 1; j <= my numberOfObservationSymbols; j ++)
- my emissionProbs [i] [j] = ( my notHidden ? ( i == j ? 1.0 : 0.0 ) : p );
- }
- }
- void HMM_setDefaultMixingProbabilities (HMM me) {
- double mp = 1.0 / my numberOfMixtureComponents;
- for (integer is = 1; is <= my numberOfObservationSymbols; is ++) {
- HMMObservation hmmo = my observationSymbols->at [is];
- for (integer im = 1; im <= my numberOfMixtureComponents; im ++)
- hmmo -> gm -> mixingProbabilities [im] = mp;
- }
- }
- void HMM_setStartProbabilities (HMM me, conststring32 probs) {
- try {
- autoVEC p = NUMwstring_to_probs (probs, my numberOfStates);
- for (integer i = 1; i <= my numberOfStates; i ++)
- my transitionProbs [0] [i] = p [i];
- } catch (MelderError) {
- Melder_throw (me, U": no start probabilities set.");
- }
- }
- void HMM_setTransitionProbabilities (HMM me, integer state_number, conststring32 state_probs) {
- try {
- Melder_require (state_number <= my states->size,
- U"State number should not exceed ", my states->size, U".");
- autoVEC p = NUMwstring_to_probs (state_probs, my numberOfStates);
- for (integer i = 1; i <= my numberOfStates + 1; i ++)
- my transitionProbs [state_number] [i] = p [i];
- } catch (MelderError) {
- Melder_throw (me, U": no transition probabilities set.");
- }
- }
- void HMM_setEmissionProbabilities (HMM me, integer state_number, conststring32 emission_probs) {
- try {
- Melder_require (state_number <= my states->size,
- U"State number should not exceed ", my states->size, U".");
- Melder_require (! my notHidden,
- U"The emission probabilities of this model are fixed.");
- autoVEC p = NUMwstring_to_probs (emission_probs, my numberOfObservationSymbols);
- for (integer i = 1; i <= my numberOfObservationSymbols; i ++) {
- my emissionProbs [state_number] [i] = p [i];
- }
- } catch (MelderError) {
- Melder_throw (me, U": no emission probabilities set.");
- }
- }
- void HMM_addObservation_move (HMM me, autoHMMObservation thee) {
- integer ns = my observationSymbols->size + 1;
- Melder_require (ns <= my numberOfObservationSymbols, U"Observation list is full.");
- my observationSymbols -> addItemAtPosition_move (thee.move(), ns);
- }
- void HMM_addState_move (HMM me, autoHMMState thee) {
- integer ns = my states->size + 1;
- Melder_require (ns <= my numberOfStates, U"States list is full.");
- my states -> addItemAtPosition_move (thee.move(), ns);
- }
- autoTableOfReal HMM_extractTransitionProbabilities (HMM me) {
- try {
- autoTableOfReal thee = TableOfReal_create (my numberOfStates + 1, my numberOfStates + 1);
- for (integer is = 1; is <= my numberOfStates; is ++) {
- HMMState hmms = my states->at [is];
- TableOfReal_setRowLabel (thee.get(), is + 1, hmms -> label.get());
- TableOfReal_setColumnLabel (thee.get(), is, hmms -> label.get());
- for (integer js = 1; js <= my numberOfStates; js ++)
- thy data [is + 1] [js] = my transitionProbs [is] [js];
- }
- TableOfReal_setRowLabel (thee.get(), 1, U"START");
- TableOfReal_setColumnLabel (thee.get(), my numberOfStates + 1, U"END");
- for (integer is = 1; is <= my numberOfStates; is ++) {
- thy data [1] [is] = my transitionProbs [0] [is];
- thy data [is + 1] [my numberOfStates + 1] = my transitionProbs [is] [my numberOfStates + 1];
- }
- return thee;
- } catch (MelderError) {
- Melder_throw (me, U": no transition probabilities extracted.");
- }
- }
- autoTableOfReal HMM_extractEmissionProbabilities (HMM me) {
- try {
- autoTableOfReal thee = TableOfReal_create (my numberOfStates, my numberOfObservationSymbols);
- for (integer js = 1; js <= my numberOfObservationSymbols; js ++) {
- HMMObservation hmms = my observationSymbols->at [js];
- TableOfReal_setColumnLabel (thee.get(), js, hmms -> label.get());
- }
- for (integer is = 1; is <= my numberOfStates; is ++) {
- HMMState hmms = my states->at [is];
- TableOfReal_setRowLabel (thee.get(), is, hmms -> label.get());
- for (integer js = 1; js <= my numberOfObservationSymbols; js ++)
- thy data [is] [js] = my emissionProbs [is] [js];
- }
- return thee;
- } catch (MelderError) {
- Melder_throw (me, U": no emission probabilities extracted.");
- };
- }
- double HMM_getExpectedValueOfDurationInState (HMM me, integer istate) {
- if (istate < 0 || istate > my numberOfStates)
- return undefined;
- return 1.0 / (1.0 - my transitionProbs [istate] [istate]);
- }
- double HMM_getProbabilityOfStayingInState (HMM me, integer istate, integer numberOfTimeUnits) {
- if (istate < 0 || istate > my numberOfStates)
- return undefined;
- return pow (my transitionProbs [istate] [istate], numberOfTimeUnits - 1.0) * (1.0 - my transitionProbs [istate] [istate]);
- }
- double HMM_HMM_getCrossEntropy (HMM me, HMM thee, integer observationLength, int symmetric) {
- double ce1 = HMM_HMM_getCrossEntropy_asym (me, thee, observationLength);
- if (! symmetric || isundef (ce1))
- return ce1;
- double ce2 = HMM_HMM_getCrossEntropy_asym (thee, me, observationLength);
- if (isundef (ce2))
- return ce2;
- return (ce1 + ce2) / 2.0;
- }
- double HMM_HMM_HMMObservationSequence_getCrossEntropy (HMM me, HMM thee, HMMObservationSequence him) {
- double ce1 = HMM_HMMObservationSequence_getCrossEntropy (me, him);
- if (isundef (ce1))
- return ce1;
- double ce2 = HMM_HMMObservationSequence_getCrossEntropy (thee, him);
- if (isundef (ce2))
- return ce2;
- return (ce1 + ce2) / 2.0;
- }
- void HMM_draw (HMM me, Graphics g, int garnish) {
- double xwidth = sqrt (my numberOfStates);
- double rstate = 0.3 / xwidth, r = xwidth / 3.0;
- double xmax = 1.2 * xwidth / 2.0, xmin = -xmax, ymin = xmin, ymax = xmax;
- autoNUMvector<double> xs (1, my numberOfStates);
- autoNUMvector<double> ys (1, my numberOfStates);
- Graphics_setWindow (g, xmin, xmax, ymin, ymax);
- // heuristic: all states on a circle until we have a better graph drawing algorithm.
- xs [1] = ys [1] = 0;
- if (my numberOfStates > 1) {
- for (integer is = 1; is <= my numberOfStates; is ++) {
- double alpha = - NUMpi + NUMpi * 2.0 * (is - 1) / my numberOfStates;
- xs [is] = r * cos (alpha); ys [is] = r * sin (alpha);
- }
- }
- // reorder the positions such that state number 1 is most left and last state number is right.
- // if > 5 may be one in the middle with the most connections
- // ...
- // find fontsize
- int fontSize = Graphics_inqFontSize (g);
- conststring32 widest_label = U"";
- double max_width = 0.0;
- for (integer is = 1; is <= my numberOfStates; is ++) {
- HMMState hmms = my states->at [is];
- double w = ( hmms -> label ? Graphics_textWidth (g, hmms -> label.get()) : 0.0 );
- if (w > max_width) {
- widest_label = hmms -> label.get();
- max_width = w;
- }
- }
- int new_fontSize = fontSize;
- while (max_width > 2.0 * rstate && new_fontSize > 4) {
- new_fontSize --;
- Graphics_setFontSize (g, new_fontSize);
- max_width = Graphics_textWidth (g, widest_label);
- }
- Graphics_setFontSize (g, new_fontSize);
- Graphics_setTextAlignment (g, Graphics_CENTRE, Graphics_HALF);
- for (integer is = 1; is <= my numberOfStates; is ++) {
- HMMState hmms = my states->at [is];
- Graphics_circle (g, xs [is], ys [is], rstate);
- Graphics_text (g, xs [is], ys [is], hmms -> label.get());
- }
- // draw connections from is to js
- // 1 -> 2 / 2-> : increase / decrease angle between 1 and 2 with pi /10
- // use cos(a+b) and cos(a -b) rules
- double cosb = cos (NUMpi / 10.0), sinb = sin (NUMpi / 10.0);
- for (integer is = 1; is <= my numberOfStates; is ++) {
- double x1 = xs [is], y1 = ys [is];
- for (integer js = 1; js <= my numberOfStates; js ++) {
- if (my transitionProbs [is] [js] > 0.0 && is != js) {
- double x2 = xs [js], y2 = ys [js];
- double dx = x2 - x1, dy = y2 - y1, h = sqrt (dx * dx + dy * dy), cosa = dx / h, sina = dy / h;
- double cosabp = cosa * cosb - sina * sinb, cosabm = cosa * cosb + sina * sinb;
- double sinabp = cosa * sinb + sina * cosb, sinabm = -cosa * sinb + sina * cosb;
- Graphics_arrow (g, x1 + rstate * cosabp, y1 + rstate * sinabp, x2 - rstate * cosabm, y2 - rstate * sinabm);
- }
- if (is == js) {
- double dx = - x1, dy = - y1, h = sqrt (dx * dx + dy * dy), cosa = dx / h, sina = dy / h;
- Graphics_doubleArrow (g, x1 - rstate * cosa, y1 - rstate * sina, x1 - 1.4 * rstate * cosa, y1 - 1.4 * rstate * sina);
- }
- }
- }
- if (garnish) {
- Graphics_drawInnerBox (g);
- }
- }
- void HMM_unExpandPCA (HMM me) {
- if (my componentDimension <= 0)
- return; // nothing to do
- for (integer is = 1; is <= my numberOfObservationSymbols; is ++) {
- HMMObservation s = my observationSymbols->at [is];
- GaussianMixture_unExpandPCA (s -> gm.get());
- }
- }
- autoHMMObservationSequence HMM_to_HMMObservationSequence (HMM me, integer startState, integer numberOfItems) {
- try {
- autoHMMObservationSequence thee = HMMObservationSequence_create (numberOfItems, my componentDimension);
- integer istate = startState == 0 ? NUMgetIndexFromProbability (my transitionProbs [0], my numberOfStates, NUMrandomUniform (0.0, 1.0)) : startState;
- if (my componentDimension > 0) {
- autoVEC obs (my componentDimension, kTensorInitializationType::RAW);
- autoVEC buf (my componentDimension, kTensorInitializationType::RAW);
- for (integer i = 1; i <= numberOfItems; i ++) {
- // Emit a symbol from istate
- integer isymbol = NUMgetIndexFromProbability (my emissionProbs [istate], my numberOfObservationSymbols, NUMrandomUniform (0.0, 1.0));
- HMMObservation s = my observationSymbols->at [isymbol];
- char32 *name;
- GaussianMixture_generateOneVector_inline (s -> gm.get(), obs.get(), &name, buf.get());
- for (integer j = 1; j <= my componentDimension; j ++)
- Table_setNumericValue (thee.get(), i, 1 + j, obs [j]);
- Table_setStringValue (thee.get(), i, 1, s -> label.get());
- // get next state
- istate = NUMgetIndexFromProbability (my transitionProbs [istate], my numberOfStates + 1, NUMrandomUniform (0.0, 1.0));
- if (istate == my numberOfStates + 1) { // final state
- for (integer j = numberOfItems; j > i; j --)
- HMMObservationSequence_removeObservation (thee.get(), j);
- break;
- }
- }
-
- } else {
- for (integer i = 1; i <= numberOfItems; i ++) {
- // Emit a symbol from istate
- integer isymbol = NUMgetIndexFromProbability (my emissionProbs [istate], my numberOfObservationSymbols, NUMrandomUniform (0.0, 1.0));
- HMMObservation s = my observationSymbols->at [isymbol];
- Table_setStringValue (thee.get(), i, 1, s -> label.get());
- // get next state
- istate = NUMgetIndexFromProbability (my transitionProbs [istate], my numberOfStates + 1, NUMrandomUniform (0.0, 1.0));
- if (istate == my numberOfStates + 1) { // final state
- for (integer j = numberOfItems; j > i; j --)
- HMMObservationSequence_removeObservation (thee.get(), j);
- break;
- }
- }
- }
- HMM_unExpandPCA (me);
- return thee;
- } catch (MelderError) {
- HMM_unExpandPCA (me);
- Melder_throw (me, U":no HMMObservationSequence created.");
- }
- }
- autoHMMBaumWelch HMM_forward (HMM me, integer *obs, integer nt) {
- try {
- autoHMMBaumWelch thee = HMMBaumWelch_create (my numberOfStates, my numberOfObservationSymbols, nt);
- HMM_HMMBaumWelch_forward (me, thee.get(), obs);
- return thee;
- } catch (MelderError) {
- Melder_throw (me, U": no HMMBaumWelch created.");
- }
- }
- autoHMMViterbi HMM_to_HMMViterbi (HMM me, integer *obs, integer ntimes) {
- try {
- autoHMMViterbi thee = HMMViterbi_create (my numberOfStates, ntimes);
- HMM_HMMViterbi_decode (me, thee.get(), obs);
- return thee;
- } catch (MelderError) {
- Melder_throw (me, U": no HMMViterbi created.");
- }
- }
- void HMMBaumWelch_reInit (HMMBaumWelch me) {
- my totalNumberOfSequences = 0;
- my lnProb = 0.0;
- /*
- The _num and _denum matrices are asigned as += in the iteration loop and therefore need to be zeroed
- at the start of each new iteration.
- The elements of alpha, beta, scale, gamma & xi are always calculated directly and need not be
- initialised.
- */
- for (integer is = 0; is <= my numberOfStates; is ++) {
- for (integer js = 1; js <= my numberOfStates + 1; js ++) {
- my aij_num [is] [js] = 0.0;
- my aij_denom [is] [js] = 0.0;
- }
- }
- for (integer is = 1; is <= my numberOfStates; is ++) {
- for (integer js = 1; js <= my numberOfSymbols; js ++) {
- my bik_num [is] [js] = 0.0;
- my bik_denom [is] [js] = 0.0;
- }
- }
- }
- void HMM_HMMObservationSequenceBag_learn (HMM me, HMMObservationSequenceBag thee, double delta_lnp, double minProb, int info) {
- try {
- // act as if all observation sequences are in memory
- integer capacity = HMMObservationSequenceBag_getLongestSequence (thee);
- autoHMMBaumWelch bw = HMMBaumWelch_create (my numberOfStates, my numberOfObservationSymbols, capacity);
- bw -> minProb = minProb;
- if (info) {
- MelderInfo_open ();
- }
- integer iter = 0;
- double lnp;
- do {
- lnp = bw -> lnProb;
- HMMBaumWelch_reInit (bw.get());
- for (integer ios = 1; ios <= thy size; ios ++) {
- HMMObservationSequence hmm_os = thy at [ios];
- autoStringsIndex si = HMM_HMMObservationSequence_to_StringsIndex (me, hmm_os); // TODO outside the loop or more efficiently
- integer *obs = si -> classIndex.at;
- integer nobs = si -> numberOfItems; // convenience
- // Interpretation of unknowns: end of sequence
- integer istart = 1, iend = nobs;
- while (istart <= nobs) {
- while (istart <= nobs && obs [istart] == 0) {
- istart ++;
- };
- if (istart > nobs) {
- break;
- }
- iend = istart + 1;
- while (iend <= nobs && obs [iend] != 0) {
- iend ++;
- }
- iend --;
- bw -> numberOfTimes = iend - istart + 1;
- bw -> totalNumberOfSequences ++;
- HMM_HMMBaumWelch_forward (me, bw.get(), obs + istart - 1); // get new alphas
- HMM_HMMBaumWelch_backward (me, bw.get(), obs + istart - 1); // get new betas
- HMMBaumWelch_getGamma (bw.get());
- HMM_HMMBaumWelch_getXi (me, bw.get(), obs + istart - 1);
- HMM_HMMBaumWelch_addEstimate (me, bw.get(), obs + istart - 1);
- istart = iend + 1;
- }
- }
- // we have processed all observation sequences, now it is time to estimate new probabilities.
- iter ++;
- HMM_HMMBaumWelch_reestimate (me, bw.get());
- if (info) {
- MelderInfo_writeLine (U"Iteration: ", iter, U" ln(prob): ", bw -> lnProb);
- }
- } while (fabs ((lnp - bw -> lnProb) / bw -> lnProb) > delta_lnp);
- if (info) {
- MelderInfo_writeLine (U"******** Learning summary *********");
- MelderInfo_writeLine (U" Processed ", thy size, U" sequences,");
- MelderInfo_writeLine (U" consisting of ", bw -> totalNumberOfSequences, U" observation sequences.");
- MelderInfo_writeLine (U" Longest observation sequence had ", capacity, U" items");
- MelderInfo_close();
- }
- } catch (MelderError) {
- Melder_throw (me, U" & ", thee, U": not learned.");
- }
- }
- // xc1 < xc2
- void HMM_HMMStateSequence_drawTrellis (HMM me, HMMStateSequence thee, Graphics g, int connect, int garnish) {
- integer numberOfTimes = thy numberOfStrings;
- autoStringsIndex si = HMM_HMMStateSequence_to_StringsIndex (me, thee);
- double xmin = 0.0, xmax = numberOfTimes + 1.0, ymin = 0.5, ymax = my numberOfStates + 0.5;
- Graphics_setInner (g);
- Graphics_setWindow (g, xmin, xmax, ymin, ymax);
- double r = 0.2 / (numberOfTimes > my numberOfStates ? numberOfTimes : my numberOfStates);
- for (integer it = 1; it <= numberOfTimes; it ++) {
- for (integer js = 1; js <= my numberOfStates; js ++) {
- double xc = it, yc = js, x2 = it, y2 = js;
- Graphics_circle (g, xc, yc, r);
- if (it > 1) {
- for (integer is = 1; is <= my numberOfStates; is ++) {
- bool indexedConnection = si -> classIndex [it - 1] == is && si -> classIndex [it] == js;
- Graphics_setLineWidth (g, indexedConnection ? 2.0 : 1.0);
- Graphics_setLineType (g, indexedConnection ? Graphics_DRAWN : Graphics_DOTTED);
- double x1 = it - 1, y1 = is;
- if (connect || indexedConnection) {
- double a = (y1 - y2) / (x1 - x2), b = y1 - a * x1;
- // double xs11 = x1 - r / (a * a + 1), ys11 = a * xs11 + b;
- double xs12 = x1 + r / (a * a + 1), ys12 = a * xs12 + b;
- double xs21 = x2 - r / (a * a + 1), ys21 = a * xs21 + b;
- // double xs22 = x2 + r / (a * a + 1), ys22 = a * xs22 + b;
- Graphics_line (g, xs12, ys12, xs21, ys21);
- }
- }
- }
- }
- }
- Graphics_unsetInner (g);
- Graphics_setLineWidth (g, 1.0);
- Graphics_setLineType (g, Graphics_DRAWN);
- if (garnish) {
- Graphics_drawInnerBox (g);
- for (integer js = 1; js <= my numberOfStates; js ++) {
- HMMState hmms = my states->at [js];
- Graphics_markLeft (g, js, false, false, false, hmms -> label.get());
- }
- Graphics_marksBottomEvery (g, 1.0, 1.0, true, true, false);
- Graphics_textBottom (g, true, U"Time index");
- }
- }
- void HMM_drawBackwardProbabilitiesIllustration (Graphics g, bool garnish) {
- double xmin = 0.0, xmax = 1.0, ymin = 0.0, ymax = 1.0;
- Graphics_setWindow (g, xmin, xmax, ymin, ymax);
- double xleft = 0.1, xright = 0.9, r = 0.03;
- integer np = 6;
- double dy = (1.0 - 0.3) / (np - 1);
- double x0 = xleft, y0 = 0.5;
- Graphics_setWindow (g, xmin, xmax, ymin, ymax);
- Graphics_circle (g, x0, y0, r);
- double x = xright, y = 0.9;
- for (integer i = 1; i <= np; i ++) {
- if (i < 4 or i == np) {
- Graphics_circle (g, x, y, r);
- double xx = x0 - x, yy = y - y0;
- double c = sqrt (xx * xx + yy * yy);
- double cosa = xx / c, sina = yy / c;
- Graphics_line (g, x0 - r * cosa, y0 + r * sina, x + r * cosa, y - r * sina);
- } else if (i == 4) {
- double ddy = 3*dy/4;
- Graphics_fillCircle (g, x, y + dy - ddy, 0.5 * r);
- Graphics_fillCircle (g, x, y + dy - 2 * ddy, 0.5 * r);
- Graphics_fillCircle (g, x, y + dy - 3 * ddy, 0.5 * r);
- }
- y -= dy;
- }
- if (garnish) {
- double x1 = xright + 1.5 * r, x2 = x1 - 0.2, y1 = 0.9;
- Graphics_setTextAlignment (g, Graphics_LEFT, Graphics_HALF);
- Graphics_text (g, x1, y1, U"%s__1_");
- Graphics_setTextAlignment (g, Graphics_CENTRE, Graphics_HALF);
- Graphics_text (g, x2, y1, U"%a__%i1_");
- y1 = 0.9 - dy;
- Graphics_setTextAlignment (g, Graphics_LEFT, Graphics_HALF);
- Graphics_text (g, x1, y1, U"%s__2_");
- Graphics_setTextAlignment (g, Graphics_CENTRE, Graphics_HALF);
- Graphics_text (g, x2, y1, U"%a__%i2_");
- y1 = 0.9 - (np - 1) * dy;
- Graphics_setTextAlignment (g, Graphics_LEFT, Graphics_HALF);
- Graphics_text (g, x1, y1, U"%s__%N_");
- Graphics_setTextAlignment (g, Graphics_CENTRE, Graphics_HALF);
- Graphics_text (g, x2, y1, U"%a__%%iN%_");
- Graphics_setTextAlignment (g, Graphics_RIGHT, Graphics_HALF);
- Graphics_text (g, x0 - 1.5 * r, y0, U"%s__%i_");
- Graphics_setTextAlignment (g, Graphics_CENTRE, Graphics_BOTTOM);
- Graphics_text (g, x0, 0.0, U"%t");
- Graphics_text (g, x, 0.0, U"%t+1");
- double y3 = 0.10;
- Graphics_setTextAlignment (g, Graphics_CENTRE, Graphics_HALF);
- Graphics_text (g, x0, y3, U"%\\be__%t_(%i)%");
- Graphics_text (g, x, y3, U"%\\be__%t+1_(%j)");
- }
- }
- void HMM_drawForwardProbabilitiesIllustration (Graphics g, bool garnish) {
- double xmin = 0.0, xmax = 1.0, ymin = 0.0, ymax = 1.0;
- Graphics_setWindow (g, xmin, xmax, ymin, ymax);
- double xleft = 0.1, xright = 0.9, r = 0.03;
- integer np = 6;
- double dy = (1.0 - 0.3) / (np - 1);
- double x0 = xright, y0 = 0.5;
- Graphics_setWindow (g, xmin, xmax, ymin, ymax);
- Graphics_circle (g, x0, y0, r);
- double x = xleft, y = 0.9;
- for (integer i = 1; i <= np; i ++) {
- if (i < 4 or i == np) {
- Graphics_circle (g, x, y, r);
- double xx = x0 - x, yy = y - y0;
- double c = sqrt (xx * xx + yy * yy);
- double cosa = xx / c, sina = yy / c;
- Graphics_line (g, x0 - r * cosa, y0 + r * sina, x + r * cosa, y - r * sina);
- } else if (i == 4) {
- double ddy = 3.0 * dy / 4.0;
- Graphics_fillCircle (g, x, y + dy - ddy, 0.5 * r);
- Graphics_fillCircle (g, x, y + dy - 2 * ddy, 0.5 * r);
- Graphics_fillCircle (g, x, y + dy - 3 * ddy, 0.5 * r);
- }
- y -= dy;
- }
- if (garnish) {
- double x1 = xleft - 1.5 * r, x2 = x1 + 0.2, y1 = 0.9;
- Graphics_setTextAlignment (g, Graphics_RIGHT, Graphics_HALF);
- Graphics_text (g, x1, y1, U"%s__1_");
- Graphics_setTextAlignment (g, Graphics_CENTRE, Graphics_HALF);
- Graphics_text (g, x2, y1, U"%a__1%j_");
- y1 = 0.9 - dy;
- Graphics_setTextAlignment (g, Graphics_RIGHT, Graphics_HALF);
- Graphics_text (g, x1, y1, U"%s__2_");
- Graphics_setTextAlignment (g, Graphics_CENTRE, Graphics_HALF);
- Graphics_text (g, x2, y1, U"%a__2%j_");
- y1 = 0.9 - (np - 1) * dy;
- Graphics_setTextAlignment (g, Graphics_RIGHT, Graphics_HALF);
- Graphics_text (g, x1, y1, U"%s__%N_");
- Graphics_setTextAlignment (g, Graphics_CENTRE, Graphics_HALF);
- Graphics_text (g, x2, y1, U"%a__%%Nj%_");
- Graphics_setTextAlignment (g, Graphics_LEFT, Graphics_HALF);
- Graphics_text (g, x0 + 1.5 * r, y0, U"%s__%j_");
- Graphics_setTextAlignment (g, Graphics_CENTRE, Graphics_BOTTOM);
- Graphics_text (g, x, 0.0, U"%t");
- Graphics_text (g, x0, 0.0, U"%t+1");
- double y3 = 0.10;
- Graphics_setTextAlignment (g, Graphics_CENTRE, Graphics_HALF);
- Graphics_text (g, x, y3, U"%\\al__%t_(%i)%");
- Graphics_text (g, x0, y3, U"%\\al__%t+1_(%j)");
- }
- }
- void HMM_drawForwardAndBackwardProbabilitiesIllustration (Graphics g, bool garnish) {
- double xfrac = 0.1, xs = 1.0 / (0.5 - xfrac), r = 0.03;
- Graphics_Viewport vp = Graphics_insetViewport (g, 0.0, 0.5-xfrac, 0.0, 1.0);
- HMM_drawForwardProbabilitiesIllustration (g, false);
- Graphics_resetViewport (g, vp);
- Graphics_insetViewport (g, 0.5 + xfrac, 1.0, 0.0, 1.0);
- HMM_drawBackwardProbabilitiesIllustration (g, false);
- Graphics_resetViewport (g, vp);
- Graphics_setWindow (g, 0.0, xs, 0.0, 1.0);
- if (garnish) {
- double rx1 = 1.0 + xs * 2.0 * xfrac + 0.1, rx2 = rx1 + 0.9 - 0.1, y1 = 0.1;
- Graphics_line (g, 0.9 + r, 0.5, rx1 - r, 0.5);
- Graphics_setTextAlignment (g, Graphics_CENTRE, Graphics_BOTTOM);
- Graphics_text (g, 0.9, 0.5 + r, U"%s__%i_");
- Graphics_text (g, rx1, 0.5 + r, U"%s__%j_");
- Graphics_setTextAlignment (g, Graphics_CENTRE, Graphics_TOP);
- Graphics_text (g, 1.0 + xfrac * xs, 0.5, U"%a__%%ij%_%b__%j_(O__%t+1_)");
- Graphics_setTextAlignment (g, Graphics_CENTRE, Graphics_BOTTOM);
- Graphics_text (g, 0.1, 0.0, U"%t-1");
- Graphics_text (g, 0.9, 0.0, U"%t");
- Graphics_text (g, rx1, 0.0, U"%t+1");
- Graphics_text (g, rx2, 0.0, U"%t+2");
- Graphics_setLineType (g, Graphics_DASHED);
- double x4 = rx1 - 0.06, x3 = 0.9 + 0.06;
- Graphics_line (g, x3, 0.7, x3, 0.0);
- Graphics_line (g, x4, 0.7, x4, 0.0);
- Graphics_setLineType (g, Graphics_DRAWN);
- Graphics_arrow (g, x4, y1, x4 + 0.2, y1);
- Graphics_arrow (g, x3, y1, x3 - 0.2, y1);
- Graphics_setTextAlignment (g, Graphics_RIGHT, Graphics_BOTTOM);
- Graphics_text (g, x3 - 0.01, y1, U"\\al__%t_(i)");
- Graphics_setTextAlignment (g, Graphics_LEFT, Graphics_BOTTOM);
- Graphics_text (g, x4 + 0.01, y1, U"\\be__%t+1_(j)");
- }
- }
- void HMM_HMMBaumWelch_getXi (HMM me, HMMBaumWelch thee, integer *obs) {
- for (integer it = 1; it <= thy numberOfTimes - 1; it ++) {
- double sum = 0.0;
- for (integer is = 1; is <= thy numberOfStates; is ++) {
- for (integer js = 1; js <= thy numberOfStates; js ++) {
- thy xi [it] [is] [js] = thy alpha [is] [it] * thy beta [js] [it + 1] *
- my transitionProbs [is] [js] * my emissionProbs [js] [ obs [it + 1] ];
- sum += thy xi [it] [is] [js];
- }
- }
- for (integer is = 1; is <= my numberOfStates; is ++) {
- for (integer js = 1; js <= my numberOfStates; js ++) {
- thy xi [it] [is] [js] /= sum;
- }
- }
- }
- }
- void HMM_HMMBaumWelch_addEstimate (HMM me, HMMBaumWelch thee, integer *obs) {
- for (integer is = 1; is <= my numberOfStates; is ++) {
- // only for valid start states with p > 0
- if (my transitionProbs [0] [is] > 0.0) {
- thy aij_num [0] [is] += thy gamma [is] [1];
- thy aij_denom [0] [is] += 1.0;
- }
- }
- for (integer is = 1; is <= my numberOfStates; is ++) {
- double gammasum = 0.0;
- for (integer it = 1; it <= thy numberOfTimes - 1; it ++) {
- gammasum += thy gamma [is] [it];
- }
- for (integer js = 1; js <= my numberOfStates; js ++) {
- double xisum = 0.0;
- for (integer it = 1; it <= thy numberOfTimes - 1; it ++) {
- xisum += thy xi [it] [is] [js];
- }
- // zero probs signal invalid connections, don't reestimate
- if (my transitionProbs [is] [js] > 0.0) {
- thy aij_num [is] [js] += xisum;
- thy aij_denom [is] [js] += gammasum;
- }
- }
- /*
- Only reestimate the emissionProbs for a hidden markov model.
- A not hidden model is emulated with fixed emissionProbs.
- */
- if (! my notHidden) {
- gammasum += thy gamma [is] [thy numberOfTimes]; // now sum all, add last term
- for (integer k = 1; k <= my numberOfObservationSymbols; k ++) {
- double gammasum_k = 0.0;
- for (integer it = 1; it <= thy numberOfTimes; it ++) {
- if (obs [it] == k) {
- gammasum_k += thy gamma [is] [it];
- }
- }
- // only reestimate probs > 0 !
- if (my emissionProbs [is] [k] > 0.0) {
- thy bik_num [is] [k] += gammasum_k;
- thy bik_denom [is] [k] += gammasum;
- }
- }
- }
- // For a left-to-right model the final state determines the transition prob to go to the END state
- if (my leftToRight) {
- thy aij_num [is] [my numberOfStates + 1] += thy gamma [is] [thy numberOfTimes];
- thy aij_denom [is] [my numberOfStates + 1] += 1.0;
- }
- }
- }
- void HMM_HMMBaumWelch_reestimate (HMM me, HMMBaumWelch thee) {
- double p;
- /*
- If we only have a couple of training sequences and they all happen to start with the same symbol,
- one or more starting probabilities can be zero.
- What to do with the P_k (see Rabiner formulas 109-110)?
- */
- for (integer is = 1; is <= my numberOfStates; is ++) {
- /*
- If we have not enough observation sequences it can happen that some probabilities
- become zero. This means that for some future observation sequences the probability evaluation
- returns p=0 for sequences where these transitions turn up. This makes recognition impossible and also comparisons between models are difficult.
- We can prevent this from happening by asumimg a minimal probability for valid transitions
- i.e. which have initially p > 0.
- */
- if (my transitionProbs [0] [is] > 0.0) {
- p = thy aij_num [0] [is] / thy aij_denom [0] [is];
- my transitionProbs [0] [is] = p > 0.0 ? p : thy minProb;
- }
- for (integer js = 1; js <= my numberOfStates; js ++) {
- if (my transitionProbs [is] [js] > 0.0) {
- p = thy aij_num [is] [js] / thy aij_denom [is] [js];
- my transitionProbs [is] [js] = p > 0.0 ? p : thy minProb;
- }
- }
- if (! my notHidden) {
- for (integer k = 1; k <= my numberOfObservationSymbols; k ++) {
- if (my emissionProbs [is] [k] > 0.0) {
- p = thy bik_num [is] [k] / thy bik_denom [is] [k];
- my emissionProbs [is] [k] = p > 0.0 ? p : thy minProb;
- }
- }
- }
- if (my leftToRight && my transitionProbs [is] [my numberOfStates + 1] > 0.0) {
- p = thy aij_num [is] [my numberOfStates + 1] / thy aij_denom [is] [my numberOfStates + 1];
- my transitionProbs [is] [my numberOfStates + 1] = p > 0.0 ? p : thy minProb;
- }
- }
- }
- void HMM_HMMBaumWelch_forward (HMM me, HMMBaumWelch thee, integer *obs) {
- // initialise at t = 1 & scale
- thy scale [1] = 0.0;
- for (integer js = 1; js <= my numberOfStates; js ++) {
- thy alpha [js] [1] = my transitionProbs [0] [js] * my emissionProbs [js] [obs [1]];
- thy scale [1] += thy alpha [js] [1];
- }
- for (integer js = 1; js <= my numberOfStates; js ++) {
- thy alpha [js] [1] /= thy scale [1];
- }
- // recursion
- for (integer it = 2; it <= thy numberOfTimes; it ++) {
- thy scale [it] = 0.0;
- for (integer js = 1; js <= my numberOfStates; js ++) {
- double sum = 0.0;
- for (integer is = 1; is <= my numberOfStates; is ++) {
- sum += thy alpha [is] [it - 1] * my transitionProbs [is] [js];
- }
- thy alpha [js] [it] = sum * my emissionProbs [js] [obs [it]];
- thy scale [it] += thy alpha [js] [it];
- }
- for (integer js = 1; js <= my numberOfStates; js ++) {
- thy alpha [js] [it] /= thy scale [it];
- }
- }
- for (integer it = 1; it <= thy numberOfTimes; it ++) {
- thy lnProb += log (thy scale [it]);
- }
- }
- void HMM_HMMBaumWelch_backward (HMM me, HMMBaumWelch thee, integer *obs) {
- for (integer is = 1; is <= my numberOfStates; is ++) {
- thy beta [is] [thy numberOfTimes] = 1.0 / thy scale [thy numberOfTimes];
- }
- for (integer it = thy numberOfTimes - 1; it >= 1; it --) {
- for (integer is = 1; is <= my numberOfStates; is ++) {
- longdouble sum = 0.0;
- for (integer js = 1; js <= my numberOfStates; js ++) {
- sum += thy beta [js] [it + 1] * my transitionProbs [is] [js] * my emissionProbs [js] [obs [it + 1]];
- }
- thy beta [is] [it] = sum / thy scale [it];
- }
- }
- }
- /*************************** HMM decoding ***********************************/
- // precondition: valid symbols, i.e. 1 <= o [i] <= my numberOfSymbols for i=1..nt
- void HMM_HMMViterbi_decode (HMM me, HMMViterbi thee, integer *obs) {
- integer ntimes = thy numberOfTimes;
- // initialisation
- for (integer is = 1; is <= my numberOfStates; is ++) {
- thy viterbi [is] [1] = my transitionProbs [0] [is] * my emissionProbs [is] [obs [1]];
- thy bp [is] [1] = 0;
- }
- // recursion
- for (integer it = 2; it <= ntimes; it ++) {
- for (integer is = 1; is <= my numberOfStates; is ++) {
- // all transitions isp -> is from previous time to current
- double max_score = -1; // any negative number is ok
- for (integer isp = 1; isp <= my numberOfStates; isp ++) {
- double score = thy viterbi [isp] [it - 1] * my transitionProbs [isp] [is]; // * my emissionProbs [is] [obs [it]]
- if (score > max_score) {
- max_score = score;
- thy bp [is] [it] = isp;
- }
- }
- thy viterbi [is] [it] = max_score * my emissionProbs [is] [ obs [it] ];
- }
- }
- // path starts at state with best end probability
- thy path [ntimes] = 1;
- thy prob = thy viterbi [1] [ntimes];
- for (integer is = 2; is <= my numberOfStates; is ++) {
- if (thy viterbi [is] [ntimes] > thy prob) {
- thy prob = thy viterbi [thy path [ntimes] = is] [ntimes];
- }
- }
- // trace back and get path
- for (integer it = ntimes; it > 1; it --) {
- thy path [it - 1] = thy bp [thy path [it]] [it];
- }
- }
- autoHMMStateSequence HMM_HMMObservationSequence_to_HMMStateSequence (HMM me, HMMObservationSequence thee) {
- try {
- autoStringsIndex si = HMM_HMMObservationSequence_to_StringsIndex (me, thee);
- integer *obs = si -> classIndex.at; // convenience
- integer numberOfUnknowns = StringsIndex_countItems (si.get(), 0);
-
- Melder_require (numberOfUnknowns == 0, U"Unknown observation symbol(s) (# = ", numberOfUnknowns, U").");
- integer numberOfTimes = thy rows.size;
- autoHMMViterbi v = HMM_to_HMMViterbi (me, obs, numberOfTimes);
- autoHMMStateSequence him = HMMStateSequence_create (numberOfTimes);
- // trace the path and get states
- for (integer it = 1; it <= numberOfTimes; it ++) {
- HMMState hmms = my states->at [v -> path [it]];
- his strings [it] = Melder_dup (hmms -> label.get());
- his numberOfStrings ++;
- }
- return him;
- } catch (MelderError) {
- Melder_throw (me, U": no HMMStateSequence created.");
- }
- }
- double HMM_HMMStateSequence_getProbability (HMM me, HMMStateSequence thee) {
- autoStringsIndex si = HMM_HMMStateSequence_to_StringsIndex (me, thee);
- integer numberOfUnknowns = StringsIndex_countItems (si.get(), 0);
- integer *index = si -> classIndex.at;
- if (index == 0) {
- return undefined;
- }
- if (numberOfUnknowns > 0) {
- Melder_warning (U"Unknown states (# = ", numberOfUnknowns, U").");
- return undefined;
- }
- double p0 = my transitionProbs [0] [index [1]];
- Melder_require (p0 > 0.0, U"You should not start with a zero probability state.");
-
- double lnp = log (p0);
- for (integer it = 2; it <= thy numberOfStrings; it ++) {
- lnp += log (my transitionProbs [index [it - 1]] [index [it]]);
- }
- return lnp;
- }
- double HMM_getProbabilityAtTimeBeingInState (HMM me, integer itime, integer istate) {
- if (istate < 1 || istate > my numberOfStates) {
- return undefined;
- }
- autoNUMvector<double> scale (1, itime);
- autoNUMvector<double> alpha_t (1, my numberOfStates);
- autoNUMvector<double> alpha_tm1 (1, my numberOfStates);
- for (integer js = 1; js <= my numberOfStates; js ++) {
- alpha_t [js] = my transitionProbs [0] [js];
- scale [1] += alpha_t [js];
- }
- for (integer js = 1; js <= my numberOfStates; js ++) {
- alpha_t [js] /= scale [1];
- }
- // recursion
- for (integer it = 2; it <= itime; it ++) {
- for (integer js = 1; js <= my numberOfStates; js ++) {
- alpha_tm1 [js] = alpha_t [js];
- }
- for (integer js = 1; js <= my numberOfStates; js ++) {
- longdouble sum = 0.0;
- for (integer is = 1; is <= my numberOfStates; is ++) {
- sum += alpha_tm1 [is] * my transitionProbs [is] [js];
- }
- alpha_t [js] = sum;
- scale [it] += alpha_t [js];
- }
- for (integer js = 1; js <= my numberOfStates; js ++) {
- alpha_t [js] /= scale [it];
- }
- }
- longdouble lnp = 0.0;
- for (integer it = 1; it <= itime; it ++) {
- lnp += log (scale [it]);
- }
- lnp = alpha_t [istate] > 0 ? lnp + log (alpha_t [istate]) : -INFINITY; // p = 0 -> ln(p)=-infinity // ppgb FIXME infinity is een laag getal
- return lnp;
- }
- double HMM_getProbabilityAtTimeBeingInStateEmittingSymbol (HMM me, integer itime, integer istate, integer isymbol) {
- // for a notHidden model emissionProbs may be zero!
- if (isymbol < 1 || isymbol > my numberOfObservationSymbols || my emissionProbs [istate] [isymbol] == 0) {
- return undefined;
- }
- double lnp = HMM_getProbabilityAtTimeBeingInState (me, itime, istate);
- return ( isundef (lnp) ? undefined : lnp + log (my emissionProbs [istate] [isymbol]) );
- }
- double HMM_getProbabilityOfObservations (HMM me, integer *obs, integer numberOfTimes) {
- autoNUMvector <double> scale (1, numberOfTimes);
- autoNUMvector <double> alpha_t (1, my numberOfStates);
- autoNUMvector <double> alpha_tm1 (1, my numberOfStates);
- // initialise
- for (integer js = 1; js <= my numberOfStates; js ++) {
- alpha_t [js] = my transitionProbs [0] [js] * my emissionProbs [js] [obs [1]];
- scale [1] += alpha_t [js];
- }
- Melder_require (scale [1] > 0.0, U"The observation sequence should not start with a symbol whose state has zero starting probability.");
-
- for (integer js = 1; js <= my numberOfStates; js ++) {
- alpha_t [js] /= scale [1];
- }
- // recursion
- for (integer it = 2; it <= numberOfTimes; it ++) {
- for (integer js = 1; js <= my numberOfStates; js ++) {
- alpha_tm1 [js] = alpha_t [js];
- }
- for (integer js = 1; js <= my numberOfStates; js ++) {
- longdouble sum = 0.0;
- for (integer is = 1; is <= my numberOfStates; is ++) {
- sum += alpha_tm1 [is] * my transitionProbs [is] [js];
- }
- alpha_t [js] = sum * my emissionProbs [js] [obs [it]];
- scale [it] += alpha_t [js];
- }
- if (scale [it] <= 0.0) {
- return -INFINITY;
- }
- for (integer js = 1; js <= my numberOfStates; js ++) {
- alpha_t [js] /= scale [it];
- }
- }
- double lnp = 0.0;
- for (integer it = 1; it <= numberOfTimes; it ++) {
- lnp += log (scale [it]);
- }
- return lnp;
- }
- double HMM_HMMObservationSequence_getProbability (HMM me, HMMObservationSequence thee) {
- autoStringsIndex si = HMM_HMMObservationSequence_to_StringsIndex (me, thee);
- integer *index = si -> classIndex.at;
- integer numberOfUnknowns = StringsIndex_countItems (si.get(), 0);
- Melder_require (numberOfUnknowns == 0, U"Unknown observations (# = ", numberOfUnknowns, U").");
-
- return HMM_getProbabilityOfObservations (me, index, thy rows.size);
- }
- double HMM_HMMObservationSequence_getCrossEntropy (HMM me, HMMObservationSequence thee) {
- double lnp = HMM_HMMObservationSequence_getProbability (me, thee);
- return isundef (lnp) ? undefined :
- -lnp / (NUMln10 * HMMObservationSequence_getNumberOfObservations (thee));
- }
- double HMM_HMMObservationSequence_getPerplexity (HMM me, HMMObservationSequence thee) {
- double ce = HMM_HMMObservationSequence_getCrossEntropy (me, thee);
- return isundef (ce) ? undefined : pow (2.0, ce);
- }
- autoHMM HMM_createFromHMMObservationSequence (HMMObservationSequence me, integer numberOfStates, int leftToRight) {
- try {
- autoHMM thee = Thing_new (HMM);
- autoStrings s = HMMObservationSequence_to_Strings (me);
- autoDistributions d = Strings_to_Distributions (s.get());
- integer numberOfObservationSymbols = d -> numberOfRows;
- thy notHidden = numberOfStates < 1;
- numberOfStates = numberOfStates > 0 ? numberOfStates : numberOfObservationSymbols;
- HMM_init (thee.get(), numberOfStates, numberOfObservationSymbols, leftToRight);
- for (integer i = 1; i <= numberOfObservationSymbols; i ++) {
- conststring32 label = d -> rowLabels [i].get();
- autoHMMObservation hmmo = HMMObservation_create (label, 0, 0, 0);
- HMM_addObservation_move (thee.get(), hmmo.move());
- if (thy notHidden) {
- autoHMMState hmms = HMMState_create (label);
- HMM_addState_move (thee.get(), hmms.move());
- }
- }
- if (! thy notHidden) {
- HMM_setDefaultStates (thee.get());
- }
- HMM_setDefaultTransitionProbs (thee.get());
- HMM_setDefaultStartProbs (thee.get());
- HMM_setDefaultEmissionProbs (thee.get());
- return thee;
- } catch (MelderError) {
- Melder_throw (me, U": no HMM created.");
- }
- }
- autoTableOfReal HMMObservationSequence_to_TableOfReal_transitions (HMMObservationSequence me, int probabilities) {
- try {
- autoStrings thee = HMMObservationSequence_to_Strings (me);
- autoTableOfReal him = Strings_to_TableOfReal_transitions (thee.get(), probabilities);
- return him;
- } catch (MelderError) {
- Melder_throw (me, U": no transitions created.");
- }
- }
- autoStringsIndex HMM_HMMObservationSequence_to_StringsIndex (HMM me, HMMObservationSequence thee) {
- try {
- autoStrings classes = Thing_new (Strings);
- classes -> strings = autostring32vector (my numberOfObservationSymbols);
- for (integer is = 1; is <= my numberOfObservationSymbols; is ++) {
- HMMObservation hmmo = my observationSymbols->at [is];
- classes -> strings [is] = Melder_dup (hmmo -> label.get());
- classes -> numberOfStrings ++;
- }
- autoStrings obs = HMMObservationSequence_to_Strings (thee);
- autoStringsIndex him = Stringses_to_StringsIndex (obs.get(), classes.get());
- return him;
- } catch (MelderError) {
- Melder_throw (me, U": no StringsIndex created.");
- }
- }
- autoStringsIndex HMM_HMMStateSequence_to_StringsIndex (HMM me, HMMStateSequence thee) {
- try {
- autoStrings classes = Thing_new (Strings);
- classes -> strings = autostring32vector (my numberOfObservationSymbols);
- for (integer is = 1; is <= my numberOfStates; is ++) {
- HMMState hmms = my states->at [is];
- classes -> strings [is] = Melder_dup (hmms -> label.get());
- classes -> numberOfStrings ++;
- }
- autoStrings sts = HMMStateSequence_to_Strings (thee);
- autoStringsIndex him = Stringses_to_StringsIndex (sts.get(), classes.get());
- return him;
- } catch (MelderError) {
- Melder_throw (me, U": no StringsIndex created.");
- }
- }
- autoTableOfReal HMM_HMMObservationSequence_to_TableOfReal_transitions (HMM me, HMMObservationSequence thee, int probabilities) {
- try {
- autoStringsIndex si = HMM_HMMObservationSequence_to_StringsIndex (me, thee);
- autoTableOfReal him = StringsIndex_to_TableOfReal_transitions (si.get(), probabilities);
- return him;
- } catch (MelderError) {
- Melder_throw (me, U": no transition table created for HMMObservationSequence.");
- }
- }
- autoTableOfReal HMM_HMMStateSequence_to_TableOfReal_transitions (HMM me, HMMStateSequence thee, int probabilities) {
- try {
- autoStringsIndex si = HMM_HMMStateSequence_to_StringsIndex (me, thee);
- autoTableOfReal him = StringsIndex_to_TableOfReal_transitions (si.get(), probabilities);
- return him;
- } catch (MelderError) {
- Melder_throw (me, U": no transition table created for HMMStateSequence.");
- }
- }
- autoTableOfReal StringsIndex_to_TableOfReal_transitions (StringsIndex me, int probabilities) {
- try {
- integer numberOfTypes = my classes->size;
- autoTableOfReal thee = TableOfReal_create (numberOfTypes + 1, numberOfTypes + 1);
- for (integer i = 1; i <= numberOfTypes; i ++) {
- SimpleString s = (SimpleString) my classes->at [i];
- TableOfReal_setRowLabel (thee.get(), i, s -> string.get());
- TableOfReal_setColumnLabel (thee.get(), i, s -> string.get());
- }
- for (integer i = 2; i <= my numberOfItems; i ++) {
- if (my classIndex [i - 1] > 0 && my classIndex [i] > 0) { // a zero is a restart!
- thy data [my classIndex [i-1]] [my classIndex [i]] ++;
- }
- }
- longdouble sum = 0.0;
- for (integer i = 1; i <= numberOfTypes; i ++) {
- double rowSum = 0.0, colSum = 0.0;
- for (integer j = 1; j <= numberOfTypes; j ++) {
- rowSum += thy data [i] [j];
- }
- thy data [i] [numberOfTypes + 1] = rowSum;
- for (integer j = 1; j <= numberOfTypes; j ++) {
- colSum += thy data [j] [i];
- }
- thy data [numberOfTypes + 1] [i] = colSum;
- sum += colSum;
- }
- thy data [numberOfTypes + 1] [numberOfTypes + 1] = sum;
- if (probabilities && sum > 0.0) {
- for (integer i = 1; i <= numberOfTypes; i ++) {
- if (thy data [i] [numberOfTypes + 1] > 0.0) {
- for (integer j = 1; j <= numberOfTypes; j ++) {
- thy data [i] [j] /= thy data [i] [numberOfTypes + 1];
- }
- }
- }
- for (integer i = 1; i <= numberOfTypes; i ++) {
- thy data [i] [numberOfTypes + 1] /= sum;
- thy data [numberOfTypes + 1] [i] /= sum;
- }
- }
- return thee;
- } catch (MelderError) {
- Melder_throw (me, U": no transition table created.");
- }
- }
- autoTableOfReal Strings_to_TableOfReal_transitions (Strings me, int probabilities) {
- try {
- autoStringsIndex him = Strings_to_StringsIndex (me);
- autoTableOfReal thee = StringsIndex_to_TableOfReal_transitions (him.get(), probabilities);
- return thee;
- } catch (MelderError) {
- Melder_throw (me, U": no transition table created.");
- }
- }
- /* End of file HMM.cpp */
|