 /*
  * Khoros: $Id: ldminorm.c,v 1.3 1992/03/20 23:28:08 dkhoros Exp $
  */

#if !defined(lint) && !defined(SABER)
static char rcsid[] = "Khoros: $Id: ldminorm.c,v 1.3 1992/03/20 23:28:08 dkhoros Exp $";
#endif

 /*
  * $Log: ldminorm.c,v $
 * Revision 1.3  1992/03/20  23:28:08  dkhoros
 * VirtualPatch5
 *
  */

/*
 *----------------------------------------------------------------------
 *
 * Copyright 1992, University of New Mexico.  All rights reserved.
 * Permission to copy and modify this software and its documen-
 * tation only for internal use in your organization is hereby
 * granted, provided that this notice is retained thereon and
 * on all copies.  UNM makes no representations as to the sui-
 * tability and operability of this software for any purpose.
 * It is provided "as is" without express or implied warranty.
 * 
 * UNM DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
 * INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FIT-
 * NESS.  IN NO EVENT SHALL UNM BE LIABLE FOR ANY SPECIAL,
 * INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY OTHER DAMAGES WHAT-
 * SOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
 * IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
 * ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PER-
 * FORMANCE OF THIS SOFTWARE.
 * 
 * No other rights, including, for example, the right to redis-
 * tribute this software and its documentation or the right to
 * prepare derivative works, are granted unless specifically
 * provided in a separate license agreement.
 *---------------------------------------------------------------------
 */

#include "unmcopyright.h"        /* Copyright 1992 by UNM */

/*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>  <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 >>>>
 >>>>         File Name: ldminorm.c
 >>>>
 >>>>      Program Name: dminorm
 >>>>
 >>>> Date Last Updated: Thu Mar  5 11:05:25 1992 
 >>>>
 >>>>          Routines: ldminorm - the library call for dminorm
 >>>>
 >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>   <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<*/


#include "vinclude.h"


/* -library_includes */
/* clean up code */
#define DMUS_CLEANUP \
{  /* cleanup workspace before leaving */ \
   if(matr!=NULL)free(matr); \
   if(mati!=NULL)free(mati); \
   if(fv1!=NULL)free(fv1); \
   if(fv2!=NULL)free(fv2); \
   if(fm1!=NULL)free(fm1); \
   if(crit!=NULL)free(crit); \
   if(eigvctr!=NULL)free(eigvctr); \
   if(eigvcti!=NULL)free(eigvcti); \
   if(eigval!=NULL)free(eigval); \
   if(rspectra!=NULL)free(rspectra); \
   if(ispectra!=NULL)free(ispectra); \
   if(temp!=NULL)free(temp); \
   if(spectra!=NULL){ \
      for(i=0;i<num_vects;i++) if(spectra[i]!=NULL)free(spectra[i]); \
      free(spectra); \
   } \
   if(psd!=NULL){ \
      for(i=0;i<num_vects;i++) if(psd[i]!=NULL)free(psd[i]); \
      free(psd); \
   } \
   if(data!=NULL){ \
      for(i=0;i<num_vects;i++) if(data[i]!=NULL)free(data[i]); \
      free(data); \
   } \
} /* end of ldminorm cleanup code */

#define AIC_CRIT  0
#define MDL_CRIT -1
#define CENTER    1
#define NOT_CENTER 0
#define EPSILON 0.0000001

static int lmus_center();
static void get_aic_crit(), get_mdl_crit();
static void m_eval_poly();

/* -library_includes_end */


