/*****************************************************************************
   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.
******************************************************************************/
   
/***
     header file for recursive partial correlation calculations
     -- RWCox, Feb 1994
     -- March 1995: this work is the basis for the paper
         "Real-Time Functional Magnetic Resonance Imaging"
         by RW Cox, A Jesmanowicz, and JS Hyde,
         Magnetic Resonance in Medicine, 33:230-236.
     -- June 1995: added get_variance and some comments
***/

#ifndef PCOR_HEADER
#define PCOR_HEADER

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/

/*** some macros ***/

#ifndef  SQR
#define SQR(x)  ((x)*(x))
#endif

#ifndef MAX
#  define MAX(x,y) (((x)>(y)) ? (x) : (y))
#endif

#ifndef MIN
#  define MIN(x,y) (((x)<(y)) ? (x) : (y))
#endif

/*** default is to expand the inner loop in update_voxel_corr for
     small values of nref;  if this is not desired, compile with
     the switch -DNO_EXPAND_UPDATE, or comment the next line out ***/

#define EXPAND_UPDATE

#ifdef NO_EXPAND_UPDATE
#undef EXPAND_UPDATE
#endif

/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/

/*** Define atomic data types ***/

#ifndef REF_FLOAT_SINGLE           /* compile time option to save space */
   typedef double ref_float ;      /* internal storage of reference data */
#  define REF_EPS  1.0e-13         /* a small number of this type */
#else
   typedef float ref_float ;
#  define REF_EPS  1.0e-7
#endif  /* REF_FLOAT_SINGLE */

/*-------------------------------*/

#ifndef VOX_SHORT                  /* define voxel data as shorts or floats */
   typedef float vox_data ;
#else
   typedef short vox_data ;
#endif  /* VOX_SHORT */

/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/

/*** this type holds the independent-of-voxel references information ***/

   /* note that the partial correlation coefficient of the data vectors
      is calculated with respect to the LAST reference vector,
      removing the effects of all the previous reference vectors */

typedef struct {
          int nref ;                    /* number of ref vectors */
          int nupdate ;                 /* number of times has been updated */
          ref_float **chol ,            /* nref X nref Cholesky factor */
                    *alp , *ff , *gg ;  /* alpha, f, and g factors */
          ref_float betasq ;            /* 1/beta^2 factor */
   } references ;

 /*** RCH is for access to the (ii,jj) element
      of the Cholesky factor of a given set of references;
      ii >= jj, since this is the lower triangular factor

      N.B.: compared to the paper's notation,
            L+1     = nref
            c_{i,j} = RCH(rr,i-1,j-1)  for i=1..L+1 , j=1..i,
        and c_{i,j} = VCH(vv,vox,j-1)  for i=L+2, j=1..L+2    ***/

#define RCH(rr,ii,jj)  (rr->chol[(ii)][(jj)])

#ifdef OV_DEBUG1
#define REF_DUMP(rr,str) \
  {  int iq,jq ; ref_float qsum ; \
     fprintf(stderr,"%s: reference dump, nref=%d betasq=%11.4e\n", \
             str , rr->nref , rr->betasq ) ; \
     for( iq=0 ; iq < rr->nref ; iq++){ \
       fprintf(stderr," ROW %d: ",iq) ; \
       qsum = 0.0 ; \
       for( jq=0 ; jq <= iq ; jq++ ){ \
          fprintf(stderr,"%11.4e ",RCH(rr,iq,jq)); \
          qsum += SQR(RCH(rr,iq,jq)) ; \
       } \
       fprintf(stderr,": qsum=%11.4e\n",qsum) ; \
       fprintf(stderr,"      alpha=%11.4e  ff=%11.4e  gg=%11.4e\n", \
               rr->alp[iq] , rr->ff[iq] , rr->gg[iq] ) ; \
     } \
   }
#endif

/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/

/*** this type holds all the voxel-specific data ***/

typedef struct {
          int nvox ;          /* number of voxels */
          int nref ;          /* number of references to allow for */
          int nupdate ;       /* number of times it has been updated */
          ref_float *chrow ;  /* last row (length nref+1)
                                 of Cholesky factor for each voxel
                                 (N.B.:  last element is actually squared) */
   } voxel_corr ;

   /* VCH is for access to the jj-th element of the last row of
      the Cholesky factor for the vox-th voxel in a given voxel_corr struct;

      N.B.:  the last element [VCH(vv,vox,nref)] is the SQUARE of
             the actual Cholesky diagonal element, since that is all that the
             partial correlation coefficient algorithm actually needs;
             this prevents taking a square root at each time step for each voxel */

#define VCH(vv,vox,jj) (vv->chrow[(jj)+(vox)*(vv->nref+1)])

#ifdef OV_DEBUG1
#define VD 2001
#define VOX_DUMP(vv,vox,str) \
   {  int jq ; ref_float qsum = 0.0 ; \
      fprintf(stderr,"%s: voxel_corr dump #%d\n  ",str,vox) ; \
      for( jq=0 ; jq < vv->nref ; jq++ ){ \
        fprintf(stderr,"%11.4e ",VCH(vv,vox,jq)) ; \
        qsum += SQR(VCH(vv,vox,jq)) ; \
      } \
      qsum += VCH(vv,vox,vv->nref) ; \
      fprintf(stderr,"%11.4e : qsum=%11.4e\n",VCH(vv,vox,vv->nref),qsum); \
   }
#endif

/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/

/*** this type holds the return data for thresholding results ***/

typedef struct {
          int   num_pcor_pos , num_pcor_neg , num_coef_pos , num_coef_neg ;
          float max_pcor_pos , max_pcor_neg , max_coef_pos , max_coef_neg ;
   } thresh_result ;

#ifdef OV_DEBUG1
#define THR_DUMP(thr,str) \
fprintf(stderr,"thresh_results dump for %s\n:") ; \
fprintf(stderr," num_pcor_pos=%d neg=%d  num_coef_pos=%d neg=%d\n", \
 (thr).num_pcor_pos,(thr).num_pcor_neg,(thr).num_coef_pos,(thr).num_coef_neg ) ; \
fprintf(stderr," max_pcor_pos=%11.4g neg=%11.4g  max_coef_pos=%11.4g neg=%11.4g\n", \
 (thr).max_pcor_pos,(thr).max_pcor_neg,(thr).max_coef_pos,(thr).max_coef_neg ) ;
#endif

/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/

/*** prototypes ***/

extern references * new_references(int numref) ;

extern void update_references(float *vec, references *ref) ;

extern voxel_corr * new_voxel_corr(int numvox, int numref) ;

extern void free_references(references *ref) ;

extern void free_voxel_corr(voxel_corr *vc) ;

extern void update_voxel_corr(vox_data *vdata, references *ref, voxel_corr *vc);

extern void get_pcor(references *ref, voxel_corr *vc, float *pcor) ;

extern void get_coef(references *ref, voxel_corr *vc, float *coef) ;

extern void get_variance(voxel_corr *vc, float *var) ;

extern void get_lsqfit(references *ref, voxel_corr *vc, float *fit[]) ;

extern void get_pcor_thresh_coef( references *ref, voxel_corr *vc,
                                  float pcthresh, float cothresh, float *pcor,
                                  float *coef, thresh_result *thr ) ;

#endif
