Main Page | Modules | Data Structures | File List | Data Fields | Globals | Related Pages

extract/lib/fit_max.c

Go to the documentation of this file.
00001 /* === Doxygen Comment ======================================= */
00013 /* =========================================================== */
00014 
00015 #include <IFU_io.h>
00016 #include <IFU_math.h>
00017 #include <snifs.h>
00018 #include <extract.h>
00019 #include <gsl/gsl_sf_erf.h>
00020 
00021 /* Global variables describing the geometric PSF */
00022 int glnggeo;               /* Nb of gaussian components */
00023 double *glIgeo, *glXgeo, *glSgeo;
00024 
00025 /* Global variables describing the global PSF */
00026 int glngpsf;               /* Nb of gaussian components */
00027 
00028 /* Global variables describing the core PSF = geom. PSF (x) global PSF */
00029 double **glIpsf, **glSpsf; /* Global var. needed by nllsqfit_globalPSF */
00030 
00031 /* Global variables describing the local PSF */
00032 int glngloc;             /* Nb of gaussian components */
00033 
00034 /* Variable used for fit_maxpack normalization */
00035 double *glipeaknorm, *glxpeakoffs;
00036 
00037 /* find_max: fit_maxpack */
00038 #define FMX_XPEAKNORM 0.3
00039 #define FMX_XPEAKOFFS(k) (glxpeakoffs[k])
00040 #define FMX_IPEAKNORM(k) (glipeaknorm[k])
00041 #define FMX_SPEAKNORM 0.1
00042 
00043 /* === Doxygen Comment ======================================= */
00058 /* =========================================================== */
00059 
00060 long nllsqfit_maxpack(long *mode, long *npts, long *npar, long *ldfj, 
00061                       double par[], double fvec[], double fjac[], long *nstate, 
00062                       long *iuser, double x_n_sig[])
00063 {
00064   long j, k, g, p;
00065   int ngeo, npsf, npeak;
00066   double qe, qs, Ge, Gs, xe, xs;
00067   double profil, *sig, **Xcs, ***Scs, ***Ics;
00068 
00069 #define Xpeak(i) (par[3*(i)  ]*FMX_XPEAKNORM + FMX_XPEAKOFFS(i))
00070 #define Ipeak(i) (par[3*(i)+1]*FMX_IPEAKNORM(i))
00071 #define Speak(i) (par[3*(i)+2]*FMX_SPEAKNORM)
00072 
00073   ngeo = glnggeo;
00074   npsf = MAX(1,glngpsf); /* Handle the case: glngpsf = 0 */
00075   npeak = (*npar)/3;     /* npar = npeak*(X,I,S) */
00076 
00077   sig  = x_n_sig + (*npts); /* sig in the 2nd half of x_n_sig */
00078 
00079   Xcs = (double  **)malloc(ngeo*sizeof(double *));
00080   Scs = (double ***)malloc(ngeo*sizeof(double **));
00081   Ics = (double ***)malloc(ngeo*sizeof(double **));
00082   for (g=0; g<ngeo; g++) {
00083     Xcs[g] = (double  *)malloc(npeak*sizeof(double));
00084     Scs[g] = (double **)malloc( npsf*sizeof(double *)); 
00085     Ics[g] = (double **)malloc( npsf*sizeof(double *)); 
00086     for (p=0; p<npsf; p++) {
00087       Scs[g][p] = (double *)malloc(npeak*sizeof(double));       
00088       Ics[g][p] = (double *)malloc(npeak*sizeof(double));       
00089     }
00090   }
00091 
00092   /* Convolution of the core PSF (ngeo*npsf gaussians) with the
00093      up-to-date local PSF model (1 gaussian) */
00094 
00095   for (k=0; k<npeak; k++) {    /* Loop over the pack peaks */
00096     for (g=0; g<ngeo; g++) {   /* Loop on the geometric gaussians */
00097       Xcs[g][k] = Xpeak(k) - glXgeo[g];
00098       for (p=0; p<npsf; p++) { /* Loop on PSF gaussians */
00099         Scs[g][p][k] = hypot(glSpsf[g][p],Speak(k));
00100         Ics[g][p][k] = M_SQRT2PI*Ipeak(k)*glIpsf[g][p]*
00101           Speak(k)*glSpsf[g][p]/Scs[g][p][k];
00102       }
00103     }
00104   }
00105 
00106   xs =-0.5; /* pixel boundaries */
00107   xe = 0.5;
00108 
00109   for (j=0; j<*npts; j++) {
00110     if (*mode != 1) fvec[j] = 0.;
00111     if (*mode != 0) for (p=0; p<*npar; p++) fjac[j+p*(*ldfj)] = 0.;
00112 
00113     for (k=0; k<npeak; k++) {    /* Loop over the pack peaks */
00114       for (g=0; g<ngeo; g++) {   /* Loop on the geom gaussians */
00115         for (p=0; p<npsf; p++) { /* Loop on the PSF gaussians */
00116           qe = (xe - Xcs[g][k])/Scs[g][p][k]/M_SQRT2;
00117           qs = (xs - Xcs[g][k])/Scs[g][p][k]/M_SQRT2;
00118           Gs = Ics[g][p][k]*exp(-SQ(qs));
00119           Ge = Ics[g][p][k]*exp(-SQ(qe));
00120           profil = sqrt(M_PI_2)*Ics[g][p][k]*Scs[g][p][k]* 
00121             (gsl_sf_erf(qe) - gsl_sf_erf(qs));
00122           if (*mode != 1) fvec[j] += profil;
00123           if (*mode == 0) continue;
00124           /* df/dXpeak */
00125           fjac[j+(3*k  )*(*ldfj)] -= (Ge - Gs)       * FMX_XPEAKNORM;
00126           /* df/dIpeak */
00127           fjac[j+(3*k+1)*(*ldfj)] += profil/Ipeak(k) * FMX_IPEAKNORM(k);
00128           /* df/dSpeak */
00129           fjac[j+(3*k+2)*(*ldfj)] +=                   FMX_SPEAKNORM *
00130             ( (SQ(glSpsf[g][p]/Speak(k)) + 1.)*profil -
00131               ((xe-Xcs[g][k])*Ge - (xs-Xcs[g][k])*Gs) ) * 
00132             Speak(k)/SQ(Scs[g][p][k]);
00133         }
00134       }
00135     }
00136     /* Normalization */
00137     fvec[j] /= sig[j];
00138     for (p=0; p<*npar; p++) fjac[j+p*(*ldfj)] /= sig[j];
00139 
00140     xs += 1.;
00141     xe += 1.;
00142   }
00143   for (g=0; g<ngeo; g++) {
00144     for (p=0; p<npsf; p++) { 
00145       free(Scs[g][p]); free(Ics[g][p]); 
00146     }
00147     free(Scs[g]); free(Ics[g]); free(Xcs[g]);
00148   }
00149   free(Scs); free(Ics); free(Xcs);
00150   return(OK);
00151 
00152 #undef Xpeak
00153 #undef Ipeak
00154 #undef Speak
00155 }
00156 
00157 /* Global variable used for fit_maxpack normalization */
00158 double *glipeaknorm, *glxpeakoffs;
00159 
00160 /* === Doxygen Comment ======================================= */
00176 /* =========================================================== */
00177 
00178 double fit_maxpack(double xarr[], double sarr[], int npts, int npeak, 
00179                    double packpar[], double dpackpar[], double xfit[])
00180 {
00181   int k, status;
00182   double *bl, *bu;
00183   double hw, chi2, gof, rms;
00184 
00185   if (glngloc > 1) {
00186     print_error("Multi-gaussian local PSF model not implemented yet");
00187     exit_session(ERR_NOIMPL);
00188   }
00189 
00190   /* Initial guesses & boundaries */
00191 
00192   bl = (double *)malloc(3*npeak*sizeof(double));
00193   bu = (double *)malloc(3*npeak*sizeof(double));
00194 
00195   hw = SNIFS_INTERSP/2;
00196   glipeaknorm = (double *)malloc(npeak*sizeof(double));
00197   glxpeakoffs = (double *)malloc(npeak*sizeof(double));
00198 
00199   for (k=0; k<npeak; k++) {    /* Loop over the pack peaks */
00200     /* Position */
00201     /* In input, packpar[3*k] already contains imax(k), k=0..npeak-1 */
00202     bl[3*k] = packpar[3*k] - hw;
00203     bu[3*k] = packpar[3*k] + hw;
00204 
00205     /* Dispersion */
00206     bl[3*k+2] = 0.;
00207     bu[3*k+2] = MAXDOUBLE;
00208 
00209     /* Intensity */
00210     bl[3*k+1] =-MAXDOUBLE; /* Let the fit deal w/ that alone */
00211     bu[3*k+1] = MAXDOUBLE;
00212 
00213     glxpeakoffs[k] = packpar[3*k];   /* Used for fit normalization */
00214     glipeaknorm[k] = packpar[3*k+1];
00215   }
00216 
00217   /* Normalization */
00218   for (k=0; k<npeak; k++) {
00219     packpar[3*k  ] -= FMX_XPEAKOFFS(k); 
00220     bl[3*k]        -= FMX_XPEAKOFFS(k); 
00221     bu[3*k]        -= FMX_XPEAKOFFS(k);
00222     packpar[3*k  ] /= FMX_XPEAKNORM; 
00223     bl[3*k]        /= FMX_XPEAKNORM;
00224     bu[3*k]        /= FMX_XPEAKNORM;
00225     packpar[3*k+1] /= FMX_IPEAKNORM(k);
00226     packpar[3*k+2] /= FMX_SPEAKNORM;
00227   }
00228 
00229   /* Least-squares */
00230   dpackpar[0] = 0; /* Eastern egg: Print-out: 1, Derivatives: 2 */
00231 
00232   status = nllsqfit_bnd(NULL,xarr,sarr,npts,packpar,bl,bu,3*npeak,
00233                         *nllsqfit_maxpack, xfit, NULL, dpackpar, &chi2,&gof);
00234 
00235   if (status==0 || status==1 || status==4) {
00236     /* Keep "Optimal but..." and "Too many iterations" results */
00237     rms = sqrt(chi2/npts);
00238 
00239     /* Renormalization */
00240     for (k=0; k<npeak; k++) {
00241       packpar[3*k  ] *= FMX_XPEAKNORM;    dpackpar[3*k  ] *= FMX_XPEAKNORM;
00242       packpar[3*k  ] += FMX_XPEAKOFFS(k);
00243       packpar[3*k+1] *= FMX_IPEAKNORM(k); dpackpar[3*k+1] *= FMX_IPEAKNORM(k);
00244       packpar[3*k+2] *= FMX_SPEAKNORM;    dpackpar[3*k+2] *= FMX_SPEAKNORM;
00245 
00246       if (packpar[3*k+1]<0) {
00247         print_warning("I[%d]<0, skipping this pack",k+1);
00248         rms = -1.;
00249       }
00250     }
00251   }
00252   else {
00253     for (k=0; k<3*npeak; k++) packpar[k] = 0.;
00254     rms = -status;
00255   }
00256 
00257   free(bl); free(bu);
00258   free(glipeaknorm); free(glxpeakoffs);
00259   return(rms);
00260 }
00261 
00262 #undef FMX_XPEAKNORM
00263 #undef FMX_XPEAKOFFS
00264 #undef FMX_IPEAKNORM
00265 #undef FMX_SPEAKNORM
00266 

Generated on Wed Oct 26 23:59:39 2005 for Snifs by doxygen 1.3.5