00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028 #ifdef HAVE_CONFIG_H
00029 #include <config.h>
00030 #endif
00031
00032
00033
00034
00035
00036 #include <math.h>
00037 #include <cpl.h>
00038
00039 #include "irplib_fit.h"
00040
00041
00045
00046
00049
00050
00051
00052
00053
00054 #define IRPLIB_FIT_AMOEBA_NMAX 5000
00055
00056
00057
00058
00059
00060
00061 #define IRPLIB_FIT_AMOEBA_SWAP(a,b) {swap=(a);(a)=(b);(b)=swap;}
00062
00063
00064
00065
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
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
00169
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
00184
00185
00186 if(rtol < ftol) {
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
00205
00206
00207
00208 ytry=irplib_fit_amotry(p,y,psum,ndim,funk,ihi,-1.0);
00209 if(ytry <= y[ilo]) {
00210
00211
00212
00213
00214 ytry=irplib_fit_amotry(p,y,psum,ndim,funk,ihi,2.0);
00215 } else if (ytry >= y[inhi]) {
00216
00217
00218
00219
00220
00221 ysave=y[ihi];
00222 ytry=irplib_fit_amotry(p,y,psum,ndim,funk,ihi,0.5);
00223 if(ytry >= ysave) {
00224
00225
00226
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;
00237 irplib_fit_amoeba_get_psum(ndim,mpts,p,psum);
00238 }
00239 } else {
00240 --(*nfunk);
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