1 /*
2  * Grace - GRaphing, Advanced Computation and Exploration of data
3  *
4  * Home page: http://plasma-gate.weizmann.ac.il/Grace/
5  *
6  * Copyright (c) 1991-1995 Paul J Turner, Portland, OR
7  * Copyright (c) 1996-2000 Grace Development Team
8  *
9  * Maintained by Evgeny Stambulchik
10  *
11  *
12  *                           All Rights Reserved
13  *
14  *    This program is free software; you can redistribute it and/or modify
15  *    it under the terms of the GNU General Public License as published by
16  *    the Free Software Foundation; either version 2 of the License, or
17  *    (at your option) any later version.
18  *
19  *    This program is distributed in the hope that it will be useful,
20  *    but WITHOUT ANY WARRANTY; without even the implied warranty of
21  *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
22  *    GNU General Public License for more details.
23  *
24  *    You should have received a copy of the GNU General Public License
25  *    along with this program; if not, write to the Free Software
26  *    Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
27  */
28 
29 /*
30  *
31  * procedures for performing transformations from the command
32  * line interpreter and the GUI.
33  *
34  */
35 
36 #include <config.h>
37 #include <cmath.h>
38 
39 #include <stdio.h>
40 #include <stdlib.h>
41 #include <string.h>
42 
43 #include "globals.h"
44 #include "utils.h"
45 #include "draw.h"
46 #include "ssdata.h"
47 #include "graphs.h"
48 #include "parser.h"
49 #include "protos.h"
50 
51 static void forwarddiff(double *x, double *y, double *resx, double *resy, int n);
52 static void backwarddiff(double *x, double *y, double *resx, double *resy, int n);
53 static void centereddiff(double *x, double *y, double *resx, double *resy, int n);
54 int get_points_inregion(int rno, int invr, int len, double *x, double *y, int *cnt, double **xt, double **yt);
55 
56 static char buf[256];
57 
58 
do_fourier_command(int gno,int setno,int ftype,int ltype)59 void do_fourier_command(int gno, int setno, int ftype, int ltype)
60 {
61     switch (ftype) {
62     case FFT_DFT:
63 	do_fourier(gno, setno, 0, 0, ltype, 0, 0, 0);
64 	break;
65     case FFT_INVDFT       :
66 	do_fourier(gno, setno, 0, 0, ltype, 1, 0, 0);
67 	break;
68     case FFT_FFT:
69 	do_fourier(gno, setno, 1, 0, ltype, 0, 0, 0);
70 	break;
71     case FFT_INVFFT       :
72 	do_fourier(gno, setno, 1, 0, ltype, 1, 0, 0);
73 	break;
74     }
75 }
76 
77 
78 /*
79  * evaluate a formula
80  */
do_compute(int gno,int setno,int graphto,int loadto,char * rarray,char * fstr)81 int do_compute(int gno, int setno, int graphto, int loadto, char *rarray, char *fstr)
82 {
83     if (is_set_active(gno, setno)) {
84 	if (gno != graphto || setno != loadto) {
85 	    if (copysetdata(gno, setno, graphto, loadto) != RETURN_SUCCESS) {
86 	        return RETURN_FAILURE;
87             }
88         }
89 	filter_set(graphto, loadto, rarray);
90         set_parser_setno(graphto, loadto);
91         if (scanner(fstr) != RETURN_SUCCESS) {
92 	    if (graphto != gno || loadto != setno) {
93 		killset(graphto, loadto);
94 	    }
95 	    return RETURN_FAILURE;
96 	} else {
97 	    set_dirtystate();
98             return RETURN_SUCCESS;
99         }
100     } else {
101         return RETURN_FAILURE;
102     }
103 }
104 
105 /*
106  * forward, backward and centered differences
107  */
forwarddiff(double * x,double * y,double * resx,double * resy,int n)108 static void forwarddiff(double *x, double *y, double *resx, double *resy, int n)
109 {
110     int i, eflag = 0;
111     double h;
112 
113     for (i = 1; i < n; i++) {
114 	resx[i - 1] = x[i - 1];
115 	h = x[i - 1] - x[i];
116 	if (h == 0.0) {
117 	    resy[i - 1] = - MAXNUM;
118 	    eflag = 1;
119 	} else {
120 	    resy[i - 1] = (y[i - 1] - y[i]) / h;
121 	}
122     }
123     if (eflag) {
124 	errmsg("Warning: infinite slope, check set status before proceeding");
125     }
126 }
127 
backwarddiff(double * x,double * y,double * resx,double * resy,int n)128 static void backwarddiff(double *x, double *y, double *resx, double *resy, int n)
129 {
130     int i, eflag = 0;
131     double h;
132 
133     for (i = 0; i < n - 1; i++) {
134 	resx[i] = x[i];
135 	h = x[i + 1] - x[i];
136 	if (h == 0.0) {
137 	    resy[i] = - MAXNUM;
138 	    eflag = 1;
139 	} else {
140 	    resy[i] = (y[i + 1] - y[i]) / h;
141 	}
142     }
143     if (eflag) {
144 	errmsg("Warning: infinite slope, check set status before proceeding");
145     }
146 }
147 
centereddiff(double * x,double * y,double * resx,double * resy,int n)148 static void centereddiff(double *x, double *y, double *resx, double *resy, int n)
149 {
150     int i, eflag = 0;
151     double h1, h2;
152 
153     for (i = 1; i < n - 1; i++) {
154 	resx[i - 1] = x[i];
155 	h1 = x[i] - x[i - 1];
156 	h2 = x[i + 1] - x[i];
157 	if (h1 + h2 == 0.0) {
158 	    resy[i - 1] = - MAXNUM;
159 	    eflag = 1;
160 	} else {
161 	    resy[i - 1] = (y[i + 1] - y[i - 1]) / (h1 + h2);
162 	}
163     }
164     if (eflag) {
165 	errmsg("Warning: infinite slope, check set status before proceeding");
166     }
167 }
168 
seasonaldiff(double * x,double * y,double * resx,double * resy,int n,int period)169 static void seasonaldiff(double *x, double *y,
170 			 double *resx, double *resy, int n, int period)
171 {
172     int i;
173 
174     for (i = 0; i < n - period; i++) {
175 	resx[i] = x[i];
176 	resy[i] = y[i] - y[i + period];
177     }
178 }
179 
180 /*
181  * trapezoidal rule
182  */
trapint(double * x,double * y,double * resx,double * resy,int n)183 double trapint(double *x, double *y, double *resx, double *resy, int n)
184 {
185     int i;
186     double sum = 0.0;
187     double h;
188 
189     if (n < 2) {
190         return 0.0;
191     }
192 
193     if (resx != NULL) {
194         resx[0] = x[0];
195     }
196     if (resy != NULL) {
197         resy[0] = 0.0;
198     }
199     for (i = 1; i < n; i++) {
200 	h = (x[i] - x[i - 1]);
201 	if (resx != NULL) {
202 	    resx[i] = x[i];
203 	}
204 	sum = sum + h * (y[i - 1] + y[i]) * 0.5;
205 	if (resy != NULL) {
206 	    resy[i] = sum;
207 	}
208     }
209     return sum;
210 }
211 
212 /*
213  * apply a digital filter
214  */
do_digfilter(int set1,int set2)215 void do_digfilter(int set1, int set2)
216 {
217     int digfiltset;
218 
219     if (!(is_set_active(get_cg(), set1) && is_set_active(get_cg(), set2))) {
220 	errmsg("Set not active");
221 	return;
222     }
223     if ((getsetlength(get_cg(), set1) < 3) || (getsetlength(get_cg(), set2) < 3)) {
224 	errmsg("Set length < 3");
225 	return;
226     }
227     digfiltset = nextset(get_cg());
228     if (digfiltset != (-1)) {
229 	activateset(get_cg(), digfiltset);
230 	setlength(get_cg(), digfiltset, getsetlength(get_cg(), set1) - getsetlength(get_cg(), set2) + 1);
231 	sprintf(buf, "Digital filter from set %d applied to set %d", set2, set1);
232 	filterser(getsetlength(get_cg(), set1),
233 		  getx(get_cg(), set1),
234 		  gety(get_cg(), set1),
235 		  getx(get_cg(), digfiltset),
236 		  gety(get_cg(), digfiltset),
237 		  gety(get_cg(), set2),
238 		  getsetlength(get_cg(), set2));
239 	setcomment(get_cg(), digfiltset, buf);
240     }
241 }
242 
243 /*
244  * linear convolution
245  */
do_linearc(int gno1,int set1,int gno2,int set2)246 void do_linearc(int gno1, int set1, int gno2, int set2)
247 {
248     int linearcset, i, itmp, cg = get_cg();
249     double *xtmp;
250 
251     if (!(is_set_active(gno1, set1) && is_set_active(gno2, set2))) {
252 	errmsg("Set not active");
253 	return;
254     }
255     if ((getsetlength(gno1, set1) < 3) || (getsetlength(gno2, set2) < 3)) {
256 	errmsg("Set length < 3");
257 	return;
258     }
259     linearcset = nextset(cg);
260     if (linearcset != (-1)) {
261 	activateset(cg, linearcset);
262 	setlength(cg, linearcset, (itmp = getsetlength(gno1, set1) + getsetlength(gno2, set2) - 1));
263 	linearconv(gety(gno2, set2), gety(gno1, set1), gety(cg, linearcset), getsetlength(gno2, set2), getsetlength(gno1, set1));
264 	xtmp = getx(cg, linearcset);
265 	for (i = 0; i < itmp; i++) {
266 	    xtmp[i] = i;
267 	}
268 	sprintf(buf, "Linear convolution of set %d with set %d", set1, set2);
269 	setcomment(cg, linearcset, buf);
270     }
271 }
272 
273 /*
274  * cross correlation/covariance
275  */
do_xcor(int gno1,int set1,int gno2,int set2,int maxlag,int covar)276 void do_xcor(int gno1, int set1, int gno2, int set2, int maxlag, int covar)
277 {
278     int xcorset, i, ierr, len, cg = get_cg();
279     double *xtmp;
280 
281     if (!(is_set_active(gno1, set1) && is_set_active(gno2, set2))) {
282 	errmsg("Set not active");
283 	return;
284     }
285 
286     len = getsetlength(gno1, set1);
287     if (getsetlength(gno2, set2) != len) {
288 	errmsg("Sets must be of the same length");
289     }
290     if (len < 2) {
291 	errmsg("Set length < 2");
292 	return;
293     }
294 
295     if (maxlag < 1 || maxlag > len) {
296 	errmsg("Lag incorrectly specified");
297 	return;
298     }
299 
300     xcorset = nextset(cg);
301 
302     if (xcorset != (-1)) {
303 	char *fname;
304         activateset(cg, xcorset);
305 	setlength(cg, xcorset, maxlag);
306 	if (covar) {
307             fname = "covariance";
308         } else {
309             fname = "correlation";
310         }
311         if (set1 != set2) {
312 	    sprintf(buf, "X-%s of G%d.S%d and G%d.S%d at maximum lag %d",
313                     fname, gno1, set1, gno2, set2, maxlag);
314 	} else {
315 	    sprintf(buf, "Auto-%s of G%d.S%d at maximum lag %d",
316                     fname, gno1, set1, maxlag);
317 	}
318 	ierr = crosscorr(gety(gno1, set1), gety(gno2, set2), len,
319                          maxlag, covar, gety(cg, xcorset));
320 	xtmp = getx(cg, xcorset);
321 	for (i = 0; i < maxlag; i++) {
322 	    xtmp[i] = i;
323 	}
324 	setcomment(cg, xcorset, buf);
325     }
326 }
327 
328 
329 /*
330  * numerical integration
331  */
do_int(int gno,int setno,int itype)332 double do_int(int gno, int setno, int itype)
333 {
334     int intset;
335     double sum = 0;
336 
337     if (!is_set_active(gno, setno)) {
338 	errmsg("Set not active");
339 	return 0.0;
340     }
341     if (getsetlength(gno, setno) < 3) {
342 	errmsg("Set length < 3");
343 	return 0.0;
344     }
345     if (itype == 0) {
346 	intset = nextset(gno);
347 	if (intset != (-1)) {
348 	    activateset(gno, intset);
349 	    setlength(gno, intset, getsetlength(gno, setno));
350 	    sprintf(buf, "Cumulative sum of set %d", setno);
351 	    sum = trapint(getx(gno, setno), gety(gno, setno), getx(gno, intset), gety(gno, intset), getsetlength(gno, setno));
352 	    setcomment(gno, intset, buf);
353 	}
354     } else {
355 	sum = trapint(getx(gno, setno), gety(gno, setno), NULL, NULL, getsetlength(gno, setno));
356     }
357     return sum;
358 }
359 
360 /*
361  * difference a set
362  * itype means
363  *  0 - forward
364  *  1 - backward
365  *  2 - centered difference
366  */
do_differ(int gno,int setno,int itype)367 void do_differ(int gno, int setno, int itype)
368 {
369     int diffset;
370 
371     if (!is_set_active(gno, setno)) {
372 	errmsg("Set not active");
373 	return;
374     }
375     if (getsetlength(gno, setno) < 3) {
376 	errmsg("Set length < 3");
377 	return;
378     }
379     diffset = nextset(gno);
380     if (diffset != (-1)) {
381 	activateset(gno, diffset);
382 	switch (itype) {
383 	case 0:
384 	    sprintf(buf, "Forward difference of set %d", setno);
385 	    setlength(gno, diffset, getsetlength(gno, setno) - 1);
386 	    forwarddiff(getx(gno, setno), gety(gno, setno), getx(gno, diffset), gety(gno, diffset), getsetlength(gno, setno));
387 	    break;
388 	case 1:
389 	    sprintf(buf, "Backward difference of set %d", setno);
390 	    setlength(gno, diffset, getsetlength(gno, setno) - 1);
391 	    backwarddiff(getx(gno, setno), gety(gno, setno), getx(gno, diffset), gety(gno, diffset), getsetlength(gno, setno));
392 	    break;
393 	case 2:
394 	    sprintf(buf, "Centered difference of set %d", setno);
395 	    setlength(gno, diffset, getsetlength(gno, setno) - 2);
396 	    centereddiff(getx(gno, setno), gety(gno, setno), getx(gno, diffset), gety(gno, diffset), getsetlength(gno, setno));
397 	    break;
398 	}
399 	setcomment(gno, diffset, buf);
400     }
401 }
402 
403 /*
404  * seasonally difference a set
405  */
do_seasonal_diff(int setno,int period)406 void do_seasonal_diff(int setno, int period)
407 {
408     int diffset;
409 
410     if (!is_set_active(get_cg(), setno)) {
411 	errmsg("Set not active");
412 	return;
413     }
414     if (getsetlength(get_cg(), setno) < 2) {
415 	errmsg("Set length < 2");
416 	return;
417     }
418     diffset = nextset(get_cg());
419     if (diffset != (-1)) {
420 	activateset(get_cg(), diffset);
421 	setlength(get_cg(), diffset, getsetlength(get_cg(), setno) - period);
422 	seasonaldiff(getx(get_cg(), setno), gety(get_cg(), setno),
423 		     getx(get_cg(), diffset), gety(get_cg(), diffset),
424 		     getsetlength(get_cg(), setno), period);
425 	sprintf(buf, "Seasonal difference of set %d, period %d", setno, period);
426 	setcomment(get_cg(), diffset, buf);
427     }
428 }
429 
430 /*
431  * regression with restrictions to region rno if rno >= 0
432  */
do_regress(int gno,int setno,int ideg,int iresid,int rno,int invr,int fitset)433 void do_regress(int gno, int setno, int ideg, int iresid, int rno, int invr, int fitset)
434 	/*
435 	 * gno, setno  - set to perform fit on
436 	 * ideg   - degree of fit
437 	 * irisid - 0 -> whole set, 1-> subset of setno
438 	 * rno    - region number of subset
439 	 * invr   - 1->invert region, 0 -> do nothing
440 	 * fitset - set to which fitted function will be loaded
441 	 * 			Y values are computed at the x values in the set
442 	 *          if -1 is specified, a set with the same x values as the
443 	 *          one being fitted will be created and used
444 	 */
445 {
446     int len, i, sdeg = ideg;
447     int cnt = 0, fitlen = 0;
448     double *x, *y, *xt = NULL, *yt = NULL, *xr = NULL, *yr = NULL;
449     char buf[256];
450     double c[20];   /* coefficients of fitted polynomial */
451 
452     if (!is_set_active(gno, setno)) {
453 		errmsg("Set not active");
454 		return;
455     }
456     len = getsetlength(gno, setno);
457     x = getx(gno, setno);
458     y = gety(gno, setno);
459     if (rno == -1) {
460 		xt = x;
461 		yt = y;
462     } else if (isactive_region(rno)) {
463 		if (!get_points_inregion(rno, invr, len, x, y, &cnt, &xt, &yt)) {
464 			if (cnt == 0) {
465 			errmsg("No points found in region, operation cancelled");
466 			}
467 			return;
468 		}
469 		len = cnt;
470     } else {
471 		errmsg("Selected region is not active");
472 		return;
473     }
474     /*
475      * first part for polynomials, second part for linear fits to transformed
476      * data
477      */
478     if ((len < ideg && ideg <= 10) || (len < 2 && ideg > 10)) {
479 		errmsg("Too few points in set, operation cancelled");
480 		return;
481     }
482 	/* determine is set provided or use abscissa from fitted set */
483     if( fitset == -1 ) {		        /* create set */
484     	if( (fitset = nextset(gno)) != -1 ) {
485 			activateset(gno, fitset);
486 			setlength(gno, fitset, len);
487 			fitlen = len;
488 			xr = getx(gno, fitset);
489 			for (i = 0; i < len; i++) {
490 	    		xr[i] = xt[i];
491 			}
492 			yr = gety(gno, fitset);
493 		}
494     } else {									/* set has been provided */
495 		fitlen = getsetlength( gno, fitset );
496 		xr = getx(gno, fitset);
497 		yr = gety(gno, fitset);
498     }
499 
500 	/* transform data so that system has the form y' = A + B * x' */
501     if (fitset != -1) {
502 	if (ideg == 12) {	/* y=A*x^B -> ln(y) = ln(A) + B * ln(x) */
503 	    ideg = 1;
504 	    for (i = 0; i < len; i++) {
505 			if (xt[i] <= 0.0) {
506 				errmsg("One of X[i] <= 0.0");
507 				return;
508 			}
509 			if (yt[i] <= 0.0) {
510 				errmsg("One of Y[i] <= 0.0");
511 				return;
512 			}
513 	    }
514 	    for (i = 0; i < len; i++) {
515 			xt[i] = log(xt[i]);
516 			yt[i] = log(yt[i]);
517 	    }
518 		for( i=0; i<fitlen; i++ )
519 			if( xr[i] <= 0.0 ) {
520 				errmsg("One of X[i] <= 0.0");
521 				return;
522 			}
523 		for( i=0; i<fitlen; i++ )
524 			xr[i] = log( xr[i] );
525 	} else if (ideg == 13) {   /*y=Aexp(B*x) -> ln(y) = ln(A) + B * x */
526 	    ideg = 1;
527 	    for (i = 0; i < len; i++) {
528 			if (yt[i] <= 0.0) {
529 				errmsg("One of Y[i] <= 0.0");
530 				return;
531 		 	}
532 	    }
533 	    for (i = 0; i < len; i++) {
534 			yt[i] = log(yt[i]);
535 	    }
536 	} else if (ideg == 14) {	/* y = A + B * ln(x) */
537 	    ideg = 1;
538 	    for (i = 0; i < len; i++) {
539 			if (xt[i] <= 0.0) {
540 				errmsg("One of X[i] <= 0.0");
541 				return;
542 			}
543 	    }
544 		for (i = 0; i < len; i++) {
545 			xt[i] = log(xt[i]);
546 		}
547 		for( i=0; i<fitlen; i++ )
548 			if( xr[i] <= 0.0 ) {
549 				errmsg("One of X[i] <= 0.0");
550 				return;
551 			}
552 		for( i=0; i<fitlen; i++ ){
553 			xr[i] = log( xr[i] );
554 		}
555 	} else if (ideg == 15) {	/* y = 1/( A + B*x ) -> 1/y = a + B*x */
556 	    ideg = 1;
557 	    for (i = 0; i < len; i++) {
558 			if (yt[i] == 0.0) {
559 			    errmsg("One of Y[i] = 0.0");
560 			    return;
561 			}
562 	    }
563 	    for (i = 0; i < len; i++) {
564 			yt[i] = 1.0 / yt[i];
565 	    }
566 	}
567 
568 	if (fitcurve(xt, yt, len, ideg, c)) {
569 	    killset(gno, fitset);
570 	    goto bustout;
571 	}
572 
573 	/* compute function at requested x ordinates */
574 	for( i=0; i<fitlen; i++ )
575 		yr[i] = leasev( c, ideg, xr[i] );
576 
577 	sprintf(buf, "\nN.B. Statistics refer to the transformed model\n");
578 	/* apply inverse transform, output fitted function in human readable form */
579 	if( sdeg<11 ) {
580 		sprintf(buf, "\ny = %.5g %c %.5g * x", c[0], c[1]>=0?'+':'-', fabs(c[1]));
581 		for( i=2; i<=ideg; i++ )
582 			sprintf( buf+strlen(buf), " %c %.5g * x^%d", c[i]>=0?'+':'-',
583 															fabs(c[i]), i );
584 		strcat( buf, "\n" );
585 	} else if (sdeg == 12) {	/* ln(y) = ln(A) + b * ln(x) */
586 		sprintf( buf, "\ny = %.5g * x^%.5g\n", exp(c[0]), c[1] );
587 	    for (i = 0; i < len; i++) {
588 			xt[i] = exp(xt[i]);
589 			yt[i] = exp(yt[i]);
590 	    }
591 	    for (i = 0; i < fitlen; i++){
592 			yr[i] = exp(yr[i]);
593 			xr[i] = exp(xr[i]);
594 		}
595 	} else if (sdeg == 13) { 	/* ln(y) = ln(A) + B * x */
596 		sprintf( buf, "\ny = %.5g * exp( %.5g * x )\n", exp(c[0]), c[1] );
597 	    for (i = 0; i < len; i++) {
598 			yt[i] = exp(yt[i]);
599 	    }
600 	    for (i = 0; i < fitlen; i++)
601 			yr[i] = exp(yr[i]);
602 	} else if (sdeg == 14) {	/* y = A + B * ln(x) */
603 		sprintf(buf, "\ny = %.5g %c %.5g * ln(x)\n", c[0], c[1]>=0?'+':'-',
604 											fabs(c[1]) );
605 	    for (i = 0; i < len; i++) {
606 			xt[i] = exp(xt[i]);
607 	    }
608 	    for (i = 0; i < fitlen; i++)
609 			xr[i] = exp(xr[i]);
610 	} else if (sdeg == 15) {	/* y = 1/( A + B*x ) */
611 		sprintf( buf, "\ny = 1/(%.5g %c %.5g * x)\n", c[0], c[1]>=0?'+':'-',
612 													fabs(c[1]) );
613 	    for (i = 0; i < len; i++) {
614 			yt[i] = 1.0 / yt[i];
615 	    }
616 	    for (i = 0; i < fitlen; i++)
617 			yr[i] = 1.0 / yr[i];
618 	}
619 	stufftext(buf);
620 	sprintf(buf, "\nRegression of set %d results to set %d\n", setno, fitset);
621 	stufftext(buf);
622 
623 	switch (iresid) {
624 	case 1:
625 		/* determine residual */
626 	    for (i = 0; i < len; i++) {
627 			yr[i] = yt[i] - yr[i];
628 	    }
629 	    break;
630 	case 2:
631 	    break;
632 	}
633 	sprintf(buf, "%d deg fit of set %d", ideg, setno);
634 	setcomment(gno, fitset, buf);
635     }
636     bustout:;
637     if (rno >= 0 && cnt != 0) {	/* had a region and allocated memory there */
638 	xfree(xt);
639 	xfree(yt);
640     }
641 }
642 
643 /*
644  * running averages, medians, min, max, std. deviation
645  */
do_runavg(int gno,int setno,int runlen,int runtype,int rno,int invr)646 void do_runavg(int gno, int setno, int runlen, int runtype, int rno, int invr)
647 {
648     int runset;
649     int len, cnt = 0;
650     double *x, *y, *xt = NULL, *yt = NULL, *xr, *yr;
651 
652     if (!is_set_active(gno, setno)) {
653 	errmsg("Set not active");
654 	return;
655     }
656     if (runlen < 2) {
657 	errmsg("Length of running average < 2");
658 	return;
659     }
660     len = getsetlength(gno, setno);
661     x = getx(gno, setno);
662     y = gety(gno, setno);
663     if (rno == -1) {
664 	xt = x;
665 	yt = y;
666     } else if (isactive_region(rno)) {
667 	if (!get_points_inregion(rno, invr, len, x, y, &cnt, &xt, &yt)) {
668 	    if (cnt == 0) {
669 		errmsg("No points found in region, operation cancelled");
670 	    }
671 	    return;
672 	}
673 	len = cnt;
674     } else {
675 	errmsg("Selected region is not active");
676 	return;
677     }
678     if (runlen >= len) {
679 	errmsg("Length of running average > set length");
680 	goto bustout;
681     }
682     runset = nextset(gno);
683     if (runset != (-1)) {
684 	activateset(gno, runset);
685 	setlength(gno, runset, len - runlen + 1);
686 	xr = getx(gno, runset);
687 	yr = gety(gno, runset);
688 	switch (runtype) {
689 	case 0:
690 	    runavg(xt, yt, xr, yr, len, runlen);
691 	    sprintf(buf, "%d-pt. avg. on set %d ", runlen, setno);
692 	    break;
693 	case 1:
694 	    runmedian(xt, yt, xr, yr, len, runlen);
695 	    sprintf(buf, "%d-pt. median on set %d ", runlen, setno);
696 	    break;
697 	case 2:
698 	    runminmax(xt, yt, xr, yr, len, runlen, 0);
699 	    sprintf(buf, "%d-pt. min on set %d ", runlen, setno);
700 	    break;
701 	case 3:
702 	    runminmax(xt, yt, xr, yr, len, runlen, 1);
703 	    sprintf(buf, "%d-pt. max on set %d ", runlen, setno);
704 	    break;
705 	case 4:
706 	    runstddev(xt, yt, xr, yr, len, runlen);
707 	    sprintf(buf, "%d-pt. std dev., set %d ", runlen, setno);
708 	    break;
709 	}
710 	setcomment(gno, runset, buf);
711     }
712   bustout:;
713     if (rno >= 0 && cnt != 0) {	/* had a region and allocated memory there */
714 	xfree(xt);
715 	xfree(yt);
716     }
717 }
718 
719 /*
720  * DFT by FFT or definition
721  */
do_fourier(int gno,int setno,int fftflag,int load,int loadx,int invflag,int type,int wind)722 void do_fourier(int gno, int setno, int fftflag, int load, int loadx, int invflag, int type, int wind)
723 {
724     int i, ilen;
725     double *x, *y, *xx, *yy, delt, T;
726     int i2 = 0, specset;
727 
728     if (!is_set_active(get_cg(), setno)) {
729 	errmsg("Set not active");
730 	return;
731     }
732     ilen = getsetlength(get_cg(), setno);
733     if (ilen < 2) {
734 	errmsg("Set length < 2");
735 	return;
736     }
737     if (fftflag) {
738 	if ((i2 = ilog2(ilen)) <= 0) {
739 	    errmsg("Set length not a power of 2");
740 	    return;
741 	}
742     }
743     specset = nextset(get_cg());
744     if (specset != -1) {
745 	activateset(get_cg(), specset);
746 	setlength(get_cg(), specset, ilen);
747 	xx = getx(get_cg(), specset);
748 	yy = gety(get_cg(), specset);
749 	x = getx(get_cg(), setno);
750 	y = gety(get_cg(), setno);
751 	copyx(get_cg(), setno, specset);
752 	copyy(get_cg(), setno, specset);
753 	if (wind != 0) {	/* apply data window if needed */
754 	    apply_window(xx, yy, ilen, type, wind);
755 	}
756 	if (type == 0) {	/* real data */
757 	    for (i = 0; i < ilen; i++) {
758 		xx[i] = yy[i];
759 		yy[i] = 0.0;
760 	    }
761 	}
762 	if (fftflag) {
763 	    fft(xx, yy, ilen, i2, invflag);
764 	} else {
765 	    dft(xx, yy, ilen, invflag);
766 	}
767 	switch (load) {
768 	case 0:
769 	    delt = (x[ilen-1] - x[0])/(ilen -1.0);
770 	    T = (x[ilen - 1] - x[0]);
771 	    xx = getx(get_cg(), specset);
772 	    yy = gety(get_cg(), specset);
773 	    for (i = 0; i < ilen / 2; i++) {
774 	      /* carefully get amplitude of complex xform:
775 		 use abs(a[i]) + abs(a[-i]) except for zero term */
776 	      if(i) yy[i] = hypot(xx[i], yy[i])+hypot(xx[ilen-i], yy[ilen-i]);
777 	      else yy[i]=2*fabs(xx[i]);
778 		switch (loadx) {
779 		case 0:
780 		    xx[i] = i;
781 		    break;
782 		case 1:
783 		    /* xx[i] = 2.0 * M_PI * i / ilen; */
784 		    xx[i] = i / T;
785 		    break;
786 		case 2:
787 		    if (i == 0) {
788 			xx[i] = T + delt;	/* the mean */
789 		    } else {
790 			/* xx[i] = (double) ilen / (double) i; */
791 			xx[i] = T / i;
792 		    }
793 		    break;
794 		}
795 	    }
796 	    setlength(get_cg(), specset, ilen / 2);
797 	    break;
798 	case 1:
799 	    delt = (x[ilen-1] - x[0])/(ilen -1.0);
800 	    T = (x[ilen - 1] - x[0]);
801 	    setlength(get_cg(), specset, ilen / 2);
802 	    xx = getx(get_cg(), specset);
803 	    yy = gety(get_cg(), specset);
804 	    for (i = 0; i < ilen / 2; i++) {
805 		yy[i] = -atan2(yy[i], xx[i]);
806 		switch (loadx) {
807 		case 0:
808 		    xx[i] = i;
809 		    break;
810 		case 1:
811 		    /* xx[i] = 2.0 * M_PI * i / ilen; */
812 		    xx[i] = i / T;
813 		    break;
814 		case 2:
815 		    if (i == 0) {
816 			xx[i] = T + delt;
817 		    } else {
818 			/* xx[i] = (double) ilen / (double) i; */
819 			xx[i] = T / i;
820 		    }
821 		    break;
822 		}
823 	    }
824 	    break;
825 	}
826 	if (fftflag) {
827 	    sprintf(buf, "FFT of set %d", setno);
828 	} else {
829 	    sprintf(buf, "DFT of set %d", setno);
830 	}
831 	setcomment(get_cg(), specset, buf);
832     }
833 }
834 
835 /*
836  * Apply a window to a set, result goes to a new set.
837  */
do_window(int setno,int type,int wind)838 void do_window(int setno, int type, int wind)
839 {
840     int ilen;
841     double *xx, *yy;
842     int specset;
843 
844     if (!is_set_active(get_cg(), setno)) {
845 	errmsg("Set not active");
846 	return;
847     }
848     ilen = getsetlength(get_cg(), setno);
849     if (ilen < 2) {
850 	errmsg("Set length < 2");
851 	return;
852     }
853     specset = nextset(get_cg());
854     if (specset != -1) {
855 	char *wtype[6];
856 	wtype[0] = "Triangular";
857 	wtype[1] = "Hanning";
858 	wtype[2] = "Welch";
859 	wtype[3] = "Hamming";
860 	wtype[4] = "Blackman";
861 	wtype[5] = "Parzen";
862 
863 	activateset(get_cg(), specset);
864 	setlength(get_cg(), specset, ilen);
865 	xx = getx(get_cg(), specset);
866 	yy = gety(get_cg(), specset);
867 	copyx(get_cg(), setno, specset);
868 	copyy(get_cg(), setno, specset);
869 	if (wind != 0) {
870 	    apply_window(xx, yy, ilen, type, wind);
871 	    sprintf(buf, "%s windowed set %d", wtype[wind - 1], setno);
872 	} else {		/* shouldn't happen */
873 	}
874 	setcomment(get_cg(), specset, buf);
875     }
876 }
877 
apply_window(double * xx,double * yy,int ilen,int type,int wind)878 void apply_window(double *xx, double *yy, int ilen, int type, int wind)
879 {
880     int i;
881 
882     for (i = 0; i < ilen; i++) {
883 	switch (wind) {
884 	case 1:		/* triangular */
885 	    if (type != 0) {
886 		xx[i] *= 1.0 - fabs((i - 0.5 * (ilen - 1.0)) / (0.5 * (ilen - 1.0)));
887 	    }
888 	    yy[i] *= 1.0 - fabs((i - 0.5 * (ilen - 1.0)) / (0.5 * (ilen - 1.0)));
889 
890 	    break;
891 	case 2:		/* Hanning */
892 	    if (type != 0) {
893 		xx[i] = xx[i] * (0.5 - 0.5 * cos(2.0 * M_PI * i / (ilen - 1.0)));
894 	    }
895 	    yy[i] = yy[i] * (0.5 - 0.5 * cos(2.0 * M_PI * i / (ilen - 1.0)));
896 	    break;
897 	case 3:		/* Welch (from Numerical Recipes) */
898 	    if (type != 0) {
899 		xx[i] *= 1.0 - pow((i - 0.5 * (ilen - 1.0)) / (0.5 * (ilen + 1.0)), 2.0);
900 	    }
901 	    yy[i] *= 1.0 - pow((i - 0.5 * (ilen - 1.0)) / (0.5 * (ilen + 1.0)), 2.0);
902 	    break;
903 	case 4:		/* Hamming */
904 	    if (type != 0) {
905 		xx[i] = xx[i] * (0.54 - 0.46 * cos(2.0 * M_PI * i / (ilen - 1.0)));
906 	    }
907 	    yy[i] = yy[i] * (0.54 - 0.46 * cos(2.0 * M_PI * i / (ilen - 1.0)));
908 	    break;
909 	case 5:		/* Blackman */
910 	    if (type != 0) {
911 		xx[i] = xx[i] * (0.42 - 0.5 * cos(2.0 * M_PI * i / (ilen - 1.0)) + 0.08 * cos(4.0 * M_PI * i / (ilen - 1.0)));
912 	    }
913 	    yy[i] = yy[i] * (0.42 - 0.5 * cos(2.0 * M_PI * i / (ilen - 1.0)) + 0.08 * cos(4.0 * M_PI * i / (ilen - 1.0)));
914 	    break;
915 	case 6:		/* Parzen (from Numerical Recipes) */
916 	    if (type != 0) {
917 		xx[i] *= 1.0 - fabs((i - 0.5 * (ilen - 1)) / (0.5 * (ilen + 1)));
918 	    }
919 	    yy[i] *= 1.0 - fabs((i - 0.5 * (ilen - 1)) / (0.5 * (ilen + 1)));
920 	    break;
921 	}
922     }
923 }
924 
925 
926 /*
927  * histograms
928  */
do_histo(int fromgraph,int fromset,int tograph,int toset,double * bins,int nbins,int cumulative,int normalize)929 int do_histo(int fromgraph, int fromset, int tograph, int toset,
930 	      double *bins, int nbins, int cumulative, int normalize)
931 {
932     int i, ndata;
933     int *hist;
934     double *x, *y, *data;
935     plotarr p;
936 
937     if (!is_set_active(fromgraph, fromset)) {
938 	errmsg("Set not active");
939 	return RETURN_FAILURE;
940     }
941     if (nbins <= 0) {
942 	errmsg("Number of bins <= 0");
943 	return RETURN_FAILURE;
944     }
945     if (toset == SET_SELECT_NEXT) {
946 	toset = nextset(tograph);
947     }
948     if (!is_valid_setno(tograph, toset)) {
949 	errmsg("Can't activate destination set");
950         return RETURN_FAILURE;
951     }
952 
953     ndata = getsetlength(fromgraph, fromset);
954     data = gety(fromgraph, fromset);
955 
956     hist = xmalloc(nbins*SIZEOF_INT);
957     if (hist == NULL) {
958         errmsg("xmalloc failed in do_histo()");
959         return RETURN_FAILURE;
960     }
961 
962     if (histogram(ndata, data, nbins, bins, hist) == RETURN_FAILURE){
963         xfree(hist);
964         return RETURN_FAILURE;
965     }
966 
967     activateset(tograph, toset);
968     setlength(tograph, toset, nbins + 1);
969     x = getx(tograph, toset);
970     y = gety(tograph, toset);
971 
972     x[0] = bins[0];
973     y[0] = 0.0;
974     for (i = 1; i < nbins + 1; i++) {
975         x[i] = bins[i];
976         y[i] = hist[i - 1];
977         if (cumulative) {
978             y[i] += y[i - 1];
979         }
980     }
981 
982     if (normalize) {
983         for (i = 1; i < nbins + 1; i++) {
984             double factor;
985             if (cumulative) {
986                 factor = 1.0/ndata;
987             } else {
988                 factor = 1.0/((bins[i] - bins[i - 1])*ndata);
989             }
990             y[i] *= factor;
991         }
992     }
993 
994     xfree(hist);
995 
996     get_graph_plotarr(tograph, toset, &p);
997     p.sym = SYM_NONE;
998     p.linet = LINE_TYPE_LEFTSTAIR;
999     p.dropline = TRUE;
1000     p.baseline = FALSE;
1001     p.baseline_type = BASELINE_TYPE_0;
1002     p.lines = 1;
1003     p.symlines = 1;
1004     sprintf(p.comments, "Histogram from G%d.S%d", fromgraph, fromset);
1005     set_graph_plotarr(tograph, toset, &p);
1006 
1007     return RETURN_SUCCESS;
1008 }
1009 
histogram(int ndata,double * data,int nbins,double * bins,int * hist)1010 int histogram(int ndata, double *data, int nbins, double *bins, int *hist)
1011 {
1012     int i, j, bsign;
1013     double minval, maxval;
1014 
1015     if (nbins < 1) {
1016         errmsg("Number of bins < 1");
1017         return RETURN_FAILURE;
1018     }
1019 
1020     bsign = monotonicity(bins, nbins + 1, TRUE);
1021     if (bsign == 0) {
1022         errmsg("Non-monotonic bins");
1023         return RETURN_FAILURE;
1024     }
1025 
1026     for (i = 0; i < nbins; i++) {
1027         hist[i] = 0;
1028     }
1029 
1030     /* TODO: binary search */
1031     for (i = 0; i < ndata; i++) {
1032         for (j = 0; j < nbins; j++) {
1033             if (bsign > 0) {
1034                 minval = bins[j];
1035                 maxval = bins[j + 1];
1036             } else {
1037                 minval = bins[j + 1];
1038                 maxval = bins[j];
1039             }
1040             if (data[i] >= minval && data[i] <= maxval) {
1041                 hist[j] += 1;
1042                 break;
1043             }
1044         }
1045     }
1046     return RETURN_SUCCESS;
1047 }
1048 
1049 
1050 /*
1051  * sample a set, by start/step or logical expression
1052  */
do_sample(int setno,int typeno,char * exprstr,int startno,int stepno)1053 void do_sample(int setno, int typeno, char *exprstr, int startno, int stepno)
1054 {
1055     int len, npts = 0, i, resset;
1056     double *x, *y;
1057     int reslen;
1058     double *result;
1059     int gno = get_cg();
1060 
1061     if (!is_set_active(gno, setno)) {
1062 	errmsg("Set not active");
1063 	return;
1064     }
1065 
1066     len = getsetlength(gno, setno);
1067 
1068     resset = nextset(gno);
1069     if (resset < 0) {
1070 	return;
1071     }
1072 
1073     x = getx(gno, setno);
1074     y = gety(gno, setno);
1075 
1076     if (typeno == 0) {
1077 	if (len <= 2) {
1078 	    errmsg("Set has <= 2 points");
1079 	    return;
1080 	}
1081 	if (startno < 1) {
1082 	    errmsg("Start point < 1 (locations in sets are numbered starting from 1)");
1083 	    return;
1084 	}
1085 	if (stepno < 1) {
1086 	    errmsg("Step < 1");
1087 	    return;
1088 	}
1089 	for (i = startno - 1; i < len; i += stepno) {
1090 	    add_point(gno, resset, x[i], y[i]);
1091 	    npts++;
1092 	}
1093 	sprintf(buf, "Sample, %d, %d set #%d", startno, stepno, setno);
1094     } else {
1095         if (set_parser_setno(gno, setno) != RETURN_SUCCESS) {
1096 	    errmsg("Bad set");
1097             killset(gno, resset);
1098 	    return;
1099         }
1100         if (v_scanner(exprstr, &reslen, &result) != RETURN_SUCCESS) {
1101             killset(gno, resset);
1102 	    return;
1103         }
1104         if (reslen != len) {
1105 	    errmsg("Internal error");
1106             killset(gno, resset);
1107 	    return;
1108         }
1109 
1110         npts = 0;
1111 	sprintf(buf, "Sample from %d, using '%s'", setno, exprstr);
1112 	for (i = 0; i < len; i++) {
1113 	    if ((int) rint(result[i])) {
1114 		add_point(gno, resset, x[i], y[i]);
1115 		npts++;
1116 	    }
1117 	}
1118         xfree(result);
1119     }
1120     if (npts > 0) {
1121 	setcomment(gno, resset, buf);
1122     }
1123 }
1124 
1125 #define prune_xconv(res,x,xtype)	\
1126     switch (deltatypeno) {		\
1127     case PRUNE_VIEWPORT:		\
1128 	res = xy_xconv(x);			\
1129 	break;				\
1130     case PRUNE_WORLD:			\
1131 	switch (xtype) {		\
1132 	case PRUNE_LIN:			\
1133 	    res = x;			\
1134 	    break;			\
1135 	case PRUNE_LOG:			\
1136 	    res = log(x);		\
1137 	    break;			\
1138 	}				\
1139     }
1140 
1141 #define prune_yconv(res,y,ytype)	\
1142     switch (deltatypeno) {		\
1143     case PRUNE_VIEWPORT:		\
1144 	res = xy_yconv(y);			\
1145 	break;				\
1146     case PRUNE_WORLD:			\
1147 	switch (ytype) {		\
1148 	case PRUNE_LIN:			\
1149 	    res = y;			\
1150 	    break;			\
1151 	case PRUNE_LOG:			\
1152 	    res = log(y);		\
1153 	    break;			\
1154 	}				\
1155     }
1156 
1157 /*
1158  * Prune data
1159  */
do_prune(int setno,int typeno,int deltatypeno,float deltax,float deltay,int dxtype,int dytype)1160 void do_prune(int setno, int typeno, int deltatypeno, float deltax, float deltay, int dxtype, int dytype)
1161 {
1162     int len, npts = 0, d, i, j, k, drop, resset;
1163     double *x, *y, *resx, *resy, xtmp, ytmp, ddx = 0.0, ddy = 0.0;
1164     double xj = 0.0, xjm = 0.0, xjp = 0.0, yj = 0.0, yjm = 0.0, yjp = 0.0;
1165 
1166     if (!is_set_active(get_cg(), setno)) {
1167         errmsg("Set not active");
1168         return;
1169     }
1170     len = getsetlength(get_cg(), setno);
1171     if (len <= 2) {
1172 	errmsg("Set has <= 2 points");
1173 	return;
1174     }
1175     x = getx(get_cg(), setno);
1176     y = gety(get_cg(), setno);
1177     switch (typeno) {
1178     case PRUNE_CIRCLE:
1179     case PRUNE_ELLIPSE:
1180     case PRUNE_RECTANGLE:
1181 	deltax = fabs(deltax);
1182 	if (deltax == 0)
1183 	    return;
1184 	break;
1185     }
1186     switch (typeno) {
1187     case PRUNE_CIRCLE:
1188 	deltay = deltax;
1189 	break;
1190     case PRUNE_ELLIPSE:
1191     case PRUNE_RECTANGLE:
1192     case PRUNE_INTERPOLATION:
1193 	deltay = fabs(deltay);
1194 	if (deltay == 0)
1195 	    return;
1196 	break;
1197     }
1198     if (deltatypeno == PRUNE_WORLD) {
1199 	if (dxtype == PRUNE_LOG && deltax < 1.0) {
1200 	    deltax = 1.0 / deltax;
1201 	}
1202 	if (dytype == PRUNE_LOG && deltay < 1.0) {
1203 	    deltay = 1.0 / deltay;
1204 	}
1205     }
1206     resset = nextset(get_cg());
1207     if (resset < 0) {
1208         return;
1209     }
1210     add_point(get_cg(), resset, x[0], y[0]);
1211     npts++;
1212     resx = getx(get_cg(), resset);
1213     resy = gety(get_cg(), resset);
1214     switch (typeno) {
1215     case PRUNE_CIRCLE:
1216     case PRUNE_ELLIPSE:
1217 	for (i = 1; i < len; i++) {
1218 	    xtmp = x[i];
1219 	    ytmp = y[i];
1220 	    drop = FALSE;
1221 	    for (j = npts - 1; j >= 0 && drop == FALSE; j--) {
1222 		switch (deltatypeno) {
1223 		case PRUNE_VIEWPORT:
1224 		    ddx = (xy_xconv(xtmp) - xy_xconv(resx[j])) / deltax;
1225 		    if (fabs(ddx) < 1.0) {
1226 			ddy = (xy_yconv(ytmp) - xy_yconv(resy[j])) / deltay;
1227 			if (ddx * ddx + ddy * ddy < 1.0) {
1228 			    drop = TRUE;
1229 			}
1230 		    }
1231 		    break;
1232 		case PRUNE_WORLD:
1233 		    switch (dxtype) {
1234 		    case PRUNE_LIN:
1235 			ddx = (xtmp - resx[j]) / deltax;
1236 			break;
1237 		    case PRUNE_LOG:
1238 			ddx = (xtmp / resx[j]);
1239 			if (ddx < 1.0) {
1240 			    ddx = 1.0 / ddx;
1241 			}
1242 			ddx /= deltax;
1243 			break;
1244 		    }
1245 		    if (fabs(ddx) < 1.0) {
1246 			switch (dytype) {
1247 			case PRUNE_LIN:
1248 			    ddy = (ytmp - resy[j]) / deltay;
1249 			    break;
1250 			case PRUNE_LOG:
1251 			    ddy = (ytmp / resy[j]);
1252 			    if (ddy < 1.0) {
1253 				ddy = 1.0 / ddy;
1254 			    }
1255 			    ddy /= deltay;
1256 			    break;
1257 			}
1258 			if (ddx * ddx + ddy * ddy < 1.0) {
1259 			    drop = TRUE;
1260 			}
1261 		    }
1262 		    break;
1263 		}
1264 	    }
1265 	    if (drop == FALSE) {
1266 		add_point(get_cg(), resset, xtmp, ytmp);
1267 		npts++;
1268 		resx = getx(get_cg(), resset);
1269 		resy = gety(get_cg(), resset);
1270 	    }
1271 	}
1272 	sprintf(buf, "Prune from %d, %s dx = %g dy = %g", setno,
1273 	    (typeno == 0) ? "Circle" : "Ellipse", deltax, deltay);
1274 	break;
1275     case PRUNE_RECTANGLE:
1276 	for (i = 1; i < len; i++) {
1277 	    xtmp = x[i];
1278 	    ytmp = y[i];
1279 	    drop = FALSE;
1280 	    for (j = npts - 1; j >= 0 && drop == FALSE; j--) {
1281 		switch (deltatypeno) {
1282 		case PRUNE_VIEWPORT:
1283 		    ddx = fabs(xy_xconv(xtmp) - xy_xconv(resx[j]));
1284 		    if (ddx < deltax) {
1285 			ddy = fabs(xy_yconv(ytmp) - xy_yconv(resy[j]));
1286 			if (ddy < deltay) {
1287 			    drop = TRUE;
1288 			}
1289 		    }
1290 		    break;
1291 		case PRUNE_WORLD:
1292 		    switch (dxtype) {
1293 		    case PRUNE_LIN:
1294 			ddx = fabs(xtmp - resx[j]);
1295 			break;
1296 		    case PRUNE_LOG:
1297 			ddx = (xtmp / resx[j]);
1298 			if (ddx < 1.0) {
1299 			    ddx = 1.0 / ddx;
1300 			}
1301 			break;
1302 		    }
1303 		    if (ddx < deltax) {
1304 			switch (dytype) {
1305 			case PRUNE_LIN:
1306 			    ddy = fabs(ytmp - resy[j]);
1307 			    break;
1308 			case PRUNE_LOG:
1309 			    ddy = (ytmp / resy[j]);
1310 			    if (ddy < 1.0) {
1311 				ddy = 1.0 / ddy;
1312 			    }
1313 			    break;
1314 			}
1315 			if (ddy < deltay) {
1316 			    drop = TRUE;
1317 			}
1318 		    }
1319 		    break;
1320 		}
1321 	    }
1322 	    if (drop == FALSE) {
1323 		add_point(get_cg(), resset, xtmp, ytmp);
1324 		npts++;
1325 		resx = getx(get_cg(), resset);
1326 		resy = gety(get_cg(), resset);
1327 	    }
1328 	}
1329 	sprintf(buf, "Prune from %d, %s dx = %g dy = %g", setno,
1330 	    "Rectangle", deltax, deltay);
1331 	break;
1332     case PRUNE_INTERPOLATION:
1333 	k = 0;
1334 	prune_xconv(xjm, x[0], dxtype);
1335 	prune_yconv(yjm, y[0], dytype);
1336 	while (k < len - 2) {
1337 	    d = 1;
1338 	    i = k + d + 1;
1339 	    drop = TRUE;
1340 	    while (TRUE) {
1341 		prune_xconv(xjp, x[i], dxtype);
1342 		prune_yconv(yjp, y[i], dytype);
1343 		for (j = k + 1; j < i && drop == TRUE; j++) {
1344 		    prune_xconv(xj, x[j], dxtype);
1345 		    prune_yconv(yj, y[j], dytype);
1346 		    if (xjp == xjm) {
1347 			ytmp = 0.5 * (yjp + yjm);
1348 		    } else {
1349 			ytmp = (yjp*(xj-xjm)+yjm*(xjp-xj))/(xjp-xjm);
1350 		    }
1351 		    switch (deltatypeno) {
1352 		    case PRUNE_VIEWPORT:
1353 			ddy = fabs(ytmp - yj);
1354 			break;
1355 		    case PRUNE_WORLD:
1356 			switch (dytype) {
1357 			case PRUNE_LIN:
1358 			    ddy = fabs(ytmp - yj);
1359 			    break;
1360 			case PRUNE_LOG:
1361 			    ddy = exp(fabs(ytmp - yj));
1362 			    break;
1363 			}
1364 		    }
1365 		    if (ddy > deltay) {
1366 			drop = FALSE;
1367 		    }
1368 		}
1369 		if (drop == FALSE || i == len - 1) {
1370 		    break;
1371 		}
1372 		d *= 2;
1373 		i = k + d + 1;
1374 		if (i >= len) {
1375 		    i = len - 1;
1376 		}
1377 	    }
1378 	    if (drop == FALSE) {
1379 		i = k + 1;
1380 		drop = TRUE;
1381 		while (d > 1) {
1382 		    d /= 2;
1383 		    i += d;
1384 		    prune_xconv(xjp, x[i], dxtype);
1385 		    prune_yconv(yjp, y[i], dytype);
1386 		    drop = TRUE;
1387 		    for (j = k + 1; j < i && drop == TRUE; j++) {
1388 			prune_xconv(xj, x[j], dxtype);
1389 			prune_yconv(yj, y[j], dytype);
1390 			ytmp = (yjp*(xj-xjm)+yjm*(xjp-xj))/(xjp-xjm);
1391 			switch (deltatypeno) {
1392 			case PRUNE_VIEWPORT:
1393 			    ddy = fabs(ytmp - yj);
1394 			    break;
1395 			case PRUNE_WORLD:
1396 			    switch (dytype) {
1397 			    case PRUNE_LIN:
1398 				ddy = fabs(ytmp - yj);
1399 				break;
1400 			    case PRUNE_LOG:
1401 				ddy = exp(fabs(ytmp - yj));
1402 				break;
1403 			    }
1404 			}
1405 			if (ddy > deltay) {
1406 			    drop = FALSE;
1407 			}
1408 		    }
1409 		    if (drop == FALSE) {
1410 			i -= d;
1411 		    }
1412 		}
1413 	    }
1414 	    k = i;
1415 	    prune_xconv(xjm, x[k], dxtype);
1416 	    prune_yconv(yjm, y[k], dytype);
1417 	    add_point(get_cg(), resset, x[k], y[k]);
1418 	    npts++;
1419 	    resx = getx(get_cg(), resset);
1420 	    resy = gety(get_cg(), resset);
1421 	}
1422 	if (k == len - 2) {
1423 	    add_point(get_cg(), resset, x[len-1], y[len-1]);
1424 	    npts++;
1425 	}
1426 	sprintf(buf, "Prune from %d, %s dy = %g", setno,
1427 	    "Interpolation", deltay);
1428 	break;
1429     }
1430     setcomment(get_cg(), resset, buf);
1431 }
1432 
get_points_inregion(int rno,int invr,int len,double * x,double * y,int * cnt,double ** xt,double ** yt)1433 int get_points_inregion(int rno, int invr, int len, double *x, double *y, int *cnt, double **xt, double **yt)
1434 {
1435     int i, clen = 0;
1436     double *xtmp, *ytmp;
1437     *cnt = 0;
1438     if (isactive_region(rno)) {
1439 	for (i = 0; i < len; i++) {
1440 	    if (invr) {
1441 		if (!inregion(rno, x[i], y[i])) {
1442 		    clen++;
1443 		}
1444 	    } else {
1445 		if (inregion(rno, x[i], y[i])) {
1446 		    clen++;
1447 		}
1448 	    }
1449 	}
1450 	if (clen == 0) {
1451 	    return 0;
1452 	}
1453 	xtmp = (double *) xcalloc(clen, SIZEOF_DOUBLE);
1454 	if (xtmp == NULL) {
1455 	    return 0;
1456 	}
1457 	ytmp = (double *) xcalloc(clen, SIZEOF_DOUBLE);
1458 	if (ytmp == NULL) {
1459 	    xfree(xtmp);
1460 	    return 0;
1461 	}
1462 	clen = 0;
1463 	for (i = 0; i < len; i++) {
1464 	    if (invr) {
1465 		if (!inregion(rno, x[i], y[i])) {
1466 		    xtmp[clen] = x[i];
1467 		    ytmp[clen] = y[i];
1468 		    clen++;
1469 		}
1470 	    } else {
1471 		if (inregion(rno, x[i], y[i])) {
1472 		    xtmp[clen] = x[i];
1473 		    ytmp[clen] = y[i];
1474 		    clen++;
1475 		}
1476 	    }
1477 	}
1478     } else {
1479 	return 0;
1480     }
1481     *cnt = clen;
1482     *xt = xtmp;
1483     *yt = ytmp;
1484     return 1;
1485 }
1486 
interpolate(double * mesh,double * yint,int meshlen,double * x,double * y,int len,int method)1487 int interpolate(double *mesh, double *yint, int meshlen,
1488     double *x, double *y, int len, int method)
1489 {
1490     double *b, *c, *d;
1491     double dx;
1492     int i, ifound;
1493     int m;
1494 
1495     /* For linear interpolation, non-strict monotonicity is fine */
1496     m = monotonicity(x, len, method == INTERP_LINEAR ? FALSE:TRUE);
1497     if (m == 0) {
1498         errmsg("Can't interpolate a set with non-monotonic abscissas");
1499         return RETURN_FAILURE;
1500     }
1501 
1502     switch (method) {
1503     case INTERP_SPLINE:
1504     case INTERP_ASPLINE:
1505         b = xcalloc(len, SIZEOF_DOUBLE);
1506         c = xcalloc(len, SIZEOF_DOUBLE);
1507         d = xcalloc(len, SIZEOF_DOUBLE);
1508         if (b == NULL || c == NULL || d == NULL) {
1509             xfree(b);
1510             xfree(c);
1511             xfree(d);
1512             return RETURN_FAILURE;
1513         }
1514         if (method == INTERP_ASPLINE){
1515             /* Akima spline */
1516             aspline(len, x, y, b, c, d);
1517         } else {
1518             /* Plain cubic spline */
1519             spline(len, x, y, b, c, d);
1520         }
1521 
1522 	seval(mesh, yint, meshlen, x, y, b, c, d, len);
1523 
1524         xfree(b);
1525         xfree(c);
1526         xfree(d);
1527         break;
1528     default:
1529         /* linear interpolation */
1530 
1531         for (i = 0; i < meshlen; i++) {
1532             ifound = find_span_index(x, len, m, mesh[i]);
1533             if (ifound < 0) {
1534                 ifound = 0;
1535             } else if (ifound > len - 2) {
1536                 ifound = len - 2;
1537             }
1538             dx = x[ifound + 1] - x[ifound];
1539             if (dx != 0.0) {
1540                 yint[i] = y[ifound] + (mesh[i] - x[ifound])*
1541                     ((y[ifound + 1] - y[ifound])/dx);
1542             } else {
1543                 yint[i] = (y[ifound] + y[ifound + 1])/2;
1544             }
1545         }
1546         break;
1547     }
1548 
1549     return RETURN_SUCCESS;
1550 }
1551 
1552 /* interpolate a set at abscissas from mesh
1553  * method - type of spline (or linear interpolation)
1554  * if strict is set, perform interpolation only within source set bounds
1555  * (i.e., no extrapolation)
1556  */
do_interp(int gno_src,int setno_src,int gno_dest,int setno_dest,double * mesh,int meshlen,int method,int strict)1557 int do_interp(int gno_src, int setno_src, int gno_dest, int setno_dest,
1558     double *mesh, int meshlen, int method, int strict)
1559 {
1560     int len, n, ncols;
1561     double *x, *xint;
1562     char *s;
1563 
1564     if (!is_valid_setno(gno_src, setno_src)) {
1565 	errmsg("Interpolated set not active");
1566 	return RETURN_FAILURE;
1567     }
1568     if (mesh == NULL || meshlen < 1) {
1569         errmsg("NULL sampling mesh");
1570         return RETURN_FAILURE;
1571     }
1572 
1573     len = getsetlength(gno_src, setno_src);
1574     ncols = dataset_cols(gno_src, setno_src);
1575 
1576     if (setno_dest == SET_SELECT_NEXT) {
1577 	setno_dest = nextset(gno_dest);
1578     }
1579     if (!is_valid_setno(gno_dest, setno_dest)) {
1580 	errmsg("Can't activate destination set");
1581 	return RETURN_FAILURE;
1582     }
1583 
1584     if (dataset_cols(gno_dest, setno_dest) != ncols) {
1585         copyset(gno_src, setno_src, gno_dest, setno_dest);
1586     }
1587 
1588     setlength(gno_dest, setno_dest, meshlen);
1589     activateset(gno_dest, setno_dest);
1590 
1591     x = getcol(gno_src, setno_src, DATA_X);
1592     for (n = 1; n < ncols; n++) {
1593         int res;
1594         double *y, *yint;
1595 
1596         y    = getcol(gno_src, setno_src, n);
1597         yint = getcol(gno_dest, setno_dest, n);
1598 
1599         res = interpolate(mesh, yint, meshlen, x, y, len, method);
1600         if (res != RETURN_SUCCESS) {
1601             killset(gno_dest, setno_dest);
1602             return RETURN_FAILURE;
1603         }
1604     }
1605 
1606     xint = getcol(gno_dest, setno_dest, DATA_X);
1607     memcpy(xint, mesh, meshlen*SIZEOF_DOUBLE);
1608 
1609     if (strict) {
1610         double xmin, xmax;
1611         int i, imin, imax;
1612         minmax(x, len, &xmin, &xmax, &imin, &imax);
1613         for (i = meshlen - 1; i >= 0; i--) {
1614             if (xint[i] < xmin || xint[i] > xmax) {
1615                 del_point(gno_dest, setno_dest, i);
1616             }
1617         }
1618     }
1619 
1620     switch (method) {
1621     case INTERP_SPLINE:
1622         s = "cubic spline";
1623         break;
1624     case INTERP_ASPLINE:
1625         s = "Akima spline";
1626         break;
1627     default:
1628         s = "linear interpolation";
1629         break;
1630     }
1631     sprintf(buf, "Interpolated from G%d.S%d using %s", gno_src, setno_src, s);
1632     setcomment(gno_dest, setno_dest, buf);
1633 
1634     return RETURN_SUCCESS;
1635 }
1636 
get_restriction_array(int gno,int setno,int rtype,int negate,char ** rarray)1637 int get_restriction_array(int gno, int setno,
1638     int rtype, int negate, char **rarray)
1639 {
1640     int i, n, regno;
1641     double *x, *y;
1642     world w;
1643     WPoint wp;
1644 
1645     if (rtype == RESTRICT_NONE) {
1646         *rarray = NULL;
1647         return RETURN_SUCCESS;
1648     }
1649 
1650     n = getsetlength(gno, setno);
1651     if (n <= 0) {
1652         *rarray = NULL;
1653         return RETURN_FAILURE;
1654     }
1655 
1656     *rarray = xmalloc(n*SIZEOF_CHAR);
1657     if (*rarray == NULL) {
1658         return RETURN_FAILURE;
1659     }
1660 
1661     x = getcol(gno, setno, DATA_X);
1662     y = getcol(gno, setno, DATA_Y);
1663 
1664     switch (rtype) {
1665     case RESTRICT_REG0:
1666     case RESTRICT_REG1:
1667     case RESTRICT_REG2:
1668     case RESTRICT_REG3:
1669     case RESTRICT_REG4:
1670         regno = rtype - RESTRICT_REG0;
1671         for (i = 0; i < n; i++) {
1672             (*rarray)[i] = inregion(regno, x[i], y[i]) ? !negate : negate;
1673         }
1674         break;
1675     case RESTRICT_WORLD:
1676         get_graph_world(gno, &w);
1677         for (i = 0; i < n; i++) {
1678             wp.x = x[i];
1679             wp.y = y[i];
1680             (*rarray)[i] = is_wpoint_inside(&wp, &w) ? !negate : negate;
1681         }
1682         break;
1683     default:
1684         errmsg("Internal error in get_restriction_array()");
1685         XCFREE(*rarray);
1686         return RETURN_FAILURE;
1687         break;
1688     }
1689     return RETURN_SUCCESS;
1690 }
1691 
monotonicity(double * array,int len,int strict)1692 int monotonicity(double *array, int len, int strict)
1693 {
1694     int i;
1695     int s0, s1;
1696 
1697     if (len < 2) {
1698         errmsg("Monotonicity of an array of length < 2 is meaningless");
1699         return 0;
1700     }
1701 
1702     s0 = sign(array[1] - array[0]);
1703     for (i = 2; i < len; i++) {
1704         s1 = sign(array[i] - array[i - 1]);
1705         if (s1 != s0) {
1706             if (strict) {
1707                 return 0;
1708             } else if (s0 == 0) {
1709                 s0 = s1;
1710             } else if (s1 != 0) {
1711                 return 0;
1712             }
1713         }
1714     }
1715 
1716     return s0;
1717 }
1718 
find_span_index(double * array,int len,int m,double x)1719 int find_span_index(double *array, int len, int m, double x)
1720 {
1721     int ind, low = 0, high = len - 1;
1722 
1723     if (len < 2 || m == 0) {
1724         errmsg("find_span_index() called with a non-monotonic array");
1725         return -2;
1726     } else if (m > 0) {
1727         /* ascending order */
1728         if (x < array[0]) {
1729             return -1;
1730         } else if (x > array[len - 1]) {
1731             return len - 1;
1732         } else {
1733             while (low <= high) {
1734 	        ind = (low + high) / 2;
1735 	        if (x < array[ind]) {
1736 	            high = ind - 1;
1737 	        } else if (x > array[ind + 1]) {
1738 	            low = ind + 1;
1739 	        } else {
1740 	            return ind;
1741 	        }
1742             }
1743         }
1744     } else {
1745         /* descending order */
1746         if (x > array[0]) {
1747             return -1;
1748         } else if (x < array[len - 1]) {
1749             return len - 1;
1750         } else {
1751             while (low <= high) {
1752 	        ind = (low + high) / 2;
1753 	        if (x > array[ind]) {
1754 	            high = ind - 1;
1755 	        } else if (x < array[ind + 1]) {
1756 	            low = ind + 1;
1757 	        } else {
1758 	            return ind;
1759 	        }
1760             }
1761         }
1762     }
1763 
1764     /* should never happen */
1765     errmsg("internal error in find_span_index()");
1766     return -2;
1767 }
1768