/****************************************************************
*
* Routine Name: ldminorm - library call for dminorm
*
* Purpose:
*    
*    1D Minimum-Norm Spectral Estimation
*    
*    

* Input:
*    
*    image          pointer to VIFF structure containing image data to
*                   be processed.
*    
*    order          order of the linear combiner.
*    
*    freq           number of computed spectral values.
*    
*    kdim           used to determine the dimension of the noise  sub-
*                   space.   If  it  is  a positive, non-zero integer,
*                   then it represents the actual value of  the  noise
*                   subspace  dimension.  If it is zero, then ldminorm
*                   will use the Akaike (AIC)  criteria  to  determine
*                   the  noise  subspace dimension.  If it is -1, then
*                   the minimum description length (MDL) criteria will
*                   be used.
*    
*    auto_type      type of auto correlation to be used.  0  specifies
*                   a  biased  estimate, 1 specifies an unbiased esti-
*                   mate, 2 specifies an FFT based estimate.
*    
*    arith_type     type of arithmetic to be used on complex data.   0
*                   specifies scalar arithmetic and 1 specifies vector
*                   arithmetic.
*    
*    psd_type       specifies the centering of the power spectral den-
*                   sity estimate.  0 specifies not centered, 1 speci-
*                   fies centered.
*    
*    procdir        process direction:  0  indicated  vector  oriented
*                   processing, 1 indicates band oriented processing.
*    
*    

* Output:
*    
*    image          pointer to VIFF structure  containing  image  data
*                   after processing.
*    
*    Return Value:  1 on success, 0 on failure.
*    
*    

*
* Written By: Jeremy Worley, Ramiro Jordan, Glen Abousleman
*    
*    Jeremy Worley 30 Jul 1990 08:38 MST
*              Changed internal routines to  static  to  prevent  con-
*              flicts with other routines in the library.
*    
*    Jeremy Worley 12 Mar 1991 13:51 MST
*              Fixed bug in static lminorm  that  generated  backwards
*              spectrum estimate.
*    
*    Jeremy Worley 05 Mar 1992 11:04 MST
*              Added many  explicit  function  declarations  in  ldmi-
*              norm(),  lrminorm()  and  lcminorm().  Fixed an fprintf
*              call in ldminorm().  Added explicit  argument  declara-
*              tions for kdim in lrminorm() and lcminorm().
*    
*    

****************************************************************/


/* -library_def */
int ldminorm(image,order,freq,kdim,auto_type,arith_type,psd_type,proc_dir)
  struct xvimage *image;
  int order,freq,kdim,auto_type,arith_type,psd_type,proc_dir;
/* -library_def_end */

