/*
 * Grace - GRaphing, Advanced Computation and Exploration of data
 * 
 * Home page: http://plasma-gate.weizmann.ac.il/Grace/
 * 
 * Copyright (c) 1991-1995 Paul J Turner, Portland, OR
 * Copyright (c) 1996-2000 Grace Development Team
 * 
 * Maintained by Evgeny Stambulchik
 * 
 * 
 *                           All Rights Reserved
 * 
 *    This program is free software; you can redistribute it and/or modify
 *    it under the terms of the GNU General Public License as published by
 *    the Free Software Foundation; either version 2 of the License, or
 *    (at your option) any later version.
 * 
 *    This program is distributed in the hope that it will be useful,
 *    but WITHOUT ANY WARRANTY; without even the implied warranty of
 *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *    GNU General Public License for more details.
 * 
 *    You should have received a copy of the GNU General Public License
 *    along with this program; if not, write to the Free Software
 *    Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/*
 *
 * procedures for performing transformations from the command
 * line interpreter and the GUI.
 *
 */

#include <config.h>
#include <cmath.h>

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include "globals.h"
#include "utils.h"
#include "draw.h"
#include "ssdata.h"
#include "graphs.h"
#include "parser.h"
#include "protos.h"

static void forwarddiff(double *x, double *y, double *resx, double *resy, int n);
static void backwarddiff(double *x, double *y, double *resx, double *resy, int n);
static void centereddiff(double *x, double *y, double *resx, double *resy, int n);
int get_points_inregion(int rno, int invr, int len, double *x, double *y, int *cnt, double **xt, double **yt);

static char buf[256];


void do_fourier_command(int gno, int setno, int ftype, int ltype)
{
    switch (ftype) {
    case FFT_DFT:
	do_fourier(gno, setno, 0, 0, ltype, 0, 0, 0);
	break;
    case FFT_INVDFT       :
	do_fourier(gno, setno, 0, 0, ltype, 1, 0, 0);
	break;
    case FFT_FFT:
	do_fourier(gno, setno, 1, 0, ltype, 0, 0, 0);
	break;
    case FFT_INVFFT       :
	do_fourier(gno, setno, 1, 0, ltype, 1, 0, 0);
	break;
    }
}


/*
 * evaluate a formula
 */
int do_compute(int gno, int setno, int graphto, int loadto, char *rarray, char *fstr)
{
    if (is_set_active(gno, setno)) {
	if (gno != graphto || setno != loadto) {
	    if (copysetdata(gno, setno, graphto, loadto) != RETURN_SUCCESS) {
	        return RETURN_FAILURE;
            }
        }
	filter_set(graphto, loadto, rarray);
        set_parser_setno(graphto, loadto);
        if (scanner(fstr) != RETURN_SUCCESS) {
	    if (graphto != gno || loadto != setno) {
		killset(graphto, loadto);
	    }
	    return RETURN_FAILURE;
	} else {
	    set_dirtystate();
            return RETURN_SUCCESS;
        }
    } else {
        return RETURN_FAILURE;
    }
}

/*
 * forward, backward and centered differences
 */
static void forwarddiff(double *x, double *y, double *resx, double *resy, int n)
{
    int i, eflag = 0;
    double h;

    for (i = 1; i < n; i++) {
	resx[i - 1] = x[i - 1];
	h = x[i - 1] - x[i];
	if (h == 0.0) {
	    resy[i - 1] = - MAXNUM;
	    eflag = 1;
	} else {
	    resy[i - 1] = (y[i - 1] - y[i]) / h;
	}
    }
    if (eflag) {
	errmsg("Warning: infinite slope, check set status before proceeding");
    }
}

static void backwarddiff(double *x, double *y, double *resx, double *resy, int n)
{
    int i, eflag = 0;
    double h;

    for (i = 0; i < n - 1; i++) {
	resx[i] = x[i];
	h = x[i + 1] - x[i];
	if (h == 0.0) {
	    resy[i] = - MAXNUM;
	    eflag = 1;
	} else {
	    resy[i] = (y[i + 1] - y[i]) / h;
	}
    }
    if (eflag) {
	errmsg("Warning: infinite slope, check set status before proceeding");
    }
}

static void centereddiff(double *x, double *y, double *resx, double *resy, int n)
{
    int i, eflag = 0;
    double h1, h2;

    for (i = 1; i < n - 1; i++) {
	resx[i - 1] = x[i];
	h1 = x[i] - x[i - 1];
	h2 = x[i + 1] - x[i];
	if (h1 + h2 == 0.0) {
	    resy[i - 1] = - MAXNUM;
	    eflag = 1;
	} else {
	    resy[i - 1] = (y[i + 1] - y[i - 1]) / (h1 + h2);
	}
    }
    if (eflag) {
	errmsg("Warning: infinite slope, check set status before proceeding");
    }
}

static void seasonaldiff(double *x, double *y,
			 double *resx, double *resy, int n, int period)
{
    int i;

    for (i = 0; i < n - period; i++) {
	resx[i] = x[i];
	resy[i] = y[i] - y[i + period];
    }
}

/*
 * trapezoidal rule
 */
double trapint(double *x, double *y, double *resx, double *resy, int n)
{
    int i;
    double sum = 0.0;
    double h;

    if (n < 2) {
        return 0.0;
    }
    
    if (resx != NULL) {
        resx[0] = x[0];
    }
    if (resy != NULL) {
        resy[0] = 0.0;
    }
    for (i = 1; i < n; i++) {
	h = (x[i] - x[i - 1]);
	if (resx != NULL) {
	    resx[i] = x[i];
	}
	sum = sum + h * (y[i - 1] + y[i]) * 0.5;
	if (resy != NULL) {
	    resy[i] = sum;
	}
    }
    return sum;
}

/*
 * apply a digital filter
 */
void do_digfilter(int set1, int set2)
{
    int digfiltset;

    if (!(is_set_active(get_cg(), set1) && is_set_active(get_cg(), set2))) {
	errmsg("Set not active");
	return;
    }
    if ((getsetlength(get_cg(), set1) < 3) || (getsetlength(get_cg(), set2) < 3)) {
	errmsg("Set length < 3");
	return;
    }
    digfiltset = nextset(get_cg());
    if (digfiltset != (-1)) {
	activateset(get_cg(), digfiltset);
	setlength(get_cg(), digfiltset, getsetlength(get_cg(), set1) - getsetlength(get_cg(), set2) + 1);
	sprintf(buf, "Digital filter from set %d applied to set %d", set2, set1);
	filterser(getsetlength(get_cg(), set1),
		  getx(get_cg(), set1),
		  gety(get_cg(), set1),
		  getx(get_cg(), digfiltset),
		  gety(get_cg(), digfiltset),
		  gety(get_cg(), set2),
		  getsetlength(get_cg(), set2));
	setcomment(get_cg(), digfiltset, buf);
    }
}

