/*
* 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.
*/
/*
* featext.c - routines to perform feature extraction on a set of curves
*/
#include <config.h>
#include <cmath.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <Xm/Xm.h>
#include <Xm/DialogS.h>
#include <Xm/Form.h>
#include <Xm/Label.h>
#include <Xm/PushB.h>
#include <Xm/RowColumn.h>
#include <Xm/ToggleB.h>
#include <Xm/Text.h>
#include "globals.h"
#include "utils.h"
#include "graphs.h"
#include "protos.h"
#include "motifinc.h"
typedef struct _Featext_ui {
Widget top;
ListStructure *tograph;
OptionStructure *feature_item;
Widget *xval_item;
ListStructure *absic_graph;
SetChoiceItem absic_set;
Widget legload_rc;
} Featext_ui;
static Featext_ui feui;
static Widget but2[3];
void do_fext_proc( Widget, XtPointer, XtPointer );
double linear_interp( double, double, double, double, double );
int dbl_comp( const void *, const void * );
int getmedian( int grno, int setno, int sorton, double *median );
int get_zero_crossing( int setl, double *xv, double *yv, double *crossing );
int get_rise_time( int setl, double *xv, double *yv, double min, double max,
double *width );
int get_fall_time( int setl, double *xv, double *yv, double min, double max,
double *width );
int mute_linear_regression(int n, double *x, double *y, double *slope,
double *intercept);
int get_half_max_width(int len, double *x, double *y, double *width);
int get_barycenter( int n, double *x, double *y, double *barycenter );
void fext_routine( int gto, int feature, int abs_src, int abs_set, int abs_graph );
void get_max_pos( double *a, double *b, int n, double max, double *d );
void do_fext_toggle (Widget w, XtPointer client_data, XtPointer call_data)
{
int value = (int) client_data;
if (value == 2 || value == 3) {
SetSensitive(feui.legload_rc, True);
} else {
SetSensitive(feui.legload_rc, False);
}
}
void create_featext_frame(void *data)
{
set_wait_cursor();
if (feui.top == NULL) {
Widget dialog;
char *label2[2];
int i;
OptionItem feature_option_items[] = {
{ 0, "Y minimum" },
{ 1, "Y maximum" },
{ 2, "Y average" },
{ 3, "Y std. dev." },
{ 4, "Y median" },
{ 5, "X minimum" },
{ 6, "X maximum" },
{ 7, "X average" },
{ 8, "X std. dev." },
{ 9, "X median" },
{10, "Frequency" },
{11, "Period" },
{12, "Zero crossing" },
{13, "Rise time" },
{14, "Fall time" },
{15, "Slope" },
{16, "Y intercept" },
{17, "Set length" },
{18, "Half maximal width"},
{19, "Barycenter X" },
{20, "Barycenter Y" },
{21, "X(Y max)" },
{22, "Y(X max)" },
{23, "Integral" }
};
label2[0] = "Accept";
label2[1] = "Close";
feui.top = XmCreateDialogShell(app_shell, "Feature Extraction", NULL, 0);
handle_close(feui.top);
dialog = XmCreateRowColumn(feui.top, "dialog_rc", NULL, 0);
feui.tograph = CreateGraphChoice(dialog, "Results to graph:", LIST_TYPE_SINGLE);
feui.feature_item =
CreateOptionChoice(dialog, "Feature:", 3, 24, feature_option_items);
feui.xval_item = CreatePanelChoice(dialog,
"X values from:", 5,
"Index", "Legends", "X from Set", "Y from set",
NULL);
for (i = 0; i < 4; i++) {
XtAddCallback(feui.xval_item[2 + i], XmNactivateCallback,
(XtCallbackProc) do_fext_toggle, (XtPointer) i);
}
CreateSeparator(dialog);
feui.legload_rc= XmCreateRowColumn(dialog, "fext_legload_rc", NULL, 0);
feui.absic_graph = CreateGraphChoice(feui.legload_rc,
"Abscissa from graph:", LIST_TYPE_SINGLE);
feui.absic_set = CreateSetSelector(feui.legload_rc, "set:",
SET_SELECT_ACTIVE,
FILTER_SELECT_NONE,
0,
SELECTION_TYPE_SINGLE);
update_save_set_list( feui.absic_set, 0 );
ManageChild(feui.legload_rc);
SetSensitive(feui.legload_rc, False);
CreateSeparator(dialog);
CreateCommandButtons(dialog, 2, but2, label2);
XtAddCallback(but2[0], XmNactivateCallback, (XtCallbackProc)
do_fext_proc,(XtPointer) & feui);
XtAddCallback(but2[1], XmNactivateCallback,
(XtCallbackProc)destroy_dialog,(XtPointer)feui.top);
ManageChild(dialog);
}
RaiseWindow(feui.top);
unset_wait_cursor();
}
void do_fext_proc( Widget w, XtPointer client_data, XtPointer call_data )
{
int gto, feature, abs_graph = -1, abs_set = -1, abs_src;
Featext_ui *ui = (Featext_ui *) client_data;
feature = GetOptionChoice(ui->feature_item);
GetSingleListChoice(ui->tograph, >o);
if( gto == -1 )
gto = get_cg();
abs_src = (int) GetChoice(ui->xval_item);
if( abs_src ==2 || abs_src==3 ) {
abs_set = GetSelectedSet(ui->absic_set);
GetSingleListChoice(ui->absic_graph, &abs_graph);
}
fext_routine( gto, feature, abs_src, abs_set, abs_graph );
update_set_lists(gto);
xdrawgraph();
}
void fext_routine( int gto, int feature, int abs_src, int abs_set, int abs_graph )
{
int i, cs, ns, fts, ncurves, extract_err;
double datum, dummy, *absy;
double y1, y2;
int iy1, iy2;
char tbuf[1024];
float *abscissa;
double xmin, xmax, ymin, ymax;
int cg = get_cg();
double *x;
abscissa = xmalloc(number_of_sets(cg)*SIZEOF_FLOAT);
if( !is_graph_active( gto ) ){
errwin("Graph for results must be active");
return;
}
if( (ns=nextset( gto ) )== -1 ) {
errwin("Choose a new graph or kill sets!");
return;
}
ncurves = nactive(get_cg());
switch( abs_src ) {
case 0: /* use index */
for( i=0; i<ncurves; i++ )
abscissa[i] = i;
break;
case 1: /* use legend label */
cs = 0;
for( i=0; i<ncurves; i++ ){
while( !is_set_active( get_cg(), cs ) )
cs++;
if(!sscanf( get_legend_string(get_cg(), cs), "%f", &abscissa[i]))
break;
cs++;
}
if( i != ncurves ) {
errwin("Bad legend label");
return;
}
break;
case 2: /* use X from set */
if( !is_set_active( abs_graph, abs_set ) ){
errwin("Abscissa set not active");
return;
}
if( getsetlength( abs_graph, abs_set ) < ncurves ) {
errwin("Not enough points in set");
return;
}
absy = getx( abs_graph, abs_set );
for( i=0; i<ncurves; i++ )
abscissa[i] = absy[i];
break;
case 3: /* use Y from set */
if( !is_set_active( abs_graph, abs_set ) ){
errwin("Abscissa set not active");
return;
}
if( getsetlength( abs_graph, abs_set ) < ncurves ) {
errwin("Not enough points in set");
return;
}
absy = gety( abs_graph, abs_set );
for( i=0; i<ncurves; i++ )
abscissa[i] = absy[i];
break;
}
cs = 0;
tbuf[0] = '\0';
for( i=0; i<ncurves; i++ ) {
while( !is_set_active( get_cg(), cs ) )
cs++;
extract_err = 0;
getsetminmax(get_cg(), cs, &xmin, &xmax, &ymin, &ymax);
switch( feature ) {
case 0: /* Y minimum */
datum = ymin;
break;
case 1: /* Y maximum */
datum = ymax;
break;
case 2: /* Y mean */
stasum(gety(cg, cs), getsetlength(cg, cs), &datum, &dummy);
break;
case 3: /* Y std dev */
stasum(gety(cg, cs), getsetlength(cg, cs), &dummy, &datum);
break;
case 4: /* Y median */
getmedian( cg, cs, DATA_Y, &datum );
break;
case 5: /* X minimum */
datum = xmin;
break;
case 6: /* X maximum */
datum = xmax;
break;
case 7: /* X mean */
stasum(getx(cg, cs), getsetlength(cg, cs), &datum, &dummy);
break;
case 8: /* X std dev */
stasum(getx(cg, cs), getsetlength(cg, cs), &dummy, &datum);
break;
case 9: /* X median */
getmedian( cg, cs, DATA_X, &datum );
break;
case 10: /* frequency and period */
case 11:
if ( ilog2(getsetlength(cg, cs)) <= 0) /* only DFT */
do_fourier(cg, cs, 0, 0, 1, 0, 0, 0);
else /* FFT */
do_fourier(cg, cs, 1, 0, 1, 0, 0, 0);
sprintf( tbuf, "FT of set %d", cs );
fts = 0;
while( strcmp( tbuf, getcomment(cg, fts)+1 ) )
fts++;
minmax(gety(cg, fts), getsetlength(cg, fts),&y1,&y2,&iy1,&iy2);
x = getx(cg, fts);
if( feature == 10 )
datum = x[iy2];
else
datum = 1./x[iy2];
killset( cg, fts ); /* get rid of Fourier set */
break;
case 12: /* first zero crossing */
if( get_zero_crossing( getsetlength( cg, cs ),
getx( cg, cs ),gety( cg, cs ), &datum ) ){
sprintf( tbuf+strlen(tbuf),
"Unable to find zero crossing of set %d\n", cs );
errwin( tbuf );
extract_err = 1;
}
break;
case 13: /* rise time */
getsetminmax(cg, cs, &xmin, &xmax, &ymin, &ymax);
if( get_rise_time( getsetlength(cg,cs), getx(cg,cs),
gety(cg,cs), ymin, ymax, &datum ) ){
sprintf( tbuf+strlen(tbuf),
"Unable to find rise time of set %d\n", cs );
errwin( tbuf );
extract_err = 1;
}
break;
case 14: /* fall time */
getsetminmax(cg, cs, &xmin, &xmax, &ymin, &ymax);
if( get_fall_time( getsetlength(cg,cs), getx(cg,cs),
gety(cg,cs), ymin, ymax, &datum ) ){
sprintf( tbuf+strlen(tbuf),
"Unable to find fall time of set %d\n", cs );
extract_err = 1;
errwin( tbuf );
}
break;
case 15: /* slope */
if( mute_linear_regression( getsetlength( cg, cs ),
getx( cg, cs ),gety( cg, cs ), &datum, &dummy ) ) {
sprintf( tbuf+strlen(tbuf),
"Unable to find slope of set %d\n", cs );
errwin( tbuf );
extract_err = 1;
}
break;
case 16: /* Y intercept */
if( mute_linear_regression( getsetlength( cg, cs ),
getx( cg, cs ), gety( cg, cs ), &dummy, &datum ) ) {
sprintf( tbuf+strlen(tbuf),
"Unable to find y-intercept of set %d\n", cs );
errwin( tbuf );
extract_err = 1;
}
break;
case 17: /* set length */
datum = getsetlength( cg, cs );
break;
case 18: /* half maximal widths */
if (get_half_max_width(getsetlength(cg, cs), getx(cg,cs),
gety(cg,cs), &datum) != RETURN_SUCCESS) {
sprintf( tbuf+strlen(tbuf),
"Unable to find half maximal width of set %d\n", cs );
extract_err = 1;
errwin( tbuf );
}
break;
case 19: /* Barycenter X */
get_barycenter( getsetlength( cg, cs ), gety(cg,cs),
getx(cg,cs), &datum );
break;
case 20: /* Barycenter Y */
get_barycenter( getsetlength( cg, cs ), getx(cg,cs),
gety(cg,cs), &datum );
break;
case 21: /* X of Maximum Y */
getsetminmax(cg, cs, &xmin, &xmax, &ymin, &ymax);
get_max_pos( gety(cg, cs), getx( cg, cs ),
getsetlength( cg, cs ), ymax, &datum );
break;
case 22: /* Y of Maximum X */
getsetminmax(cg, cs, &xmin, &xmax, &ymin, &ymax);
get_max_pos( getx(cg, cs), gety( cg, cs ),
getsetlength( cg, cs ), xmax, &datum );
break;
case 23: /* cumulative sum */
datum = do_int(cg, cs, 1);
break;
}
if( !extract_err )
add_point(gto, ns, abscissa[i], datum);
cs++;
}
/* set comment */
switch( feature ) {
case 0: /* Y minimum */
sprintf(tbuf,"Y minima of graph %d",cg);
break;
case 1: /* Y maximum */
sprintf(tbuf,"Y maxima of graph %d",cg);
break;
case 2: /* Y mean */
sprintf(tbuf,"Y means of graph %d",cg);
break;
case 3: /* Y std dev */
sprintf(tbuf,"Y std. dev.'s of graph %d",cg);
break;
case 4: /* Y median */
sprintf(tbuf,"Y medians of graph %d",cg);
break;
case 5: /* X minimum */
sprintf(tbuf,"X minima of graph %d",cg);
break;
case 6: /* X maximum */
sprintf(tbuf,"X maxima of graph %d",cg);
break;
case 7: /* X mean */
sprintf(tbuf,"X means of graph %d",cg);
break;
case 8: /* X std dev */
sprintf(tbuf,"X std. dev.'s of graph %d",cg);
break;
case 9: /* X median */
sprintf(tbuf,"X medians of graph %d",cg);
break;
case 10: /* frequency and period */
sprintf(tbuf,"frequencies of graph %d",cg);
break;
case 11:
sprintf(tbuf,"periods of graph %d",cg);
break;
case 12: /* first zero crossing */
sprintf(tbuf,"zero crossings of graph %d",cg);
break;
case 13: /* rise time */
sprintf(tbuf,"rise times of graph %d",cg);
break;
case 14: /* fall time */
sprintf(tbuf,"fall times of graph %d",cg);
break;
case 15: /* slopes */
sprintf(tbuf,"slopes of graph %d",cg);
break;
case 16: /* Y intercepts */
sprintf(tbuf,"Y intercepts of graph %d",cg);
break;
case 17: /* set lengths */
sprintf(tbuf,"set lengths of graph %d",cg);
break;
case 18: /* 1/2 maximal widths */
sprintf(tbuf,"half maximal widths of graph %d",cg);
break;
case 19: /* barycenter X */
sprintf(tbuf,"X barycenters of graph %d",cg);
break;
case 20: /* barycenter Y */
sprintf(tbuf,"Y barycenters of graph %d",cg);
break;
case 21: /* X of maximum Y */
sprintf(tbuf,"X positions of maximum Y's of graph %d",cg);
break;
case 22: /* Y of maximum X */
sprintf(tbuf,"Y positions of maximum X's of graph %d",cg);
break;
case 23: /* integral */
sprintf(tbuf,"integrals of sets of graph %d",cg);
break;
}
set_set_hidden(gto, ns, FALSE);
setcomment( gto, ns, tbuf );
xfree( abscissa );
}
/* linear regression without posting results to log */
int mute_linear_regression(int n, double *x, double *y, double *slope,
double *intercept)
{
double xbar, ybar; /* sample means */
double sdx, sdy; /* sample standard deviations */
double SXX, SYY, SXY; /* sums of squares */
int i;
if (n <= 3) {
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;
}
*slope = SXY / SXX;
*intercept = ybar - *slope * xbar;
return 0;
}
/*
* assume graph starts off at ymin and rises to ymax
* Determine time to go from 10% to 90% of rise
*/
int get_rise_time( int setl, double *xv, double *yv, double min, double max,
double *width )
{
int x10=0, x90;
double amp10, amp90;
amp10 = min + (max-min)*0.1;
amp90 = min + (max-min)*0.9;
while( x10<setl && yv[x10]<amp10 )
x10++;
if( x10==setl || x10==0)
return 1;
x90 = x10+1;
while( x90<setl && yv[x90]<amp90 )
x90++;
*width = linear_interp( yv[x90-1], xv[x90-1], yv[x90], xv[x90], amp90 ) -
linear_interp( yv[x10-1], xv[x10-1], yv[x10], xv[x10], amp10 );
return 0;
}
/* assume graph starts off at ymax and drops to ymin
Determine time to go from 90% to 10% of drop */
int get_fall_time( int setl, double *xv, double *yv, double min, double max,
double *width )
{
int x10, x90=0;
double amp10, amp90;
amp10 = min + (max-min)*0.1;
amp90 = min + (max-min)*0.9;
while( x90<setl && yv[x90]>amp90 )
x90++;
if( x90==setl || x90==0)
return 1;
x10= x90+1;
while( x10<setl && yv[x10]>amp10 )
x10++;
if( x10==setl )
return 1;
*width = linear_interp( yv[x10-1], xv[x10-1], yv[x10], xv[x10], amp10 )-
linear_interp( yv[x90-1], xv[x90-1], yv[x90], xv[x90], amp90 );
return 0;
}
int get_zero_crossing( int setl, double *xv, double *yv, double *crossing )
{
int i=0;
while( i<setl && yv[i] != 0. && yv[i]*yv[i+1]>0. )
i++;
if( i==setl )
return 1;
if( yv[i] == 0 )
*crossing = xv[i];
else
*crossing = linear_interp( yv[i], xv[i], yv[i+1], xv[i+1], 0 );
return 0;
}
/* Get FWHM of the highest peak */
int get_half_max_width(int len, double *x, double *y, double *width)
{
int i, imin, imax;
double ymin, ymax, yavg;
double x_u, x_d;
minmax(y, len, &ymin, &ymax, &imin, &imax);
yavg = (ymin + ymax)/2.0;
i = imax;
while (i >= 0 && y[i] > yavg) {
i--;
}
if (i < 0) {
return RETURN_FAILURE;
}
x_d = linear_interp(y[i], x[i], y[i + 1], x[i + 1], yavg);
i = imax;
while (i < len && y[i] > yavg) {
i++;
}
if (i == len) {
return RETURN_FAILURE;
}
x_u = linear_interp(y[i - 1], x[i - 1], y[i], x[i], yavg);
*width = fabs(x_u - x_d);
return RETURN_SUCCESS;
}
/* linear interpolate between two points, return a y value given an x */
double linear_interp( double x1, double y1, double x2, double y2, double x )
{
return y1 + ( x-x1 )*(y2-y1)/(x2-x1);
}
/* get the median of the X or Y portion of a set */
int getmedian( int grno, int setno, int sorton, double *median )
{
int setlen;
double *setdata;
setlen = getsetlength( get_cg(), setno );
setdata = (double *)xmalloc( setlen*sizeof(double) );
if( sorton == DATA_X )
memcpy( setdata, getx( grno, setno ), setlen*sizeof(double) );
else
memcpy( setdata, gety( grno, setno ), setlen*sizeof(double) );
qsort( setdata, setlen, sizeof(double), dbl_comp );
if( setlen%2 ) /* odd set */
*median = setdata[(setlen+1)/2-1];
else
*median = ( setdata[setlen/2-1] + setdata[setlen/2] )/2.;
xfree( setdata );
return 0;
}
int dbl_comp( const void *a, const void *b )
{
if( *(double *)a < *(double *)b )
return -1;
else if( *(double *)a > *(double *)b )
return 1;
else
return 0;
}
int get_barycenter( int n, double *x, double *y, double *barycenter )
{
double wsum=0;
*barycenter = 0;
for( n--; n>=0; n-- ) {
wsum += x[n]*y[n];
*barycenter += x[n];
}
*barycenter = wsum/(*barycenter);
return 0;
}
/* given maximum of set a, find the corresponding entry in set b */
void get_max_pos( double *a, double *b, int n, double max, double *d )
{
int i=-1;
while( ++i<n && a[i] != max );
if( i==n )
return;
else
*d = b[i];
}
syntax highlighted by Code2HTML, v. 0.9.1