/* -library_code */
{
  float **data, **psd, *temp=NULL;  
  double *matr=NULL,*mati=NULL,*fv1=NULL,*fv2=NULL,*fm1=NULL,
         *eigvctr=NULL,*eigvcti=NULL,*eigval=NULL, *crit=NULL,
         **spectra=NULL,*rspectra=NULL,*ispectra=NULL;
  char  **dload_vector();
  int   dtype,j,i,n,num_vects,dimension;
  char *program = "ldminorm";
  int propertype(), lrminorm(), lcminorm(), dunload_vector();

/*
** check for proper data types.
*/
 
    if( ! (propertype(program,image,VFF_TYP_FLOAT,FALSE)) && 
        ! (propertype(program,image,VFF_TYP_COMPLEX,FALSE)) ) { 
          fprintf(stderr,"%s: data storage type must be float or complex.", 
                  program); 
          return(0); 
    }

/*
** keep track of data type for later.
*/

   dtype = 1;
   if(image->data_storage_type==VFF_TYP_COMPLEX)dtype = 2;
      
/*
** reorganize the data into our internal representation
*/

  if((data = (float **)dload_vector(image, &num_vects, &dimension, proc_dir))
              == NULL)
  {
      (void) fprintf(stderr,"%s: dload_vector failed \n",program);
      return(0);
  }

/*
** allocate psd arrays.  the spectra array is actually a double precision
** temporary for psd.  This data space is used when unload_vector is called
** because it accepts floats, not doubles.
*/

  psd = (float **)malloc((unsigned)num_vects * sizeof(float *));
  if(psd==NULL){
     fprintf(stderr,"%s:  [1]memory allocation failure.\n",program);
     DMUS_CLEANUP; 
     return(0);
  }
  for(i=0;i<num_vects;i++){
      psd[i] = (float *)malloc((unsigned)freq *sizeof(float));
      if(psd[i]==NULL){
         fprintf(stderr,"%s:  [2]memory allocation failure.\n",program);
         DMUS_CLEANUP;
         return(0);
      }
  }

/*
** spectra will hold the spectral estimate when returning from lrminorm and
** lcminorm.  Since all operations are done in double precision, this work
** space is also in double precision.
*/
  
  spectra = (double **)malloc((unsigned)num_vects * sizeof(double *));
  if(spectra==NULL){
     fprintf(stderr,"%s:  [3]memory allocation failure.\n",program);
     DMUS_CLEANUP;
     return(0);
  }
  for(i=0;i<num_vects;i++){
      spectra[i] = (double *)malloc((unsigned)freq *sizeof(double));
      if(spectra[i]==NULL){
         fprintf(stderr,"%s:  [4]memory allocation failure.\n",program);
         DMUS_CLEANUP;
         return(0);
      }
  }

/*
** the next three work arrays are used by Eispack routines is some
** mysterious way.   
*/

  fv1 = (double *)malloc((unsigned)order * sizeof(double));
  if(fv1==NULL){
     fprintf(stderr,"%s:  [5]memory allocation failure.\n",program);
     DMUS_CLEANUP;
     return(0);
  }

  fv2 = (double *)malloc((unsigned)order * sizeof(double));
  if(fv2==NULL){
     fprintf(stderr,"%s:  [6]memory allocation failure.\n",program);
     DMUS_CLEANUP;
     return(0);
  }

  fm1 = (double *)malloc((unsigned)order * 2 * sizeof(double));
  if(fm1==NULL){
     fprintf(stderr,"%s:  [7]memory allocation failure.\n",program);
     DMUS_CLEANUP;
     return(0);
  }

/*
** crit is used as work space for the AIC and MDL criterion calculation
*/

  crit = (double *)malloc((unsigned)order * dtype * sizeof(double));
  if(crit==NULL){
     fprintf(stderr,"%s:  [8]memory allocation failure.\n",program);
     DMUS_CLEANUP;
     return(0);
  }

/*
** allocate memory for covariance matrix
*/

  matr = (double *)malloc((unsigned)order * order * sizeof(double));
  if(matr==NULL){
     fprintf(stderr,"%s:  [9]memory allocation failure.\n",program);
     DMUS_CLEANUP;
     return(0);
  }

  mati = (double *)malloc((unsigned)order * order * sizeof(double));
  if(mati==NULL){
     fprintf(stderr,"%s:  [10]memory allocation failure.\n",program);
     DMUS_CLEANUP;
     return(0);
  }

/*
** allocate memory for output of vectors and values
*/

  eigvctr = (double *)malloc((unsigned)order * order * sizeof(double));
  if(eigvctr==NULL){
     fprintf(stderr,"%s:  [11]memory allocation failure.\n",program);
     DMUS_CLEANUP;
     return(0);
  }

  eigvcti = (double *)malloc((unsigned)order * order * sizeof(double));
  if(eigvcti==NULL){
     fprintf(stderr,"%s:  [12]memory allocation failure.\n",program);
     DMUS_CLEANUP;
     return(0);
  }

  eigval = (double *)malloc((unsigned)order * sizeof(double));
  if(eigval==NULL){
     fprintf(stderr,"%s:  [13]memory allocation failure.\n",program);
     DMUS_CLEANUP;
     return(0);
  }

/*
** allocate memory for scalar arith type processing
*/

  if(!arith_type && dtype == 2){
     temp = (float *)malloc((unsigned)dimension * sizeof(float));
     if(temp==NULL){
        fprintf(stderr,"%s:  [14]memory allocation failure.\n",program);
        DMUS_CLEANUP;
        return(0);
     }
     rspectra = (double *)malloc((unsigned)freq * sizeof(double));
     if(rspectra==NULL){
        fprintf(stderr,"%s:  [15]memory allocation failure.\n",program);
        DMUS_CLEANUP;
        return(0);
     }
     ispectra = (double *)malloc((unsigned)freq * sizeof(double));
     if(ispectra==NULL){
        fprintf(stderr,"%s:  [16]memory allocation failure.\n",program);
        DMUS_CLEANUP;
        return(0);
     }
  }

/*
** Now that all of the memory allocation crap is out of the way, lets do some
** real work. 
**
** now use big loop to perform minorm algorithm on all vectors
*/

  for(n=0; n<num_vects;n++){
      if(dtype==1){
         if(!lrminorm(data[n],spectra[n],dimension,order,freq,kdim,crit,
                     matr,mati,fv1,fv2,fm1,eigvctr,eigvcti,eigval,auto_type)){
            fprintf(stderr,"%s:  [1]Failure in lrminorm.\n",program);
            DMUS_CLEANUP;
            return(0);
         }
      }else if(arith_type==1){
         if(!lcminorm(data[n],spectra[n],dimension,order,freq,kdim,crit,
                     matr,mati,fv1,fv2,fm1,eigvctr,eigvcti,eigval,auto_type)){
            fprintf(stderr,"%s:  Failure in lcminorm.\n",program);
            DMUS_CLEANUP;
            return(0);
         }
      }else{
         for(i=0,j=0;i<dimension;i++,j+=2){
             temp[i] = data[n][j];
         }
         if(!lrminorm(temp,rspectra,dimension,order,freq,kdim,crit,
                     matr,mati,fv1,fv2,fm1,eigvctr,eigvcti,eigval,auto_type)){
            fprintf(stderr,"%s:  [2]Failure in lrminorm.\n",program);
            DMUS_CLEANUP;
            return(0);
         }
         for(i=0,j=1;i<dimension;i++,j+=2){
             temp[i] = data[n][j];
         }
         if(!lrminorm(temp,ispectra,dimension,order,freq,kdim,crit,
                     matr,mati,fv1,fv2,fm1,eigvctr,eigvcti,eigval,auto_type)){
            fprintf(stderr,"%s:  [3]Failure in lrminorm.\n",program);
            DMUS_CLEANUP;
            return(0);
         }
         for(i=0;i<freq;i++){
             spectra[n][2*i] = rspectra[i];
             spectra[n][2*i+1] = ispectra[i];
         }
      }
  } /* end for statement */

/*
** Center the data if that was requested
*/

  if(arith_type == 0)freq *= 2;

  if(psd_type==CENTER){
     if(!lmus_center(spectra,num_vects,freq)){
        fprintf(stderr,"%s:  Something went wrong while centering data.\n",
                program);
        return(0);
     }
  }

/*
** return data to its original format and exit program
*/

  for(i=0;i<num_vects;i++){
      for(j=0;j<freq;j++){
          psd[i][j] = spectra[i][j];
      }
  }

  if(arith_type || dtype == 1){
     if (!dunload_vector((char **)psd, image, (unsigned long)VFF_TYP_FLOAT, 
        num_vects, freq, proc_dir))
       {
          (void) fprintf (stderr,"%s: dunload_vector failed \n",program);
          return(0);
       }
  }else{
     freq /= 2;
     if(!dunload_vector((char **)psd, image, (unsigned long)VFF_TYP_COMPLEX, 
        num_vects, freq, proc_dir))
       {
          fprintf(stderr,"%s:  dunload_vector failed\n",program);
          return(0);
       }
  }

    DMUS_CLEANUP;
    return(1);
}