/*
 * linear convolution
 */
void do_linearc(int set1, int set2)
{
    int linearcset, i, itmp;
    double *xtmp;

    if (!(is_set_active(get_cg(), set1) && is_set_active(get_cg(), set2))) {
	errmsg("Set not active");
	return;
    }
    if ((getsetlength(get_cg(), set1) < 3) || (getsetlength(get_cg(), set2) < 3)) {
	errmsg("Set length < 3");
	return;
    }
    linearcset = nextset(get_cg());
    if (linearcset != (-1)) {
	activateset(get_cg(), linearcset);
	setlength(get_cg(), linearcset, (itmp = getsetlength(get_cg(), set1) + getsetlength(get_cg(), set2) - 1));
	linearconv(gety(get_cg(), set2), gety(get_cg(), set1), gety(get_cg(), linearcset), getsetlength(get_cg(), set2), getsetlength(get_cg(), set1));
	xtmp = getx(get_cg(), linearcset);
	for (i = 0; i < itmp; i++) {
	    xtmp[i] = i;
	}
	sprintf(buf, "Linear convolution of set %d with set %d", set1, set2);
	setcomment(get_cg(), linearcset, buf);
    }
}

/*
 * cross correlation/covariance
 */
void do_xcor(int gno1, int set1, int gno2, int set2, int maxlag, int covar)
{
    int xcorset, i, ierr, len, cg = get_cg();
    double *xtmp;

    if (!(is_set_active(gno1, set1) && is_set_active(gno2, set2))) {
	errmsg("Set not active");
	return;
    }
    
    len = getsetlength(gno1, set1);
    if (getsetlength(gno2, set2) != len) {
	errmsg("Sets must be of the same length");
    }
    if (len < 2) {
	errmsg("Set length < 2");
	return;
    }
    
    if (maxlag < 1 || maxlag > len) {
	errmsg("Lag incorrectly specified");
	return;
    }
    
    xcorset = nextset(cg);
    
    if (xcorset != (-1)) {
	char *fname;
        activateset(cg, xcorset);
	setlength(cg, xcorset, maxlag);
	if (covar) {
            fname = "covariance";
        } else {
            fname = "correlation";
        }
        if (set1 != set2) {
	    sprintf(buf, "X-%s of G%d.S%d and G%d.S%d at maximum lag %d",
                    fname, gno1, set1, gno2, set2, maxlag);
	} else {
	    sprintf(buf, "Auto-%s of G%d.S%d at maximum lag %d",
                    fname, gno1, set1, maxlag);
	}
	ierr = crosscorr(gety(gno1, set1), gety(gno2, set2), len,
                         maxlag, covar, gety(cg, xcorset));
	xtmp = getx(cg, xcorset);
	for (i = 0; i < maxlag; i++) {
	    xtmp[i] = i;
	}
	setcomment(cg, xcorset, buf);
    }
}


/*
 * numerical integration
 */
double do_int(int gno, int setno, int itype)
{
    int intset;
    double sum = 0;

    if (!is_set_active(gno, setno)) {
	errmsg("Set not active");
	return 0.0;
    }
    if (getsetlength(gno, setno) < 3) {
	errmsg("Set length < 3");
	return 0.0;
    }
    if (itype == 0) {
	intset = nextset(gno);
	if (intset != (-1)) {
	    activateset(gno, intset);
	    setlength(gno, intset, getsetlength(gno, setno));
	    sprintf(buf, "Cumulative sum of set %d", setno);
	    sum = trapint(getx(gno, setno), gety(gno, setno), getx(gno, intset), gety(gno, intset), getsetlength(gno, setno));
	    setcomment(gno, intset, buf);
	}
    } else {
	sum = trapint(getx(gno, setno), gety(gno, setno), NULL, NULL, getsetlength(gno, setno));
    }
    return sum;
}

/*
 * difference a set
 * itype means
 *  0 - forward
 *  1 - backward
 *  2 - centered difference
 */
void do_differ(int gno, int setno, int itype)
{
    int diffset;

    if (!is_set_active(gno, setno)) {
	errmsg("Set not active");
	return;
    }
    if (getsetlength(gno, setno) < 3) {
	errmsg("Set length < 3");
	return;
    }
    diffset = nextset(gno);
    if (diffset != (-1)) {
	activateset(gno, diffset);
	switch (itype) {
	case 0:
	    sprintf(buf, "Forward difference of set %d", setno);
	    setlength(gno, diffset, getsetlength(gno, setno) - 1);
	    forwarddiff(getx(gno, setno), gety(gno, setno), getx(gno, diffset), gety(gno, diffset), getsetlength(gno, setno));
	    break;
	case 1:
	    sprintf(buf, "Backward difference of set %d", setno);
	    setlength(gno, diffset, getsetlength(gno, setno) - 1);
	    backwarddiff(getx(gno, setno), gety(gno, setno), getx(gno, diffset), gety(gno, diffset), getsetlength(gno, setno));
	    break;
	case 2:
	    sprintf(buf, "Centered difference of set %d", setno);
	    setlength(gno, diffset, getsetlength(gno, setno) - 2);
	    centereddiff(getx(gno, setno), gety(gno, setno), getx(gno, diffset), gety(gno, diffset), getsetlength(gno, setno));
	    break;
	}
	setcomment(gno, diffset, buf);
    }
}

/*
 * seasonally difference a set
 */
void do_seasonal_diff(int setno, int period)
{
    int diffset;

    if (!is_set_active(get_cg(), setno)) {
	errmsg("Set not active");
	return;
    }
    if (getsetlength(get_cg(), setno) < 2) {
	errmsg("Set length < 2");
	return;
    }
    diffset = nextset(get_cg());
    if (diffset != (-1)) {
	activateset(get_cg(), diffset);
	setlength(get_cg(), diffset, getsetlength(get_cg(), setno) - period);
	seasonaldiff(getx(get_cg(), setno), gety(get_cg(), setno),
		     getx(get_cg(), diffset), gety(get_cg(), diffset),
		     getsetlength(get_cg(), setno), period);
	sprintf(buf, "Seasonal difference of set %d, period %d", setno, period);
	setcomment(get_cg(), diffset, buf);
    }
}

/*
 * regression with restrictions to region rno if rno >= 0
 */
