/*
 * 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.
 * See also do_fourier in fourier.c
 */

#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"

#include "computils.h"

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 comment_and_legend_set (int gno ,int setno ,char *buf)
{
  setcomment (gno ,setno ,buf);
  if (legend_for_new_set) set_legend_string (gno ,setno ,buf);
}

/*
 * 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_z_scaled (graphto ,loadto ,FALSE);
	    set_dirtystate();
            return RETURN_SUCCESS;
        }
    } else {
        return RETURN_FAILURE;
    }
}

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));
	comment_and_legend_set (get_cg(), digfiltset, buf);
    }
}

/*
 * linear convolution
 */
void do_linearc(int gno1, int set1, int gno2, int set2)
{
    int linearcset, i, itmp, cg = get_cg();
    double *xtmp;

    if (!(is_set_active(gno1, set1) && is_set_active(gno2, set2))) {
	errmsg("Set not active");
	return;
    }
    if ((getsetlength(gno1, set1) < 3) || (getsetlength(gno2, set2) < 3)) {
	errmsg("Set length < 3");
	return;
    }
    linearcset = nextset(cg);
    if (linearcset != (-1)) {
  	activateset(cg, linearcset);
	setlength(cg, linearcset, (itmp = getsetlength(gno1, set1) + getsetlength(gno2, set2) - 1));
	linearconv(gety(gno2, set2), gety(gno1, set1), gety(cg, linearcset), getsetlength(gno2, set2), getsetlength(gno1, set1));
	xtmp = getx(cg, linearcset);
	for (i = 0; i < itmp; i++) {
	    xtmp[i] = i;
	}
	sprintf(buf, "Linear convolution of set %d with set %d", set1, set2);
	comment_and_legend_set (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;
	}
	comment_and_legend_set (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));
	    comment_and_legend_set (gno, intset, buf);
	}
    } else {
	sum = trapint(getx(gno, setno), gety(gno, setno), NULL, NULL, getsetlength(gno, setno));
    }
    return sum;
}

/*
 * 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);
	comment_and_legend_set (get_cg(), diffset, buf);
    }
}


/****** Finite differencing/derivating *******/

static DIFFParams diff;

void diff_set_defaults (void)
{
  diff.derivative = DIFF_TYPE_DERIVATIVE;
  diff.xplace     = DIFF_XPLACE_CENTER;
  diff.period     = 1;
}

int diff_params_get (DIFFParams *p)
{
  if (p != NULL) {
    memcpy (p ,&diff ,sizeof(DIFFParams));
    return RETURN_SUCCESS;
  } else {
    return RETURN_FAILURE;
  }
}
int diff_params_set (DIFFParams *p)
{
  if (p != NULL) {
    memcpy (&diff ,p ,sizeof(DIFFParams));
    return RETURN_SUCCESS;
  } else {
    return RETURN_FAILURE;
  }
}

/**
 *   replace do_differ
 */