/***********************************************************************
*
*  Routine Name: lrminorm
*
*          Date:
*        
*       Purpose: Sets up matrices, etc and calls lminorm on real data set 
*
*         Input: 
*
*        Output: 
*
*    Written By: Jeremy Worley 
*
* Modifications:
*
***********************************************************************/

int lrminorm(data,psd,dimension,order,freq,kdim,crit,matr,mati,fv1,fv2,fm1,
            eigvctr,eigvcti,eigval,auto_type)
   float *data;
   double *psd;
   double *matr,*mati,*fv1,*fv2,*fm1,*eigvctr,*eigvcti,*eigval,*crit;
   int dimension,order,freq,auto_type,kdim;
{
   int i,j,k,idx, opt, ierr;
   float minimum;
   void  ch_();
   void  lminorm(),swap();
   int ldfacorr(), ldtacorr();
   char *program = "lrminorm";

/*
** first calculate the autocorrelation of the data.
*/

   if(auto_type==2){
      if(!ldfacorr(data,dimension,order,0,0)) {
         fprintf ( stderr,"%s:  ldfacorr failed\n",program);
         return(0);
      }
   }else{
      if(!ldtacorr(data,dimension,order,0,auto_type,0)) {
         fprintf ( stderr,"%s:  ldtacorr failed\n",program);
         return(0);
      }
   }

   for(k=0;k<order;k++) {
       matr[k] = data[k];
   }

/*
** now generate the autocorrelation matrix
*/

   for(i=0;i<order;i++){
     double temp;
     idx = 0;
     for(j=i;j<order;j++){
       temp = matr[idx];
       matr[order*j+i] = temp;
       matr[order*i+j] = temp;
       mati[order*j+i] = 0.0 ; 
       mati[order*i+j] = 0.0 ; /* zero out imaginary part also */
       idx += 1;
     }
   }

/*
** get eigenvectors and values of autocorrelation matrix
*/

   opt = 1; /* used when you want both eigenvalues and eigenvectors */
   (void)ch_(&order,&order,matr,mati,eigval,&opt,eigvctr,eigvcti,
             fv1,fv2,fm1,&ierr);
   if(ierr!=0){
      fprintf(stderr,"%s:  Eispack routine ch was unable e-values\n",
              program);
      return(0);
   }

/*
** determine noise subspace dimension
*/

   if(kdim==AIC_CRIT){
      get_aic_crit(order-1,eigval,(int)(dimension - order), crit);
      minimum = crit[0];
      k = 1;
      for (i=1;i<order;i++)
        if (crit[i] < minimum)  {
          minimum = crit[i];
          k = i + 1;
        }
   }else if(kdim==MDL_CRIT){
      get_mdl_crit(order-1,eigval,(int)(dimension - order), crit);
      minimum = crit[0];
      k = 1;
      for (i=1;i<order;i++)
        if (crit[i] < minimum)  {
          minimum = crit[i];
          k = i + 1;
        }
   }else if(kdim<0){
      fprintf(stderr,"%s:  illegal noise subspace dimension type.\n",program);
      return(0);
   }else{
      k=kdim;
   }

   lminorm(order,eigvctr,eigvcti,k,freq,psd,fv1,fv2);
   /* now psd is the output */

/*
** normal return
*/
   return(1);
}

