irplib_fit.c

00001 /* $Id: irplib_fit.c,v 1.12 2006/11/21 07:54:07 amodigli Exp $
00002  *
00003  * This file is part of the irplib package
00004  * Copyright (C) 2002,2003 European Southern Observatory
00005  *
00006  * This program is free software; you can redistribute it and/or modify
00007  * it under the terms of the GNU General Public License as published by
00008  * the Free Software Foundation; either version 2 of the License, or
00009  * (at your option) any later version.
00010  *
00011  * This program is distributed in the hope that it will be useful,
00012  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00013  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00014  * GNU General Public License for more details.
00015  *
00016  * You should have received a copy of the GNU General Public License
00017  * along with this program; if not, write to the Free Software
00018  * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02111-1307  USA
00019  */
00020 
00021 /*
00022  * $Author: amodigli $
00023  * $Date: 2006/11/21 07:54:07 $
00024  * $Revision: 1.12 $
00025  * $Name:  $
00026  */
00027 
00028 #ifdef HAVE_CONFIG_H
00029 #include <config.h>
00030 #endif
00031 
00032 /*-----------------------------------------------------------------------------
00033                                    Includes
00034  -----------------------------------------------------------------------------*/
00035 
00036 #include <math.h>
00037 #include <cpl.h>
00038 
00039 #include "irplib_fit.h"
00040 
00041 /*---------------------------------------------------------------------------*/
00045 /*---------------------------------------------------------------------------*/
00046 
00049 /*---------------------------------------------------------------------------
00050    Defines
00051    --------------------------------------------------------------------------*/
00052 
00053 
00054 #define IRPLIB_FIT_AMOEBA_NMAX 5000
00055 
00056 
00057 /*---------------------------------------------------------------------------
00058    Macros
00059    --------------------------------------------------------------------------*/
00060 
00061 #define IRPLIB_FIT_AMOEBA_SWAP(a,b) {swap=(a);(a)=(b);(b)=swap;}
00062 
00063 
00064 /*---------------------------------------------------------------------------
00065    Function prototypes
00066    --------------------------------------------------------------------------*/
00067 
00068 static inline void 
00069 irplib_fit_amoeba_get_psum(int ndim, int mpts, double** p, double* psum);
00070 
00071 static  double 
00072 irplib_fit_amotry(double** p, 
00073              double y[], 
00074              double psum[], 
00075              int ndim, 
00076              double (*funk)(double[]),
00077              int ihi, 
00078              double fac);
00079 
00080 /*---------------------------------------------------------------------------
00081    Function implementation
00082    --------------------------------------------------------------------------*/
00083 
00095 static inline void 
00096 irplib_fit_amoeba_get_psum(int ndim, int mpts, double** p, double* psum)
00097 {
00098   int i=0;
00099   int j=0;
00100   double sum=0;
00101   for (j=0;j<ndim;j++) {
00102     for (sum=0.0,i=0;i<mpts;i++) {
00103       sum += p[i][j];
00104     }
00105     psum[j]=sum;
00106   }
00107 
00108 } 
00109 
00110 
00135 void 
00136 irplib_fit_amoeba(
00137          double**p, 
00138          double y[], 
00139          int ndim,
00140          double ftol,
00141          double (*funk)(double[]),
00142          int* nfunk)
00143 {
00144 
00145 
00146   int i=0;
00147   int ihi=0;
00148   int ilo=0;
00149   int inhi=0;
00150   int j=0;
00151   int mpts=ndim+1;
00152   double rtol=0;
00153   double swap=0;
00154   double ysave=0;
00155   double ytry=0;
00156   cpl_vector* sum=NULL;
00157   double* psum=NULL;
00158 
00159   sum=cpl_vector_new(ndim);
00160   psum=cpl_vector_get_data(sum);
00161   *nfunk=0;
00162 
00163   irplib_fit_amoeba_get_psum(ndim,mpts,p,psum);
00164 
00165   for(;;) {
00166     ilo=0;
00167     /*
00168     First we must determine which point is the highest (worst), 
00169     next-highest, and lowest (best), by looping over the points in the simplex
00170     */
00171     ihi=y[0]>y[1] ? (inhi=1,0) : (inhi=0,1);
00172 
00173     for (i=0;i< mpts;i++) {
00174       if (y[i] <= y[ilo]) ilo=i;
00175       if (y[i] > y[ihi]) {
00176     inhi=ihi;
00177     ihi=i;
00178       } else if (y[i] > y[inhi] && i != ihi) inhi=i;
00179     }
00180     rtol=2.0*fabs(y[ihi]-y[ilo])/(fabs(y[ihi])+fabs(y[ilo]));
00181 
00182     /*
00183     compute the fractional range from highest to lowest and return if 
00184     satisfactory
00185     */
00186     if(rtol < ftol) { // if returning, but best point and value is in slot 1
00187       IRPLIB_FIT_AMOEBA_SWAP(y[0],y[ilo]);
00188       for (i=0;i<ndim;i++) IRPLIB_FIT_AMOEBA_SWAP(p[1][i],p[ilo][i]);
00189       break;
00190     }
00191     if (*nfunk >= IRPLIB_FIT_AMOEBA_NMAX) {
00192          cpl_msg_error(cpl_func,"NMAX exceeded");
00193          IRPLIB_FIT_AMOEBA_SWAP(y[0],y[ilo]);
00194          for (i=0;i<ndim;i++) IRPLIB_FIT_AMOEBA_SWAP(p[1][i],p[ilo][i]);
00195      for (i=0;i<ndim;i++) {
00196            cpl_msg_info(cpl_func,"p[1][i]=%g p[ilo][i]=%g ilo=%d",
00197               p[1][i],p[ilo][i],ilo);
00198      }
00199          break;
00200 
00201     }
00202     *nfunk +=2;
00203     /*
00204     Begin a new iteration. First extrapolate by a Factor -1 through the face
00205     of the simplex across the high point, i.e. reflect the simplex from the 
00206     high point
00207     */
00208     ytry=irplib_fit_amotry(p,y,psum,ndim,funk,ihi,-1.0);
00209     if(ytry <= y[ilo]) {
00210       /* 
00211       Gives a result better than the best point, so try an additional 
00212       extrapolation by a factor 2 
00213       */
00214       ytry=irplib_fit_amotry(p,y,psum,ndim,funk,ihi,2.0);
00215     } else if (ytry >= y[inhi]) {
00216 
00217       /*
00218       The reflected point is worse than the second highest, so look for an 
00219       intermediate lower point, i.e. do a one-dimensional contraction
00220       */
00221       ysave=y[ihi];
00222       ytry=irplib_fit_amotry(p,y,psum,ndim,funk,ihi,0.5);
00223       if(ytry >= ysave) { 
00224        /* 
00225           Can't seem to get rid of that high point.
00226           Better contract around the lowest (best) point
00227        */
00228     for(i=0;i<mpts;i++) {
00229       if(i != ilo) {
00230         for( j=0;j<ndim;j++) {
00231           p[i][j]=psum[j]=0.5*(p[i][j]+p[ilo][j]);
00232         }
00233         y[i]=(*funk)(psum);
00234       }
00235     }
00236     *nfunk += ndim; /* Keep track of function evaluations */
00237         irplib_fit_amoeba_get_psum(ndim,mpts,p,psum);/* Recomputes psum */
00238       }
00239     } else {
00240       --(*nfunk); /* Go back for the test of doneness and the next iteration */
00241     }
00242   }
00243 
00244   cpl_vector_delete(sum);
00245 }
00246 
00247 
00263 static  double 
00264 irplib_fit_amotry(double** p, double y[], double psum[], int ndim, 
00265              double (*funk)(double[]),int ihi, double fac)
00266 {
00267   int j;
00268   double fac1=0;
00269   double fac2=0;
00270   double ytry=0;
00271   cpl_vector * vtry=NULL;
00272   double *ptry=NULL;
00273 
00274   vtry=cpl_vector_new(ndim);
00275   ptry=cpl_vector_get_data(vtry);
00276 
00277   fac1=(1.0-fac)/ndim;
00278   fac2=fac1-fac;
00279         
00280   for (j=0;j<ndim;j++) {
00281     ptry[j]=psum[j]*fac1-p[ihi][j]*fac2;
00282   }
00283   ytry=(*funk)(ptry);
00284   if (ytry < y[ihi]) {
00285     y[ihi]=ytry;
00286     for (j=0;j<ndim;j++) {
00287       psum[j] += ptry[j]-p[ihi][j];
00288       p[ihi][j]=ptry[j];
00289     }
00290   }
00291   cpl_vector_delete(vtry);
00292   return ytry;
00293 }
00294 
00295 
00296 

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