1 /* NUMinterpol.cpp
2 *
3 * Copyright (C) 1992-2008,2011,2012,2014,2015,2017,2018,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 /*
20 * pb 2002/03/07 GPL
21 * pb 2003/06/19 ridders3 replaced with ridders
22 * pb 2003/07/09 gsl
23 * pb 2007/01/27 use #defines for value interpolation
24 * pb 2007/08/20 built a "weird value" check into NUMviterbi (bug report by Adam Jacks)
25 * pb 2011/03/29 C++
26 */
27
28 #include "melder.h"
29 #include "../dwsys/NUM2.h"
30
31 #if defined (__POWERPC__)||1
32 #define RECOMPUTE_SINES 0
33 #else
34 #define RECOMPUTE_SINES 1
35 #endif
NUM_interpolate_sinc(constVEC const & y,double x,integer maxDepth)36 double NUM_interpolate_sinc (constVEC const& y, double x, integer maxDepth) {
37 const integer midleft = (integer) floor (x), midright = midleft + 1;
38 double result = 0.0;
39 if (y.size < 1)
40 return undefined; // there exists no best guess
41 if (x < 1)
42 return y [1]; // offleft: constant extrapolation
43 if (x > y.size)
44 return y [y.size]; // offright: constant extrapolation
45 if (x == midleft)
46 return y [midleft]; // the interpolated curve goes through the points
47 /*
48 1 < x < y.size && x not integer: interpolate.
49 */
50 Melder_clipRight (& maxDepth, midright - 1);
51 Melder_clipRight (& maxDepth, y.size - midleft);
52 if (maxDepth <= NUM_VALUE_INTERPOLATE_NEAREST)
53 return y [(integer) floor (x + 0.5)];
54 if (maxDepth == NUM_VALUE_INTERPOLATE_LINEAR)
55 return y [midleft] + (x - midleft) * (y [midright] - y [midleft]);
56 if (maxDepth == NUM_VALUE_INTERPOLATE_CUBIC) {
57 const double yl = y [midleft], yr = y [midright];
58 const double dyl = 0.5 * (yr - y [midleft - 1]), dyr = 0.5 * (y [midright + 1] - yl);
59 const double fil = x - midleft, fir = midright - x;
60 return yl * fir + yr * fil - fil * fir * (0.5 * (dyr - dyl) + (fil - 0.5) * (dyl + dyr - 2 * (yr - yl)));
61 }
62 /*
63 maxDepth >= 3: sinc interpolation
64 */
65 const integer left = midright - maxDepth;
66 const integer right = midleft + maxDepth;
67 double a = NUMpi * (x - midleft);
68 double halfsina = 0.5 * sin (a);
69 double aa = a / (x - left + 1.0);
70 double daa = NUMpi / (x - left + 1.0);
71 #if ! RECOMPUTE_SINES
72 double cosaa = cos (aa);
73 double sinaa = sin (aa);
74 double cosdaa = cos (daa);
75 double sindaa = sin (daa);
76 #endif
77 for (integer ix = midleft; ix >= left; ix --) {
78 #if RECOMPUTE_SINES
79 const double d = halfsina / a * (1.0 + cos (aa));
80 #else
81 const double d = halfsina / a * (1.0 + cosaa);
82 #endif
83 result += y [ix] * d;
84 a += NUMpi;
85 #if RECOMPUTE_SINES
86 aa += daa;
87 #else
88 const double help = cosaa * cosdaa - sinaa * sindaa;
89 sinaa = cosaa * sindaa + sinaa * cosdaa;
90 cosaa = help;
91 #endif
92 halfsina = - halfsina;
93 }
94 a = NUMpi * (midright - x);
95 halfsina = 0.5 * sin (a);
96 aa = a / (right - x + 1.0);
97 daa = NUMpi / (right - x + 1.0);
98 #if ! RECOMPUTE_SINES
99 cosaa = cos (aa);
100 sinaa = sin (aa);
101 cosdaa = cos (daa);
102 sindaa = sin (daa);
103 #endif
104 for (integer ix = midright; ix <= right; ix ++) {
105 #if RECOMPUTE_SINES
106 const double d = halfsina / a * (1.0 + cos (aa));
107 #else
108 const double d = halfsina / a * (1.0 + cosaa);
109 #endif
110 result += y [ix] * d;
111 a += NUMpi;
112 #if RECOMPUTE_SINES
113 aa += daa;
114 #else
115 const double help = cosaa * cosdaa - sinaa * sindaa;
116 sinaa = cosaa * sindaa + sinaa * cosdaa;
117 cosaa = help;
118 #endif
119 halfsina = - halfsina;
120 }
121 return result;
122 }
123
124 /********** Improving extrema **********/
125 #pragma mark Improving extrema
126
127 struct improve_params {
128 integer depth;
129 constVEC y;
130 bool isMaximum;
131 };
132
improve_evaluate(double x,void * closure)133 static double improve_evaluate (double x, void *closure) {
134 struct improve_params *me = (struct improve_params *) closure;
135 const double y = NUM_interpolate_sinc (my y, x, my depth);
136 return my isMaximum ? - y : y;
137 }
138
NUMimproveExtremum(constVEC const & y,integer ixmid,integer interpolationDepth,double * ixmid_real,bool isMaximum)139 double NUMimproveExtremum (constVEC const& y, integer ixmid, integer interpolationDepth, double *ixmid_real, bool isMaximum) {
140 struct improve_params params;
141 double result;
142 if (ixmid <= 1) {
143 *ixmid_real = double (1);
144 return y [1];
145 }
146 if (ixmid >= y.size) {
147 *ixmid_real = double (y.size);
148 return y [y.size];
149 }
150 if (interpolationDepth <= NUM_PEAK_INTERPOLATE_NONE) {
151 *ixmid_real = double (ixmid);
152 return y [ixmid];
153 }
154 if (interpolationDepth == NUM_PEAK_INTERPOLATE_PARABOLIC) {
155 const double dy = 0.5 * (y [ixmid + 1] - y [ixmid - 1]);
156 const double d2y = 2 * y [ixmid] - y [ixmid - 1] - y [ixmid + 1];
157 *ixmid_real = ixmid + dy / d2y;
158 return y [ixmid] + 0.5 * dy * dy / d2y;
159 }
160 /*
161 Cubic or sinc interpolation.
162 */
163 params. depth = (
164 interpolationDepth == NUM_PEAK_INTERPOLATE_CUBIC ? NUM_VALUE_INTERPOLATE_CUBIC :
165 interpolationDepth == NUM_PEAK_INTERPOLATE_SINC70 ? NUM_VALUE_INTERPOLATE_SINC70 :
166 NUM_VALUE_INTERPOLATE_SINC700
167 );
168 params. y = y;
169 params. isMaximum = isMaximum;
170 *ixmid_real = NUMminimize_brent (improve_evaluate, ixmid - 1, ixmid + 1, & params, 1e-10, & result);
171 return isMaximum ? - result : result;
172 }
173
NUMimproveMinimum(constVEC const & y,integer ixmid,integer interpolationDepth,double * ixmid_real)174 double NUMimproveMinimum (constVEC const& y, integer ixmid, integer interpolationDepth, double *ixmid_real) {
175 return NUMimproveExtremum (y, ixmid, interpolationDepth, ixmid_real, false);
176 }
NUMimproveMaximum(constVEC const & y,integer ixmid,integer interpolationDepth,double * ixmid_real)177 double NUMimproveMaximum (constVEC const& y, integer ixmid, integer interpolationDepth, double *ixmid_real) {
178 return NUMimproveExtremum (y, ixmid, interpolationDepth, ixmid_real, true);
179 }
180
181 /********** Viterbi **********/
182
NUM_viterbi(integer numberOfFrames,integer maxnCandidates,integer (* getNumberOfCandidates)(integer iframe,void * closure),double (* getLocalCost)(integer iframe,integer icand,void * closure),double (* getTransitionCost)(integer iframe,integer icand1,integer icand2,void * closure),void (* putResult)(integer iframe,integer place,void * closure),void * closure)183 void NUM_viterbi (
184 integer numberOfFrames, integer maxnCandidates,
185 integer (*getNumberOfCandidates) (integer iframe, void *closure),
186 double (*getLocalCost) (integer iframe, integer icand, void *closure),
187 double (*getTransitionCost) (integer iframe, integer icand1, integer icand2, void *closure),
188 void (*putResult) (integer iframe, integer place, void *closure),
189 void *closure)
190 {
191 autoMAT delta = raw_MAT (numberOfFrames, maxnCandidates);
192 autoINTMAT psi = raw_INTMAT (numberOfFrames, maxnCandidates);
193 autoINTVEC numberOfCandidates = raw_INTVEC (numberOfFrames);
194 for (integer iframe = 1; iframe <= numberOfFrames; iframe ++) {
195 numberOfCandidates [iframe] = getNumberOfCandidates (iframe, closure);
196 for (integer icand = 1; icand <= numberOfCandidates [iframe]; icand ++)
197 delta [iframe] [icand] = - getLocalCost (iframe, icand, closure);
198 }
199 for (integer iframe = 2; iframe <= numberOfFrames; iframe ++) {
200 for (integer icand2 = 1; icand2 <= numberOfCandidates [iframe]; icand2 ++) {
201 double maximum = -1e308;
202 integer place = 0;
203 for (integer icand1 = 1; icand1 <= numberOfCandidates [iframe - 1]; icand1 ++) {
204 const double value = delta [iframe - 1] [icand1] + delta [iframe] [icand2]
205 - getTransitionCost (iframe, icand1, icand2, closure);
206 if (value > maximum) {
207 maximum = value;
208 place = icand1;
209 }
210 }
211 if (place == 0)
212 Melder_throw (U"Viterbi algorithm cannot compute a track because of weird values.");
213 delta [iframe] [icand2] = maximum;
214 psi [iframe] [icand2] = place;
215 }
216 }
217 /*
218 Find the end of the most probable path.
219 */
220 integer place;
221 double maximum = delta [numberOfFrames] [place = 1];
222 for (integer icand = 2; icand <= numberOfCandidates [numberOfFrames]; icand ++) {
223 if (delta [numberOfFrames] [icand] > maximum)
224 maximum = delta [numberOfFrames] [place = icand];
225 }
226 /*
227 Backtrack.
228 */
229 for (integer iframe = numberOfFrames; iframe >= 1; iframe --) {
230 putResult (iframe, place, closure);
231 place = psi [iframe] [place];
232 }
233 }
234
235 /******************/
236
237 struct parm2 {
238 integer ntrack;
239 integer ncomb;
240 INTMAT indices;
241 double (*getLocalCost) (integer iframe, integer icand, integer itrack, void *closure);
242 double (*getTransitionCost) (integer iframe, integer icand1, integer icand2, integer itrack, void *closure);
243 void (*putResult) (integer iframe, integer place, integer itrack, void *closure);
244 void *closure;
245 };
246
getNumberOfCandidates_n(integer iframe,void * closure)247 static integer getNumberOfCandidates_n (integer iframe, void *closure) {
248 struct parm2 *me = (struct parm2 *) closure;
249 (void) iframe;
250 return my ncomb;
251 }
getLocalCost_n(integer iframe,integer jcand,void * closure)252 static double getLocalCost_n (integer iframe, integer jcand, void *closure) {
253 struct parm2 *me = (struct parm2 *) closure;
254 double localCost = 0.0;
255 for (integer itrack = 1; itrack <= my ntrack; itrack ++)
256 localCost += my getLocalCost (iframe, my indices [jcand] [itrack], itrack, my closure);
257 return localCost;
258 }
getTransitionCost_n(integer iframe,integer jcand1,integer jcand2,void * closure)259 static double getTransitionCost_n (integer iframe, integer jcand1, integer jcand2, void *closure) {
260 struct parm2 *me = (struct parm2 *) closure;
261 double transitionCost = 0.0;
262 for (integer itrack = 1; itrack <= my ntrack; itrack ++)
263 transitionCost += my getTransitionCost (iframe,
264 my indices [jcand1] [itrack], my indices [jcand2] [itrack], itrack, my closure);
265 return transitionCost;
266 }
putResult_n(integer iframe,integer jplace,void * closure)267 static void putResult_n (integer iframe, integer jplace, void *closure) {
268 struct parm2 *me = (struct parm2 *) closure;
269 for (integer itrack = 1; itrack <= my ntrack; itrack ++)
270 my putResult (iframe, my indices [jplace] [itrack], itrack, my closure);
271 }
272
NUM_viterbi_multi(integer nframe,integer ncand,integer ntrack,double (* getLocalCost)(integer iframe,integer icand,integer itrack,void * closure),double (* getTransitionCost)(integer iframe,integer icand1,integer icand2,integer itrack,void * closure),void (* putResult)(integer iframe,integer place,integer itrack,void * closure),void * closure)273 void NUM_viterbi_multi (
274 integer nframe, integer ncand, integer ntrack,
275 double (*getLocalCost) (integer iframe, integer icand, integer itrack, void *closure),
276 double (*getTransitionCost) (integer iframe, integer icand1, integer icand2, integer itrack, void *closure),
277 void (*putResult) (integer iframe, integer place, integer itrack, void *closure),
278 void *closure)
279 {
280 struct parm2 parm;
281
282 if (ntrack > ncand) Melder_throw (U"(NUM_viterbi_multi:) "
283 U"Number of tracks (", ntrack, U") should not exceed number of candidates (", ncand, U").");
284 const integer ncomb = Melder_iround (NUMcombinations (ncand, ntrack));
285 if (ncomb > 10'000'000) Melder_throw (U"(NUM_viterbi_multi:) "
286 U"Unrealistically high number of combinations (", ncomb, U").");
287 parm. ntrack = ntrack;
288 parm. ncomb = ncomb;
289
290 /*
291 For ncand == 5 and ntrack == 3, parm.indices is going to contain:
292 1 2 3
293 1 2 4
294 1 2 5
295 1 3 4
296 1 3 5
297 1 4 5
298 2 3 4
299 2 3 5
300 2 4 5
301 3 4 5
302 */
303 autoINTMAT indices = zero_INTMAT (ncomb, ntrack);
304 autoINTVEC icand = to_INTVEC (ntrack); // start out with "1 2 3"
305 integer jcomb = 0;
306 for (;;) {
307 jcomb ++;
308 for (integer itrack = 1; itrack <= ntrack; itrack ++)
309 indices [jcomb] [itrack] = icand [itrack];
310 integer itrack = ntrack;
311 for (; itrack >= 1; itrack --) {
312 if (++ icand [itrack] <= ncand - (ntrack - itrack)) {
313 for (integer jtrack = itrack + 1; jtrack <= ntrack; jtrack ++)
314 icand [jtrack] = icand [itrack] + jtrack - itrack;
315 break;
316 }
317 }
318 if (itrack == 0) break;
319 }
320 Melder_assert (jcomb == ncomb);
321 parm. indices = indices.get();
322 parm. getLocalCost = getLocalCost;
323 parm. getTransitionCost = getTransitionCost;
324 parm. putResult = putResult;
325 parm. closure = closure;
326 NUM_viterbi (nframe, ncomb, getNumberOfCandidates_n, getLocalCost_n, getTransitionCost_n, putResult_n, & parm);
327 }
328
329 /* End of file NUMinterpol.cpp */
330