/***********************************************************************
*
*  Routine Name: lcminorm
*
*          Date:
*        
*       Purpose: Sets up matrices, etc and calls lminorm on complex data set 
*
*         Input: 
*
*        Output: 
*
*    Written By: Jeremy Worley 
*
* Modifications:
*
***********************************************************************/

int lcminorm(data,psd,dimension,order,freq,kdim,crit,matr,mati,fv1,fv2,fm1,
            eigvctr,eigvcti,eigval,auto_type)
   float *data;
   double *psd;
   double *matr,*mati,*fv1,*fv2,*fm1,*eigvctr,*eigvcti,*eigval,*crit;
   int kdim,dimension,order,freq,auto_type;
{
   int i,j,k,idx, opt, ierr;
   float minimum;
   void  ch_();
   void lminorm(),swap();
   char *program = "lcminorm";
   int ldfacorr(), ldtacorr();

/*
** first calculate the autocorrelation of the data.
*/

   if(auto_type==2){
      if(!ldfacorr(data,dimension,order,1,auto_type)) {
         fprintf ( stderr,"%s:  ldfacorr failed\n",program);
         return(0);
      }
   }else{
      if(!ldtacorr(data,dimension,order,1,auto_type,1)) {
         fprintf ( stderr,"%s:  ldtacorr failed\n",program);
         return(0);
      }
   }

   for(k=0,i=0;k<2*order;k+=2,i++) {
       matr[i] = data[k];
       mati[i] = data[k+1];
   }

/*
** now generate the autocorrelation matrix
*/

   for(i=0;i<order;i++){
     double temp;
     idx = 0;
     for(j=i;j<order;j++){
       temp = matr[idx];
       matr[order*j+i] = temp;
       matr[order*i+j] = temp;
       temp = mati[idx];
       mati[order*i+j] = temp ; 
       mati[order*j+i] = (-temp) ; /* complex matrix must be hermitian */
       idx += 1;
     }
   }

/*
** get eigenvectors and values of autocorrelation matrix
*/

   opt = 1; /* used when you want both eigenvalues and eigenvectors */
   (void)ch_(&order,&order,matr,mati,eigval,&opt,eigvctr,eigvcti,
             fv1,fv2,fm1,&ierr);
   if(ierr!=0){
      fprintf(stderr,"%s:  Eispack routine ch was unable e-values\n",
              program);
      return(0);
   }

/*
** determine noise subspace dimension
*/

   if(kdim==AIC_CRIT){
      get_aic_crit(order-1,eigval,(int)(dimension - order), crit);
      minimum = crit[0];
      k = 1;
      for (i=1;i<order;i++)
        if (crit[i] < minimum)  {
          minimum = crit[i];
          k = i + 1;
        }
   }else if(kdim==MDL_CRIT){
      get_mdl_crit(order-1,eigval,(int)(dimension - order), crit);
      minimum = crit[0];
      k = 1;
      for (i=1;i<order;i++)
        if (crit[i] < minimum)  {
          minimum = crit[i];
          k = i + 1;
        }
   }else if(kdim<0){
      fprintf(stderr,"%s:  illegal noise subspace dimension type.\n",program);
      return(0);
   }else{
      k=kdim;
   }

   lminorm(order,eigvctr,eigvcti,k,freq,psd,fv1,fv2);
   /* now psd is the output */

/*
** normal return
*/
   return(1);
}

/***********************************************************************
*
*  Routine Name: lminorm
*
*          Date:
*        
*       Purpose:  
*
*         Input: 
*
*        Output: 
*
*    Written By:  Glen Abousleman
*
* Modifications:
*
***********************************************************************/