void do_regress(int gno, int setno, int ideg, int iresid, int rno, int invr, int fitset)
	/* 
	 * gno, setno  - set to perform fit on
	 * ideg   - degree of fit
	 * irisid - 0 -> whole set, 1-> subset of setno
	 * rno    - region number of subset
	 * invr   - 1->invert region, 0 -> do nothing
	 * fitset - set to which fitted function will be loaded
	 * 			Y values are computed at the x values in the set
	 *          if -1 is specified, a set with the same x values as the
	 *          one being fitted will be created and used
	 */
{
    int len, i, sdeg = ideg;
    int cnt = 0, fitlen = 0;
    double *x, *y, *xt = NULL, *yt = NULL, *xr = NULL, *yr = NULL;
    char buf[256];
    double c[20];   /* coefficients of fitted polynomial */

    if (!is_set_active(gno, setno)) {
		errmsg("Set not active");
		return;
    }
    len = getsetlength(gno, setno);
    x = getx(gno, setno);
    y = gety(gno, setno);
    if (rno == -1) {
		xt = x;
		yt = y;
    } else if (isactive_region(rno)) {
		if (!get_points_inregion(rno, invr, len, x, y, &cnt, &xt, &yt)) {
			if (cnt == 0) {
			errmsg("No points found in region, operation cancelled");
			}
			return;
		}
		len = cnt;
    } else {
		errmsg("Selected region is not active");
		return;
    }
    /*
     * first part for polynomials, second part for linear fits to transformed
     * data
     */
    if ((len < ideg && ideg <= 10) || (len < 2 && ideg > 10)) {
		errmsg("Too few points in set, operation cancelled");
		return;
    }
	/* determine is set provided or use abscissa from fitted set */
    if( fitset == -1 ) {		        /* create set */
    	if( (fitset = nextset(gno)) != -1 ) {
			activateset(gno, fitset);
			setlength(gno, fitset, len);
			fitlen = len;
			xr = getx(gno, fitset);
			for (i = 0; i < len; i++) {
	    		xr[i] = xt[i];
			}	
			yr = gety(gno, fitset);
		}
    } else {									/* set has been provided */
		fitlen = getsetlength( gno, fitset );
		xr = getx(gno, fitset);
		yr = gety(gno, fitset);
    }

	/* transform data so that system has the form y' = A + B * x' */
    if (fitset != -1) {
	if (ideg == 12) {	/* y=A*x^B -> ln(y) = ln(A) + B * ln(x) */
	    ideg = 1;
	    for (i = 0; i < len; i++) {
			if (xt[i] <= 0.0) {
				errmsg("One of X[i] <= 0.0");
				return;
			}
			if (yt[i] <= 0.0) {
				errmsg("One of Y[i] <= 0.0");
				return;
			}
	    }
	    for (i = 0; i < len; i++) {
			xt[i] = log(xt[i]);
			yt[i] = log(yt[i]);
	    }
		for( i=0; i<fitlen; i++ )
			if( xr[i] <= 0.0 ) {
				errmsg("One of X[i] <= 0.0");
				return;
			} 
		for( i=0; i<fitlen; i++ )
			xr[i] = log( xr[i] );
	} else if (ideg == 13) {   /*y=Aexp(B*x) -> ln(y) = ln(A) + B * x */
	    ideg = 1;
	    for (i = 0; i < len; i++) {
			if (yt[i] <= 0.0) {
				errmsg("One of Y[i] <= 0.0");
				return;
		 	}
	    }
	    for (i = 0; i < len; i++) {
			yt[i] = log(yt[i]);
	    }
	} else if (ideg == 14) {	/* y = A + B * ln(x) */
	    ideg = 1;
	    for (i = 0; i < len; i++) {
			if (xt[i] <= 0.0) {
				errmsg("One of X[i] <= 0.0");
				return;
			}
	    }
		for (i = 0; i < len; i++) {
			xt[i] = log(xt[i]);	
		}
		for( i=0; i<fitlen; i++ )
			if( xr[i] <= 0.0 ) {
				errmsg("One of X[i] <= 0.0");
				return;
			} 
		for( i=0; i<fitlen; i++ ){
			xr[i] = log( xr[i] );
		}	  
	} else if (ideg == 15) {	/* y = 1/( A + B*x ) -> 1/y = a + B*x */
	    ideg = 1;
	    for (i = 0; i < len; i++) {
			if (yt[i] == 0.0) {
			    errmsg("One of Y[i] = 0.0");
			    return;
			}
	    }
	    for (i = 0; i < len; i++) {
			yt[i] = 1.0 / yt[i];
	    }
	}

	if (fitcurve(xt, yt, len, ideg, c)) {
	    killset(gno, fitset);
	    goto bustout;
	}

	/* compute function at requested x ordinates */
	for( i=0; i<fitlen; i++ )
		yr[i] = leasev( c, ideg, xr[i] );

	sprintf(buf, "\nN.B. Statistics refer to the transformed model\n");
	/* apply inverse transform, output fitted function in human readable form */
	if( sdeg<11 ) {
		sprintf(buf, "\ny = %.5g %c %.5g * x", c[0], c[1]>=0?'+':'-', fabs(c[1]));
		for( i=2; i<=ideg; i++ )
			sprintf( buf+strlen(buf), " %c %.5g * x^%d", c[i]>=0?'+':'-', 
															fabs(c[i]), i );
		strcat( buf, "\n" );
	} else if (sdeg == 12) {	/* ln(y) = ln(A) + b * ln(x) */
		sprintf( buf, "\ny = %.5g * x^%.5g\n", exp(c[0]), c[1] );
	    for (i = 0; i < len; i++) {
			xt[i] = exp(xt[i]);
			yt[i] = exp(yt[i]);
	    }
	    for (i = 0; i < fitlen; i++){
			yr[i] = exp(yr[i]);
			xr[i] = exp(xr[i]);
		}
	} else if (sdeg == 13) { 	/* ln(y) = ln(A) + B * x */
		sprintf( buf, "\ny = %.5g * exp( %.5g * x )\n", exp(c[0]), c[1] );
	    for (i = 0; i < len; i++) {
			yt[i] = exp(yt[i]);
	    }
	    for (i = 0; i < fitlen; i++)
			yr[i] = exp(yr[i]);
	} else if (sdeg == 14) {	/* y = A + B * ln(x) */
		sprintf(buf, "\ny = %.5g %c %.5g * ln(x)\n", c[0], c[1]>=0?'+':'-',
											fabs(c[1]) );
	    for (i = 0; i < len; i++) {
			xt[i] = exp(xt[i]);
	    }
	    for (i = 0; i < fitlen; i++)
			xr[i] = exp(xr[i]);
	} else if (sdeg == 15) {	/* y = 1/( A + B*x ) */
		sprintf( buf, "\ny = 1/(%.5g %c %.5g * x)\n", c[0], c[1]>=0?'+':'-',
													fabs(c[1]) );
	    for (i = 0; i < len; i++) {
			yt[i] = 1.0 / yt[i];
	    }
	    for (i = 0; i < fitlen; i++)
			yr[i] = 1.0 / yr[i];
	}
	stufftext(buf);
	sprintf(buf, "\nRegression of set %d results to set %d\n", setno, fitset);
	stufftext(buf);
	
	switch (iresid) {
	case 1:
		/* determine residual */
	    for (i = 0; i < len; i++) {
			yr[i] = yt[i] - yr[i];
	    }
	    break;
	case 2:
	    break;
	}
	sprintf(buf, "%d deg fit of set %d", ideg, setno);
	setcomment(gno, fitset, buf);
    }
    bustout:;
    if (rno >= 0 && cnt != 0) {	/* had a region and allocated memory there */
	xfree(xt);
	xfree(yt);
    }
}

