/*
* 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.
*/
/*
*
* curve fitting, and other numerical routines used in compose.
*
* Contents:
*
* void stasum() - compute mean and variance
* double leasev() - evaluate least squares polynomial
* int fitcurve() - compute coefficients for a polynomial fit of degree >1
* void runavg() - compute a running average
* void runstddev() - compute a running standard deviation
* void runmedian() - compute a running median
* void runminmax() - compute a running minimum or maximum
* void filterser() - apply a digital filter
* void linearconv() - convolve one set with another
* int crosscorr() - cross/auto correlation
* int linear_regression() - linear regression
* void spline() - compute a spline fit
* int seval() - evaluate the spline computed in spline()
*/
#include <config.h>
#include <cmath.h>
#include <stdio.h>
#include <stdlib.h>
#include "defines.h"
#include "utils.h"
#include "protos.h"
#include "as274c.h"
static char buf[256];
/*
compute mean and standard dev
*/
void stasum(double *x, int n, double *xbar, double *sd)
{
int i;
*xbar = 0;
*sd = 0;
if (x == NULL) {
return;
}
if (n < 1) {
return;
}
for (i = 0; i < n; i++) {
*xbar = (*xbar) + x[i];
}
*xbar = (*xbar)/n;
if (n > 1) {
for (i = 0; i < n; i++) {
*sd = (*sd) + (x[i] - *xbar) * (x[i] - *xbar);
}
*sd = sqrt(*sd/(n - 1));
}
}
/*
evaluate least squares polynomial
*/
double leasev(double *c, int degree, double x)
{
double temp;
int i;
/*
* Commented out - refer to "Numerical Methods in C" for the
* reason, section 5.2! The revolution has begun.
temp = 0.0;
for (i = 0; i <= degree; i++) {
if ((i == 0) && (x == 0.0)) {
temp = temp + c[i];
} else {
temp = temp + c[i] * pow(x, (double) (i));
}
} */
temp = c[degree];
for( i=degree-1; i>=0; i-- )
temp = c[i] + temp*x;
return (temp);
}
/*
polynomial curve fitting
*/
int fitcurve(double *x, double *y, int n, int ideg, double *coeff)
/*
* x - x values
* y - y values
* n - number of points
* ideg - degree of fit
* coeff - coefficients of fit
*/
{
int i, ifail;
ifail = 1;
if (ideg > 1) {
dofitcurve(n, x, y, ideg, coeff);
} else {
ifail = linear_regression(n, x, y, coeff);
if (ifail == 1) {
errmsg("Linear_regression entered with N < 2");
return ifail;
} else if (ifail == 2) {
errmsg("Linear_regression - all values of x or y are the same");
return ifail;
}
}
/* check coefficients */
for (i = 0; i <= ideg; i++) {
if (!finite(coeff[i])) {
errmsg("Linear_regression - all values of x or y are the same");
ifail = 3;
return ifail;
}
}
ifail = 0;
return ifail;
}
/*
compute a running average
*/
void runavg(double *x, double *y, double *ax, double *ay, int n, int ilen)
{
int i;
double sumy = 0.0;
double sumx = 0.0;
for (i = 0; i < ilen; i++) {
sumx = sumx + x[i];
sumy = sumy + y[i];
}
ax[0] = sumx / ilen;
ay[0] = sumy / ilen;
for (i = 1; i < (n - ilen + 1); i++) {
sumx = x[i + ilen - 1] - x[i - 1] + sumx;
ax[i] = sumx / ilen;
sumy = y[i + ilen - 1] - y[i - 1] + sumy;
ay[i] = sumy / ilen;
}
}
/*
compute a running standard deviation
*/
void runstddev(double *x, double *y, double *ax, double *ay, int n, int ilen)
{
int i;
double ybar, ysd;
double sumx = 0.0;
for (i = 0; i < ilen; i++) {
sumx = sumx + x[i];
}
ax[0] = sumx / ilen;
stasum(y, ilen, &ybar, &ysd);
ay[0] = ysd;
for (i = 1; i < (n - ilen + 1); i++) {
stasum(y + i, ilen, &ybar, &ysd);
sumx = x[i + ilen - 1] - x[i - 1] + sumx;
ax[i] = sumx / ilen;
ay[i] = ysd;
}
}
/*
compute a running median
*/
void runmedian(double *x, double *y, double *ax, double *ay, int n, int ilen)
{
int i, j, nlen = n - ilen + 1;
double *tmpx, *tmpy;
tmpx = (double *) xcalloc(ilen, sizeof(double));
if (tmpx == NULL) {
errmsg("Can't xcalloc tmpx in runmedian");
return;
}
tmpy = (double *) xcalloc(ilen, sizeof(double));
if (tmpy == NULL) {
errmsg("Can't xcalloc tmpy in runmedian");
XCFREE(tmpx);
return;
}
for (i = 0; i < nlen; i++) {
for (j = 0; j < ilen; j++) {
tmpx[j] = x[j + i];
tmpy[j] = y[j + i];
}
sort_xy(tmpx, tmpy, ilen, 1, 0);
if (ilen % 2) {
ax[i] = x[i + (ilen / 2)];
ay[i] = tmpy[ilen / 2];
} else {
ax[i] = (x[i + ilen / 2] + x[i + (ilen - 1) / 2]) * 0.5;
ay[i] = (tmpy[ilen / 2] + tmpy[(ilen - 1) / 2]) * 0.5;
}
}
XCFREE(tmpx);
XCFREE(tmpy);
}
/*
compute a running minimum or maximum
*/
void runminmax(double *x, double *y, double *ax, double *ay, int n, int ilen, int type)
{
int i, j;
double min, max;
double sumx = 0.0;
min = max = y[0];
for (i = 0; i < ilen; i++) {
sumx = sumx + x[i];
if (min > y[i])
min = y[i];
if (max < y[i])
max = y[i];
}
ax[0] = sumx / ilen;
if (type == 0) {
ay[0] = min;
} else if (type == 1) {
ay[0] = max;
} else {
errmsg("Unknown type in runminmax, setting type = min");
type = 0;
}
for (i = 1; i < (n - ilen + 1); i++) {
sumx = x[i + ilen - 1] - x[i - 1] + sumx;
ax[i] = sumx / ilen;
min = y[i];
max = y[i];
for (j = 0; j < ilen; j++) {
if (min > y[i + j])
min = y[i + j];
if (max < y[i + j])
max = y[i + j];
}
if (type == 0) {
ay[i] = min;
} else if (type == 1) {
ay[i] = max;
}
}
}
/*
Apply a digital filter of length len to a set in x, y,
of length n with the results going to resx, resy.
the length of the result is set by the caller
*/
void filterser(int n, double *x, double *y, double *resx, double *resy, double *h, int len)
{
int i, j, outlen, eo, ld2;
double sum;
outlen = n - len + 1;
eo = len % 2;
ld2 = len / 2;
for (i = 0; i < outlen; i++) {
sum = 0.0;
for (j = 0; j < len; j++) {
sum = sum + h[j] * y[j + i];
}
resy[i] = sum;
if (eo)
resx[i] = x[i + ld2];
else
resx[i] = (x[i + ld2] + x[i + ld2 - 1]) / 2.0;
}
}
/*
linear convolution of set x (length n) with h (length m) and
result to y. the length of y is set by the caller
*/
void linearconv(double *x, double *h, double *y, int n, int m)
{
int i, j, itmp;
for (i = 0; i < n + m - 1; i++) {
for (j = 0; j < m; j++) {
itmp = i - j;
if ((itmp >= 0) && (itmp < n)) {
y[i] = y[i] + h[j] * x[itmp];
}
}
}
}
/*
* cross correlation/covariance
*/
int crosscorr(double *x, double *y, int n, int maxlag, int covar, double *xres)
{
double cnorm = 1.0;
double xbar, ybar, sd;
int i, j;
if (!x || !y || !xres || n < 2 || maxlag > n) {
return RETURN_FAILURE;
}
if (covar) {
stasum(x, n, &xbar, &sd);
stasum(y, n, &ybar, &sd);
} else {
xbar = 0.0;
ybar = 0.0;
}
for (i = 0; i < maxlag; i++) {
xres[i] = 0.0;
for (j = 0; j < n - i; j++) {
xres[i] += (y[j] - ybar)*(x[j + i] - xbar);
}
if (i == 0) {
if (xres[0] != 0.0) {
cnorm = fabs(xres[0]);
}
}
xres[i] /= cnorm;
}
return RETURN_SUCCESS;
}
/*
References,
_Aplied Linear Regression_, Weisberg
_Elements of Statistical Computing_, Thisted
Fits y = coef*x + intercept + e
uses a 2 pass method for means and variances
*/
int linear_regression(int n, double *x, double *y, double *coeff)
{
double xbar, ybar; /* sample means */
double sdx, sdy; /* sample standard deviations */
double sxy, rxy; /* sample covariance and sample correlation */
double SXX, SYY, SXY; /* sums of squares */
double RSS; /* residual sum of squares */
double rms; /* residual mean square */
double sereg; /* standard error of regression */
double seslope, seintercept;
double slope, intercept; /* */
double SSreg, F, R2;
int i;
if (n < 2) {
return 1;
}
xbar = ybar = 0.0;
SXX = SYY = SXY = 0.0;
for (i = 0; i < n; i++) {
xbar = xbar + x[i];
ybar = ybar + y[i];
}
xbar = xbar / n;
ybar = ybar / n;
for (i = 0; i < n; i++) {
SXX = SXX + (x[i] - xbar) * (x[i] - xbar);
SYY = SYY + (y[i] - ybar) * (y[i] - ybar);
SXY = SXY + (x[i] - xbar) * (y[i] - ybar);
}
sdx = sqrt(SXX / (n - 1));
sdy = sqrt(SYY / (n - 1));
if (sdx == 0.0) {
return 2;
}
if (sdy == 0.0) {
return 2;
}
sxy = SXY / (n - 1);
rxy = sxy / (sdx * sdy);
slope = SXY / SXX;
intercept = ybar - slope * xbar;
RSS = SYY - slope * SXY;
sprintf(buf, "Number of observations\t\t\t = %d\n", n);
stufftext(buf);
sprintf(buf, "Mean of independent variable\t\t = %.7g\n", xbar);
stufftext(buf);
sprintf(buf, "Mean of dependent variable\t\t = %.7g\n", ybar);
stufftext(buf);
sprintf(buf, "Standard dev. of ind. variable\t\t = %.7g\n", sdx);
stufftext(buf);
sprintf(buf, "Standard dev. of dep. variable\t\t = %.7g\n", sdy);
stufftext(buf);
sprintf(buf, "Correlation coefficient\t\t\t = %.7g\n", rxy);
stufftext(buf);
sprintf(buf, "Regression coefficient (SLOPE)\t\t = %.7g\n", slope);
stufftext(buf);
if (n == 2) {
coeff[1] = (y[1] - y[0])/(x[1] - x[0]);
coeff[0] = y[0] - coeff[1]*x[0];
sprintf(buf, "Regression constant (INTERCEPT)\t\t = %.7g\n", intercept);
stufftext(buf);
return 0;
}
rms = RSS / (n - 2);
sereg = sqrt(rms);
seintercept = sqrt(rms * (1.0 / n + xbar * xbar / SXX));
seslope = sqrt(rms / SXX);
SSreg = SYY - RSS;
F = SSreg / rms;
R2 = SSreg / SYY;
sprintf(buf, "Standard error of coefficient\t\t = %.7g\n", seslope);
stufftext(buf);
sprintf(buf, "t - value for coefficient\t\t = %.7g\n", slope / seslope);
stufftext(buf);
sprintf(buf, "Regression constant (INTERCEPT)\t\t = %.7g\n", intercept);
stufftext(buf);
sprintf(buf, "Standard error of constant\t\t = %.7g\n", seintercept);
stufftext(buf);
sprintf(buf, "t - value for constant\t\t\t = %.7g\n", intercept / seintercept);
stufftext(buf);
sprintf(buf, "\nAnalysis of variance\n");
stufftext(buf);
sprintf(buf, "Source\t\t d.f\t Sum of squares\t Mean Square\t F\n");
stufftext(buf);
sprintf(buf, "Regression\t 1\t%.7g\t%.7g\t%.7g\n", SSreg, SSreg, F);
stufftext(buf);
sprintf(buf, "Residual\t%5d\t%.7g\t%.7g\n", n - 2, RSS, rms);
stufftext(buf);
sprintf(buf, "Total\t\t%5d\t%.7g\n\n", n - 1, SYY);
stufftext(buf);
for (i = 0; i < n; i++) {
coeff[0] = intercept;
coeff[1] = slope;
}
return 0;
}
/*
a literal translation of the spline routine in
Forsyth, Malcolm, and Moler
*/
void spline(int n, double *x, double *y, double *b, double *c, double *d)
{
/*
c
c the coefficients b(i), c(i), and d(i), i=1,2,...,n are computed
c for a cubic interpolating spline
c
c s(x) = y(i) + b(i)*(x-x(i)) + c(i)*(x-x(i))**2 + d(i)*(x-x(i))**3
c
c for x(i) .le. x .le. x(i+1)
c
c input..
c
c n = the number of data points or knots (n.ge.2)
c x = the abscissas of the knots in strictly increasing order
c y = the ordinates of the knots
c
c output..
c
c b, c, d = arrays of spline coefficients as defined above.
c
c using p to denote differentiation,
c
c y(i) = s(x(i))
c b(i) = sp(x(i))
c c(i) = spp(x(i))/2
c d(i) = sppp(x(i))/6 (derivative from the right)
c
c the accompanying function subprogram seval can be used
c to evaluate the spline.
c
c
*/
int nm1, ib, i;
double t;
/*
Gack!
*/
x--;
y--;
b--;
c--;
d--;
/*
Fortran 66
*/
nm1 = n - 1;
if (n < 2)
return;
if (n < 3)
goto l50;
/*
c
c set up tridiagonal system
c
c b = diagonal, d = offdiagonal, c = right hand side.
c
*/
d[1] = x[2] - x[1];
c[2] = (y[2] - y[1]) / d[1];
for (i = 2; i <= nm1; i++) {
d[i] = x[i + 1] - x[i];
b[i] = 2.0 * (d[i - 1] + d[i]);
c[i + 1] = (y[i + 1] - y[i]) / d[i];
c[i] = c[i + 1] - c[i];
}
/*
c
c end conditions. third derivatives at x(1) and x(n)
c obtained from divided differences
c
*/
b[1] = -d[1];
b[n] = -d[n - 1];
c[1] = 0.0;
c[n] = 0.0;
if (n == 3)
goto l15;
c[1] = c[3] / (x[4] - x[2]) - c[2] / (x[3] - x[1]);
c[n] = c[n - 1] / (x[n] - x[n - 2]) - c[n - 2] / (x[n - 1] - x[n - 3]);
c[1] = c[1] * d[1] * d[1] / (x[4] - x[1]);
c[n] = -c[n] * d[n - 1] * d[n - 1] / (x[n] - x[n - 3]);
/*
c
c forward elimination
c
*/
l15:;
for (i = 2; i <= n; i++) {
t = d[i - 1] / b[i - 1];
b[i] = b[i] - t * d[i - 1];
c[i] = c[i] - t * c[i - 1];
}
/*
c
c back substitution
c
*/
c[n] = c[n] / b[n];
for (ib = 1; ib <= nm1; ib++) {
i = n - ib;
c[i] = (c[i] - d[i] * c[i + 1]) / b[i];
}
/*
c
c c(i) is now the sigma(i) of the text
c
c compute polynomial coefficients
c
*/
b[n] = (y[n] - y[nm1]) / d[nm1] + d[nm1] * (c[nm1] + 2.0 * c[n]);
for (i = 1; i <= nm1; i++) {
b[i] = (y[i + 1] - y[i]) / d[i] - d[i] * (c[i + 1] + 2.0 * c[i]);
d[i] = (c[i + 1] - c[i]) / d[i];
c[i] = 3.0 * c[i];
}
c[n] = 3.0 * c[n];
d[n] = d[n - 1];
return;
l50:;
b[1] = (y[2] - y[1]) / (x[2] - x[1]);
c[1] = 0.0;
d[1] = 0.0;
b[2] = b[1];
c[2] = 0.0;
d[2] = 0.0;
return;
}
/***************************************************************************
* aspline - modified version of David Frey's spline.c *
* *
* aspline does an Akima spline interpolation. *
***************************************************************************/
void aspline(int n, double *x, double *y, double *b, double *c, double *d)
{
int i;
double num, den;
double m_m1, m_m2, m_p1, m_p2;
double x_m1, x_m2, x_p1, x_p2;
double y_m1, y_m2, y_p1, y_p2;
#define dx(i) (x[i+1]-x[i])
#define dy(i) (y[i+1]-y[i])
#define m(i) (dy(i)/dx(i))
if (n > 0) /* we have data to process */
{
/*
* calculate the coefficients of the spline
* (the Akima interpolation itself)
*/
/* b) interpolate the missing points: */
x_m1 = x[0] + x[1] - x[2];
y_m1 = (x[0]-x_m1) * (m(1) - 2 * m(0)) + y[0];
m_m1 = (y[0]-y_m1)/(x[0]-x_m1);
x_m2 = 2 * x[0] - x[2];
y_m2 = (x_m1-x_m2) * (m(0) - 2 * m_m1) + y_m1;
m_m2 = (y_m1-y_m2)/(x_m1-x_m2);
x_p1 = x[n-1] + x[n-2] - x[n-3];
y_p1 = (2 * m(n-2) - m(n-3)) * (x_p1 - x[n-1]) + y[n-1];
m_p1 = (y_p1-y[n-1])/(x_p1-x[n-1]);
x_p2 = 2 * x[n-1] - x[n-3];
y_p2 = (2 * m_p1 - m(n-2)) * (x_p2 - x_p1) + y_p1;
m_p2 = (y_p2-y_p1)/(x_p2-x_p1);
/* i = 0 */
num=fabs(m(1) - m(0)) * m_m1 + fabs(m_m1 - m_m2) * m(0);
den=fabs(m(1) - m(0)) + fabs(m_m1 - m_m2);
if (den != 0.0) b[0]=num / den;
else b[0]=0.0;
/* i = 1 */
num=fabs(m(2) - m(1)) * m(0) + fabs(m(0) - m_m1) * m(1);
den=fabs(m(2) - m(1)) + fabs(m(0) - m_m1);
if (den != 0.0) b[1]=num / den;
else b[1]=0.0;
for (i=2; i < n-2; i++)
{
num=fabs(m(i+1) - m(i)) * m(i-1) + fabs(m(i-1) - m(i-2)) * m(i);
den=fabs(m(i+1) - m(i)) + fabs(m(i-1) - m(i-2));
if (den != 0.0) b[i]=num / den;
else b[i]=0.0;
}
/* i = n - 2 */
num=fabs(m_p1 - m(n-2)) * m(n-3) + fabs(m(n-3) - m(n-4)) * m(n-2);
den=fabs(m_p1 - m(n-2)) + fabs(m(n-3) - m(n-4));
if (den != 0.0) b[n-2]=num / den;
else b[n-2]=0.0;
/* i = n - 1 */
num=fabs(m_p2 - m_p1) * m(n-2) + fabs(m(n-2) - m(n-3)) * m_p1;
den=fabs(m_p2 - m_p1) + fabs(m(n-2) - m(n-3));
if (den != 0.0) b[n-1]=num / den;
else b[n-1]=0.0;
for (i=0; i < n-1; i++)
{
double dxv = dx(i);
c[i]=(3 * m(i) - 2 * b[i] - b[i+1]) / dxv;
d[i]=(b[i] + b[i+1] - 2 * m(i)) / (dxv * dxv);
}
}
#undef dx
#undef dy
#undef m
}
int seval(double *u, double *v, int ulen,
double *x, double *y, double *b, double *c, double *d, int n)
{
/*
*
* this subroutine evaluates the cubic spline function on a mesh
*
* seval = y(i) + b(i)*(u-x(i)) + c(i)*(u-x(i))**2 + d(i)*(u-x(i))**3
*
* where x(i) .lt. u .lt. x(i+1), using horner's rule
*
* if u .lt. x(1) then i = 1 is used.
* if u .ge. x(n) then i = n is used.
*
* input..
*
* u = the array of abscissas at which the spline is to be evaluated
* ulen = length of the mesh
*
* x,y = the arrays of data abscissas and ordinates
* b,c,d = arrays of spline coefficients computed by spline
* n = the number of data points
*
* output..
*
* v = the array of evaluated values
*/
int j, m;
m = monotonicity(x, n, FALSE);
if (m == 0) {
errmsg("seval() called with a non-monotonic array");
return RETURN_FAILURE;
}
for (j = 0; j < ulen; j++) {
double dx;
int ifound;
ifound = find_span_index(x, n, m, u[j]);
if (ifound < 0) {
ifound = 0;
} else if (ifound > n - 2) {
ifound = n - 1;
}
dx = u[j] - x[ifound];
v[j] = y[ifound] + dx*(b[ifound] + dx*(c[ifound] + dx*d[ifound]));
}
return RETURN_SUCCESS;
}
syntax highlighted by Code2HTML, v. 0.9.1