void lminorm(order,eigvctr,eigvcti,k,nf,spectra,ar,ai)
  int order, k, nf;
  double *eigvctr,*eigvcti;
  double *spectra,*ar,*ai;
{
  int j,i,idx;
  double zreal,zimag,treal,timag, omega;

/*
** zero out data first
*/

  bzero(spectra,nf*sizeof(double));
  bzero(ar,order*sizeof(double));
  bzero(ai,order*sizeof(double));

  for(j=0;j<k;j++){
      idx = j * order; /* idx to start of jth e-vector */

/*
** get frequency response of an eigenvector
*/

      for(i=0;i<order;i++){
          ar[j] = eigvctr[i]*eigvctr[idx+i] + eigvcti[i]*eigvcti[idx+i];
          ai[j] = eigvctr[i]*eigvcti[idx+i] - eigvcti[i]*eigvctr[idx+i];
      } /* end for i */

      for(i=0;i<nf;i++){
          omega = 2*XV_PI*i/nf;
          zreal = cos(omega);
          zimag = sin(omega);
          m_eval_poly(&treal,&timag,order,&(eigvctr[idx]),&(eigvcti[idx]),
                    zreal,zimag);

          spectra[i] += treal*treal + timag*timag;
      } /* end for i */
 
  } /* end for j */

/* 
** invert the data to get inverse response
*/

  for(j=0;j<nf;j++){
      if(spectra[j] < EPSILON)
         spectra[j] = EPSILON;
      spectra[j] = 1/spectra[j];
  }
}

/***********************************************************************
*
*  Routine Name: m_eval_poly
*
*          Date:
*        
*       Purpose:  
*
*         Input: 
*
*        Output: 
*
*    Written By:  Glen Abousleman, Jeremy Worley 
*
* Modifications:
*
***********************************************************************/

static void m_eval_poly(treal,timag,order,eigvctr,eigvcti,zreal,zimag)
   int order;
   double *treal, *timag, *eigvctr, *eigvcti, zreal, zimag;
{
   int i;
   double preal,pimag,tempr,tempi;
 
   preal = eigvctr[order-1];
   pimag = eigvcti[order-1];

   for(i=order-2;i>=0;i--){
       tempr = zreal*preal - zimag*pimag;
       tempi = zreal*pimag + zimag*preal;
       preal = tempr + eigvctr[i];
       pimag = tempi + eigvcti[i];
   }

   *treal = preal;
   *timag = pimag; 
}

/***********************************************************************
*
*  Routine Name: get_aic_crit
*
*          Date:
*        
*       Purpose:  
*
*         Input: 
*
*        Output: 
*
*    Written By:  Glen Abousleman, Jeremy Worley 
*
* Modifications:
*
***********************************************************************/

static void get_aic_crit(M,lambda,N,AIC)
double *lambda;
double *AIC;
int M, N;
{
  int k;
  double L, avlog = 0, av = 0;
  for (k=1;k<=M+1;k++)  {
    av += (lambda[k-1] - av)/k;
    avlog += (log(lambda[k-1]) - avlog)/k;
    L = avlog - log(av);
    AIC[k-1] = -2*N*k*L + 2*(M+1-k)*(M+1+k);
  }
}


/***********************************************************************
*
*  Routine Name: 
*
*          Date:
*        
*       Purpose:  
*
*         Input: 
*
*        Output: 
*
*    Written By:  
*
* Modifications:
*
***********************************************************************/

static void get_mdl_crit(M,lambda,N,MDL)
double *lambda;
double *MDL;
int M, N;
{
  int k;
  double L, avlog = 0, av = 0, logN = log((double)N);
  for (k=1;k<=M+1;k++)  {
    av += (lambda[k-1] - av)/k;
    avlog += (log(lambda[k-1]) - avlog)/k;
    L = avlog - log(av);
    MDL[k-1] = -N*k*L + 0.5*(M+1-k)*(M+1+k)*logN;
  }
}

/***********************************************************************
*
*  Routine Name: lmus_center
*
*          Date:
*        
*       Purpose: center data if necessary 
*
*         Input: 
*
*        Output: 
*
*    Written By: Jeremy Worley 
*
* Modifications:
*
***********************************************************************/
  
static int lmus_center(data,num_vects,dimension)
  double **data;
  int num_vects,dimension;
{
    int v,i,dim2;
    float temp;
    
    dim2 = dimension/2;

    for (v = 0; v < num_vects; v++)
    {
      for (i = 0; i < dim2; i ++ )
      {
         temp = data[v][i];
         data[v][i] = data[v][i+dim2];
         data[v][i+dim2] = temp;
      }
    }

  return(1);
}

/* -library_code_end */
