1 /* RealTier.cpp
2  *
3  * Copyright (C) 1992-2012,2014-2021 Paul Boersma
4  *
5  * This code is free software; you can redistribute it and/or modify
6  * it under the terms of the GNU General Public License as published by
7  * the Free Software Foundation; either version 2 of the License, or (at
8  * your option) any later version.
9  *
10  * This code is distributed in the hope that it will be useful, but
11  * WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
13  * See the GNU General Public License for more details.
14  *
15  * You should have received a copy of the GNU General Public License
16  * along with this work. If not, see <http://www.gnu.org/licenses/>.
17  */
18 
19 #include "RealTier.h"
20 #include "Formula.h"
21 
22 #include "oo_DESTROY.h"
23 #include "RealTier_def.h"
24 #include "oo_COPY.h"
25 #include "RealTier_def.h"
26 #include "oo_EQUAL.h"
27 #include "RealTier_def.h"
28 #include "oo_CAN_WRITE_AS_ENCODING.h"
29 #include "RealTier_def.h"
30 #include "oo_WRITE_TEXT.h"
31 #include "RealTier_def.h"
32 #include "oo_READ_TEXT.h"
33 #include "RealTier_def.h"
34 #include "oo_WRITE_BINARY.h"
35 #include "RealTier_def.h"
36 #include "oo_READ_BINARY.h"
37 #include "RealTier_def.h"
38 #include "oo_DESCRIPTION.h"
39 #include "RealTier_def.h"
40 
41 /********** class RealPoint **********/
42 
43 Thing_implement (RealPoint, AnyPoint, 0);
44 
RealPoint_create(double time,double value)45 autoRealPoint RealPoint_create (double time, double value) {
46 	autoRealPoint me = Thing_new (RealPoint);
47 	my number = time;
48 	my value = value;
49 	return me;
50 }
51 
52 /********** class RealTier **********/
53 
v_info()54 void structRealTier :: v_info () {
55 	structFunction :: v_info ();
56 	MelderInfo_writeLine (U"Number of points: ", our points.size);
57 	MelderInfo_writeLine (U"Minimum value: ", RealTier_getMinimumValue (this));
58 	MelderInfo_writeLine (U"Maximum value: ", RealTier_getMaximumValue (this));
59 }
60 
v_getVector(integer irow,integer icol)61 double structRealTier :: v_getVector (integer irow, integer icol) {
62 	(void) irow;
63 	return RealTier_getValueAtIndex (this, icol);
64 }
65 
v_getFunction1(integer irow,double x)66 double structRealTier :: v_getFunction1 (integer irow, double x) {
67 	(void) irow;
68 	return RealTier_getValueAtTime (this, x);
69 }
70 
71 Thing_implement (RealTier, AnyTier, 0);   // the semantic superclass (see RealTier_def.h for explanation)
72 
RealTier_init(RealTier me,double tmin,double tmax)73 void RealTier_init (RealTier me, double tmin, double tmax) {
74 	my xmin = tmin;
75 	my xmax = tmax;
76 }
77 
RealTier_create(double tmin,double tmax)78 autoRealTier RealTier_create (double tmin, double tmax) {
79 	try {
80 		autoRealTier me = Thing_new (RealTier);
81 		RealTier_init (me.get(), tmin, tmax);
82 		return me;
83 	} catch (MelderError) {
84 		Melder_throw (U"RealTier not created.");
85 	}
86 }
87 
RealTier_createWithClass(double tmin,double tmax,ClassInfo klas)88 autoRealTier RealTier_createWithClass (double tmin, double tmax, ClassInfo klas) {
89 	try {
90 		autoRealTier me = Thing_newFromClass (klas).static_cast_move <structRealTier> ();
91 		RealTier_init (me.get(), tmin, tmax);
92 		return me;
93 	} catch (MelderError) {
94 		Melder_throw (klas -> className, U" not created.");
95 	}
96 }
97 
Thing_create()98 template <typename T> autoSomeThing <T> Thing_create () {
99 	return Thing_newFromClass (nullptr);
100 }
101 
102 template <>
Thing_create()103 autoSomeThing <structRealTier> Thing_create <structRealTier> () {
104 	return Thing_newFromClass (classRealTier). static_cast_move<structRealTier>();
105 }
106 
107 template <typename structSomeRealTier>
SomeRealTier_create(double tmin,double tmax)108 autoSomeThing <structSomeRealTier> SomeRealTier_create (double tmin, double tmax) {
109 	try {
110 		autoSomeThing <structSomeRealTier> me = Thing_create <structSomeRealTier> ();
111 		RealTier_init (me.get(), tmin, tmax);
112 		return me;
113 	} catch (MelderError) {
114 		Melder_throw (U"RealTier not created.");
115 	}
116 }
117 
RealTier_addPoint(RealTier me,double t,double value)118 void RealTier_addPoint (RealTier me, double t, double value) {
119 	try {
120 		autoRealPoint point = RealPoint_create (t, value);
121 		my points. addItem_move (point.move());
122 	} catch (MelderError) {
123 		Melder_throw (me, U": point not added.");
124 	}
125 }
126 
RealTier_getValueAtIndex(RealTier me,integer i)127 double RealTier_getValueAtIndex (RealTier me, integer i) {
128 	if (i < 1 || i > my points.size) return undefined;
129 	return my points.at [i] -> value;
130 }
131 
RealTier_getValueAtTime(RealTier me,double t)132 double RealTier_getValueAtTime (RealTier me, double t) {
133 	const integer n = my points.size;
134 	if (n == 0)
135 		return undefined;
136 	const RealPoint firstPoint = my points.at [1];
137 	if (t <= firstPoint -> number)
138 		return firstPoint -> value;   // constant extrapolation
139 	const RealPoint lastPoint = my points.at [n];
140 	if (t >= lastPoint -> number)
141 		return lastPoint -> value;   // constant extrapolation
142 	Melder_assert (n >= 2);
143 	const integer ileft = AnyTier_timeToLowIndex (me->asAnyTier(), t), iright = ileft + 1;
144 	Melder_assert (ileft >= 1 && iright <= n);
145 	const RealPoint pointLeft = my points.at [ileft];
146 	const RealPoint pointRight = my points.at [iright];
147 	const double tleft = pointLeft -> number, fleft = pointLeft -> value;
148 	const double tright = pointRight -> number, fright = pointRight -> value;
149 	return t == tright ? fright   // be very accurate
150 		: tleft == tright ? 0.5 * (fleft + fright)   // unusual, but possible; no preference
151 		: fleft + (t - tleft) * (fright - fleft) / (tright - tleft);   // linear interpolation
152 }
153 
RealTier_getMaximumValue(RealTier me)154 double RealTier_getMaximumValue (RealTier me) {
155 	double result = undefined;
156 	integer n = my points.size;
157 	for (integer i = 1; i <= n; i ++) {
158 		RealPoint point = my points.at [i];
159 		if (isundef (result) || point -> value > result)
160 			result = point -> value;
161 	}
162 	return result;
163 }
164 
RealTier_getMinimumValue(RealTier me)165 double RealTier_getMinimumValue (RealTier me) {
166 	double result = undefined;
167 	integer n = my points.size;
168 	for (integer i = 1; i <= n; i ++) {
169 		RealPoint point = my points.at [i];
170 		if (isundef (result) || point -> value < result)
171 			result = point -> value;
172 	}
173 	return result;
174 }
175 
RealTier_getArea(RealTier me,double tmin,double tmax)176 double RealTier_getArea (RealTier me, double tmin, double tmax) {
177 	integer n = my points.size, imin, imax;
178 	if (n == 0) return undefined;
179 	if (n == 1) return (tmax - tmin) * my points.at [1] -> value;
180 	imin = AnyTier_timeToLowIndex (me->asAnyTier(), tmin);
181 	if (imin == n) return (tmax - tmin) * my points.at [n] -> value;
182 	imax = AnyTier_timeToHighIndex (me->asAnyTier(), tmax);
183 	if (imax == 1) return (tmax - tmin) * my points.at [1] -> value;
184 	Melder_assert (imin < n);
185 	Melder_assert (imax > 1);
186 	/*
187 	 * Sum the areas between the points.
188 	 * This works even if imin is 0 (offleft) and/or imax is n + 1 (offright).
189 	 */
190 	longdouble area = 0.0;
191 	for (integer i = imin; i < imax; i ++) {
192 		double tleft, fleft, tright, fright;
193 		if (i == imin) {
194 			tleft = tmin;
195 			fleft = RealTier_getValueAtTime (me, tmin);
196 		} else {
197 			tleft = my points.at [i] -> number;
198 			fleft = my points.at [i] -> value;
199 		}
200 		if (i + 1 == imax) {
201 			tright = tmax;
202 			fright = RealTier_getValueAtTime (me, tmax);
203 		} else {
204 			tright = my points.at [i + 1] -> number;
205 			fright = my points.at [i + 1] -> value;
206 		}
207 		area += 0.5 * (fleft + fright) * (tright - tleft);
208 	}
209 	return (double) area;
210 }
211 
RealTier_getMean_curve(RealTier me,double tmin,double tmax)212 double RealTier_getMean_curve (RealTier me, double tmin, double tmax) {
213 	Function_unidirectionalAutowindow (me, & tmin, & tmax);
214 	const double area = RealTier_getArea (me, tmin, tmax);
215 	if (isundef (area)) return undefined;
216 	return area / (tmax - tmin);
217 }
218 
RealTier_getStandardDeviation_curve(RealTier me,double tmin,double tmax)219 double RealTier_getStandardDeviation_curve (RealTier me, double tmin, double tmax) {
220 	Function_unidirectionalAutowindow (me, & tmin, & tmax);
221 	const integer n = my points.size;
222 	if (n == 0) return undefined;
223 	if (n == 1) return 0.0;
224 	integer imin = AnyTier_timeToLowIndex (me->asAnyTier(), tmin);
225 	if (imin == n) return 0.0;
226 	integer imax = AnyTier_timeToHighIndex (me->asAnyTier(), tmax);
227 	if (imax == 1) return 0.0;
228 	Melder_assert (imin < n);
229 	Melder_assert (imax > 1);
230 	/*
231 	 * Add the areas between the points.
232 	 * This works even if imin is 0 (offleft) and/or imax is n + 1 (offright).
233 	 */
234 	double mean = RealTier_getMean_curve (me, tmin, tmax);
235 	longdouble integral = 0.0;
236 	for (integer i = imin; i < imax; i ++) {
237 		double tleft, fleft, tright, fright;
238 		if (i == imin) {
239 			tleft = tmin;
240 			fleft = RealTier_getValueAtTime (me, tmin);
241 		} else {
242 			tleft = my points.at [i] -> number;
243 			fleft = my points.at [i] -> value - mean;
244 		}
245 		if (i + 1 == imax) {
246 			tright = tmax;
247 			fright = RealTier_getValueAtTime (me, tmax);
248 		} else {
249 			tright = my points.at [i + 1] -> number;
250 			fright = my points.at [i + 1] -> value - mean;
251 		}
252 		/*
253 		 * The area is integral dt f^2
254 		 *   = integral dt [f1 + (f2-f1)/(t2-t1) (t-t1)]^2
255 		 *   = int dt f1^2 + int dt 2 f1 (f2-f1)/(t2-t1) (t-t1) + int dt [(f2-f1)/(t2-t1)]^2 (t-t1)^2
256 		 *   = f1^2 (t2-t1) + f1 (f2-f1)/(t2-t1) (t2-t1)^2 + 1/3 [(f2-f1)/(t2-t1)]^2 (t2-t1)^3
257 		 *   = (t2-t1) [f1 f2 + 1/3 (f2-f1)^2]
258 		 *   = (t2-t1) (f1^2 + f2^2 + 1/3 f1 f2)
259 		 *   = (t2-t1) [1/4 (f1+f2)^2 + 1/12 (f1-f2)^2]
260 		 * In the last expression, we have a sum of squares, which is computationally best.
261 		 */
262 		double sum = fleft + fright;
263 		double diff = fleft - fright;
264 		integral += (sum * sum + (1.0/3.0) * diff * diff) * (tright - tleft);
265 	}
266 	return sqrt (0.25 * (double) integral / (tmax - tmin));
267 }
268 
RealTier_getMean_points(RealTier me,double tmin,double tmax)269 double RealTier_getMean_points (RealTier me, double tmin, double tmax) {
270 	Function_unidirectionalAutowindow (me, & tmin, & tmax);
271 	integer imin, imax;
272 	integer n = AnyTier_getWindowPoints (me->asAnyTier(), tmin, tmax, & imin, & imax);
273 	if (n == 0) return undefined;
274 	longdouble sum = 0.0;
275 	for (integer i = imin; i <= imax; i ++)
276 		sum += my points.at [i] -> value;
277 	return (double) sum / n;
278 }
279 
RealTier_getStandardDeviation_points(RealTier me,double tmin,double tmax)280 double RealTier_getStandardDeviation_points (RealTier me, double tmin, double tmax) {
281 	Function_unidirectionalAutowindow (me, & tmin, & tmax);
282 	integer imin, imax;
283 	integer n = AnyTier_getWindowPoints (me->asAnyTier(), tmin, tmax, & imin, & imax);
284 	if (n < 2) return undefined;
285 	double mean = RealTier_getMean_points (me, tmin, tmax);
286 	longdouble sum = 0.0;
287 	for (integer i = imin; i <= imax; i ++) {
288 		double diff = my points.at [i] -> value - mean;
289 		sum += diff * diff;
290 	}
291 	return sqrt ((double) sum / (n - 1));
292 }
293 
RealTier_multiplyPart(RealTier me,double tmin,double tmax,double factor)294 void RealTier_multiplyPart (RealTier me, double tmin, double tmax, double factor) {
295 	for (integer ipoint = 1; ipoint <= my points.size; ipoint ++) {
296 		RealPoint point = my points.at [ipoint];
297 		double t = point -> number;
298 		if (t >= tmin && t <= tmax) {
299 			point -> value *= factor;
300 		}
301 	}
302 }
303 
RealTier_draw(RealTier me,Graphics g,double tmin,double tmax,double fmin,double fmax,bool garnish,conststring32 method,conststring32 quantity)304 void RealTier_draw (RealTier me, Graphics g, double tmin, double tmax, double fmin, double fmax,
305 	bool garnish, conststring32 method, conststring32 quantity)
306 {
307 	Function_unidirectionalAutowindow (me, & tmin, & tmax);
308 	const bool drawLines = str32str (method, U"lines") || str32str (method, U"Lines");
309 	const bool drawSpeckles = str32str (method, U"speckles") || str32str (method, U"Speckles");
310 	const integer n = my points.size;
311 	Graphics_setWindow (g, tmin, tmax, fmin, fmax);
312 	Graphics_setInner (g);
313 	integer imin = AnyTier_timeToHighIndex (me->asAnyTier(), tmin);
314 	integer imax = AnyTier_timeToLowIndex (me->asAnyTier(), tmax);
315 	if (n == 0) {
316 	} else if (imax < imin) {
317 		double fleft = RealTier_getValueAtTime (me, tmin);
318 		double fright = RealTier_getValueAtTime (me, tmax);
319 		if (drawLines)
320 			Graphics_line (g, tmin, fleft, tmax, fright);
321 	} else for (integer i = imin; i <= imax; i ++) {
322 		RealPoint point = my points.at [i];
323 		const double t = point -> number, f = point -> value;
324 		if (drawSpeckles)
325 			Graphics_speckle (g, t, f);
326 		if (drawLines) {
327 			if (i == 1)
328 				Graphics_line (g, tmin, f, t, f);
329 			else if (i == imin)
330 				Graphics_line (g, t, f, tmin, RealTier_getValueAtTime (me, tmin));
331 			if (i == n)
332 				Graphics_line (g, t, f, tmax, f);
333 			else if (i == imax)
334 				Graphics_line (g, t, f, tmax, RealTier_getValueAtTime (me, tmax));
335 			else {
336 				RealPoint pointRight = my points.at [i + 1];
337 				Graphics_line (g, t, f, pointRight -> number, pointRight -> value);
338 			}
339 		}
340 	}
341 	Graphics_unsetInner (g);
342 	if (garnish) {
343 		Graphics_drawInnerBox (g);
344 		Graphics_textBottom (g, true, my v_getUnitText (0, 0, 0));
345 		Graphics_marksBottom (g, 2, true, true, false);
346 		Graphics_marksLeft (g, 2, true, true, false);
347 		if (quantity)
348 			Graphics_textLeft (g, true, quantity);
349 	}
350 }
351 
RealTier_downto_TableOfReal(RealTier me,conststring32 timeLabel,conststring32 valueLabel)352 autoTableOfReal RealTier_downto_TableOfReal (RealTier me, conststring32 timeLabel, conststring32 valueLabel) {
353 	try {
354 		autoTableOfReal thee = TableOfReal_create (my points.size, 2);
355 		TableOfReal_setColumnLabel (thee.get(), 1, timeLabel);
356 		TableOfReal_setColumnLabel (thee.get(), 2, valueLabel);
357 		for (integer i = 1; i <= my points.size; i ++) {
358 			RealPoint point = my points.at [i];
359 			thy data [i] [1] = point -> number;
360 			thy data [i] [2] = point -> value;
361 		}
362 		return thee;
363 	} catch (MelderError) {
364 		Melder_throw (me, U": not converted to TableOfReal.");
365 	}
366 }
367 
RealTier_interpolateQuadratically(RealTier me,integer numberOfPointsPerParabola,int logarithmically)368 void RealTier_interpolateQuadratically (RealTier me, integer numberOfPointsPerParabola, int logarithmically) {
369 	try {
370 		autoRealTier thee = Data_copy (me);
371 		for (integer ipoint = 1; ipoint < my points.size; ipoint ++) {
372 			RealPoint point1 = my points.at [ipoint], point2 = my points.at [ipoint + 1];
373 			double time1 = point1 -> number, time2 = point2 -> number, tmid = 0.5 * (time1 + time2);
374 			double value1 = point1 -> value, value2 = point2 -> value, valuemid;
375 			double timeStep = (tmid - time1) / (numberOfPointsPerParabola + 1);
376 			if (logarithmically) {
377 				value1 = log (value1);
378 				value2 = log (value2);
379 			}
380 			valuemid = 0.5 * (value1 + value2);
381 			/*
382 				Left from the midpoint.
383 			*/
384 			for (integer inewpoint = 1; inewpoint <= numberOfPointsPerParabola; inewpoint ++) {
385 				double newTime = time1 + inewpoint * timeStep;
386 				double phase = (newTime - time1) / (tmid - time1);
387 				double newValue = value1 + (valuemid - value1) * phase * phase;
388 				if (logarithmically)
389 					newValue = exp (newValue);
390 				RealTier_addPoint (thee.get(), newTime, newValue);
391 			}
392 			/*
393 				The midpoint.
394 			*/
395 			RealTier_addPoint (thee.get(), tmid, logarithmically ? exp (valuemid) : valuemid);
396 			/*
397 				Right from the midpoint.
398 			*/
399 			for (integer inewpoint = 1; inewpoint <= numberOfPointsPerParabola; inewpoint ++) {
400 				double newTime = tmid + inewpoint * timeStep;
401 				double phase = (time2 - newTime) / (time2 - tmid);
402 				double newValue = value2 + (valuemid - value2) * phase * phase;
403 				if (logarithmically)
404 					newValue = exp (newValue);
405 				RealTier_addPoint (thee.get(), newTime, newValue);
406 			}
407 		}
408 		Thing_swap (me, thee.get());
409 	} catch (MelderError) {
410 		Melder_throw (me, U": not interpolated quadratically.");
411 	}
412 }
413 
RealTier_downto_Table(RealTier me,conststring32 indexText,conststring32 timeText,conststring32 valueText)414 autoTable RealTier_downto_Table (RealTier me, conststring32 indexText, conststring32 timeText, conststring32 valueText) {
415 	try {
416 		autoTable thee = Table_createWithoutColumnNames (my points.size,
417 			(!! indexText) + (!! timeText) + (!! valueText));
418 		integer icol = 0;
419 		if (indexText) Table_setColumnLabel (thee.get(), ++ icol, indexText);
420 		if (timeText ) Table_setColumnLabel (thee.get(), ++ icol, timeText);
421 		if (valueText) Table_setColumnLabel (thee.get(), ++ icol, valueText);
422 		for (integer ipoint = 1; ipoint <= my points.size; ipoint ++) {
423 			RealPoint point = my points.at [ipoint];
424 			icol = 0;
425 			if (indexText) Table_setNumericValue (thee.get(), ipoint, ++ icol, ipoint);
426 			if (timeText)  Table_setNumericValue (thee.get(), ipoint, ++ icol, point -> number);
427 			if (valueText) Table_setNumericValue (thee.get(), ipoint, ++ icol, point -> value);
428 		}
429 		return thee;
430 	} catch (MelderError) {
431 		Melder_throw (me, U": not converted to Table.");
432 	}
433 }
434 
Vector_to_RealTier(Vector me,integer channel,ClassInfo klas)435 autoRealTier Vector_to_RealTier (Vector me, integer channel, ClassInfo klas) {
436 	try {
437 		autoRealTier thee = RealTier_createWithClass (my xmin, my xmax, klas);
438 		for (integer i = 1; i <= my nx; i ++)
439 			RealTier_addPoint (thee.get(), Sampled_indexToX (me, i), my z [channel] [i]);
440 		return thee;
441 	} catch (MelderError) {
442 		Melder_throw (me, U": not converted to ", klas -> className, U".");
443 	}
444 }
445 
Vector_to_RealTier_peaks(Vector me,integer channel,ClassInfo klas)446 autoRealTier Vector_to_RealTier_peaks (Vector me, integer channel, ClassInfo klas) {
447 	try {
448 		autoRealTier thee = RealTier_createWithClass (my xmin, my xmax, klas);
449 		for (integer i = 2; i < my nx; i ++) {
450 			double left = my z [channel] [i - 1], centre = my z [channel] [i], right = my z [channel] [i + 1];
451 			if (left <= centre && right < centre) {
452 				double x, maximum;
453 				Vector_getMaximumAndX (me, my x1 + (i - 2.5) * my dx, my x1 + (i + 0.5) * my dx,
454 						channel, kVector_peakInterpolation :: PARABOLIC, & maximum, & x);
455 				RealTier_addPoint (thee.get(), x, maximum);
456 			}
457 		}
458 		return thee;
459 	} catch (MelderError) {
460 		Melder_throw (me, U": not converted to ", klas -> className, U" (peaks).");
461 	}
462 }
463 
Vector_to_RealTier_valleys(Vector me,integer channel,ClassInfo klas)464 autoRealTier Vector_to_RealTier_valleys (Vector me, integer channel, ClassInfo klas) {
465 	try {
466 		autoRealTier thee = RealTier_createWithClass (my xmin, my xmax, klas);
467 		for (integer i = 2; i < my nx; i ++) {
468 			double left = my z [channel] [i - 1], centre = my z [channel] [i], right = my z [channel] [i + 1];
469 			if (left >= centre && right > centre) {
470 				double x, minimum;
471 				Vector_getMinimumAndX (me, my x1 + (i - 2.5) * my dx, my x1 + (i + 0.5) * my dx,
472 						channel, kVector_peakInterpolation :: PARABOLIC, & minimum, & x);
473 				RealTier_addPoint (thee.get(), x, minimum);
474 			}
475 		}
476 		return thee;
477 	} catch (MelderError) {
478 		Melder_throw (me, U": not converted to ", klas -> className, U" (valleys).");
479 	}
480 }
481 
PointProcess_upto_RealTier(PointProcess me,double value,ClassInfo klas)482 autoRealTier PointProcess_upto_RealTier (PointProcess me, double value, ClassInfo klas) {
483 	try {
484 		autoRealTier thee = RealTier_createWithClass (my xmin, my xmax, klas);
485 		for (integer i = 1; i <= my nt; i ++)
486 			RealTier_addPoint (thee.get(), my t [i], value);
487 		return thee;
488 	} catch (MelderError) {
489 		Melder_throw (me, U": not converted to RealTier.");
490 	}
491 }
492 
RealTier_formula(RealTier me,conststring32 expression,Interpreter interpreter,RealTier thee)493 void RealTier_formula (RealTier me, conststring32 expression, Interpreter interpreter, RealTier thee) {
494 	try {
495 		Formula_compile (interpreter, me, expression, kFormula_EXPRESSION_TYPE_NUMERIC, true);
496 		Formula_Result result;
497 		if (! thee)
498 			thee = me;
499 		for (integer icol = 1; icol <= my points.size; icol ++) {
500 			Formula_run (0, icol, & result);
501 			if (isundef (result. numericResult))
502 				Melder_throw (U"Cannot put an undefined value into the tier.");
503 			thy points.at [icol] -> value = result. numericResult;
504 		}
505 	} catch (MelderError) {
506 		Melder_throw (me, U": formula not completed.");
507 	}
508 }
509 
RealTier_removePointsBelow(RealTier me,double level)510 void RealTier_removePointsBelow (RealTier me, double level) {
511 	for (integer ipoint = my points.size; ipoint > 0; ipoint --) {
512 		RealPoint point = my points.at [ipoint];
513 		if (point -> value < level)
514 			AnyTier_removePoint (me->asAnyTier(), ipoint);
515 	}
516 }
517 
RealTier_PointProcess_into_RealTier(RealTier me,PointProcess pp,RealTier thee)518 void RealTier_PointProcess_into_RealTier (RealTier me, PointProcess pp, RealTier thee) {
519 	for (integer i = 1; i <= pp -> nt; i ++) {
520 		const double time = pp -> t [i];
521 		const double value = RealTier_getValueAtTime (me, time);
522 		RealTier_addPoint (thee, time, value);
523 	}
524 }
525 
RealTier_PointProcess_to_RealTier(RealTier me,PointProcess pp)526 autoRealTier RealTier_PointProcess_to_RealTier (RealTier me, PointProcess pp) {
527 	try {
528 		if (my points.size == 0)
529 			Melder_throw (U"No points.");
530 		autoRealTier thee = RealTier_create (pp -> xmin, pp -> xmax);
531 		RealTier_PointProcess_into_RealTier (me, pp, thee.get());
532 		return thee;
533 	} catch (MelderError) {
534 		Melder_throw (me, U" & ", pp, U": not converted to RealTier.");
535 	}
536 }
537 
AnyRealTier_downto_RealTier(RealTier me)538 autoRealTier AnyRealTier_downto_RealTier (RealTier me) {
539 	try {
540 		autoRealTier thee = Thing_new (RealTier);
541 		my structRealTier :: v_copy (thee.get());
542 		return thee;
543 	} catch (MelderError) {
544 		Melder_throw (me, U": not converted to RealTier.");
545 	}
546 }
547 
RealTier_checkThatNoPointFallsOutsideDefinedTimeDomain(RealTier me)548 static void RealTier_checkThatNoPointFallsOutsideDefinedTimeDomain (RealTier me) {
549 	if (my points.size == 0) {
550 		// nothing to check
551 	} else if (my points.size == 1) {
552 		const double onlyTime = my points.at [1] -> number;
553 		if (isdefined (my xmin))
554 			Melder_require (my xmin <= onlyTime,
555 				U"The only point (at time ", onlyTime, U" seconds) falls outside the time domain, i.e. before ", my xmin, U" seconds.");
556 		if (isdefined (my xmax))
557 			Melder_require (my xmax >= onlyTime,
558 				U"The only point (at time ", onlyTime, U" seconds) falls outside the time domain, i.e. after ", my xmax, U" seconds.");
559 	} else {
560 		const double firstTime = my points.at [1] -> number;
561 		const double lastTime = my points.at [my points.size] -> number;
562 		if (isdefined (my xmin))
563 			Melder_require (my xmin <= firstTime,
564 				U"The first point (at time ", firstTime, U" seconds) falls outside the time domain, i.e. before ", my xmin, U" seconds.");
565 		if (isdefined (my xmax))
566 			Melder_require (my xmax >= lastTime,
567 				U"The last point (at time ", lastTime, U" seconds) falls outside the time domain, i.e. after ", my xmax, U" seconds.");
568 	}
569 }
570 
RealTier_fitUndefinedTimeDomainToData(RealTier me)571 static void RealTier_fitUndefinedTimeDomainToData (RealTier me) {
572 	if (my points.size == 0) {
573 		if (isundef (my xmin) && isundef (my xmax)) {
574 			my xmin = 0.0;
575 			my xmax = 1.0;
576 		} else if (isundef (my xmin)) {
577 			my xmin = my xmax - 1.0;
578 		} else if (isundef (my xmax)) {
579 			my xmax = my xmin + 1.0;
580 		}
581 	} else if (my points.size == 1) {
582 		const double onlyTime = my points.at [1] -> number;
583 		if (isundef (my xmin) && isundef (my xmax)) {
584 			my xmin = onlyTime - 1.0;
585 			my xmax = onlyTime + 1.0;
586 		} else if (isundef (my xmin)) {
587 			my xmin = onlyTime - 1.0 * ( my xmax == onlyTime );
588 		} else if (isundef (my xmax)) {
589 			my xmax = onlyTime + 1.0 * ( my xmin == onlyTime );
590 		}
591 	} else {
592 		const double firstTime = my points.at [1] -> number;
593 		const double lastTime = my points.at [my points.size] -> number;
594 		if (isundef (my xmin))
595 			my xmin = firstTime;
596 		if (isundef (my xmax))
597 			my xmax = lastTime;
598 	}
599 }
600 
Table_to_RealTier(Table me,integer timeColumn,integer valueColumn,double tmin,double tmax)601 autoRealTier Table_to_RealTier (Table me, integer timeColumn, integer valueColumn, double tmin, double tmax) {
602 	try {
603 		Melder_require (timeColumn >= 1 && timeColumn <= my numberOfColumns,
604 			U"The column number (for the times) should be between 1 and ", my numberOfColumns);
605 		Melder_require (valueColumn >= 1 && valueColumn <= my numberOfColumns,
606 			U"The column number (for the values) should be between 1 and ", my numberOfColumns);
607 		Melder_require (! (tmax <= tmin),   // NaN-safe
608 			U"The end of the time domain (", tmax, U") should be greater than the start of the time domain (", tmin, U").");
609 		autoRealTier thee = RealTier_create (tmin, tmax);
610 		Table_numericize_Assert (me, timeColumn);
611 		Table_numericize_Assert (me, valueColumn);
612 		for (integer irow = 1; irow <= my rows.size; irow ++) {
613 			TableRow row = my rows.at [irow];
614 			RealTier_addPoint (thee.get(), row -> cells [timeColumn]. number, row -> cells [valueColumn]. number);
615 		}
616 		/*
617 			At this point, all times are in sorted order and unique,
618 			because RealTier_addPoint inserts its time in order and complains if the time already exists.
619 			The data-dependent tests therefore need to be only about the time domain.
620 		*/
621 		RealTier_checkThatNoPointFallsOutsideDefinedTimeDomain (thee.get());
622 		RealTier_fitUndefinedTimeDomainToData (thee.get());
623 		return thee;
624 	} catch (MelderError) {
625 		Melder_throw (me, U": not converted to RealTier.");
626 	}
627 }
628 
Matrix_to_RealTier(Matrix me,integer timeColumn,integer valueColumn,double tmin,double tmax)629 autoRealTier Matrix_to_RealTier (Matrix me, integer timeColumn, integer valueColumn, double tmin, double tmax) {
630 	try {
631 		Melder_require (timeColumn >= 1 && timeColumn <= my nx,
632 			U"The column number (for the times) should be between 1 and ", my nx);
633 		Melder_require (valueColumn >= 1 && valueColumn <= my nx,
634 			U"The column number (for the values) should be between 1 and ", my nx);
635 		Melder_require (! (tmax <= tmin),   // NaN-safe
636 			U"The end of the time domain (", tmax, U") should be greater than the start of the time domain (", tmin, U").");
637 		autoRealTier thee = RealTier_create (tmin, tmax);
638 		for (integer irow = 1; irow <= my ny; irow ++)
639 			RealTier_addPoint (thee.get(), my z [irow] [timeColumn], my z [irow] [valueColumn]);
640 		/*
641 			At this point, all times are in sorted order and unique,
642 			because RealTier_addPoint inserts its time in order and complains if the time already exists.
643 			The data-dependent tests therefore need to be only about the time domain.
644 		*/
645 		RealTier_checkThatNoPointFallsOutsideDefinedTimeDomain (thee.get());
646 		RealTier_fitUndefinedTimeDomainToData (thee.get());
647 		return thee;
648 	} catch (MelderError) {
649 		Melder_throw (me, U": not converted to RealTier.");
650 	}
651 }
652 
653 /* End of file RealTier.cpp */
654