/*
 * running averages, medians, min, max, std. deviation
 */
void do_runavg(int gno, int setno, int runlen, int runtype, int rno, int invr)
{
    int runset;
    int len, cnt = 0;
    double *x, *y, *xt = NULL, *yt = NULL, *xr, *yr;

    if (!is_set_active(gno, setno)) {
	errmsg("Set not active");
	return;
    }
    if (runlen < 2) {
	errmsg("Length of running average < 2");
	return;
    }
    len = getsetlength(gno, setno);
    x = getx(gno, setno);
    y = gety(gno, setno);
    if (rno == -1) {
	xt = x;
	yt = y;
    } else if (isactive_region(rno)) {
	if (!get_points_inregion(rno, invr, len, x, y, &cnt, &xt, &yt)) {
	    if (cnt == 0) {
		errmsg("No points found in region, operation cancelled");
	    }
	    return;
	}
	len = cnt;
    } else {
	errmsg("Selected region is not active");
	return;
    }
    if (runlen >= len) {
	errmsg("Length of running average > set length");
	goto bustout;
    }
    runset = nextset(gno);
    if (runset != (-1)) {
	activateset(gno, runset);
	setlength(gno, runset, len - runlen + 1);
	xr = getx(gno, runset);
	yr = gety(gno, runset);
	switch (runtype) {
	case 0:
	    runavg(xt, yt, xr, yr, len, runlen);
	    sprintf(buf, "%d-pt. avg. on set %d ", runlen, setno);
	    break;
	case 1:
	    runmedian(xt, yt, xr, yr, len, runlen);
	    sprintf(buf, "%d-pt. median on set %d ", runlen, setno);
	    break;
	case 2:
	    runminmax(xt, yt, xr, yr, len, runlen, 0);
	    sprintf(buf, "%d-pt. min on set %d ", runlen, setno);
	    break;
	case 3:
	    runminmax(xt, yt, xr, yr, len, runlen, 1);
	    sprintf(buf, "%d-pt. max on set %d ", runlen, setno);
	    break;
	case 4:
	    runstddev(xt, yt, xr, yr, len, runlen);
	    sprintf(buf, "%d-pt. std dev., set %d ", runlen, setno);
	    break;
	}
	setcomment(gno, runset, buf);
    }
  bustout:;
    if (rno >= 0 && cnt != 0) {	/* had a region and allocated memory there */
	xfree(xt);
	xfree(yt);
    }
}

/*
 * DFT by FFT or definition
 */
void do_fourier(int gno, int setno, int fftflag, int load, int loadx, int invflag, int type, int wind)
{
    int i, ilen;
    double *x, *y, *xx, *yy, delt, T;
    int i2 = 0, specset;

    if (!is_set_active(get_cg(), setno)) {
	errmsg("Set not active");
	return;
    }
    ilen = getsetlength(get_cg(), setno);
    if (ilen < 2) {
	errmsg("Set length < 2");
	return;
    }
    if (fftflag) {
	if ((i2 = ilog2(ilen)) <= 0) {
	    errmsg("Set length not a power of 2");
	    return;
	}
    }
    specset = nextset(get_cg());
    if (specset != -1) {
	activateset(get_cg(), specset);
	setlength(get_cg(), specset, ilen);
	xx = getx(get_cg(), specset);
	yy = gety(get_cg(), specset);
	x = getx(get_cg(), setno);
	y = gety(get_cg(), setno);
	copyx(get_cg(), setno, specset);
	copyy(get_cg(), setno, specset);
	if (wind != 0) {	/* apply data window if needed */
	    apply_window(xx, yy, ilen, type, wind);
	}
	if (type == 0) {	/* real data */
	    for (i = 0; i < ilen; i++) {
		xx[i] = yy[i];
		yy[i] = 0.0;
	    }
	}
	if (fftflag) {
	    fft(xx, yy, ilen, i2, invflag);
	} else {
	    dft(xx, yy, ilen, invflag);
	}
	switch (load) {
	case 0:
	    delt = (x[ilen-1] - x[0])/(ilen -1.0);
	    T = (x[ilen - 1] - x[0]);
	    xx = getx(get_cg(), specset);
	    yy = gety(get_cg(), specset);
	    for (i = 0; i < ilen / 2; i++) {
	      /* carefully get amplitude of complex xform: 
		 use abs(a[i]) + abs(a[-i]) except for zero term */
	      if(i) yy[i] = hypot(xx[i], yy[i])+hypot(xx[ilen-i], yy[ilen-i]);
	      else yy[i]=fabs(xx[i]);
		switch (loadx) {
		case 0:
		    xx[i] = i;
		    break;
		case 1:
		    /* xx[i] = 2.0 * M_PI * i / ilen; */
		    xx[i] = i / T;
		    break;
		case 2:
		    if (i == 0) {
			xx[i] = T + delt;	/* the mean */
		    } else {
			/* xx[i] = (double) ilen / (double) i; */
			xx[i] = T / i;
		    }
		    break;
		}
	    }
	    setlength(get_cg(), specset, ilen / 2);
	    break;
	case 1:
	    delt = (x[ilen-1] - x[0])/(ilen -1.0);
	    T = (x[ilen - 1] - x[0]);
	    setlength(get_cg(), specset, ilen / 2);
	    xx = getx(get_cg(), specset);
	    yy = gety(get_cg(), specset);
	    for (i = 0; i < ilen / 2; i++) {
		yy[i] = -atan2(yy[i], xx[i]);
		switch (loadx) {
		case 0:
		    xx[i] = i;
		    break;
		case 1:
		    /* xx[i] = 2.0 * M_PI * i / ilen; */
		    xx[i] = i / T;
		    break;
		case 2:
		    if (i == 0) {
			xx[i] = T + delt;
		    } else {
			/* xx[i] = (double) ilen / (double) i; */
			xx[i] = T / i;
		    }
		    break;
		}
	    }
	    break;
	}
	if (fftflag) {
	    sprintf(buf, "FFT of set %d", setno);
	} else {
	    sprintf(buf, "DFT of set %d", setno);
	}
	setcomment(get_cg(), specset, buf);
    }
}

/*
 * Apply a window to a set, result goes to a new set.
 */
