/*---------------------------------------------------------------------------*/
/*
  This file contains routines for implementing the 3ddelay functions.
	Based on fim+.c  by B. Douglas Ward

  This software is copyrighted and owned by the Medical College of Wisconsin.
  See the file README.Copyright for details.

*/

/*---------------------------------------------------------------------------*/
/*
  Include software for linear regression analysis and sorting numbers.
*/

#include "RegAna.c"
#include "ranks.c"

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

/*---------------------------------------------------------------------------*/
/* 
   Initialize independent variable X matrix 
*/

int init_indep_var_matrix 
(
  int q,                      /* number of parameters in the baseline model */
  int p,                      /* number of parameters in the baseline model 
			         plus number of ideals */
  int NFirst,                 /* index of first image to use in the analysis */
  int N,                      /* total number of images used in the analysis */
  int polort,                 /* degree of polynomial ort */
  int num_ort_files,          /* number of ort time series files */
  int num_ideal_files,        /* number of ideal time series files */
  MRI_IMAGE ** ort_array,     /* ort time series arrays */
  int ** ort_list,            /* list of ort time series */
  MRI_IMAGE ** ideal_array,   /* ideal time series arrays */
  int ** ideal_list,          /* list of ideal time series */
  float * x_bot,              /* minimum of stimulus time series */
  float * x_ave,              /* average of stimulus time series */
  float * x_top,              /* maximum of stimulus time series */
  int * good_list,            /* list of good time points */
  matrix * x                  /* independent variable matrix */
)

{
  const int BIG_NUMBER = 33333;
  int i=-1;                    /* time index */
  int m=-1;                    /* X matrix column index */
  int n=-1;                    /* X matrix row index */
  int is=-1;                   /* input ideal index */
  float * far = NULL;
  int nx=-1, ny=-1, iq=-1, nq=-1;
  int Ngood=-1;
  int maxcols=-1;
  matrix xgood;

   if (num_ort_files) {
      fprintf(stderr,"Not ready to deal with ort files at all!\n");
      exit(1);
   }
  /*----- Initialize X matrix -----*/
  maxcols = 0;
  for (is = 0;  is < num_ort_files;  is++) {/* add orts */
      if (ort_list[is] == NULL) maxcols += ort_array[is]->ny; 
      else maxcols += ort_list[is][0]+1;
  } 
  for (is = 0;  is < num_ideal_files;  is++) {/* add ideals */
      if (ideal_list[is] == NULL) maxcols += ideal_array[is]->ny;
      else maxcols += ideal_list[is][0]+1;
  }
  maxcols = maxcols + polort + 1;  /* add polorts */
  /* fprintf(stderr,"N=%d,p=%d, polort=%d, maxcols = %d\n", N, p, polort, maxcols); */
  matrix_create (N+NFirst, maxcols, x);
  matrix_initialize (&xgood);

  /*----- Set up columns of X matrix corresponding to polynomial orts -----*/
  for (m = 0;  m <= polort;  m++)
    for (n = 0;  n < N;  n++)
      {
	i = n + NFirst;
	(*x).elts[n][m] = pow ((double)i, (double)m);
      }

  /*----- Set up columns of X matrix corresponding to ort time series -----*/
  for (is = 0;  is < num_ort_files;  is++)
    {
      far = MRI_FLOAT_PTR (ort_array[is]);
      nx = ort_array[is]->nx;
      ny = ort_array[is]->ny;

      if (ort_list[is] == NULL)
	for (iq = 0;  iq < ny;  iq++)
	  {
	    for (n = 0;  n < N;  n++)
	      {
		i = n + NFirst;
		(*x).elts[n][m] = *(far + iq*nx + i);
	      }
	    m++;
	  }
      else
	{
	  nq = ort_list[is][0];
	  for (iq = 1;  iq <= nq;  iq++)
	    {
	      for (n = 0;  n < N;  n++)
		{
		  i = n + NFirst;
		  (*x).elts[n][m] = *(far + ort_list[is][iq]*nx + i);
		}
	      m++;
	    }
	}
    }

  /*----- Set up columns of X matrix corresponding to ideal time series -----*/
  for (is = 0;  is < num_ideal_files;  is++)
    {
      far = MRI_FLOAT_PTR (ideal_array[is]);
      nx = ideal_array[is]->nx;
      ny = ideal_array[is]->ny;

      if (ideal_list[is] == NULL)
	for (iq = 0;  iq < ny;  iq++)
	  {
	    for (n = 0;  n < N;  n++)
	      {
		i = n + NFirst;
		(*x).elts[n][m] = *(far + iq*nx + i);
	      }
	    
	    m++;
	  }
      else
	{
	  nq = ideal_list[is][0];
	  for (iq = 1;  iq <= nq;  iq++)
	    {
	      for (n = 0;  n < N;  n++)
		{
		  i = n + NFirst;
		  (*x).elts[n][m] = *(far + ideal_list[is][iq]*nx + i);
		}

	      m++;
	    }
	}
    }

  /*----- Remove row if ideal contains value over 33333 -----*/
  Ngood = N;
  m = 0;
  for (n = 0;  n < N;  n++)
    {
      for (is = q;  is < p;  is++)
	{
	  if ((*x).elts[n][is] >= BIG_NUMBER)  break;
	}
      if (is < p)
	{
	  Ngood--;
	}
      else
	{
	  good_list[m] = n;
	  m++;
	}
    }
  matrix_extract_rows ((*x), Ngood, good_list, &xgood);
  matrix_equate (xgood, x);

  /*----- Find min, max, and ave for each column of the X matrix -----*/
  for (is = 0;  is < p;  is++)
    {      
      x_bot[is] = x_top[is] = (*x).elts[0][is];
      x_ave[is] = 0.0;
      for (n = 0;  n < Ngood;  n++)
	{
	  if ((*x).elts[n][is] < x_bot[is])  x_bot[is] = (*x).elts[n][is];  
	  if ((*x).elts[n][is] > x_top[is])  x_top[is] = (*x).elts[n][is];
	  x_ave[is] += (*x).elts[n][is] / Ngood;
	}
    }
  
  
  matrix_destroy (&xgood);

  return (Ngood);

}

