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_pack normalization */
00035 double *glipeaknorm, *glxpeakoffs;
00036 
00037 /* find_max: Fit_pack */
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_pack normalization */
00158 double *glipeaknorm, *glxpeakoffs;
00159 
00160 /* === Doxygen Comment ======================================= */
00173 /* =========================================================== */
00174 
00175 double Fit_pack(double xarr[], int npts, int npeak, 
00176                 double packpar[], double dpackpar[], double xfit[])
00177 {
00178   int k, status;
00179   double *bl, *bu;
00180   double hw, chi2, gof, rms;
00181 
00182   if (glngloc > 1) {
00183     print_error("Multi-gaussian local PSF model not implemented yet");
00184     exit_session(ERR_NOIMPL);
00185   }
00186 
00187   /* Initial guesses & boundaries */
00188 
00189   bl = (double *)malloc(3*npeak*sizeof(double));
00190   bu = (double *)malloc(3*npeak*sizeof(double));
00191 
00192   hw = SNIFS_INTERSP/2;
00193   glipeaknorm = (double *)malloc(npeak*sizeof(double));
00194   glxpeakoffs = (double *)malloc(npeak*sizeof(double));
00195 
00196   for (k=0; k<npeak; k++) {    /* Loop over the pack peaks */
00197     /* Position */
00198     /* In input, packpar[3*k] already contains imax(k), k=0..npeak-1 */
00199     bl[3*k] = packpar[3*k] - hw;
00200     bu[3*k] = packpar[3*k] + hw;
00201 
00202     /* Dispersion */
00203     bl[3*k+2] = 0.;
00204     bu[3*k+2] = MAXDOUBLE;
00205 
00206     /* Intensity */
00207     bl[3*k+1] =-MAXDOUBLE; /* Let the fit deal w/ that alone */
00208     bu[3*k+1] = MAXDOUBLE;
00209 
00210     glxpeakoffs[k] = packpar[3*k];   /* Used for fit normalization */
00211     glipeaknorm[k] = packpar[3*k+1];
00212   }
00213 
00214   /* Normalization */
00215   for (k=0; k<npeak; k++) {
00216     packpar[3*k  ] -= FMX_XPEAKOFFS(k); 
00217     bl[3*k] -= FMX_XPEAKOFFS(k); bu[3*k] -= FMX_XPEAKOFFS(k);
00218     packpar[3*k  ] /= FMX_XPEAKNORM; 
00219     bl[3*k] /= FMX_XPEAKNORM;    bu[3*k] /= FMX_XPEAKNORM;
00220     packpar[3*k+1] /= FMX_IPEAKNORM(k);
00221     packpar[3*k+2] /= FMX_SPEAKNORM;
00222   }
00223 
00224   /* Least-squares */
00225   dpackpar[0] = 0; /* Eastern egg: Print-out: 1, Derivatives: 2 */
00226 
00227   status = nllsqfit_bnd(NULL,xarr,NULL,npts,packpar,bl,bu,3*npeak,
00228                         *nllsqfit_maxpack, xfit, NULL, dpackpar, &chi2,&gof);
00229 
00230   if (status==0 || status==1 || status==4) {
00231     /* Keep "Optimal but..." and "Too many iterations" results */
00232     rms = sqrt(chi2/npts);
00233 
00234     /* Renormalization */
00235     for (k=0; k<npeak; k++) {
00236       packpar[3*k  ] *= FMX_XPEAKNORM;    dpackpar[3*k  ] *= FMX_XPEAKNORM;
00237       packpar[3*k  ] += FMX_XPEAKOFFS(k);
00238       packpar[3*k+1] *= FMX_IPEAKNORM(k); dpackpar[3*k+1] *= FMX_IPEAKNORM(k);
00239       packpar[3*k+2] *= FMX_SPEAKNORM;    dpackpar[3*k+2] *= FMX_SPEAKNORM;
00240 
00241       if (packpar[3*k+1]<0) {
00242         print_warning("I[%d]<0, skipping this pack",k+1);
00243         rms = -1.;
00244       }
00245     }
00246   }
00247   else {
00248     for (k=0; k<3*npeak; k++) packpar[k] = 0.;
00249     rms = -status;
00250   }
00251 
00252   free(bl); free(bu);
00253   free(glipeaknorm); free(glxpeakoffs);
00254   return(rms);
00255 }
00256 
00257 #undef FMX_XPEAKNORM
00258 #undef FMX_XPEAKOFFS
00259 #undef FMX_IPEAKNORM
00260 #undef FMX_SPEAKNORM
00261 

Generated on Tue Nov 23 18:04:19 2004 for Snifs by doxygen 1.3.3