/*
Following are all the files & routines that compute the extrema
(minima & maxima) of a data set, as described in the April 1996
issue of Dr. Dobbs Journal in "A Cubic Spline Algorithm for Extrema
Computation". An electronic version of the test/validation data used
in the article are also available from the author (me) electronically
upon request.
Mike J. Courtney 1996
/************************** cubic_io.c ****************************/
/*
File: cubic_io.c
Purpose: Input/output program for the FindCubicExtrema() routine.
Reads an ascii file of x,y data, calls the routine, and lists the
returned extrema. Naturally, this is not required if
FindCubicExtrema() is called directly from an application.
Author: Mike Courtney 1996
*/
#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <math.h>
#include "cubic.h"
main() {
/* the following are required for main() */
unsigned int i /* array indices */
unsigned int num_extr; /* number of extrema found */
BYTE len; /* input data parsing variable */
char data_file[80]; /* ascii data file name */
char str[40], line[80]; /* input data parsing variables */
FILE *fp; /* ascii data file pointer */
struct point *first_extr; /* ptr to the first extreme structure */
/* the following are required for FindCubicExtrema() */
unsigned int num_pnts; /* number of x,y pairs */
float *x_in, *y_in; /* arrays of x and y input data values */
struct point *extr; /* ptr to current extreme structure */
/* open input data file and determine number of input points */
printf("\n Enter name of input file: ");
gets(data_file);
if((fp = fopen(data_file, "r")) != NULL) {
num_pnts = 0;
while(!feof(fp)) {
fgets(line, 80, fp);
num_pnts++;
}
printf("\n The number of input points was %d.",num_pnts);
/* allocate the input data arrays */
x_in = malloc(num_pnts * sizeof(float));
y_in = malloc(num_pnts * sizeof(float));
/* read in the each data line and parse into x and y arrays */
rewind(fp);
for (i = 0; i < num_pnts; i++) {
/* get a line */
fgets(line, 80, fp);
len = strcspn(line,",");
/* get the x value */
strncpy(str, line, len);
/* NULL terminate the x value */
str[len] = '\0';
x_in[i] = atof(str);
/* get the y value */
strcpy(str, line+len+1);
y_in[i] = atof(str);
}
fclose(fp);
/* allocate first structure of linked list of output extrema */
extr = (struct point *)malloc(sizeof(struct point));
/* save the address of the first structure in the list */
first_extr = extr;
/* call the routine that computes the extrema */
if (FindCubicExtrema(num_pnts, x_in, y_in, extr) == FAILURE)
printf("\n\7 No extrema found !\n");
else {
/* print the linked list of extrema */
printf("\n\n Relative extrema computed:");
extr = first_extr;
num_extr = 0;
while (extr) {
printf("\n %d. y = %f at x = %f",
num_extr+1, extr->y, extr->x);
extr = extr->next;
num_extr++;
}
}
free(x_in);
free(y_in);
/* FREE THE LINKED LIST OF EXTREME STRUCTURES */
do {
/* point to first structure */
extr = first_extr;
/* save address of next extreme structure so that when we free
the current one, we still have a pointer to the next one */
first_extr = extr->next;
free(extr);
} while (first_extr != NULL);
}
else
printf("\n Couldn't open %s !", data_file);
}
/*********************** End of cubic_io.c ************************/
/************************* cubic_extrema.c ************************/
/*
File: cubic_extrema.c
Purpose: Contains FindCubicExtrema() and supporting routines that
implement the cubic spline extrema algorithm. Given a set of x,y data
points, determine in the cubic spline sense the relative extrema of
the function describing the data.
Author: Mike Courtney 1996
*/
#include <stdlib.h>
#include <math.h>
#include <stdio.h>
#include <assert.h>
#include "cubic.h"
/*
Primary routine that implements the cubic spline extrema algorithm.
Calls ComputeSecDerivs() to compute the second derivatives, computes
the quadratic coefficients, calls FindQuadRoots() to solve the
quadratic roots, determines if the roots are valid abscissa, and
calls ComputeY() to compute ordinates.
*/
BOOL FindCubicExtrema (
unsigned int num_pnts, /* input - number of x & y points */
float *x, /* input - array of x values */
float *y, /* input - array of y values */
struct point *extr) /* output - singly linked list of extrema */
{
float a, b, c; /* coefficients of quadratic equation */
float x1, x2; /* roots of quadratic equation to be computed */
unsigned int i; /* array index */
float *sec_deriv; /* second derivatives of each data interval */
BOOL root_stat; /* computation status returned by FindRoots() */
BOOL valid_flag; /* TRUE if at least one valid root found */
BOOL first_flag; /* TRUE if current root is the first one */
/* allocate array for second derivatives */
sec_deriv = malloc(num_pnts * sizeof(float));
/* compute the second derivatives */
ComputeSecDerivs(num_pnts, x, y, sec_deriv);
/* initialize extrema flags */
valid_flag = FALSE;
first_flag = TRUE;
/* loop through all the input points and find the extrema */
for (i = 0; i < num_pnts - 1; i++) {
/* compute the quadratic coefficients */
a = 3.0 * (sec_deriv[i+1] - sec_deriv[i]);
b = 6.0 * (x[i+1] * sec_deriv[i] - x[i] * sec_deriv[i+1]);
c = 6.0 * (y[i+1] - y[i]) - (2.0 * x[i+1] * x[i+1] - x[i] *
x[i] + 2 * x[i] * x[i+1]) * sec_deriv[i];
c -= (x[i+1] * x[i+1] - 2.0 * x[i] * x[i+1] - 2.0 * x[i] * x[i])
* sec_deriv[i+1];
/* determine the roots of the cubic extrema quadratic equation */
root_stat = FindQuadRoots(a, b, c, &x1, &x2);
if (root_stat != FAILURE) {
/* if root x1 was calculated successfully */
if (root_stat != FAILURE1) {
/* Determine if root is within the interval */
if ((x1 > x[i]) && (x1 < x[i+1])) {
/* first root (extremum) */
if (first_flag == TRUE)
first_flag = FALSE;
/* beyond first valid root so allocate next extremum structure */
else {
extr->next =
(struct point *)malloc(sizeof(struct point));
assert(extr->next != NULL);
extr = extr->next;
}
extr->next = NULL;
extr->x = x1;
/* compute corresponding value of y at extreme x value */
ComputeY(i, x, extr->x, y, sec_deriv, &extr->y);
valid_flag = TRUE;
}
}
/* if root x2 was calculated successfully */
if (root_stat != FAILURE2) {
/* Determine if root is within the current interval */
if ((x2 > x[i]) && (x2 < x[i+1])) {
/* first root (extremum) */
if (first_flag == TRUE)
first_flag = FALSE;
/* beyond 1st valid root so allocate next extremum structure */
else {
extr->next =
(struct point *)malloc(sizeof(struct point));
extr = extr->next;
}
extr->next = NULL;
extr->x = x2;
/* compute corresponding value of y at extreme x value */
ComputeY(i, x, extr->x, y, sec_deriv, &extr->y);
valid_flag = TRUE;
}
}
} /* end of if(root_stat != FAILURE) */
} /* end of for(i) */
free(sec_deriv);
if (valid_flag == TRUE)
return SUCCESS;
else {
/* set next to NULL just in case it was not set in the loop -
this is so that free loop will operate properly upon return */
extr->next = NULL;
return FAILURE;
}
}
/*
Use input x,y data to form tridiagonal matrix and compute second
derivatives of function in the cubic spline sense.
*/
void ComputeSecDerivs (
unsigned int num_pnts, /* input - number of x & y points */
float *x, /* input - x values */
float *y, /* input - y values */
float *sec_deriv) /* output - 2nd derivatives of intervals */
{
unsigned int i; /* index */
float ftemp; /* temporary float */
float *main_diag; /* matrix main diagonal array */
float *diag; /* matrix diagonal array */
float *right; /* equation right hand side values */
main_diag = malloc((num_pnts - 2) * sizeof(float));
diag = malloc((num_pnts - 1) * sizeof(float));
right = malloc((num_pnts - 2) * sizeof(float));
/* compute the matrix main and off-diagonal values */
/* even though the calling program is suppose to have guaranteed
that the x values are increasing, assert that neither of the
diagonal differences are zero to avoid a divide by zero
condition */
for (i = 1; i < num_pnts - 1; i++) {
main_diag[i-1] = 2.0 * (x[i+1] - x[i-1]);
assert(main_diag[i-1] > 0);
}
for (i = 0; i < num_pnts - 1; i++) {
diag[i-1] = x[i+1] - x[i];
assert(diag[i-1] > 0);
}
/* compute right hand side of equation */
for (i = 1; i < num_pnts - 1; i++)
right[i-1] = 6.0 * ((y[i+1] - y[i]) / diag[i-1] -
(y[i] - y[i-1]) / diag[i-2]);
/* forward eliminate tridiagonal */
sec_deriv[0] = 0.0;
sec_deriv[num_pnts - 1] = 0.0;
for (i = 1; i < num_pnts - 2; i++) {
ftemp = diag[i] / main_diag[i];
right[i] -= (right[i-1] * ftemp);
main_diag[i] -= (diag[i-1] * ftemp);
}
/* backward substitution to solve for 2nd derivative at each knot */
for (i = num_pnts - 2; i > 0; i--)
sec_deriv[i] = (right[i-1] - diag[i-1] * sec_deriv[i+1]) / main_diag[i-1];
free(main_diag);
free(diag);
free(right);
return;
}
/*
Solve for roots x1 and x2 of a quadratic equation of the form
a * (x * x) + b * x + c = 0
using
d = -0.5 * [b + sgn(b) * sqrt(b*b - 4ac)]
where
x1 = d / a and
x2 = c / d.
This root-finder yields particulary accurate results when a
and/or c are small values.
*/
BOOL FindQuadRoots (
float a, /* input - coefficient a of quadratic equation */
float b, /* input - coefficient b of quadratic equation */
float c, /* input - coefficient c of quadratic equation */
float *x1, /* output - first root computed */
float *x2) /* output - second root computed */
{
float d; /* intermediate value */
BOOL root_stat; /* status of root computations */
d = b * b - 4 * a * c;
if (d < 0)
return FAILURE;
else {
d = (float)sqrt((double)d);
/* make the result of sqrt the sign of b */
if (b < 0 )
d = -d;
d = -0.5 * (b + d);
/* solve for the roots of the quadratic */
/* if both root computations will yield divide by 0 ... forget it! */
if ( (a == 0) && (d == 0) )
return FAILURE;
root_stat = SUCCESS;
/* compute first root */
if (a == 0)
root_stat = FAILURE1;
else
*x1 = d / a;
/* compute second root */
if (d == 0)
root_stat = FAILURE2;
else
*x2 = c / d;
return root_stat;
}
}
/*
Given an abscissa (x) location, computes the corresponding cubic
spline ordinate (y) value.
*/
void ComputeY (
unsigned int i, /* input - array index */
float *x, /* input - array of x values */
float x_value, /* input - x value at which to solve for y */
float *y, /* input - array of y values */
float *sec_deriv, /* input - array of second derivatives of each data interval */
float *y_value) /* output - y extreme value at x */
{
float A, B, C, D; /* cubic spline coefficients */
float ftemp; /* temporary float value */
/* compute the standard cubic spline coefficients */
A = (x[i + 1] - x_value) / (x[i + 1] - x[i]);
B = 1 - A;
ftemp = (float) pow((double)(x[i + 1] - x[i]), 2.0) / 6.0;
C = (A * A * A - A) * ftemp;
D = (B * B * B - B) * ftemp;
/* compute the ordinate value at the abscissa location */
*y_value = A * y[i] + B * y[i + 1] + C * sec_deriv[i] + D * sec_deriv[i + 1];
return;
}
/********************* End of cubic_extrema.c **********************/
/**************************** cubic.h ******************************/
/*
File: cubic.h
Purpose: #defines, typedefs, and prototypes for cubic_io.c and cubic_extrema.c
Author: Mike Courtney 1996
*/
/* status defines */
#define SUCCESS 1
#define FAILURE 0
#define FAILURE1 -1
#define FAILURE2 -2
/* flag defines */
#define TRUE 1
#define FALSE 0
/* typedefs */
typedef int BOOL;
typedef unsigned char BYTE;
/* structures */
struct point {
struct point *next;
float x;
float y;
};
/* prototypes */
void ComputeSecDerivs (unsigned int, float *, float *, float *);
BOOL FindCubicExtrema (unsigned int, float *,float *, struct point *);
BOOL FindQuadRoots (float, float, float, float *, float *);
void ComputeY (unsigned int, float *, float, float *, float *, float *);
/********************* End of cubic.h ***************************/