/*---------------------------------------------------------------------------*/
/*
  Initialization for the delay analysis.
*/

int init_delay 
(
  int q,                      /* number of parameters in the baseline model */
  int p,                      /* number of parameters in the baseline model 
			         plus number of ideals */
  int N,                      /* total number of images used in the analysis */
  int num_idealts,            /* number of ideal time series */
  matrix xdata,               /* independent variable matrix */
  matrix * x_base,            /* extracted X matrix    for baseline model */
  matrix * xtxinvxt_base,     /* matrix:  (1/(X'X))X'  for baseline model */
  matrix * x_ideal,           /* extracted X matrices  for ideal models */
  matrix * xtxinvxt_ideal,    /* matrix:  (1/(X'X))X'  for ideal models */
  float ** rarray             /* ranked arrays of ideal time series */
)

{
  int * plist = NULL;         /* list of model parameters */
  int ip, it;                 /* parameter indices */
  int is, js;                 /* ideal indices */ 
  int jm;                     /* lag index */
  int ok;                     /* flag for successful matrix calculation */
  matrix xtxinv_temp;         /* intermediate results */
  vector ideal;               /* ideal vector */
  vector coef_temp;           /* intermediate results */
  vector xres;                /* vector of residuals */
  float sse_base;             /* baseline model error sum of squares */ 
        

  /*----- Initialize matrix -----*/
  matrix_initialize (&xtxinv_temp);
  vector_initialize (&ideal);
  vector_initialize (&coef_temp);
  vector_initialize (&xres);

  /*----- Allocate memory -----*/
  plist = (int *) malloc (sizeof(int) * p);   MTEST (plist);

  /*----- Initialize matrices for the baseline model -----*/
  for (ip = 0;  ip < q;  ip++)
    plist[ip] = ip;
  ok = calc_matrices (xdata, q, plist, x_base, &xtxinv_temp, xtxinvxt_base);
  if (!ok)  { matrix_destroy (&xtxinv_temp);  return (0); };

  /*----- Initialize matrices for ideal functions -----*/
  for (is = 0;  is < num_idealts;  is++)
    {
      for (ip = 0;  ip < q;  ip++)
	{
	  plist[ip] = ip;
	}

      plist[q] = q + is;

      ok = calc_matrices (xdata, q+1, plist, 
			  &(x_ideal[is]), &xtxinv_temp, &(xtxinvxt_ideal[is]));
      if (!ok)  { matrix_destroy (&xtxinv_temp);  return (0); };
    }

  /*----- Set up the ranked array for each ideal -----*/
  for (is = 0;  is < num_idealts;  is++)
    {
      /*----- Convert column of matrix to vector -----*/
      column_to_vector (xdata, q+is, &ideal);

      /*----- Calculate regression coefficients for baseline model -----*/
      calc_coef (*xtxinvxt_base, ideal, &coef_temp);

      /*----- Calculate the error sum of squares for the baseline model -----*/
      sse_base = calc_resids (*x_base, coef_temp, ideal, &xres);
    
      /*----- Form rank array from residual array -----*/
      rarray[is] = rank_darray (N, xres.elts);

    }

  /*----- Destroy matrix -----*/
  matrix_destroy (&xtxinv_temp);
  vector_destroy (&ideal);
  vector_destroy (&coef_temp);
  vector_destroy (&xres);

  /*----- Deallocate memory -----*/
  free (plist);   plist = NULL;

  return (1);
}

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

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

