/*****************************************************************************
Major portions of this software are copyrighted by the Medical College
of Wisconsin, 1994-2000, and are released under the Gnu General Public
License, Version 2. See the file README.Copyright for details.
******************************************************************************/
#include "afni.h"
#ifndef ALLOW_PLUGINS
# error "Plugins not properly set up -- see machdep.h"
#endif
/***********************************************************************
Plugin to provide a least squares fitting 1D function for graphs
************************************************************************/
/*------------- string to 'help' the user -------------*/
static char helpstring[] =
" Purpose: Control the 'LSqFit' and 'LSqDtr 1D functions.\n"
"\n"
" Parameters: Baseline = 'Constant' or 'Linear'\n"
" Is the baseline 'a' or 'a+b*t'?\n"
" Ignore = Number of points to ignore at\n"
" start of each timeseries.\n"
" \n"
" Sinusoids: Period = Fundamental period to use.\n"
" Harmonics = Number of overtones to use.\n"
" \n"
" Timeseries: File = Input timeseries file to use.\n"
;
static char plehstring[] =
" Purpose: Generate a timeseries and store in AFNI list\n"
"\n"
" Parameters: Delta = time step between points\n"
" Length = number of points\n"
"\n"
" Output: Label = String to label timeseries by\n"
" in AFNI choosers\n"
"\n"
" Polynomial: Order = Maximum power to include\n"
" (Chebyshev polynomials are used)\n"
"\n"
" Sinusoid: Period = Fundamental period to use.\n"
" Harmonics = Number of overtones to use.\n"
;
/*------- Strings for baseline control ------*/
#define NBASE 3
static char * baseline_strings[NBASE] = { "Constant" , "Linear" , "Quadratic" } ;
/*--------------- prototypes for internal routines ---------------*/
char * LSQ_main( PLUGIN_interface * ) ; /* the entry point */
void LSQ_fitter( int nt, double to, double dt, float * vec, char ** label ) ;
void LSQ_detrend( int nt, double to, double dt, float * vec, char ** label ) ;
void LSQ_worker( int nt, double dt, float * vec, int dofit, char ** label ) ;
PLUGIN_interface * TSGEN_init(void) ;
char * TSGEN_main( PLUGIN_interface * ) ;
PLUGIN_interface * EXP0D_init(void) ;
char * EXP0D_main( PLUGIN_interface * ) ;
void EXP0D_worker( int num , float * vec ) ;
#ifdef ALLOW_LOMO
PLUGIN_interface * LOMOR_init(void) ;
char * LOMOR_main( PLUGIN_interface * ) ;
void LOMOR_fitter() ;
#endif
/*---------------- global data -------------------*/
static PLUGIN_interface * global_plint = NULL ;
#define NRMAX_SIN 2
#define NRMAX_TS 2
#define HARM_MAX 22
static int polort=1 , ignore=3 , nrsin=0 , nrts=0 , initialize=1 ;
static float sinper[NRMAX_SIN] ;
static int sinharm[NRMAX_SIN] ;
static MRI_IMAGE * tsim[NRMAX_TS] ;
/***********************************************************************
Set up the interface to the user:
1) Create a new interface using "PLUTO_new_interface";
2) For each line of inputs, create the line with "PLUTO_add_option"
(this line of inputs can be optional or mandatory);
3) For each item on the line, create the item with
"PLUTO_add_dataset" for a dataset chooser,
"PLUTO_add_string" for a string chooser,
"PLUTO_add_number" for a number chooser,
"PLUTO_add_timeseries" for a timeseries chooser.
************************************************************************/
DEFINE_PLUGIN_PROTOTYPE
PLUGIN_interface * PLUGIN_init( int ncall )
{
int ii ;
PLUGIN_interface * plint ; /* will be the output of this routine */
if( ncall > 3 ) return NULL ; /* generate interfaces for ncall 0-3 */
if( ncall == 1 ) return TSGEN_init() ; /* interface # 1 */
if( ncall == 2 ) return EXP0D_init() ; /* interface # 2 */
#ifdef ALLOW_LOMO
if( ncall == 3 ) return LOMOR_init() ; /* interface # 3 */
#else
if( ncall == 3 ) return NULL ;
#endif
/***** otherwise, do interface # 0 *****/
/*---------------- set titles and call point ----------------*/
plint = PLUTO_new_interface( "LSqFit & Dtr" ,
"Control LSqFit and LSqDtr Functions" ,
helpstring ,
PLUGIN_CALL_VIA_MENU , LSQ_main ) ;
global_plint = plint ; /* make global copy */
PLUTO_set_sequence( plint , "A:funcs:fitting" ) ;
PLUTO_add_hint( plint , "Control LSqFit and LSqDtr Functions" ) ;
PLUTO_set_runlabels( plint , "Set+Keep" , "Set+Close" ) ; /* 04 Nov 2003 */
/*----- Parameters -----*/
PLUTO_add_option( plint , "Parameters" , "Parameters" , TRUE ) ;
PLUTO_add_string( plint , "Baseline" , NBASE , baseline_strings , 1 ) ;
PLUTO_add_number( plint , "Ignore" , 0,20,0,3 , FALSE ) ;
/*----- Sinusoid -----*/
for( ii=0 ; ii < NRMAX_SIN ; ii++ ){
PLUTO_add_option( plint , "Sinusoid" , "Sinusoid" , FALSE ) ;
PLUTO_add_number( plint , "Period" , 0,99999,0,20, TRUE ) ;
PLUTO_add_number( plint , "Harmonics" , 1,HARM_MAX,0,1 , FALSE ) ;
}
/*----- Timeseries -----*/
for( ii=0 ; ii < NRMAX_TS ; ii++ ){
PLUTO_add_option( plint , "Timeseries" , "Timeseries" , FALSE ) ;
PLUTO_add_timeseries( plint , "File" ) ;
}
/*--------- done with interface setup ---------*/
PLUTO_register_1D_funcstr( "LSqFit" , LSQ_fitter ) ;
PLUTO_register_1D_funcstr( "LSqDtr" , LSQ_detrend ) ;
return plint ;
}
/***************************************************************************
Main routine for this plugin (will be called from AFNI).
If the return string is not NULL, some error transpired, and
AFNI will popup the return string in a message box.
****************************************************************************/
char * LSQ_main( PLUGIN_interface * plint )
{
char * str ;
int ii ;
float * tsar ;
/*--------- go to first input line ---------*/
PLUTO_next_option(plint) ;
str = PLUTO_get_string(plint) ;
polort = PLUTO_string_index( str , NBASE , baseline_strings ) ;
ignore = PLUTO_get_number(plint) ;
/*------ loop over remaining options, check their tags, process them -----*/
nrsin = nrts = 0 ;
do {
str = PLUTO_get_optiontag(plint) ; if( str == NULL ) break ;
if( strcmp(str,"Sinusoid") == 0 ){
sinper[nrsin] = PLUTO_get_number(plint) ;
sinharm[nrsin] = PLUTO_get_number(plint) - 1.0 ;
if( sinper[nrsin] <= 0.0 )
return "************************\n"
"Illegal Sinusoid Period!\n"
"************************" ;
nrsin++ ;
} else if( strcmp(str,"Timeseries") == 0 ){
tsim[nrts] = PLUTO_get_timeseries(plint) ;
if( tsim[nrts] == NULL || tsim[nrts]->nx < 3 || tsim[nrts]->kind != MRI_float )
return "*************************\n"
"Illegal Timeseries Input!\n"
"*************************" ;
tsar = MRI_FLOAT_PTR(tsim[nrts]) ;
for( ii=ignore ; ii < tsim[nrts]->nx && tsar[ii] >= WAY_BIG ; ii++ ) ; /* nada */
ignore = ii ;
nrts++ ;
} else {
return "************************\n"
"Illegal optiontag found!\n"
"************************" ;
}
} while(1) ;
/*--- nothing left to do until data arrives ---*/
initialize = 1 ; /* force re-initialization */
/*** compute how many ref functions are ordered ***/
{ int nref , ks ;
char str[64] ;
nref = (polort+1) + nrts ;
for( ks=0 ; ks < nrsin ; ks++ ) nref += 2*(sinharm[ks]+1) ;
sprintf(str," \nNumber of fit parameters = %d\n",nref) ;
PLUTO_popup_transient( plint , str ) ;
}
return NULL ;
}
/*---------------------------------------------------------------*/
/** 22 Apr 1997: added label that will go to graphs **/
void LSQ_fitter( int nt , double to , double dt , float * vec , char ** label )
{
LSQ_worker( nt , dt , vec , TRUE , label ) ;
return ;
}
void LSQ_detrend( int nt , double to , double dt , float * vec , char ** label )
{
LSQ_worker( nt , dt , vec , FALSE , label ) ;
return ;
}
static char lbuf[4096] ; /* 22 Apr 1997: will hold label for graphs */
static char sbuf[256] ;
void LSQ_worker( int nt , double dt , float * vec , int dofit , char ** label )
{
int nlen , nref ;
static int nlen_old = -666 , nref_old = -666 ;
static double dt_old = -666.666 ;
static float ** ref = NULL ;
static double * dch = NULL ;
int ir , ii , ks,jh , j;
float fac , tm , val ;
float * fit , * tsar ;
/*** compute how many ref functions are ordered ***/
nref = (polort+1) + nrts ;
for( ks=0 ; ks < nrsin ; ks++ ) nref += 2*(sinharm[ks]+1) ;
/*** do nothing if not enough data to fit ***/
nlen = nt - ignore ;
if( nlen <= nref ) return ; /* do nothing if not enough data to fit */
/** if data vectors are new length,
or have a new number of reference vectors,
or have a new time step and need sinusoids,
or the initialize flag is set,
then reinitialize reference vectors and Choleski factor **/
if( nlen != nlen_old || nref != nref_old ||
initialize || (dt != dt_old && nrsin > 0) ){
/* free old storage */
if( ref != NULL ){
for( ir=0 ; ir < nref_old ; ir++ ) if( ref[ir] != NULL ) free(ref[ir]) ;
free(ref) ;
}
if( dch != NULL ) free(dch) ;
/* make space for ref vectors */
ref = (float **) malloc( sizeof(float *) * nref ) ;
if( ref == NULL ){fprintf(stderr,"\nmalloc error in plug_lsqfit\n\a");
return;
/* EXIT(1); */
}
for( ir=0 ; ir < nref ; ir++ ){
ref[ir] = (float *) malloc( sizeof(float) * nlen ) ;
if( ref[ir] == NULL ){
fprintf(stderr,"\nmalloc error in plug_lsqfit\n\a");
for(j=0;j<=ir;j++)
free(ref[j]);
free(ref);
ref = NULL;
return;
/* EXIT(1); */
}
}
nlen_old = nlen ;
nref_old = nref ;
dt_old = dt ;
/**** fill ref vectors ****/
/* r(t) = 1 */
for( ii=0 ; ii < nlen ; ii++ ) ref[0][ii] = 1.0 ;
ir = 1 ;
if( polort > 0 ){
/* r(t) = t - tmid */
tm = 0.5 * (nlen-1.0) ; fac = 2.0 / nlen ;
for( ii=0 ; ii < nlen ; ii++ ) ref[1][ii] = (ii-tm)*fac ;
ir = 2 ;
/* r(t) = (t-tmid)**ir */
for( ; ir <= polort ; ir++ )
for( ii=0 ; ii < nlen ; ii++ )
ref[ir][ii] = pow( (ii-tm)*fac , (double)ir ) ;
}
if( dt == 0.0 ) dt = 1.0 ;
/* r(t) = sinusoids */
for( ks=0 ; ks < nrsin ; ks++ ){
for( jh=0 ; jh <= sinharm[ks] ; jh++ ){
fac = (2.0*PI) * dt * (jh+1) / sinper[ks] ;
for( ii=0 ; ii < nlen ; ii++ ){
ref[ir] [ii] = cos( fac * ii ) ;
ref[ir+1][ii] = sin( fac * ii ) ;
}
ir += 2 ;
}
}
/* r(t) = timeseries files */
for( ks=0 ; ks < nrts ; ks++ ){
if( tsim[ks] == NULL || tsim[ks]->nx - ignore < nlen ){
initialize = 1 ;
fprintf(stderr,"Inadequate time series #%d in LSQ plugin\n\a",ks+1) ;
return ;
}
tsar = MRI_FLOAT_PTR(tsim[ks]) ;
for( ii=0 ; ii < nlen ; ii++ ) ref[ir][ii] = tsar[ii+ignore] ;
ir++ ;
}
/* Cholesky-ize */
dch = startup_lsqfit( nlen , NULL , nref , ref ) ;
if( dch == NULL ){
initialize = 1 ;
fprintf(stderr,"Choleski error in LSQ plugin\n\a") ;
return ;
}
initialize = 0 ;
}
/** find least squares fit coefficients **/
fit = delayed_lsqfit( nlen , vec+ignore , nref , ref , dch ) ;
for( ii=0 ; ii < nlen ; ii++ ){
val = 0.0 ;
for( ir=0 ; ir < nref ; ir++ ) val += fit[ir] * ref[ir][ii] ;
vec[ii+ignore] = (dofit) ? val : vec[ii+ignore] - val ;
}
/** 22 Apr 1997: create label if desired by AFNI **/
/** [This is in static storage, since AFNI will copy it] **/
if( label != NULL ){ /* assemble this 1 line at a time from sbuf */
lbuf[0] = '\0' ; /* make this a 0 length string to start */
/** for each reference, make a string into sbuf **/
ir = 0 ;
sprintf(sbuf,"Coef of 1 = %g\n" , fit[ir++] ) ;
strcat(lbuf,sbuf) ;
for( ; ir <= polort ; ){
sprintf(sbuf,"Coef of t**%d = %g\n" , ir,fit[ir++] ) ;
strcat(lbuf,sbuf) ;
}
for( ks=0 ; ks < nrsin ; ks++ ){
for( jh=0 ; jh <= sinharm[ks] ; jh++ ){
fac = sinper[ks] / (jh+1) ;
sprintf(sbuf,"Coef of cos(2*Pi*t/%g) = %g\n" , fac , fit[ir++] ) ;
strcat(lbuf,sbuf) ;
sprintf(sbuf,"Coef of sin(2*Pi*t/%g) = %g\n" , fac , fit[ir++] ) ;
strcat(lbuf,sbuf) ;
}
}
for( ks=0 ; ks < nrts ; ks++ ){
sprintf(sbuf,"Coef of %s = %g\n" , tsim[ks]->name , fit[ir++] ) ;
strcat(lbuf,sbuf) ;
}
*label = lbuf ; /* send address of lbuf back in what label points to */
}
free(fit) ;
return ;
}
/*********************************************************************************/
PLUGIN_interface * TSGEN_init(void)
{
PLUGIN_interface * plint ;
/*---------------- set titles and call point ----------------*/
plint = PLUTO_new_interface( "TS Generate" ,
"Generate a Timeseries" ,
plehstring ,
PLUGIN_CALL_VIA_MENU , TSGEN_main ) ;
PLUTO_add_hint( plint , "Generate a 1D Timeseries" ) ;
/*----- Parameters -----*/
PLUTO_add_option( plint , "Parameters" , "Parameters" , TRUE ) ;
PLUTO_add_number( plint , "Delta" , 0,99999,1, 0 , TRUE ) ;
PLUTO_add_number( plint , "Length" , 3,9999,0,3 , TRUE ) ;
/*----- Output -----*/
PLUTO_add_option( plint , "Output" , "Output" , TRUE ) ;
PLUTO_add_string( plint , "Label" , 0,NULL , 19 ) ;
/*----- Polynomial -----*/
PLUTO_add_option( plint , "Polynomial" , "Polynomial" , FALSE ) ;
PLUTO_add_number( plint , "Order" , 2,20,0,2 , FALSE ) ;
/*----- Sinusoid -----*/
PLUTO_add_option( plint , "Sinusoid" , "Sinusoid" , FALSE ) ;
PLUTO_add_number( plint , "Period" , 0,99999,0,20, TRUE ) ;
PLUTO_add_number( plint , "Harmonics" , 1,HARM_MAX,0,1 , FALSE ) ;
return plint ;
}
char * TSGEN_main( PLUGIN_interface * plint )
{
char * label , * str ;
int ii , jj ;
float * tsar ;
MRI_IMAGE * tsim ;
float delta , period=0.0 ;
int nx , ny=0 , npol=0 , nharm=-1 ;
int pp ;
double fac , val ;
/*--------- go to first input line ---------*/
PLUTO_next_option(plint) ;
delta = PLUTO_get_number(plint) ;
if( delta <= 0.0 ) return "**********************\n"
"Illegal value of Delta\n"
"**********************" ;
nx = PLUTO_get_number(plint) ;
/*----- next input line -----*/
PLUTO_next_option(plint) ;
label = PLUTO_get_string(plint) ;
if( label == NULL || strlen(label) == 0 ) return "**********************\n"
"Illegal value of Label\n"
"**********************" ;
/*----- rest of input lines -----*/
do {
str = PLUTO_get_optiontag(plint) ; if( str == NULL ) break ;
if( strcmp(str,"Sinusoid") == 0 ){
period = PLUTO_get_number(plint) ;
nharm = PLUTO_get_number(plint) - 1.0 ;
if( period <= 0.0 ) return "***********************\n"
"Illegal Sinusoid Period\n"
"***********************" ;
} else if( strcmp(str,"Polynomial") == 0 ){
npol = PLUTO_get_number(plint) ;
} else {
return "***********************\n"
"Illegal optiontag found\n"
"***********************" ;
}
} while(1) ;
/********** Make the timeseries ***********/
ny = 0 ;
if( npol > 0 ) ny = npol-1 ;
if( period > 0.0 ) ny += 2*(nharm+1) ;
if( ny < 1 ) return "***********************\n"
"No timeseries specified\n"
"***********************" ;
tsim = mri_new( nx , ny , MRI_float ) ;
jj = 0 ;
fac = 1.99999 / (nx-1) ;
for( pp=2 ; pp <= npol ; pp++,jj++ ){
tsar = MRI_FLOAT_PTR(tsim) + (jj*nx) ;
for( ii=0 ; ii < nx ; ii++ ){
val = fac * ii - 0.999995 ;
tsar[ii] = cos( pp * acos(val) ) ;
}
}
for( pp=0 ; pp <= nharm ; pp++ , jj+=2 ){
fac = (2.0*PI) * delta * (pp+1) / period ;
tsar = MRI_FLOAT_PTR(tsim) + (jj*nx) ;
for( ii=0 ; ii < nx ; ii++ ){
tsar[ii] = cos( fac * ii ) ;
tsar[ii+nx] = sin( fac * ii ) ;
}
}
PLUTO_register_timeseries( label , tsim ) ;
mri_free(tsim) ;
return NULL ;
}
/***************************************************************************/
#include "parser.h"
static char fredstring[] =
" Purpose: Control the Expr 0D Transformation function\n"
"\n"
" Variable = letter used as input to expression\n"
" Expression = arithmetic expression to evaluate\n"
;
#define NALPHA 26
static char * vstring[NALPHA] = {
" A ", " B ", " C ", " D ", " E ",
" F ", " G ", " H ", " I ", " J ",
" K ", " L ", " M ", " N ", " O ",
" P ", " Q ", " R ", " S ", " T ",
" U ", " V ", " W ", " X ", " Y ", " Z " } ;
static int exp0d_var = 23 ;
static PARSER_code * exp0d_pc = NULL ;
static PLUGIN_interface *plint_EXP0D=NULL ;
static void EXP0D_func_init(void) /* 21 Jul 2003 */
{
PLUG_startup_plugin_CB( NULL , (XtPointer)plint_EXP0D , NULL ) ;
}
PLUGIN_interface * EXP0D_init(void)
{
PLUGIN_interface * plint ;
/*---------------- set titles and call point ----------------*/
plint = PLUTO_new_interface( "Expr 0D" ,
"Control the Expr 0D transformation" ,
fredstring ,
PLUGIN_CALL_VIA_MENU , EXP0D_main ) ;
PLUTO_add_option( plint , "Variable" , "Variable" , TRUE ) ;
PLUTO_add_string( plint , NULL , NALPHA,vstring , exp0d_var ) ;
PLUTO_add_option( plint , "Expression" , "Expression" , TRUE ) ;
PLUTO_add_string( plint , NULL , 0,NULL , 50 ) ;
PLUTO_register_0D_function( "Expr 0D" , EXP0D_worker ) ;
plint_EXP0D = plint ;
AFNI_register_nD_func_init( 0 , (generic_func *)EXP0D_func_init ) ; /* 21 Jul 2003 */
return plint ;
}
char * EXP0D_main( PLUGIN_interface * plint )
{
char * str ;
PLUTO_next_option(plint) ;
str = PLUTO_get_string(plint) ;
exp0d_var = PLUTO_string_index( str , NALPHA , vstring ) ;
if( exp0d_pc != NULL ){ free(exp0d_pc) ; exp0d_pc = NULL ; }
PLUTO_next_option(plint) ;
str = PLUTO_get_string(plint) ;
exp0d_pc = PARSER_generate_code( str ) ;
if( exp0d_pc == NULL ) return "*******************************\n"
"Error when compiling expression\n"
"*******************************" ;
return NULL ;
}
#define VSIZE 64
void EXP0D_worker( int num , float * vec )
{
int ii , jj , jbot,jtop ;
static int first = 1 ;
static double * atoz[26] ;
static double tvec[VSIZE] ;
if( num <= 0 || vec == NULL || exp0d_pc == NULL ) return ;
#if 0
fprintf(stderr,"Enter EXP0D_worker\n") ;
#endif
if( first ){
for( ii=0 ; ii < 26 ; ii++)
atoz[ii] = (double *) malloc(sizeof(double) * VSIZE ) ;
first = 0 ;
#if 0
fprintf(stderr,"Allocated atoz\n") ;
#endif
}
for( ii=0 ; ii < 26 ; ii++ )
for (jj=0; jj<VSIZE; jj++) atoz[ii][jj] = 0.0 ;
#if 0
fprintf(stderr,"Zeroed atoz\n") ;
#endif
for( ii=0 ; ii < num ; ii+=VSIZE ){
jbot = ii ;
jtop = MIN( ii + VSIZE , num ) ;
for( jj=jbot ; jj < jtop ; jj ++ ) atoz[exp0d_var][jj-ii] = vec[jj] ;
PARSER_evaluate_vector( exp0d_pc , atoz , jtop-jbot , tvec ) ;
for( jj=jbot ; jj < jtop ; jj ++ ) vec[jj] = tvec[jj-ii] ;
}
#if 0
fprintf(stderr,"Exit EXP0D_worker\n") ;
#endif
return ;
}
/*------------------------------------------------------------------*/
#ifdef ALLOW_LOMO
static int lomo_order = 3 ;
static int lomo_levels = 32 ;
int lomo_regress( int , int , int * , int * ) ;
PLUGIN_interface * LOMOR_init(void)
{
PLUGIN_interface * plint ;
/*---------------- set titles and call point ----------------*/
plint = PLUTO_new_interface( "Lomo Regression" ,
"Control the Local Monotone Regression" ,
NULL ,
PLUGIN_CALL_VIA_MENU , LOMOR_main ) ;
PLUTO_add_option( plint , "Parameters" , "Parameters" , TRUE ) ;
PLUTO_add_number( plint , "Order" , 3, 20,0,lomo_order , FALSE ) ;
PLUTO_add_number( plint , "Levels", 4,128,0,lomo_levels , FALSE ) ;
PLUTO_register_1D_function( "Lomo 1D" , LOMOR_fitter ) ;
return plint ;
}
char * LOMOR_main( PLUGIN_interface * plint )
{
PLUTO_next_option(plint) ;
lomo_order = PLUTO_get_number(plint) ;
lomo_levels = PLUTO_get_number(plint) ;
return NULL ;
}
static int lnum = 0 ;
static int * xin = NULL ;
static int * yout = NULL ;
void LOMOR_fitter( int num , double to , double dt , float * vec )
{
int ii ;
float top , bot , scl ;
if( num <= 3 || vec == NULL ) return ;
if( lnum < num ){
if( xin != NULL ){ free(xin) ; free(yout) ; }
xin = (int *) malloc( sizeof(int) * num ) ;
yout = (int *) malloc( sizeof(int) * num ) ;
if( xin == NULL || yout == NULL ) return ;
lnum = num ;
}
bot = top = vec[0] ;
for( ii=1 ; ii < num ; ii++ ){
if( vec[ii] < bot ) bot = vec[ii] ;
else if( vec[ii] > top ) top = vec[ii] ;
}
if( bot >= top ) return ;
scl = (lomo_levels - 0.01) / (top-bot) ;
fprintf(stderr,"LOMO: range %f .. %f; scl=%f\n",bot,top,scl) ;
for( ii=0 ; ii < num ; ii++ ) xin[ii] = (int)((vec[ii]-bot)*scl) ;
ii = lomo_regress( num , lomo_order , xin , yout ) ;
if( ii == -1 ) return ;
scl = 1.0 / scl ;
for( ii=0 ; ii < num ; ii++ ) vec[ii] = bot + yout[ii]*scl ;
return ;
}
/*******************************************************/
/* Fast Digital Locally Monotonic Regression */
/*******************************************************/
/* Copyright (c) 1995 */
/* University of Maryland at College Park */
/* All Rights Reserved */
/*******************************************************/
/* by Nicholas Sidiropoulos, Aug. 1995 */
/*******************************************************/
/* compile using -lm option */
/*******************************************************/
/* input: from stdinput.dat: standard matlab vector */
/* (i.e., ASCII file containing */
/* a long line of N vector elements */
/* separated by spaces) */
/* output: in stdoutput.dat.*: same format as input */
/* control: in stdcontrol_switch.dat: the */
/* ``effective'' M is */
/* the first entry here, followed by `\n` */
/* Rest: fixed to 1 (future option) */
/* So stdcontrol_switch .dat should be ASCII */
/* file starting with the effective M \n */
/* and followed by N `1''s \n, e.g., for M=10 */
/* 10<newline>1<newline>...1<newline> */
/* a total of N ones */
/*******************************************************/
/* This is just an ``Academic'' piece of software - */
/* it has been checked for correctness to the best */
/* of my ability, however, no guarantees whatsoever */
/* are given - all disclaimers here! */
/* It has been optimized for minimum development effort*/
/* and NOT for optimum speed and/or minimum comput. */
/* resources. Note, in particular, that I do not take */
/* advantage of trellis path merging to reduce storage */
/* requirements. As a result, depending on your choice */
/* of parameters M,N,R, the program may require */
/* considerable amounts of static memory. Therefore, */
/* if your computer lacks it you need to rlogin to */
/* a machine (like oxygen or apollo at SIL lab) which */
/* has sufficient memory */
/*******************************************************/
/*=====================================================*/
/*== Modified Feb 1997 by RWCox, to be a C function. ==*/
/*=====================================================*/
#include <stdio.h>
#include <math.h>
#include <string.h>
typedef struct _state {
int value, length, cumcost;
struct _state *prevstate;
} state;
#ifdef ABS
#undef ABS
#endif
#define ABS(x) (((x)<0) ? (-x) : (x))
#define USE_STATIC_TRELLIS
#ifdef USE_STATIC_TRELLIS
# define MMAX 35 /* (strictly) > max lomo degree */
# define NMAX 512 /* input data length */
# define RMAX 100 /* range of input: 0 -> R-1 inc. */
static state TRELLIS[RMAX][2*MMAX][NMAX] ;
# define trellis(i,j,k) TRELLIS[i][j][k]
#else
static state * TRELLIS = NULL ;
static int TDIM1 , TDIM2 , TDIM12 , TDIM3 ;
# define trellis(i,j,k) TRELLIS[i+j*TDIM1+k*TDIM12]
#endif
/*=====================================================*/
/*== N = number of input points ==*/
/*== m = local monotonicity order to impose ==*/
/*== yin = pointer to inputs (array of length N) ==*/
/*== x = pointer to outputs (array of length N) ==*/
/*== ==*/
/*== return value is 0 if all is OK, -1 if not ==*/
/*=====================================================*/
int lomo_regress( int N , int m , int * yin , int * x )
{
int base , R , * y ;
int n,v,l,pv,pl,i,j;
int cost, maxcost;
int peakval;
state dummy_state; /* dummy initial state */
/*== check inputs for OK-ness ==*/
if( N < 3 || m < 3 || m >= N || yin == NULL || x == NULL ) return -1 ;
/*== Compute range of input into R ==*/
y = (int *) malloc(sizeof(int)*N) ; if( y == NULL ) return -1 ;
base = yin[0] ;
for( i=1 ; i < N ; i++ ) if( yin[i] < base ) base = yin[i] ;
R = y[0] = yin[0] - base ;
for( i=1 ; i < N ; i++ ){
y[i] = yin[i] - base ; if( y[i] > R ) R = y[i] ;
}
R++ ;
if( R == 1 ){
for( i=0 ; i < N ; i++ ) x[i] = yin[i] ;
free(y) ; return 0 ;
}
fprintf(stderr,"LOMO: %d points; %d levels; %d order\n",N,R,m) ;
/* compute maxcost */
maxcost = 0;
for (n = 0; n < N; n++) maxcost += ABS(y[n]);
/* now init FLOMOR : */
dummy_state.value = -1; dummy_state.length = m ;
dummy_state.cumcost = 0 ; dummy_state.prevstate = NULL;
#ifndef USE_STATIC_TRELLIS
/*== malloc trellis space ==*/
TDIM1 = R ; TDIM2 = 2*m+2 ; TDIM12 = TDIM1*TDIM2 ; TDIM3 = N ;
TRELLIS = (state *) calloc( TDIM1*TDIM2*TDIM3 , sizeof(state) ) ;
if( TRELLIS == NULL ){ free(y) ; return -1 ; }
#endif
fprintf(stderr,"LOMO: init trellis\n") ;
for (n = 0; n < N; n++)
{
for (l = 1; l <= 2*m; l++)
{
for (v = 0; v < R; v++) { trellis(v,l,n).value = v;
trellis(v,l,n).length = l;
if ((n == 0) && ((l == 1) || (l == m+1)))
{ trellis(v,l,n).cumcost = ABS((v-y[0]));
trellis(v,l,n).prevstate = &dummy_state;
}
else
{ trellis(v,l,n).cumcost = maxcost;
trellis(v,l,n).prevstate = NULL;
}
}
}
}
/************************ main FLOMOR loop ************************************/
/* note that l = -1, ..., -m is mapped to l = m+1, ..., 2m resp. */
/******************************************************************************/
fprintf(stderr,"LOMO: compute trellis") ; fflush(stderr) ;
for (n = 1; n < N; n++)
{
/* states of the first type: (v,1) */
fprintf(stderr,".") ; fflush(stderr) ;
for (v = 0; v < R; v++)
{
for (pv = 0; pv < v; pv++)
{
for (pl = 1; pl <= m; pl++)
{
if (trellis(pv,pl,n-1).cumcost != maxcost)
{
cost = trellis(pv,pl,n-1).cumcost + ABS((v-y[n]));
if (cost < trellis(v,1,n).cumcost)
{
trellis(v,1,n).cumcost = cost;
trellis(v,1,n).prevstate = &trellis(pv,pl,n-1);
}
}
}
if (trellis(pv,2*m,n-1).cumcost != maxcost)
{
cost = trellis(pv,2*m,n-1).cumcost + ABS((v-y[n]));
if (cost < trellis(v,1,n).cumcost)
{
trellis(v,1,n).cumcost = cost;
trellis(v,1,n).prevstate = &trellis(pv,2*m,n-1);
}
}
}
}
/* states of the second type: (v,-1) */
for (v = 0; v < R; v++)
{
for (pv = R - 1; pv > v; pv--)
{
for (pl = m + 1; pl <= 2*m; pl++)
{
if (trellis(pv,pl,n-1).cumcost != maxcost)
{
cost = trellis(pv,pl,n-1).cumcost + ABS((v-y[n]));
if (cost < trellis(v,m+1,n).cumcost)
{
trellis(v,m+1,n).cumcost = cost;
trellis(v,m+1,n).prevstate = &trellis(pv,pl,n-1);
}
}
}
if (trellis(pv,m,n-1).cumcost != maxcost)
{
cost = trellis(pv,m,n-1).cumcost + ABS((v-y[n]));
if (cost < trellis(v,m+1,n).cumcost)
{
trellis(v,m+1,n).cumcost = cost;
trellis(v,m+1,n).prevstate = &trellis(pv,m,n-1);
}
}
}
}
/* states of the third type: (v,l), 1 < l < m */
for (l = 2; l < m; l++)
{
for (v = 0; v < R; v++)
{
if (trellis(v,l-1,n-1).cumcost != maxcost)
{
trellis(v,l,n).cumcost = trellis(v,l-1,n-1).cumcost + ABS((v-y[n]));
trellis(v,l,n).prevstate = &trellis(v,l-1,n-1);
}
}
}
/* states of the fourth type: (v,l), -m < l < -1 */
for (l = m+2; l < 2*m; l++)
{
for (v = 0; v < R; v++)
{
if (trellis(v,l-1,n-1).cumcost != maxcost)
{
trellis(v,l,n).cumcost = trellis(v,l-1,n-1).cumcost + ABS((v-y[n]));
trellis(v,l,n).prevstate = &trellis(v,l-1,n-1);
}
}
}
/* states of the fifth type: (v,l), l = m */
for (v = 0; v < R; v++)
{
if (trellis(v,m,n-1).cumcost != maxcost)
{
cost = trellis(v,m,n-1).cumcost + ABS((v-y[n]));
if (cost < trellis(v,m,n).cumcost)
{
trellis(v,m,n).cumcost = cost;
trellis(v,m,n).prevstate = &trellis(v,m,n-1);
}
}
if (trellis(v,m-1,n-1).cumcost != maxcost)
{
cost = trellis(v,m-1,n-1).cumcost + ABS((v-y[n]));
if (cost < trellis(v,m,n).cumcost)
{
trellis(v,m,n).cumcost = cost;
trellis(v,m,n).prevstate = &trellis(v,m-1,n-1);
}
}
}
/* states of the sixth type: (v,l), l = -m */
for (v = 0; v < R; v++)
{
if (trellis(v,2*m,n-1).cumcost != maxcost)
{
cost = trellis(v,2*m,n-1).cumcost + ABS((v-y[n]));
if (cost < trellis(v,2*m,n).cumcost)
{
trellis(v,2*m,n).cumcost = cost;
trellis(v,2*m,n).prevstate = &trellis(v,2*m,n-1);
}
}
if (trellis(v,2*m-1,n-1).cumcost != maxcost)
{
cost = trellis(v,2*m-1,n-1).cumcost + ABS((v-y[n]));
if (cost < trellis(v,2*m,n).cumcost)
{
trellis(v,2*m,n).cumcost = cost;
trellis(v,2*m,n).prevstate = &trellis(v,2*m-1,n-1);
}
}
}
}
/************************ eof main FLOMOR loop *******************************/
/* now pick best path, and write it to x[n] */
fprintf(stderr,"\nLOMO: pick best path\n") ; fflush(stderr) ;
cost = maxcost;
for (l = 1; l <= m; l++)
{
for (v = 0; v < R; v++)
{
if (trellis(v,l,N-1).cumcost < cost)
{
cost = trellis(v,l,N-1).cumcost;
pl = l; pv = v;
}
}
}
/* now traverse backwards */
if (cost < maxcost)
{
x[N-1] = pv;
fprintf(stderr," [%d] = %d\n",N-1,pv) ;
for (n = N-2; n >=0; n--)
{
fprintf(stderr," [%d] = %d",n,trellis(pv,pl,n+1).prevstate->value) ; fflush(stderr) ;
x[n] = trellis(pv,pl,n+1).prevstate->value;
fprintf(stderr," pv = %d",trellis(pv,pl,n+1).prevstate->value) ; fflush(stderr) ;
pv = trellis(pv,pl,n+1).prevstate->value;
fprintf(stderr," pl = %d",trellis(pv,pl,n+1).prevstate->length) ; fflush(stderr) ;
pl = trellis(pv,pl,n+1).prevstate->length;
fprintf(stderr,"\n") ;
}
}
else { for (n=0;n<N;n++) x[n] = 0; /* as good as any... */ }
/*== add base back to output ==*/
fprintf(stderr,"\nLOMO: add base back on\n") ;
for( i=0 ; i < N ; i++ ) x[n] += base ;
fprintf(stderr,"LOMO: exit regression\n") ;
free(y) ;
#ifndef USE_STATIC_TRELLIS
free(TRELLIS) ;
#endif
return 0 ;
}
#endif /* ALLOW_LOMO */
syntax highlighted by Code2HTML, v. 0.9.1