void do_window(int setno, int type, int wind)
{
    int ilen;
    double *xx, *yy;
    int specset;

    if (!is_set_active(get_cg(), setno)) {
	errmsg("Set not active");
	return;
    }
    ilen = getsetlength(get_cg(), setno);
    if (ilen < 2) {
	errmsg("Set length < 2");
	return;
    }
    specset = nextset(get_cg());
    if (specset != -1) {
	char *wtype[6];
	wtype[0] = "Triangular";
	wtype[1] = "Hanning";
	wtype[2] = "Welch";
	wtype[3] = "Hamming";
	wtype[4] = "Blackman";
	wtype[5] = "Parzen";

	activateset(get_cg(), specset);
	setlength(get_cg(), specset, ilen);
	xx = getx(get_cg(), specset);
	yy = gety(get_cg(), specset);
	copyx(get_cg(), setno, specset);
	copyy(get_cg(), setno, specset);
	if (wind != 0) {
	    apply_window(xx, yy, ilen, type, wind);
	    sprintf(buf, "%s windowed set %d", wtype[wind - 1], setno);
	} else {		/* shouldn't happen */
	}
	setcomment(get_cg(), specset, buf);
    }
}

void apply_window(double *xx, double *yy, int ilen, int type, int wind)
{
    int i;

    for (i = 0; i < ilen; i++) {
	switch (wind) {
	case 1:		/* triangular */
	    if (type != 0) {
		xx[i] *= 1.0 - fabs((i - 0.5 * (ilen - 1.0)) / (0.5 * (ilen - 1.0)));
	    }
	    yy[i] *= 1.0 - fabs((i - 0.5 * (ilen - 1.0)) / (0.5 * (ilen - 1.0)));

	    break;
	case 2:		/* Hanning */
	    if (type != 0) {
		xx[i] = xx[i] * (0.5 - 0.5 * cos(2.0 * M_PI * i / (ilen - 1.0)));
	    }
	    yy[i] = yy[i] * (0.5 - 0.5 * cos(2.0 * M_PI * i / (ilen - 1.0)));
	    break;
	case 3:		/* Welch (from Numerical Recipes) */
	    if (type != 0) {
		xx[i] *= 1.0 - pow((i - 0.5 * (ilen - 1.0)) / (0.5 * (ilen + 1.0)), 2.0);
	    }
	    yy[i] *= 1.0 - pow((i - 0.5 * (ilen - 1.0)) / (0.5 * (ilen + 1.0)), 2.0);
	    break;
	case 4:		/* Hamming */
	    if (type != 0) {
		xx[i] = xx[i] * (0.54 - 0.46 * cos(2.0 * M_PI * i / (ilen - 1.0)));
	    }
	    yy[i] = yy[i] * (0.54 - 0.46 * cos(2.0 * M_PI * i / (ilen - 1.0)));
	    break;
	case 5:		/* Blackman */
	    if (type != 0) {
		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)));
	    }
	    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)));
	    break;
	case 6:		/* Parzen (from Numerical Recipes) */
	    if (type != 0) {
		xx[i] *= 1.0 - fabs((i - 0.5 * (ilen - 1)) / (0.5 * (ilen + 1)));
	    }
	    yy[i] *= 1.0 - fabs((i - 0.5 * (ilen - 1)) / (0.5 * (ilen + 1)));
	    break;
	}
    }
}


/*
 * histograms
 */
int do_histo(int fromgraph, int fromset, int tograph, int toset,
	      double *bins, int nbins, int cumulative, int normalize)
{
    int i, ndata;
    int *hist;
    double *x, *y, *data;
    plotarr p;
    
    if (!is_set_active(fromgraph, fromset)) {
	errmsg("Set not active");
	return RETURN_FAILURE;
    }
    if (nbins <= 0) {
	errmsg("Number of bins <= 0");
	return RETURN_FAILURE;
    }
    if (toset == SET_SELECT_NEXT) {
	toset = nextset(tograph);
    }
    if (!is_valid_setno(tograph, toset)) {
	errmsg("Can't activate destination set");
        return RETURN_FAILURE;
    }
    
    ndata = getsetlength(fromgraph, fromset);
    data = gety(fromgraph, fromset);
    
    hist = xmalloc(nbins*SIZEOF_INT);
    if (hist == NULL) {
        errmsg("xmalloc failed in do_histo()");
        return RETURN_FAILURE;
    }

    if (histogram(ndata, data, nbins, bins, hist) == RETURN_FAILURE){
        xfree(hist);
        return RETURN_FAILURE;
    }
    
    activateset(tograph, toset);
    setlength(tograph, toset, nbins + 1);
    x = getx(tograph, toset);
    y = gety(tograph, toset);
    
    x[0] = bins[0];
    y[0] = 0.0;
    for (i = 1; i < nbins + 1; i++) {
        x[i] = bins[i];
        y[i] = hist[i - 1];
        if (cumulative) {
            y[i] += y[i - 1];
        }
    }
    
    if (normalize) {
        for (i = 1; i < nbins + 1; i++) {
            double factor;
            if (cumulative) {
                factor = 1.0/ndata;
            } else {
                factor = 1.0/((bins[i] - bins[i - 1])*ndata);
            }
            y[i] *= factor;
        }
    }
    
    xfree(hist);

    get_graph_plotarr(tograph, toset, &p);
    p.sym = SYM_NONE;
    p.linet = LINE_TYPE_LEFTSTAIR;
    p.dropline = TRUE;
    p.baseline = FALSE;
    p.baseline_type = BASELINE_TYPE_0;
    p.lines = 1;
    p.symlines = 1;
    sprintf(p.comments, "Histogram from G%d.S%d", fromgraph, fromset);
    set_graph_plotarr(tograph, toset, &p);

    return RETURN_SUCCESS;
}

int histogram(int ndata, double *data, int nbins, double *bins, int *hist)
{
    int i, j, bsign;
    double minval, maxval;
    
    if (nbins < 1) {
        errmsg("Number of bins < 1");
        return RETURN_FAILURE;
    }
    
    bsign = monotonicity(bins, nbins + 1, TRUE);
    if (bsign == 0) {
        errmsg("Non-monotonic bins");
        return RETURN_FAILURE;
    }
    
    for (i = 0; i < nbins; i++) {
        hist[i] = 0;
    }
    
    /* TODO: binary search */
    for (i = 0; i < ndata; i++) {
        for (j = 0; j < nbins; j++) {
            if (bsign > 0) {
                minval = bins[j];
                maxval = bins[j + 1];
            } else {
                minval = bins[j + 1];
                maxval = bins[j];
            }
            if (data[i] >= minval && data[i] <= maxval) {
                hist[j] += 1;
                break;
            }
        }
    }
    return RETURN_SUCCESS;
}


/*
 * sample a set, by start/step or logical expression
 */
