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