/*---------------------------------------------------------------------------*/
/* 
   Calculate the sign function.
*/

float sign (float x)

{
  if (x > 0.0)  return (1.0);
  if (x < 0.0)  return (-1.0);
  return (0.0);
}

/*---------------------------------------------------------------------------*/
int Read_part_file_delay(float *x,
					char *f_name,
					int a,
					int b)
   
   { 
      float *far=NULL;
      int i=0, cnt=0;
      MRI_IMAGE *im_data=NULL;
      
      im_data = mri_read_1D(f_name);
      if (!im_data) { fprintf(stderr,"Failed to read 1D file."); return(-1);}
      far = MRI_FLOAT_PTR(im_data);
      if (im_data->ny != 1) { fprintf(stderr,"More than one column in 1D file."); return(-1);}
      if (b>=im_data->nx) { fprintf(stderr,"Nlast >= number of values in 1D file!"); return(-1);}
      cnt = 0;
      for (i=a;i<=b;++i) {
         x[cnt] = far[i];
         #ifdef ZDBG
            if (cnt < 10) { fprintf(stderr,"%d: %.3f\t", i, x[cnt]);}
         #endif
         ++cnt;
      }
      mri_free(im_data);  im_data = NULL;
      
      return(cnt);
   
   }
   
int Read_part_file_delay_OLD (float *x,
					char *f_name,
					int a,
					int b)
   
    { 
     
     int cnt=0,ex,line_num;
     float buf;
     FILE*file_in;
     
     file_in = fopen (f_name,"r");
     if (file_in == NULL) {
            printf ("\aCould not open %s \n",f_name);
           printf ("Exiting @ Read_file function\n");
            exit (0);
         }
     
     if (a > b || a==0) {
     				printf ("\a\n\33[1mError in (from , to) line numbers\n\33\[0m");
     				printf ("Exiting @Read_part_file function \n");
     				exit (0);
     		   }
     
     line_num = 1;	
     if (a == 1) {
     			ex = fscanf (file_in,"%f",&x[cnt]);	
     			++cnt;
     			}				   	
      else  ex = fscanf (file_in,"%f",&buf);					   	
     ++line_num;
     while (ex != EOF && line_num <= b)
      {
        if (line_num >=a && line_num <=b) 
        {
         ex = fscanf (file_in,"%f",&x[cnt]);
         ++cnt;
         if (ex == EOF) --cnt;
         }
        else 
        {
         ex = fscanf (file_in,"%f",&buf);
         }
        ++line_num;
        
      }
      
      if (ex == EOF) 
      	{
      	    --line_num;
      		printf ("\n\33[1mEOF reached before line \33[0m%d\n",b);
      		printf ("Only %d lines were read, from line %d to %d\n",cnt,a,line_num-1);
      	}
      
      fclose (file_in);
      return (cnt);  							     
   }

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

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

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