void do_sample(int setno, int typeno, char *exprstr, int startno, int stepno)
{
    int len, npts = 0, i, resset;
    double *x, *y;
    int reslen;
    double *result;
    int gno = get_cg();

    if (!is_set_active(gno, setno)) {
	errmsg("Set not active");
	return;
    }

    len = getsetlength(gno, setno);

    resset = nextset(gno);
    if (resset < 0) {
	return;
    }

    x = getx(gno, setno);
    y = gety(gno, setno);

    if (typeno == 0) {
	if (len <= 2) {
	    errmsg("Set has <= 2 points");
	    return;
	}
	if (startno < 1) {
	    errmsg("Start point < 1 (locations in sets are numbered starting from 1)");
	    return;
	}
	if (stepno < 1) {
	    errmsg("Step < 1");
	    return;
	}
	for (i = startno - 1; i < len; i += stepno) {
	    add_point(gno, resset, x[i], y[i]);
	    npts++;
	}
	sprintf(buf, "Sample, %d, %d set #%d", startno, stepno, setno);
    } else {
        if (set_parser_setno(gno, setno) != RETURN_SUCCESS) {
	    errmsg("Bad set");
            killset(gno, resset);
	    return;
        }
        if (v_scanner(exprstr, &reslen, &result) != RETURN_SUCCESS) {
            killset(gno, resset);
	    return;
        }
        if (reslen != len) {
	    errmsg("Internal error");
            killset(gno, resset);
	    return;
        }

        npts = 0;
	sprintf(buf, "Sample from %d, using '%s'", setno, exprstr);
	for (i = 0; i < len; i++) {
	    if ((int) rint(result[i])) {
		add_point(gno, resset, x[i], y[i]);
		npts++;
	    }
	}
        xfree(result);
    }
    if (npts > 0) {
	setcomment(gno, resset, buf);
    }
}

#define prune_xconv(res,x,xtype)	\
    switch (deltatypeno) {		\
    case PRUNE_VIEWPORT:		\
	res = xy_xconv(x);			\
	break;				\
    case PRUNE_WORLD:			\
	switch (xtype) {		\
	case PRUNE_LIN:			\
	    res = x;			\
	    break;			\
	case PRUNE_LOG:			\
	    res = log(x);		\
	    break;			\
	}				\
    }

#define prune_yconv(res,y,ytype)	\
    switch (deltatypeno) {		\
    case PRUNE_VIEWPORT:		\
	res = xy_yconv(y);			\
	break;				\
    case PRUNE_WORLD:			\
	switch (ytype) {		\
	case PRUNE_LIN:			\
	    res = y;			\
	    break;			\
	case PRUNE_LOG:			\
	    res = log(y);		\
	    break;			\
	}				\
    }

/*
 * Prune data
 */
void do_prune(int setno, int typeno, int deltatypeno, float deltax, float deltay, int dxtype, int dytype)
{
    int len, npts = 0, d, i, j, k, drop, resset;
    double *x, *y, *resx, *resy, xtmp, ytmp, ddx = 0.0, ddy = 0.0;
    double xj = 0.0, xjm = 0.0, xjp = 0.0, yj = 0.0, yjm = 0.0, yjp = 0.0;

    if (!is_set_active(get_cg(), setno)) {
        errmsg("Set not active");
        return;
    }
    len = getsetlength(get_cg(), setno);
    if (len <= 2) {
	errmsg("Set has <= 2 points");
	return;
    }
    x = getx(get_cg(), setno);
    y = gety(get_cg(), setno);
    switch (typeno) {
    case PRUNE_CIRCLE:
    case PRUNE_ELLIPSE:
    case PRUNE_RECTANGLE:
	deltax = fabs(deltax);
	if (deltax == 0)
	    return;
	break;
    }
    switch (typeno) {
    case PRUNE_CIRCLE:
	deltay = deltax;
	break;
    case PRUNE_ELLIPSE:
    case PRUNE_RECTANGLE:
    case PRUNE_INTERPOLATION:
	deltay = fabs(deltay);
	if (deltay == 0)
	    return;
	break;
    }
    if (deltatypeno == PRUNE_WORLD) {
	if (dxtype == PRUNE_LOG && deltax < 1.0) {
	    deltax = 1.0 / deltax;
	}
	if (dytype == PRUNE_LOG && deltay < 1.0) {
	    deltay = 1.0 / deltay;
	}
    }
    resset = nextset(get_cg());
    if (resset < 0) {
        return;
    }
    add_point(get_cg(), resset, x[0], y[0]);
    npts++;
    resx = getx(get_cg(), resset);
    resy = gety(get_cg(), resset);
    switch (typeno) {
    case PRUNE_CIRCLE:
    case PRUNE_ELLIPSE:
	for (i = 1; i < len; i++) {
	    xtmp = x[i];
	    ytmp = y[i];
	    drop = FALSE;
	    for (j = npts - 1; j >= 0 && drop == FALSE; j--) {
		switch (deltatypeno) {
		case PRUNE_VIEWPORT:
		    ddx = (xy_xconv(xtmp) - xy_xconv(resx[j])) / deltax;
		    if (fabs(ddx) < 1.0) {
			ddy = (xy_yconv(ytmp) - xy_yconv(resy[j])) / deltay;
			if (ddx * ddx + ddy * ddy < 1.0) {
			    drop = TRUE;
			}
		    }
		    break;
		case PRUNE_WORLD:
		    switch (dxtype) {
		    case PRUNE_LIN:
			ddx = (xtmp - resx[j]) / deltax;
			break;
		    case PRUNE_LOG:
			ddx = (xtmp / resx[j]);
			if (ddx < 1.0) {
			    ddx = 1.0 / ddx;
			}
			ddx /= deltax;
			break;
		    }
		    if (fabs(ddx) < 1.0) {
			switch (dytype) {
			case PRUNE_LIN:
			    ddy = (ytmp - resy[j]) / deltay;
			    break;
			case PRUNE_LOG:
			    ddy = (ytmp / resy[j]);
			    if (ddy < 1.0) {
				ddy = 1.0 / ddy;
			    }
			    ddy /= deltay;
			    break;
			}
			if (ddx * ddx + ddy * ddy < 1.0) {
			    drop = TRUE;
			}
		    }
		    break;
		}
	    }
	    if (drop == FALSE) {
		add_point(get_cg(), resset, xtmp, ytmp);
		npts++;
		resx = getx(get_cg(), resset);
		resy = gety(get_cg(), resset);
	    }
	}
	sprintf(buf, "Prune from %d, %s dx = %g dy = %g", setno, 
	    (typeno == 0) ? "Circle" : "Ellipse", deltax, deltay);
	break;
    case PRUNE_RECTANGLE:
	for (i = 1; i < len; i++) {
	    xtmp = x[i];
	    ytmp = y[i];
	    drop = FALSE;
	    for (j = npts - 1; j >= 0 && drop == FALSE; j--) {
		switch (deltatypeno) {
		case PRUNE_VIEWPORT:
		    ddx = fabs(xy_xconv(xtmp) - xy_xconv(resx[j]));
		    if (ddx < deltax) {
			ddy = fabs(xy_yconv(ytmp) - xy_yconv(resy[j]));
			if (ddy < deltay) {
			    drop = TRUE;
			}
		    }
		    break;
		case PRUNE_WORLD:
		    switch (dxtype) {
		    case PRUNE_LIN:
			ddx = fabs(xtmp - resx[j]);
			break;
		    case PRUNE_LOG:
			ddx = (xtmp / resx[j]);
			if (ddx < 1.0) {
			    ddx = 1.0 / ddx;
			}
			break;
		    }
		    if (ddx < deltax) {
			switch (dytype) {
			case PRUNE_LIN:
			    ddy = fabs(ytmp - resy[j]);
			    break;
			case PRUNE_LOG:
			    ddy = (ytmp / resy[j]);
			    if (ddy < 1.0) {
				ddy = 1.0 / ddy;
			    }
			    break;
			}
			if (ddy < deltay) {
			    drop = TRUE;
			}
		    }
		    break;
		}
	    }
	    if (drop == FALSE) {
		add_point(get_cg(), resset, xtmp, ytmp);
		npts++;
		resx = getx(get_cg(), resset);
		resy = gety(get_cg(), resset);
	    }
	}
	sprintf(buf, "Prune from %d, %s dx = %g dy = %g", setno, 
	    "Rectangle", deltax, deltay);
	break;
    case PRUNE_INTERPOLATION:
	k = 0;
	prune_xconv(xjm, x[0], dxtype);
	prune_yconv(yjm, y[0], dytype);
	while (k < len - 2) {
	    d = 1;
	    i = k + d + 1;
	    drop = TRUE;
	    while (TRUE) {
		prune_xconv(xjp, x[i], dxtype);
		prune_yconv(yjp, y[i], dytype);
		for (j = k + 1; j < i && drop == TRUE; j++) {
		    prune_xconv(xj, x[j], dxtype);
		    prune_yconv(yj, y[j], dytype);
		    if (xjp == xjm) {
			ytmp = 0.5 * (yjp + yjm);
		    } else {
			ytmp = (yjp*(xj-xjm)+yjm*(xjp-xj))/(xjp-xjm);
		    }
		    switch (deltatypeno) {
		    case PRUNE_VIEWPORT:
			ddy = fabs(ytmp - yj);
			break;
		    case PRUNE_WORLD:
			switch (dytype) {
			case PRUNE_LIN:
			    ddy = fabs(ytmp - yj);
			    break;
			case PRUNE_LOG:
			    ddy = exp(fabs(ytmp - yj));
			    break;
			}
		    }
		    if (ddy > deltay) {
			drop = FALSE;
		    }
		}
		if (drop == FALSE || i == len - 1) {
		    break;
		}
		d *= 2;
		i = k + d + 1;
		if (i >= len) {
		    i = len - 1;
		}
	    }
	    if (drop == FALSE) {
		i = k + 1;
		drop = TRUE;
		while (d > 1) {
		    d /= 2;
		    i += d;
		    prune_xconv(xjp, x[i], dxtype);
		    prune_yconv(yjp, y[i], dytype);
		    drop = TRUE;
		    for (j = k + 1; j < i && drop == TRUE; j++) {
			prune_xconv(xj, x[j], dxtype);
			prune_yconv(yj, y[j], dytype);
			ytmp = (yjp*(xj-xjm)+yjm*(xjp-xj))/(xjp-xjm);
			switch (deltatypeno) {
			case PRUNE_VIEWPORT:
			    ddy = fabs(ytmp - yj);
			    break;
			case PRUNE_WORLD:
			    switch (dytype) {
			    case PRUNE_LIN:
				ddy = fabs(ytmp - yj);
				break;
			    case PRUNE_LOG:
				ddy = exp(fabs(ytmp - yj));
				break;
			    }
			}
			if (ddy > deltay) {
			    drop = FALSE;
			}
		    }
		    if (drop == FALSE) {
			i -= d;
		    }
		}
	    }
	    k = i;
	    prune_xconv(xjm, x[k], dxtype);
	    prune_yconv(yjm, y[k], dytype);
	    add_point(get_cg(), resset, x[k], y[k]);
	    npts++;
	    resx = getx(get_cg(), resset);
	    resy = gety(get_cg(), resset);
	}
	if (k == len - 2) {
	    add_point(get_cg(), resset, x[len-1], y[len-1]);
	    npts++;
	}
	sprintf(buf, "Prune from %d, %s dy = %g", setno, 
	    "Interpolation", deltay);
	break;
    }
    setcomment(get_cg(), resset, buf);
}