int do_diff (int gno1 ,int setno1 ,int gno2 ,int setno2)
{
  int i ,inlen ,outlen;
  double *in_x ,*in_y ,h ,*out_x ,*out_y;
  char *stype ,pbuf[32] ,buf[256];

  /* get input and perform safety check on it */
  inlen = getsetlength (gno1 ,setno1);
  if (inlen < 3) {
    errmsg("Set length < 3");
    return RETURN_FAILURE;
  }
  if (diff.period < 1) {
    errmsg("Non-positive period");
    return RETURN_FAILURE;
  }
  outlen = inlen - diff.period;
  if (outlen <= 0) {
    errmsg("Source set length <= differentiation period");
    return RETURN_FAILURE;
  }
  in_x = getcol (gno1 ,setno1 ,DATA_X);
  if (diff.derivative) {
    for (i = 0; i < outlen; i++) {
      if (in_x[i + diff.period] - in_x[i] == 0.0) {
	sprintf(buf ,"Can't evaluate derivative, in_x[%d] = in_x[%d]" ,i ,i + diff.period);
	errmsg(buf);
	return RETURN_FAILURE;
      }
    }
  }
  in_y = getcol (gno1 ,setno1 ,DATA_Y);


  /******* output preparation ********/

  /* Reallocate set(s) */
  setlength (gno2 ,setno2 ,outlen);
  out_x = getx (gno2 ,setno2);
  out_y = gety (gno2 ,setno2);
  /* Compute abscissas */
  for (i = 0; i < outlen; i++) {
    switch (diff.xplace) {
    case DIFF_XPLACE_LEFT:
      out_x[i] = in_x[i];
      break;
    case DIFF_XPLACE_RIGHT:
      out_x[i] = in_x[i+diff.period];
      break;
    case DIFF_XPLACE_CENTER:
      out_x[i] = (in_x[i+diff.period] + in_x[i]) / 2.0;
      break;
    }
  }
  /* EFFECTIVE computation of difference/derivative */
  for (i = 0; i < outlen; i++) {
    out_y[i] = in_y[i + diff.period] - in_y[i];
    if (diff.derivative) {
      h = in_x[i + diff.period] - in_x[i];
      out_y[i] /= h;
    }
  }

  /* Write appropriate comment */
  switch (diff.xplace) {
  case DIFF_XPLACE_LEFT:
    stype = "Left";
    break;
  case DIFF_XPLACE_CENTER:
    stype = "Centered";
    break;
  case DIFF_XPLACE_RIGHT:
    stype = "Right";
    break;
  default:
    errmsg ("Wrong xplace passed to do_diff()");
    return RETURN_FAILURE;
    break;
  }
  if (diff.period != 1) {
    sprintf (pbuf ," (period = %d)" ,diff.period);
  } else {
    pbuf[0] = '\0';
  }
  sprintf (buf ,"%s %s%s of G%d.s%d" ,stype ,pbuf,
	   diff.derivative ? "derivative":"difference" ,gno1 ,setno1 );
  comment_and_legend_set (gno2 ,setno2 ,buf);
  return RETURN_SUCCESS;
}

void do_diff_command (int gno1 ,int setno1 ,DIFF_xplace xplace)
{
  diff.xplace = xplace;
  int setno2  = nextset (gno1);
  do_copyset (gno1 ,setno1 ,gno1 ,setno2);
  do_diff    (gno1 ,setno1 ,gno1 ,setno2);
}

void do_diff_command_1 (int gno1 ,int setno1 ,DIFF_xplace xplace
			,DIFF_type   derivative
			,int         period
			)
{
  diff.xplace = xplace;
  int setno2  = nextset (gno1);
  do_copyset (gno1 ,setno1 ,gno1 ,setno2);
  do_diff    (gno1 ,setno1 ,gno1 ,setno2);
}

/*
 * 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);
	comment_and_legend_set (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
 *     rno is the region number
 *  the output sets setno2 and setno22 must exists before the call
 */


void do_runavg (int gno, int setno, int runlen, int runtype, int rno, int invr)
{
  int setno2  = nextset (gno);
  do_copyset (gno ,setno ,gno ,setno2);
  do_runavg12 (gno ,setno ,gno ,setno2 ,runlen ,runtype ,rno ,invr);
}

void do_runavg12 (int gno1 ,int setno1 ,int gno2 ,int setno2
		  ,int runlen ,int runtype ,int rno ,int invr)
{
  int inlen ,outlen ,buflen = 0;
  double *in_x ,*in_y ,*buf_x ,*buf_y ,*out_x ,*out_y;

  /* get input and perform safety check on it */
  if (runlen < 2) {
    errmsg("Length of running average < 2");
    return;
  }
  inlen = getsetlength (gno1 ,setno1);
  if (inlen < 2) {
    errmsg("Set length < 2");
    return;
  }
  in_x = getcol (gno1 ,setno1 ,DATA_X);
  in_y = getcol (gno1 ,setno1 ,DATA_Y);
  if (rno == -1) {
    buf_x = in_x;
    buf_y = in_y;
    buflen = inlen;
  } else if (isactive_region(rno)) {
    if (!get_points_inregion (rno, invr, inlen, in_x, in_y, &buflen, &buf_x, &buf_y)) {
      if (buflen == 0) {
	errmsg ("No points found in region, operation cancelled");
      }
      return;
    }
  } else {
    errmsg ("Selected region is not active");
    return;
  }
  if (runlen >= buflen) {
    errmsg("Length of running average > length of set or region");
    goto bustout;
  }

  /* Reallocate set */
  outlen = buflen - runlen + 1;
  setlength (gno2 ,setno2 ,outlen);
  out_x = getx (gno2 ,setno2);
  out_y = gety (gno2 ,setno2);

  /* EFFECTIVE COMPUTATION of the averaged set */
  switch (runtype) {
  case 0:
    runavg (buf_x ,buf_y ,out_x ,out_y ,buflen ,runlen);
    sprintf(buf ,"%d-pt. avg. on G%d.s%d " ,runlen ,gno1 ,setno1);
    break;
  case 1:
    runmedian (buf_x ,buf_y ,out_x ,out_y ,buflen ,runlen);
    sprintf(buf ,"%d-pt. median on set %d " ,runlen ,setno1);
    break;
  case 2:
    runminmax (buf_x ,buf_y ,out_x ,out_y ,buflen ,runlen ,0);
    sprintf(buf ,"%d-pt. min on set %d " ,runlen ,setno1);
    break;
  case 3:
    runminmax (buf_x ,buf_y ,out_x ,out_y ,buflen ,runlen ,1);
    sprintf(buf ,"%d-pt. max on set %d " ,runlen ,setno1);
    break;
  case 4:
    runstddev (buf_x ,buf_y ,out_x ,out_y ,buflen ,runlen);
    sprintf(buf ,"%d-pt. std dev. ,set %d " ,runlen ,setno1);
    break;
  }
  comment_and_legend_set (gno2 ,setno2 ,buf);

bustout:;
if (rno >= 0 && buflen != 0) {	/* had a region and allocated memory there */
  xfree (buf_x);
  xfree (buf_y);
 }


}




