sinfo_recipes.c

00001 /*
00002  * This file is part of the ESO SINFONI Pipeline
00003  * Copyright (C) 2004,2005 European Southern Observatory
00004  *
00005  * This program is free software; you can redistribute it and/or modify
00006  * it under the terms of the GNU General Public License as published by
00007  * the Free Software Foundation; either version 2 of the License, or
00008  * (at your option) any later version.
00009  *
00010  * This program is distributed in the hope that it will be useful,
00011  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00012  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00013  * GNU General Public License for more details.
00014  *
00015  * You should have received a copy of the GNU General Public License
00016  * along with this program; if not, write to the Free Software
00017  * Foundation, 51 Franklin St, Fifth Floor, Boston, MA  02111-1307  USA
00018  */
00019 /***************************************************************************
00020 * E.S.O. - VLT project
00021 *
00022 *
00023 *
00024 * who       when      what
00025 * --------  --------  ----------------------------------------------
00026 * schreib  05/06/00  created
00027 */
00028 
00029 #ifdef HAVE_CONFIG_H
00030 #  include <config.h>
00031 #endif
00032 
00033 #include "sinfo_vltPort.h"
00034 
00035 /*
00036  * System Headers
00037  */
00038 
00039 /*
00040  * Local Headers
00041  */
00042 
00043 #include "sinfo_recipes.h"
00044 #include "sinfo_globals.h"
00045 
00046 /*----------------------------------------------------------------------------
00047  *                                 Local variables
00048  *--------------------------------------------------------------------------*/
00049 
00050 static float  sqrarg ;
00051 
00052 static double chi1 ;                    /* old reduced chi-squared */
00053 static double sinfo_chi2 ;                    /* new reduced chi-squared */
00054 static double labda ;                   /* mixing parameter */
00055 static double vec[MAXPAR] ;             /* correction sinfo_vector */
00056 static double matrix1[MAXPAR][MAXPAR] ; /* original sinfo_matrix */
00057 static double matrix2[MAXPAR][MAXPAR] ; /* inverse of matrix1 */
00058 static int    nfree ;                   /* number of free parameters */
00059 static int    parptr[MAXPAR] ;          /* parameter pointer */
00060 
00061 /*----------------------------------------------------------------------------
00062  *                              Defines
00063  *--------------------------------------------------------------------------*/
00064 
00065 #define SQR(a) (sqrarg = (a) , sqrarg*sqrarg)
00066 
00067 /*----------------------------------------------------------------------------
00068  *                  Functions private to this module
00069  *--------------------------------------------------------------------------*/
00070 
00071 
00072 static int new_inv_mat (void) ;
00073 
00074 static void new_get_mat ( float  * xdat,
00075                      int    * xdim,
00076                      float  * ydat,
00077                      float  * wdat,
00078                      int    * ndat,
00079                      float  * fpar,
00080                      float  * epar/*,
00081                      int    * npar */) ;
00082 
00083 static int new_get_vec ( float  * xdat,
00084                     int    * xdim,
00085                     float  * ydat,
00086                     float  * wdat,
00087                     int    * ndat,
00088                     float  * fpar,
00089                     float  * epar,
00090                     int    * npar ) ;
00091 
00092 static float 
00093 new_gaussian ( float * xdat, float * parlist/*, int * npar*/ );
00094 static void 
00095 new_gaussian_deriv( float * xdat, float * parlist, 
00096                     float * dervs/*, int * npar*/ );
00097 
00098 
00099 
00107 /*----------------------------------------------------------------------------
00108  *                          Function codes
00109  *--------------------------------------------------------------------------*/
00110 
00111 
00112 float sinfo_new_f_median(float * array, int n)
00113 {
00114  pixelvalue p_array[100];
00115  int i;
00116  
00117  for (i=0;i<n;i++)
00118     p_array[i]= (pixelvalue) array[i];
00119 
00120  return (float) sinfo_new_median(p_array, n);
00121 }
00122 
00123 
00140 float sinfo_new_clean_mean( float * array, 
00141                   int     n_elements,
00142                   float   throwaway_low,
00143                   float   throwaway_high )
00144 {
00145     int i, n ;
00146     int lo_n, hi_n ;
00147     float sum ;
00148     
00149     if ( array == NULL )
00150     {
00151         sinfo_msg_error(" no array given in sinfo_clean_mean!") ;
00152         return FLT_MAX ;
00153     }
00154   
00155     if ( n_elements <= 0 )
00156     {
00157         sinfo_msg_error("wrong number of elements given") ;
00158         return FLT_MAX ;
00159     }
00160 
00161     if ( throwaway_low < 0. || throwaway_high < 0. ||
00162          throwaway_low + throwaway_high >= 100. )
00163     {
00164         sinfo_msg_error("wrong throw away percentage given!") ;
00165         return FLT_MAX ;
00166     }
00167 
00168     lo_n = (int) (throwaway_low * (float)n_elements / 100.) ;
00169     hi_n = (int) (throwaway_high * (float)n_elements / 100.) ;
00170 
00171     /* sort the array */
00172     sinfo_pixel_qsort( array, n_elements ) ;
00173 
00174     n = 0 ;
00175     sum = 0. ;
00176     for ( i = lo_n ; i < n_elements - hi_n ; i++ )
00177     {
00178         if ( !isnan(array[i]) )
00179         {
00180             sum += array[i] ;
00181             n++ ;
00182         }
00183     }
00184     if ( n == 0 )  
00185     {
00186         return FLAG ;
00187     }
00188     else
00189     {
00190         return sum/(float)n ;
00191     }
00192 }
00193 
00194 /*--------------------------------------------------------------------------*/
00207 pixelvalue sinfo_new_median(pixelvalue * array, int n)
00208 {
00209     pixelvalue med ;
00210     
00211     if ( array == NULL || n <= 0 )
00212     {
00213         sinfo_msg_warning("nothing in the pixelvalue array, ZERO returned");
00214         return ZERO ;
00215     }
00216     
00217     if ( n == 1 )
00218     {
00219         return array[0] ;
00220     }
00221  
00222     sinfo_pixel_qsort((float*) array, n) ;
00223     if ( n % 2 == 1 )
00224     {
00225         med = array[n/2] ;
00226     }
00227     else
00228     {
00229         med = (array[n/2] + array[n/2 - 1])/2. ;
00230     }
00231     return med ;
00232 }
00233 
00234 
00235 
00236 
00237 
00285 int sinfo_new_lsqfit_c ( float * xdat,
00286                int   * xdim,
00287                float * ydat,
00288                float * wdat,
00289                int   * ndat,
00290                float * fpar,
00291                float * epar,
00292                int   * mpar,
00293                int   * npar,
00294                float * tol ,
00295                int   * its ,
00296                float * lab  )
00297 {
00298     int i, n, r ;
00299     int itc ;                      /* fate of fit */
00300     int found ;                    /* fit converged: 1, not yet converged: 0 */
00301     int  nuse ;                    /* number of useable data points */
00302     double tolerance ;             /* accuracy */
00303 
00304     itc   = 0 ;                    /* fate of fit */
00305     found = 0 ;                    /* reset */
00306     nfree = 0 ;                    /* number of free parameters */
00307     nuse  = 0 ;                    /* number of legal data points */
00308 
00309     if ( *tol < (FLT_EPSILON * 10.0 ) )
00310     {
00311         tolerance = FLT_EPSILON * 10.0 ;  /* default tolerance */
00312     }
00313     else
00314     {
00315         tolerance = *tol ;                /* tolerance */
00316     }
00317     
00318     labda = fabs( *lab ) * LABFAC ;   /* start value for mixing parameter */
00319     for ( i = 0 ; i < (*npar) ; i++ )
00320     {
00321         if ( mpar[i] )
00322         {
00323             if ( nfree > MAXPAR )         /* too many free parameters */
00324             {
00325                 return -1 ;
00326             }
00327             parptr[nfree++] = i ;         /* a free parameter */
00328         }
00329     }
00330     
00331     if (nfree == 0)                       /* no free parameters */     
00332     {
00333         return -2 ;
00334     }
00335     
00336     for ( n = 0 ; n < (*ndat) ; n++ )
00337     {
00338         if ( wdat[n] > 0.0 )              /* legal weight */
00339         {
00340             nuse ++ ;
00341         }
00342     }
00343     
00344     if ( nfree >= nuse )
00345     {
00346         return -3 ;                       /* no degrees of freedom */
00347     }
00348     if ( labda == 0.0 )                   /* linear fit */
00349     {
00350         /* initialize fpar array */
00351         for ( i = 0 ; i < nfree ; fpar[parptr[i++]] = 0.0 ) ;  
00352         new_get_mat ( xdat, xdim, ydat, wdat, ndat, fpar, epar/*, npar*/ ) ;
00353         r =  new_get_vec ( xdat, xdim, ydat, wdat, ndat, fpar, epar, npar ) ;
00354         if ( r )                         /* error */
00355         {
00356             return r ;
00357         }
00358         for ( i = 0 ; i < (*npar) ; i++ )
00359         {
00360             fpar[i] = epar[i] ;           /* save new parameters */
00361             epar[i] = 0.0 ;               /* and set errors to zero */
00362         }
00363         chi1 = sqrt( chi1 / (double) (nuse - nfree) ) ;
00364         for ( i = 0 ; i < nfree ; i++ )
00365         {
00366             if ( (matrix1[i][i] <= 0.0 ) || (matrix2[i][i] <= 0.0) )
00367             {
00368                 return -7 ;
00369             }
00370             epar[parptr[i]] = chi1 * sqrt( matrix2[i][i] ) / 
00371                                      sqrt( matrix1[i][i] ) ;
00372         }
00373     }
00374     else                                  /* non-linear fit */
00375     {
00376         /*----------------------------------------------------------------
00377          * the non-linear fit uses the steepest descent method in combination
00378          * with the Taylor method. The mixing of these methods is controlled
00379          * by labda. In the outer loop ( called the iteration loop ) we build
00380          * the matrix and calculate the correction vector. In the inner loop
00381          * (called the interpolation loop) we check whether we have obtained a
00382          * better solution than the previous one. If so, we leave the inner 
00383            loop else we increase labda (give more weight to the steepest 
00384            descent method) calculate the correction vector and check again. 
00385            After the inner loop we do a final check on the goodness of the 
00386            fit and if this satisfies
00387          * the tolerance we calculate the errors of the fitted parameters.
00388          */
00389         while ( !found )                  /* iteration loop */
00390         {      
00391             if ( itc++ == (*its) )        /* increase iteration counter */
00392             {
00393                 return -4 ;               
00394             }
00395             new_get_mat( xdat, xdim, ydat, wdat, ndat, fpar, epar/*, npar */) ;
00396             
00397             /*-------------------------------------------------------------
00398              * here we decrease labda since we may assume that each iteration
00399              * brings us closer to the answer.
00400              */
00401             if ( labda > LABMIN )
00402             {
00403                 labda = labda / LABFAC ;         /* decrease labda */
00404             }
00405             r = new_get_vec( xdat, xdim, ydat, wdat, ndat, fpar, epar, npar ) ;
00406             if ( r )                      /* error */
00407             {
00408                 return r ;
00409             }
00410 
00411             while ( chi1 >= sinfo_chi2 )        /* interpolation loop */
00412             {
00413                 /*-----------------------------------------------------------
00414                  * The next statement is based on experience, not on the 
00415                    mathematics of the problem. It is assumed that we have 
00416                    reached convergence when the pure steepest descent method 
00417                    does not produce a better solution.
00418                  */
00419                 if ( labda > LABMAX )    /* assume solution found */
00420                 {
00421                     break ;
00422                 }
00423                 labda = labda * LABFAC ;     /* increase mixing parameter */
00424                 r = new_get_vec(xdat,xdim,ydat,wdat,ndat,fpar,epar,npar) ;
00425                 if ( r )                  /* error */
00426                 {
00427                     return r ;
00428                 }
00429             }
00430 
00431             if ( labda <= LABMAX )        /* save old parameters */
00432             {
00433                 for ( i = 0 ; i < *npar ; i++ )
00434                 {
00435                     fpar[i] = epar[i] ;
00436                 }
00437             }
00438             if ( (fabs( sinfo_chi2 - chi1 ) <= (tolerance * chi1)) || 
00439                  (labda > LABMAX) )
00440             {
00441                 /*-------------------------------------------------------------
00442                  * we have a satisfying solution, so now we need to calculate 
00443                    the correct errors of the fitted parameters. This we do by 
00444                    using the pure Taylor method because we are very close to 
00445                    the real solution.
00446                  */
00447                 labda = 0.0 ;              /* for Taylor solution */
00448                 new_get_mat(xdat,xdim,ydat,wdat,ndat,fpar,epar/*, npar */) ;
00449                 r = new_get_vec(xdat,xdim,ydat,wdat,ndat,fpar,epar,npar) ;
00450 
00451                 if ( r )                    /* error */
00452                 {
00453                     return r ;
00454                 }
00455                 for ( i = 0 ; i < (*npar) ; i++ )
00456                 {
00457                     epar[i] = 0.0 ;          /* set error to zero */
00458                 }
00459                 sinfo_chi2 = sqrt ( sinfo_chi2 / (double) (nuse - nfree) ) ;
00460 
00461                 for ( i = 0 ; i < nfree ; i++ )
00462                 {
00463                     if ( (matrix1[i][i] <= 0.0) || (matrix2[i][i] <= 0.0) )
00464                     {
00465                         return -7 ;
00466                     }
00467                     epar[parptr[i]] = sinfo_chi2 * sqrt( matrix2[i][i] ) / 
00468                                                    sqrt( matrix1[i][i] ) ;
00469                 }
00470                 found = 1 ;                  /* we found a solution */
00471             }
00472         }
00473     }
00474     return itc ;                             /* return number of iterations */
00475 }
00476 
00477 
00478 
00485 void sinfo_new_convert_ZEROs_to_0_for_images(cpl_image * im)
00486 {
00487     int i ;
00488     int ilx=0;
00489     int ily=0;
00490     float* pidata=NULL;
00491 
00492     if ( im == NULL )
00493     {
00494         sinfo_msg_error("no input image given!\n") ; 
00495         return ;
00496     }
00497     ilx=cpl_image_get_size_x(im);
00498     ily=cpl_image_get_size_y(im);
00499     pidata=cpl_image_get_data(im);
00500     for ( i = 0 ; i < (int) ilx*ily ; i++ )
00501     {
00502         if( isnan(pidata[i]) )
00503         {
00504             pidata[i] = 0. ;
00505         }
00506     }
00507     return ;
00508 }
00509 
00518 void sinfo_new_convert_ZEROs_to_0_for_cubes(cpl_imagelist * cube)
00519 {
00520     int i ;
00521     int inp=0;
00522     cpl_image* i_img=NULL;  
00523 
00524     if ( cube == NULL )
00525     {
00526         sinfo_msg_error("no input cube given!") ; 
00527         return ;
00528     }
00529     inp=cpl_imagelist_get_size(cube);
00530 
00531     for ( i = 0 ; i < inp ; i++ )
00532     {
00533       i_img=cpl_imagelist_get(cube,i);
00534       sinfo_new_convert_ZEROs_to_0_for_images(i_img) ;
00535       cpl_imagelist_set(cube,i_img,i);
00536     }
00537     return ;
00538 }
00539 
00540 
00549 void 
00550 sinfo_new_convert_ZEROs_to_0_for_cubes_range(cpl_imagelist * cube,
00551                                              const int z_min,const int z_max)
00552 {
00553     int i ;
00554     cpl_image* i_img=NULL;  
00555 
00556     if ( cube == NULL )
00557     {
00558         sinfo_msg_error("no input cube given!") ; 
00559         return ;
00560     }
00561     for ( i = z_min ; i < z_max ; i++ )
00562     {
00563       i_img=cpl_imagelist_get(cube,i);
00564       sinfo_new_convert_ZEROs_to_0_for_images(i_img) ;
00565       cpl_imagelist_set(cube,i_img,i);
00566     }
00567     return ;
00568 }
00575 void sinfo_new_convert_0_to_ZEROs_for_images(cpl_image * im)
00576 {
00577     int i ;
00578     int ilx=0;
00579     int ily=0;
00580     float* pidata=NULL;
00581 
00582     if ( im == NULL )
00583     {
00584         sinfo_msg_error("no input image given!") ;
00585         return ;
00586     }
00587     ilx=cpl_image_get_size_x(im);
00588     ily=cpl_image_get_size_y(im);
00589     pidata=cpl_image_get_data(im);
00590     for ( i = 0 ; i < (int) ilx*ily ; i++ )
00591     {
00592         if( pidata[i] == 0. )
00593         {
00594             pidata[i] = ZERO ;
00595         }
00596     }
00597     return ;
00598 }
00599 
00606 void sinfo_new_convert_0_to_ZERO_for_cubes(cpl_imagelist * cube)
00607 {
00608     int i ;
00609     int inp=0;
00610     cpl_image* i_img=NULL;
00611 
00612     if ( cube == NULL )
00613     {
00614         sinfo_msg_error("no input cube given!") ;
00615         return ;
00616     }
00617     inp=cpl_imagelist_get_size(cube);
00618     for ( i = 0 ; i < inp ; i++ )
00619     {
00620       i_img=cpl_imagelist_get(cube,i);
00621       sinfo_new_convert_0_to_ZEROs_for_images(i_img) ;
00622       cpl_imagelist_set(cube,i_img,i);
00623     }
00624     return ;
00625 }
00626 
00627 
00636 void 
00637 sinfo_new_convert_0_to_ZERO_for_cubes_range(cpl_imagelist * cube,
00638                                             const int z_min,const int z_max)
00639 {
00640     int i ;
00641     int inp=0;
00642     cpl_image* i_img=NULL;
00643 
00644     if ( cube == NULL )
00645     {
00646         sinfo_msg_error("no input cube given!") ;
00647         return ;
00648     }
00649     inp=cpl_imagelist_get_size(cube);
00650     for ( i = z_min ; i < z_max ; i++ )
00651     {
00652       i_img=cpl_imagelist_get(cube,i);
00653       sinfo_new_convert_0_to_ZEROs_for_images(i_img) ;
00654       cpl_imagelist_set(cube,i_img,i);
00655     }
00656     return ;
00657 }
00664 void sinfo_new_invert(cpl_image * im)
00665 {
00666     int i ;
00667     int ilx=0;
00668     int ily=0;
00669     float* pidata=NULL;
00670 
00671     ilx=cpl_image_get_size_x(im);
00672     ily=cpl_image_get_size_y(im);
00673     pidata=cpl_image_get_data(im);
00674 
00675     for ( i = 0 ; i < (int) ilx*ily ; i++ )
00676     {
00677         pidata[i] = -pidata[i] ;
00678     }
00679     return ;
00680 }
00681 
00689 int sinfo_new_nint ( double x ) 
00690 {
00691     int k ;
00692 
00693     k = x ;
00694     if ( x >= 0. )
00695     {
00696         if ( (x - (double) k) <= 0.5 )
00697         {
00698             return k ;
00699         }
00700         else
00701         {
00702             return k + 1 ;
00703         }
00704     }
00705     else
00706     {
00707         if ( (x - (double) k) <= -0.5 )
00708         {
00709             return k - 1;
00710         }
00711         else
00712         {
00713             return k ;
00714         }
00715     }
00716 }
00717 
00718 
00732 #define STEP_MIN        (-half_search)
00733 #define STEP_MAX        (half_search)
00734 
00735 double * sinfo_new_xcorrel(
00736     float      *    line_i,
00737     int             width_i,
00738     float      *    line_t,
00739     int             width_t,
00740     int             half_search,
00741     int     *       delta,
00742     int       *     maxpos,
00743     double     *    xcorr_max 
00744     
00745 )
00746 {
00747     double  * xcorr ;
00748     double   mean_i, mean_t ;
00749     double   rms_i, rms_t ;
00750     double   sum, sqsum ;
00751     double   norm ;
00752     int      nsteps ;
00753     int      i ;
00754     int      step ;
00755     int      nval ;
00756     /*double   r;*/
00757 
00758 
00759     /* Compute normalization factors */
00760     sum = sqsum = 0.00 ;
00761     for (i=0 ; i<width_i ; i++) {
00762         sum += (double)line_i[i] ;
00763         sqsum += (double)line_i[i] * (double)line_i[i];
00764     }
00765     mean_i = sum / (double)width_i ;
00766     sqsum /= (double)width_i ;
00767     rms_i = sqsum - mean_i*mean_i ;
00768 
00769     sum = sqsum = 0.00 ;
00770     for (i=0 ; i<width_t ; i++) {
00771         sum += (double)line_t[i] ;
00772         sqsum += (double)line_t[i] * (double)line_t[i];
00773     }
00774     mean_t = sum / (double)width_t ;
00775     sqsum /= (double)width_t ;
00776     rms_t = sqsum - mean_t*mean_t ;
00777 
00778     norm = 1.00 / sqrt(rms_i * rms_t);
00779 
00780     nsteps = (STEP_MAX - STEP_MIN)  ;
00781     xcorr = cpl_malloc(nsteps * sizeof(double));
00782     for (step=STEP_MIN ; step<STEP_MAX ; step++) {
00783         xcorr[step-STEP_MIN] = 0.00 ;
00784         nval = 0 ;
00785         for (i=0 ; i<width_t ; i++) {
00786             if ((i+step >= 0) &&
00787                 (i+step < width_i)) {
00788             xcorr[step-STEP_MIN] += ((double)line_t[i] - mean_t) *
00789                                     ((double)line_i[i+step] - mean_i) *
00790                                     norm ;
00791                 nval++ ;
00792             }
00793         }
00794         xcorr[step-STEP_MIN] /= (double)nval ;
00795     }
00796     *xcorr_max = xcorr[0] ;
00797     *maxpos    = 0 ;
00798     for (i=0 ; i<nsteps ; i++) {
00799         if (xcorr[i]>*xcorr_max) {
00800             *maxpos = i ;
00801             *xcorr_max = xcorr[i];
00802         }
00803     }
00804     (*delta) = + (STEP_MIN + *maxpos);
00805     return xcorr ;
00806 }
00807 
00808 /* FILE ELEMENT: sinfo_nev_ille.c                                           */
00809 /*                                                                    */
00810 /**********************************************************************/
00811 /*                                                                    */
00812 /*                      double sinfo_nev_ille()                             */
00813 /*                                                                    */
00814 /**********************************************************************/
00815 /*                                                                    */
00816 /*  DESCRIPTION:                                                      */
00817 /*  For a given table (x , f(x )), i = 0(1)n  and a given argument z  */
00818 /*  the function computes the interpolated value for the argument z   */
00819 /*  using Neville's interpolation/ extrapolation algorithm.           */
00820 /*                                                                    */
00821 /*  FUNCTIONS CALLED:                                                 */
00822 /*  System library: <stdio.h> printf(), fabs();                       */
00823 /*  Numlib library: None                                              */
00824 /*  Local functions: nevtable();                                      */
00825 /*  User supplied: None                                               */
00826 /*                                                                    */
00827 /*  PROGRAMMED BY: T.Haavie                                           */
00828 /*  DATE/VERSION: 88-07-06/1.0                                        */
00829 /*                                                                    */
00830 /**********************************************************************/
00831 
00832 float sinfo_new_nev_ille(float x[], float f[], int n, float z, int* flag)
00833                /* PARAMETERS(input):                                  */
00834 /* float x[];     Abscissae values in the table.                      */
00835 /* float f[];     Function values in the table.                       */
00836 /* int n;         The number of elements in the table is n+1.         */
00837 /* float z;       Argument to be used in interpolation/extrapolation. */
00838 
00839 
00840 /*                PARAMETERS(input/output):                           */
00841 /* int *flag;    Flag parameter(output):                             */
00842                /* = 0, n < 0 and/or eps < 0, should be positive+.     */
00843                /* = 1, required rel.err. is not obtained.             */
00844                /* = 2, required rel. err. is obtained.                */
00845 
00846 /* the computed estimate for the interpolated/extrapolated  value re- */
00847 /* turned through function name sinfo_nev_ille.                             */
00848 
00849 {
00850         float p[11]; /* Array used for storing the new row elements */
00851                        /* in the interpolation/extrapolation table.   */
00852         float q[11]; /* Array used for storing the old row elements */
00853                        /* in the interpolation/extrapolation table    */
00854 
00855         float factor;
00856         int m, k;
00857 
00858        
00859 
00860         if (n < 0 )
00861         {
00862                 *flag = 0;
00863                 return(0.);
00864         }
00865 
00866 
00867         q[0] = f[0];               /* Set initial value in the table. */
00868 
00869         for (k = 1; k <= n; k++)   /* k counts rows in the table.     */
00870         {
00871                 p[0] = f[k];
00872                 for (m = 1; m <= k; m++) /* m counts element in row.  */
00873                 {
00874                         factor = (z - x[k]) / (x[k] - x[k-m]);
00875                         p[m] = p[m-1] + factor * (p[m-1] - q[m-1]);
00876                 }
00877 
00878 
00879                 for (m = 0; m <= k; m++) /* Shift old row to new row.  */
00880                         q[m] = p[m];
00881 
00882         } /* End of k-loop. */
00883 
00884         *flag = 1;              /* Required rel.error is not obtained. */
00885         return(p[n]);
00886 
00887 } /* End of sinfo_nev_ille(). */
00888 
00889 
00914 static int new_get_vec ( float * xdat,
00915                     int   * xdim,
00916                     float * ydat,
00917                     float * wdat,
00918                     int   * ndat,
00919                     float * fpar,
00920                     float * epar,
00921                     int   * npar )
00922 {
00923     double dj ;
00924     double dy ;
00925     double mii ;
00926     double mji ;
00927     double mjj ;
00928     double wn ;
00929     int i, j, n, r ;
00930 
00931     /* loop to modify and scale the sinfo_matrix */
00932     for ( j = 0 ; j < nfree ; j++ )
00933     {
00934         mjj = matrix1[j][j] ;
00935         if ( mjj <= 0.0 )             /* diagonal element wrong */
00936         {
00937             return -5 ;
00938         }
00939         mjj = sqrt( mjj ) ;
00940         for ( i = 0 ; i < j ; i++ )
00941         {
00942             mji = matrix1[j][i] / mjj / sqrt( matrix1[i][i] ) ;
00943             matrix2[i][j] = matrix2[j][i] = mji ;
00944         }
00945         matrix2[j][j] = 1.0 + labda ;  /* scaled value on diagonal */
00946     }    
00947     
00948     if ( (r = new_inv_mat()) )       /* sinfo_invert sinfo_matrix inlace */
00949     {
00950         return r ;
00951     }
00952     
00953     for ( i = 0 ; i < (*npar) ; i ++ )
00954     {
00955         epar[i] = fpar[i] ;
00956     }
00957     
00958     /* loop to calculate correction sinfo_vector */
00959     for ( j = 0 ; j < nfree ; j++ )
00960     {
00961         dj = 0.0 ;
00962         mjj = matrix1[j][j] ;
00963         if ( mjj <= 0.0)               /* not allowed */
00964         {
00965             return -7 ;
00966         }
00967         mjj = sqrt ( mjj ) ;
00968         for ( i = 0 ; i < nfree ; i++ )
00969         {
00970             mii = matrix1[i][i] ;
00971             if ( mii <= 0.0 )
00972             {
00973                 return -7 ;
00974             }
00975             mii = sqrt( mii ) ;
00976             dj += vec[i] * matrix2[j][i] / mjj / mii ;
00977         }
00978         epar[parptr[j]] += dj ;       /* new parameters */
00979     }    
00980     chi1 = 0.0 ;                      /* reset reduced chi-squared */
00981  
00982     /* loop through the data points */
00983     for ( n = 0 ; n < (*ndat) ; n++ )
00984     {
00985         wn = wdat[n] ;               /* get weight */
00986         if ( wn > 0.0 )              /* legal weight */
00987         {
00988             dy = ydat[n] - new_gaussian( &xdat[(*xdim) * n], epar/*, npar*/ ) ;
00989             chi1 += wdat[n] * dy * dy ;
00990         }
00991     }
00992     return 0 ;
00993 }   
00994     
00995 
01011 static void new_get_mat ( float * xdat,
01012                      int   * xdim,
01013                      float * ydat,
01014                      float * wdat,
01015                      int   * ndat,
01016                      float * fpar,
01017                      float * epar/*,
01018                      int   * npar */)
01019 {
01020     double wd ;
01021     double wn ;
01022     double yd ;
01023     int i, j, n ;
01024 
01025     for ( j = 0 ; j < nfree ; j++ )
01026     {
01027         vec[j] = 0.0 ; /* zero sinfo_vector */
01028         for ( i = 0 ; i<= j ; i++ )   
01029         /* zero sinfo_matrix only on and below diagonal */
01030         {
01031             matrix1[j][i] = 0.0 ;
01032         }
01033     }
01034     sinfo_chi2 = 0.0 ;  /* reset reduced chi-squared */
01035     
01036     /* loop through data points */
01037     for ( n = 0 ; n < (*ndat) ; n++ )
01038     {
01039         wn = wdat[n] ;
01040         if ( wn > 0.0 )  /* legal weight ? */
01041         {
01042             yd = ydat[n] - new_gaussian( &xdat[(*xdim) * n], fpar/*, npar*/ ) ;
01043             new_gaussian_deriv( &xdat[(*xdim) * n], fpar, epar/*, npar*/ ) ;
01044             sinfo_chi2 += yd * yd * wn ; /* add to chi-squared */
01045             for ( j = 0 ; j < nfree ; j++ )
01046             {
01047                 wd = epar[parptr[j]] * wn ;  /* weighted derivative */
01048                 vec[j] += yd * wd ;       /* fill sinfo_vector */
01049                 for ( i = 0 ; i <= j ; i++ ) /* fill sinfo_matrix */
01050                 {
01051                     matrix1[j][i] += epar[parptr[i]] * wd ;
01052                 }
01053             }
01054         }
01055     }                   
01056 }  
01057    
01058 
01059 
01060 
01061 
01062  
01071 static int new_inv_mat (void)
01072 {
01073     double even ;
01074     double hv[MAXPAR] ;
01075     double mjk ;
01076     double rowmax ;
01077     int evin ;
01078     int i, j, k, row ;
01079     int per[MAXPAR] ;
01080    
01081     /* set permutation array */
01082     for ( i = 0 ; i < nfree ; i++ )
01083     {
01084         per[i] = i ;
01085     }
01086     
01087     for ( j = 0 ; j < nfree ; j++ ) /* in j-th column */
01088     {
01089         /* determine largest element of a row */                               
01090         rowmax = fabs ( matrix2[j][j] ) ;
01091         row = j ;                         
01092 
01093         for ( i = j + 1 ; i < nfree ; i++ )
01094         {
01095             if ( fabs ( matrix2[i][j] ) > rowmax )
01096             {
01097                 rowmax = fabs( matrix2[i][j] ) ;
01098                 row = i ;
01099             }
01100         }
01101 
01102         /* determinant is zero! */
01103         if ( matrix2[row][j] == 0.0 )
01104         {
01105             return -6 ;
01106         }
01107         
01108         /*if the largest element is not on the diagonal, then permutate rows */
01109         if ( row > j )
01110         {
01111             for ( k = 0 ; k < nfree ; k++ )
01112             {
01113                 even = matrix2[j][k] ;
01114                 matrix2[j][k] = matrix2[row][k] ;
01115                 matrix2[row][k] = even ;
01116             }
01117             /* keep track of permutation */
01118             evin = per[j] ;
01119             per[j] = per[row] ;
01120             per[row] = evin ;
01121         }
01122         
01123         /* modify column */
01124         even = 1.0 / matrix2[j][j] ;
01125         for ( i = 0 ; i < nfree ; i++ )
01126         {
01127             matrix2[i][j] *= even ;
01128         }
01129         matrix2[j][j] = even ;
01130         
01131         for ( k = 0 ; k < j ; k++ )
01132         {
01133             mjk = matrix2[j][k] ;
01134             for ( i = 0 ; i < j ; i++ )
01135             {
01136                 matrix2[i][k] -= matrix2[i][j] * mjk ;
01137             }
01138             for ( i = j + 1 ; i < nfree ; i++ )
01139             {
01140                 matrix2[i][k] -= matrix2[i][j] * mjk ;
01141             }
01142             matrix2[j][k] = -even * mjk ;
01143         }
01144     
01145         for ( k = j + 1 ; k < nfree ; k++ )
01146         {
01147             mjk = matrix2[j][k] ;
01148             for ( i = 0 ; i < j ; i++ )
01149             {
01150                 matrix2[i][k]  -= matrix2[i][j] * mjk ;
01151             }
01152             for ( i = j + 1 ; i < nfree ; i++ )
01153             {
01154                 matrix2[i][k]  -= matrix2[i][j] * mjk ;
01155             }
01156             matrix2[j][k] = -even * mjk ;
01157         }
01158     }
01159     
01160     /* finally, repermute the columns */
01161     for ( i = 0 ; i < nfree ; i++ )
01162     {
01163         for ( k = 0 ; k < nfree ; k++ )
01164         {
01165             hv[per[k]] = matrix2[i][k] ;
01166         }
01167         for ( k = 0 ; k < nfree ; k++ )
01168         {
01169             matrix2[i][k] = hv[k] ;
01170         }
01171     }
01172         
01173     /* all is well */
01174     return 0 ;
01175 }       
01176 
01177 
01178 
01179 
01200 float new_gaussian ( float * xdat, float * parlist/*, int * npar*/ )
01201 {
01202     double  xd ;  /* FWHM's of gauss function */
01203     double   x ;  /* position */
01204 
01205     xd = fabs((double) parlist[1]) ;
01206     x  = (double) xdat[0] - (double) parlist[2] ;
01207     return (float) (
01208            (double) parlist[0] * exp( -4.0 * log(2.0) * (x/xd) * (x/xd) )
01209            + (double) parlist[3] ) ;
01210 }
01211       
01212        
01237 void 
01238 new_gaussian_deriv(float * xdat,float * parlist,float * dervs/*, int * npar*/ )
01239 {
01240     double xd ; /* FWHM of sinfo_gaussian */
01241     double x, expon ; /* position and exponent */
01242 
01243     xd = fabs( (double) parlist[1] ) ;
01244     
01245     /* offset from peak position */
01246     x = (double) xdat[0] - (double) parlist[2] ;
01247 
01248     /* determine the derivatives: */
01249     expon = -4.0 * log(2.0) * (x/xd) * (x/xd) ;
01250     expon = exp( expon ) ;
01251 
01252     /* partial derivative by the amplitude */
01253     dervs[0] = (float) expon ;
01254 
01255     /* calculate a * exp(-arg) */
01256     expon = (double) parlist[0] * expon ;
01257 
01258     /* partial derivative by FWHM */
01259     dervs[1] = (float) ( expon * 8.0 * log(2.0) * x*x / (xd*xd*xd) ) ;
01260 
01261     /* partial derivative by the x position (parlist[2]) */
01262     dervs[2] = (float) (expon * 8.0 * log(2.0) * x/(xd*xd) ) ;
01263 
01264     /* partial derivative by the zero level */
01265     dervs[3] = 1.0 ;
01266 }
01267 
01268 
01269 /*==================================================================*/
01270 
01271 
01291 void 
01292 sinfo_my_fit(float x[], float y[], int ndata, float sig[], int mwt, float *a, 
01293            float *b, float *siga, float *sigb, float *chi2, float *q)
01294 {
01295     int i ;
01296     float wt, t, sxoss, sx=0., sy=0., st2=0., ss, sigdat ;
01297 
01298     *b = 0. ;             /*accumulate sums ...*/
01299     if ( mwt )
01300     {
01301         ss = 0. ;
01302         for ( i = 0 ; i < ndata ; i++ )  /*... with weights*/
01303         {
01304             wt = 1./SQR(sig[i]) ;
01305             ss += wt ;
01306             sx += x[i]*wt ;
01307             sy += y[i]*wt ;
01308         }
01309     }
01310     else
01311     {
01312         for ( i = 0 ; i < ndata ; i++ ) /*... or without weights*/
01313         {
01314              sx += x[i] ;
01315              sy += y[i] ;
01316         }
01317         ss = ndata ;
01318     }
01319     sxoss = sx/ss ;
01320              
01321     if ( mwt )
01322     {
01323         for ( i = 0 ; i < ndata ; i ++ )
01324         {
01325             t = (x[i] - sxoss)/sig[i] ;
01326             st2 += t*t ;
01327             *b += t*y[i]/sig[i] ;
01328         }
01329     }
01330     else
01331     {
01332         for ( i = 0 ; i < ndata ; i++ )
01333         {
01334             t = x[i] - sxoss ;
01335             st2 += t*t ;
01336             *b += t*y[i] ;           
01337         }
01338     }
01339 
01340     *b /= st2 ;
01341     *a = (sy - sx*(*b))/ss ;
01342     *siga = sqrt ((1.0 + sx*sx/(ss*st2))/ss) ;
01343     *sigb = sqrt (1.0/st2) ;
01344     *chi2 = 0.0 ;  /*calculate chi-square*/
01345     if ( mwt == 0 )
01346     {
01347         for ( i = 0 ; i < ndata ; i++ )
01348         {
01349             *chi2 += SQR (y[i] - (*a) - (*b)*x[i]) ;
01350         }
01351         *q = 1. ;
01352         
01353         /*------------------------------------------------------------------
01354          * for unweighted data evaluate typical sig using chi2, and adjust
01355          * the standard deviation
01356          */
01357         sigdat = sqrt ((*chi2)/(ndata - 2)) ;
01358         *siga *= sigdat ;
01359         *sigb *= sigdat ;
01360     }
01361     else
01362     {
01363         for (i = 0 ; i < ndata ; i++)
01364         {
01365             *chi2 += SQR ((y[i] - (*a) - (*b) * x[i])/sig[i]) ;
01366         }    
01367         *q = 1. ; /* delete rest of lines. q is not a good value */
01368     }
01369 }
01370 
01385 int sinfo_new_correlation ( float * data1, float * data2, int ndata )
01386 {
01387     /*float help[3*ndata] ; 
01388     float corsum[3*ndata] ;*/
01389     float* help=NULL ; 
01390     float* corsum=NULL ;
01391     float maxval ;
01392     int i, j, k, position, shift ;
01393     int /*start,end,size,*/ndata3,limit;
01394     
01395     
01396     /*ndata3=3*ndata;*/
01397     ndata3=ndata+300;
01398 
01399     if ( NULL == data1 || NULL == data2 || ndata <= 1 )
01400     {
01401         sinfo_msg_error(" wrong input for sinfo_correlation\n") ;
01402         return INT32_MAX ;
01403     }
01404 
01405     /* initialize the help arrays with zeros */
01406     help=cpl_calloc(ndata+300,sizeof(float));
01407     for ( i = 0 ; i < ndata3 ; i++ )
01408     {
01409         help[i] = 0. ;
01410     }
01411 
01412     /* shift the second data array by ndata in the help array */
01413     for ( i = 0 ; i < ndata ; i++ )
01414     {
01415         help[(300/2) + i] = data2[i] ;
01416     }
01417 
01418     /* compute the cross sinfo_correlation sum array */
01419     corsum=cpl_calloc(ndata+300,sizeof(float));
01420     for ( j = 0 ; j < ndata3 ; j++ )
01421     {
01422         if ( ndata3-j <= ndata) 
01423         limit = ndata3-j;
01424     else
01425         limit = ndata;
01426         corsum[j] = 0. ;
01427         for ( k = 0 ; k < limit ; k++ )
01428         {
01429             /*if ( k + j >= ndata3 )
01430             {
01431                 break ;
01432             }*/
01433             corsum[j] += data1[k] * help[k + j] ;
01434         }
01435     }
01436 
01437     /* search for the maximal corsum value and determine its position */
01438     maxval = -FLT_MAX ;
01439     position = -1 ;
01440     for ( i = 0 ; i < ndata3 ; i++ )
01441     {
01442         if ( maxval < corsum[i] )
01443         {
01444             maxval = corsum[i] ;
01445             position = i ;
01446         }
01447     }
01448     
01449     /* determine shift of data2 relative to the data1 array */
01450     shift = position - 300/2 ;
01451     cpl_free(help);
01452     cpl_free(corsum);
01453  
01454     return shift ;
01455 }
01456 
01457 /*--------------------------------------------------------------------------*/

Generated on Wed Jan 17 08:33:44 2007 for SINFONI Pipeline Reference Manual by  doxygen 1.4.4