int get_points_inregion(int rno, int invr, int len, double *x, double *y, int *cnt, double **xt, double **yt)
{
    int i, clen = 0;
    double *xtmp, *ytmp;
    *cnt = 0;
    if (isactive_region(rno)) {
	for (i = 0; i < len; i++) {
	    if (invr) {
		if (!inregion(rno, x[i], y[i])) {
		    clen++;
		}
	    } else {
		if (inregion(rno, x[i], y[i])) {
		    clen++;
		}
	    }
	}
	if (clen == 0) {
	    return 0;
	}
	xtmp = (double *) xcalloc(clen, SIZEOF_DOUBLE);
	if (xtmp == NULL) {
	    return 0;
	}
	ytmp = (double *) xcalloc(clen, SIZEOF_DOUBLE);
	if (ytmp == NULL) {
	    xfree(xtmp);
	    return 0;
	}
	clen = 0;
	for (i = 0; i < len; i++) {
	    if (invr) {
		if (!inregion(rno, x[i], y[i])) {
		    xtmp[clen] = x[i];
		    ytmp[clen] = y[i];
		    clen++;
		}
	    } else {
		if (inregion(rno, x[i], y[i])) {
		    xtmp[clen] = x[i];
		    ytmp[clen] = y[i];
		    clen++;
		}
	    }
	}
    } else {
	return 0;
    }
    *cnt = clen;
    *xt = xtmp;
    *yt = ytmp;
    return 1;
}

int interpolate(double *mesh, double *yint, int meshlen,
    double *x, double *y, int len, int method)
{
    double *b, *c, *d;
    double dx;
    int i, ifound;
    int m;

    /* For linear interpolation, non-strict monotonicity is fine */
    m = monotonicity(x, len, method == INTERP_LINEAR ? FALSE:TRUE);
    if (m == 0) {
        errmsg("Can't interpolate a set with non-monotonic abscissas");
        return RETURN_FAILURE;
    }

    switch (method) {
    case INTERP_SPLINE:
    case INTERP_ASPLINE:
        b = xcalloc(len, SIZEOF_DOUBLE);
        c = xcalloc(len, SIZEOF_DOUBLE);
        d = xcalloc(len, SIZEOF_DOUBLE);
        if (b == NULL || c == NULL || d == NULL) {
            xfree(b);
            xfree(c);
            xfree(d);
            return RETURN_FAILURE;
        }
        if (method == INTERP_ASPLINE){
            /* Akima spline */
            aspline(len, x, y, b, c, d);
        } else {
            /* Plain cubic spline */
            spline(len, x, y, b, c, d);
        }

	seval(mesh, yint, meshlen, x, y, b, c, d, len);

        xfree(b);
        xfree(c);
        xfree(d);
        break;
    default:
        /* linear interpolation */

        for (i = 0; i < meshlen; i++) {
            ifound = find_span_index(x, len, m, mesh[i]);
            if (ifound < 0) {
                ifound = 0;
            } else if (ifound > len - 2) {
                ifound = len - 2;
            }
            dx = x[ifound + 1] - x[ifound];
            if (dx != 0.0) {
                yint[i] = y[ifound] + (mesh[i] - x[ifound])*
                    ((y[ifound + 1] - y[ifound])/dx);
            } else {
                yint[i] = (y[ifound] + y[ifound + 1])/2;
            }
        }
        break;
    }
    
    return RETURN_SUCCESS;
}