/** 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 */
	}
	comment_and_legend_set (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, const 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) {
	comment_and_legend_set (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;
    }
    comment_and_legend_set (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;
}

static int interpolate(double *mesh, double *yint, int meshlen,
    double *x, double *y, int len, int method ,double xspline_s)
{
    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("Only XY-spline interpolation allows 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;
    case INTERP_XSPLINE:
      xspline_eval_y (mesh ,yint ,meshlen ,x ,y ,len ,xspline_s);
      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 ,double xspline_s)
{
  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 ,xspline_s);
	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;
    case INTERP_XSPLINE:
        s = "X-spline";
        break;
    default:
        s = "linear interpolation";
        break;
    }
    sprintf(buf, "%s (G%d.S%d)",s , gno_src, setno_src);
    comment_and_legend_set(gno_dest, setno_dest, buf);

    return RETURN_SUCCESS;
}


int do_xyspline (int gno_src ,int setno_src ,double xspline_s ,double xspline_p)
{
  int setlen ,setno_dest ,n_interp;
  double *x ,*y ,*xint ,*yint;
  int gno_dest = get_cg ();

  if (!is_valid_setno (gno_src ,setno_src)) {
    errmsg ("Interpolated set not active");
    return RETURN_FAILURE;
  }

  setlen     = getsetlength (gno_src ,setno_src);
  setno_dest = nextset (gno_dest);

  if (!is_valid_setno (gno_dest ,setno_dest)) {
    errmsg ("Can't activate destination set");
    return RETURN_FAILURE;
  }
  xint = xcalloc (setlen ,SIZEOF_DOUBLE);
  yint = xcalloc (setlen ,SIZEOF_DOUBLE);
  x    = getx (gno_src ,setno_src);
  y    = gety (gno_src ,setno_src);

  if (xspline_eval_xy (&xint ,&yint ,setlen ,x ,y ,setlen ,&n_interp ,xspline_s ,xspline_p) == RETURN_FAILURE) {
    xfree (xint);
    xfree (yint);
    return RETURN_FAILURE;
  }
  copyset     (gno_src ,setno_src ,gno_dest ,setno_dest);
  killsetdata (gno_dest ,setno_dest);
  setcol      (gno_dest ,setno_dest ,DATA_X ,xint ,n_interp);
  setcol      (gno_dest ,setno_dest ,DATA_Y ,yint ,n_interp);
  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;
}

/**
 * Picked from grace-5.99.dev5
 */
int monospaced(double *array, int len, double *space)
{
    int i;
    double eps;

    if (len < 2) {
        errmsg("Monospacing of an array of length < 2 is meaningless");
        return FALSE;
    }

    *space = array[1] - array[0];
    eps = fabs((array[len - 1] - array[0]))*1.0e-6; /* FIXME */
    for (i = 2; i < len; i++) {
        if (fabs(array[i] - array[i - 1] - *space) > eps) {
            return FALSE;
        }
    }

    return TRUE;
}