/* interpolate a set at abscissas from mesh
 * method - type of spline (or linear interpolation)
 * if strict is set, perform interpolation only within source set bounds
 * (i.e., no extrapolation)
 */
int do_interp(int gno_src, int setno_src, int gno_dest, int setno_dest,
    double *mesh, int meshlen, int method, int strict)
{
    int len, n, ncols;
    double *x, *xint;
    char *s;
	
    if (!is_valid_setno(gno_src, setno_src)) {
	errmsg("Interpolated set not active");
	return RETURN_FAILURE;
    }
    if (mesh == NULL || meshlen < 1) {
        errmsg("NULL sampling mesh");
        return RETURN_FAILURE;
    }
    
    len = getsetlength(gno_src, setno_src);
    ncols = dataset_cols(gno_src, setno_src);

    if (setno_dest == SET_SELECT_NEXT) {
	setno_dest = nextset(gno_dest);
    }
    if (!is_valid_setno(gno_dest, setno_dest)) {
	errmsg("Can't activate destination set");
	return RETURN_FAILURE;
    }

    if (dataset_cols(gno_dest, setno_dest) != ncols) {
        copyset(gno_src, setno_src, gno_dest, setno_dest);
    }
    
    setlength(gno_dest, setno_dest, meshlen);
    activateset(gno_dest, setno_dest);
    
    x = getcol(gno_src, setno_src, DATA_X);
    for (n = 1; n < ncols; n++) {
        int res;
        double *y, *yint;
        
        y    = getcol(gno_src, setno_src, n);
        yint = getcol(gno_dest, setno_dest, n);
        
        res = interpolate(mesh, yint, meshlen, x, y, len, method);
        if (res != RETURN_SUCCESS) {
            killset(gno_dest, setno_dest);
            return RETURN_FAILURE;
        }
    }

    xint = getcol(gno_dest, setno_dest, DATA_X);
    memcpy(xint, mesh, meshlen*SIZEOF_DOUBLE);

    if (strict) {
        double xmin, xmax;
        int i, imin, imax;
        minmax(x, len, &xmin, &xmax, &imin, &imax);
        for (i = meshlen - 1; i >= 0; i--) {
            if (xint[i] < xmin || xint[i] > xmax) {
                del_point(gno_dest, setno_dest, i);
            }
        }
    }
    
    switch (method) {
    case INTERP_SPLINE:
        s = "cubic spline";
        break;
    case INTERP_ASPLINE:
        s = "Akima spline";
        break;
    default:
        s = "linear interpolation";
        break;
    }
    sprintf(buf, "Interpolated from G%d.S%d using %s", gno_src, setno_src, s);
    setcomment(gno_dest, setno_dest, buf);
    
    return RETURN_SUCCESS;
}

int get_restriction_array(int gno, int setno,
    int rtype, int negate, char **rarray)
{
    int i, n, regno;
    double *x, *y;
    world w;
    WPoint wp;
    
    if (rtype == RESTRICT_NONE) {
        *rarray = NULL;
        return RETURN_SUCCESS;
    }
    
    n = getsetlength(gno, setno);
    if (n <= 0) {
        *rarray = NULL;
        return RETURN_FAILURE;
    }
    
    *rarray = xmalloc(n*SIZEOF_CHAR);
    if (*rarray == NULL) {
        return RETURN_FAILURE;
    }
    
    x = getcol(gno, setno, DATA_X);
    y = getcol(gno, setno, DATA_Y);
    
    switch (rtype) {
    case RESTRICT_REG0:
    case RESTRICT_REG1:
    case RESTRICT_REG2:
    case RESTRICT_REG3:
    case RESTRICT_REG4:
        regno = rtype - RESTRICT_REG0;
        for (i = 0; i < n; i++) {
            (*rarray)[i] = inregion(regno, x[i], y[i]) ? !negate : negate;
        }
        break;
    case RESTRICT_WORLD:
        get_graph_world(gno, &w);
        for (i = 0; i < n; i++) {
            wp.x = x[i];
            wp.y = y[i];
            (*rarray)[i] = is_wpoint_inside(&wp, &w) ? !negate : negate;
        }
        break;
    default:
        errmsg("Internal error in get_restriction_array()");
        XCFREE(*rarray);
        return RETURN_FAILURE;
        break;
    }
    return RETURN_SUCCESS;
}

int monotonicity(double *array, int len, int strict)
{
    int i;
    int s0, s1;
    
    if (len < 2) {
        errmsg("Monotonicity of an array of length < 2 is meaningless");
        return 0;
    }
    
    s0 = sign(array[1] - array[0]);
    for (i = 2; i < len; i++) {
        s1 = sign(array[i] - array[i - 1]);
        if (s1 != s0) {
            if (strict) {
                return 0;
            } else if (s0 == 0) {
                s0 = s1;
            } else if (s1 != 0) {
                return 0;
            }
        }
    }
    
    return s0;
}

int find_span_index(double *array, int len, int m, double x)
{
    int ind, low = 0, high = len - 1;
    
    if (len < 2 || m == 0) {
        errmsg("find_span_index() called with a non-monotonic array");
        return -2;
    } else if (m > 0) {
        /* ascending order */
        if (x < array[0]) {
            return -1;
        } else if (x > array[len - 1]) {
            return len - 1;
        } else {
            while (low <= high) {
	        ind = (low + high) / 2;
	        if (x < array[ind]) {
	            high = ind - 1;
	        } else if (x > array[ind + 1]) {
	            low = ind + 1;
	        } else {
	            return ind;
	        }
            }
        }
    } else {
        /* descending order */
        if (x > array[0]) {
            return -1;
        } else if (x < array[len - 1]) {
            return len - 1;
        } else {
            while (low <= high) {
	        ind = (low + high) / 2;
	        if (x > array[ind]) {
	            high = ind - 1;
	        } else if (x < array[ind + 1]) {
	            low = ind + 1;
	        } else {
	            return ind;
	        }
            }
        }
    }

    /* should never happen */
    errmsg("internal error in find_span_index()");
    return -2;
}