Changeset 37642

Timestamp:

Nov 19, 2014, 7:56:39 PM (12 years ago)

Author:

eugene

Message:

working on the SH chi-sq fitting

Location:

branches/eam_branches/ipp-20140904/Ohana/src/relastro/src

Files:

• FrameCorrection.c (modified) (3 diffs)
• UpdateObjects.c (modified) (2 diffs)

branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/FrameCorrection.c

@@ -176 +176 @@
 }
 
+// I need to generate a collection of values dR,dD(R,D), where dR = average.R - ICRF.R, etc.  
+// I will then fit the set of dR,dD values to a set of spherical harmonics (up to Lmax).
+// I will then generate correction images dR,dD(R,D)  
+// I will then correct all values average.R,D by the interpolated values from the images
+// average.R' = average.R - dR(R,d), average.D' = average.D - dD(R,d)
+
+int SHfitWithMask (double *R, double *D, double *value, int *mask, int Npts, SHterms *fit);
+
 int FrameCorrectionFitSH (Catalog *catalog, int Ncatalog, SHterms *dRc, SHterms *dDc) {
 
+  int i;
+
   myAssert (dRc->lmax == dDc->lmax, "dR and dD must match\n");
-
-  // allocate an SHterms structure to hold the Ylm values 
-  SHterms *SH = SHtermsInit (dRc->lmax);
 
   int Nicrf = ICRFmax();
   fprintf (stderr, "start Fit %d SH l-modes for %d QSOs\n", dRc->lmax, Nicrf);
 
-  FILE *f = fopen ("sh.dat", "w");
-  myAssert (f, "oops");
-
-  int i, j;
+  double *Rave, *Dave, *dRoff, *dDoff;
+
+  int Npts = 0;
+  ALLOCATE (Rave,  double, Nicrf);
+  ALLOCATE (Dave,  double, Nicrf);
+  ALLOCATE (dRoff, double, Nicrf);
+  ALLOCATE (dDoff, double, Nicrf);
+
+  int *mask = NULL;
+  ALLOCATE (mask, int, Nicrf);
+  memset (mask, 0, Nicrf*sizeof(int));
+
+  // select the ICRF QSOS and save the necessary data
   for (i = 0; i < Nicrf; i++) {
 
     int cat, meas, ave;
     ICRFdata (i, &cat, &ave, &meas);
-
-    // I need to generate a collection of values dR,dD(R,D), where dR = average.R - ICRF.R, etc.  
-    // I will then fit the set of dR,dD values to a set of spherical harmonics (up to Lmax).
-    // I will then generate correction images dR,dD(R,D)  
-    // I will then correct all values average.R,D by the interpolated values from the images
-    // average.R' = average.R - dR(R,d), average.D' = average.D - dD(R,d)
 
     Average *average = &catalog[cat].average[ave];
@@ -209 +219 @@
     double dD = 3600.0*(average->D - measure->D);
 
-    fprintf (f, "%12.8f %12.8f %7.3f %7.3f\n", average->R, average->D, dR, dD);
     if (isnan(dR)) continue;
     if (isnan(dD)) continue;
 
-    // get the value of Ylm at this coordinate
-    SHtermsForRD (SH, average->R, average->D);
+
+    Rave[Npts] = average->R;
+    Dave[Npts] = average->D;
+    dRoff[Npts] = dR;
+    dDoff[Npts] = dD;
+    Npts ++;
+  }
+
+  SHfitWithMask (Rave, Dave, dRoff, mask, Npts, dRc);
+  SHfitWithMask (Rave, Dave, dDoff, mask, Npts, dDc);
+
+  // allocate an SHterms structure to hold the Ylm values 
+  SHterms *Ylm = SHtermsInit (dRc->lmax);
+
+  FILE *f = fopen ("sh.dat", "w");
+  myAssert (f, "oops");
+
+  for (i = 0; i < Nicrf; i++) {
     
-    // accumulate the dot products
-    for (j = 0; j < SH->Nterms; j++) {
-      dRc->Fr[j] += dR * SH->Fr[j];
-      dRc->Fi[j] += dR * SH->Fi[j];
-      dDc->Fr[j] += dD * SH->Fr[j];
-      dDc->Fi[j] += dD * SH->Fi[j];
-    }
+    SHtermsForRD (Ylm, Rave[i], Dave[i]);
+      
+    double dRfit = 0.0;
+    double dDfit = 0.0;
+
+    int j;
+    for (j = 0; j < Ylm->Nterms; j++) {
+      dRfit += dRc->Fr[j]*Ylm->Fr[j] + dRc->Fi[j]*Ylm->Fi[j];
+      dDfit += dDc->Fr[j]*Ylm->Fr[j] + dDc->Fi[j]*Ylm->Fi[j];
+    }
+    fprintf (f, "%12.8f %12.8f %7.3f %7.3f : %7.3f %7.3f\n", Rave[i], Dave[i], dRoff[i], dDoff[i], dRfit, dDfit);
   }
   fclose (f);
 
-  for (j = 0; j < SH->Nterms; j++) {
-    dRc->Fr[j] *= 4.0*M_PI/(float)Nicrf;
-    dRc->Fi[j] *= 4.0*M_PI/(float)Nicrf;
-    dDc->Fr[j] *= 4.0*M_PI/(float)Nicrf;
-    dDc->Fi[j] *= 4.0*M_PI/(float)Nicrf;
-  }
-  SHtermsFree (SH);
+  free (Rave);
+  free (Dave);
+  free (dRoff);
+  free (dDoff);
+  free (mask);
+
+  return TRUE;
+}
+
+int SHfitWithMask (double *R, double *D, double *value, int *mask, int Npts, SHterms *fit) {
+
+  int i, j, k;
+
+  // allocate an SHterms structure to hold the Ylm values 
+  SHterms *Ylm = SHtermsInit (fit->lmax);
+
+  // we only fit the linearly independent terms: Re(m >= 0), Im(m > 0)
+  int NtermRE = 0;
+  int NtermIM = 0;
+  for (i = 0; i < Ylm->Nterms; i++) {
+    if (Ylm->m[i] >= 0) NtermRE ++;
+    if (Ylm->m[i] >  0) NtermIM ++;
+  }
+  int *Nre = NULL;
+  int *Nim = NULL;
+  ALLOCATE (Nre, int, NtermRE);
+  ALLOCATE (Nim, int, NtermIM);
+
+  NtermRE = NtermIM = 0;
+  for (i = 0; i < Ylm->Nterms; i++) {
+    if (Ylm->m[i] >= 0) {
+      Nre[NtermRE] = i;
+      NtermRE ++;
+    }
+    if (Ylm->m[i] >  0) {
+      Nim[NtermIM] = i;
+      NtermIM ++;
+    }
+  }
+
+  double **Are, **bre, **Aim, **bim;
+  ALLOCATE (Are, double *, NtermRE);
+  ALLOCATE (bre, double *, NtermRE);
+  ALLOCATE (Aim, double *, NtermIM);
+  ALLOCATE (bim, double *, NtermIM);
+  for (i = 0; i < NtermRE; i++) {
+    ALLOCATE_ZERO (Are[i], double, NtermRE);
+    ALLOCATE_ZERO (bre[i], double, 1);
+  }
+  for (i = 0; i < NtermIM; i++) {
+    ALLOCATE_ZERO (Aim[i], double, NtermIM);
+    ALLOCATE_ZERO (bim[i], double, 1);
+  }
+  
+  // measure the dot product \sum(F_i * Ylm_i) and the cross terms (\sum(Y_lm * Y_jk))
+  int Nfit = 0;
+  for (i = 0; i < Npts; i++) {
+    if (mask[i]) continue;
+    Nfit ++;
+
+    // set the values of Ylm(R[i],D[i])
+    SHtermsForRD (Ylm, R[i], D[i]); 
+
+    double Fv = value[i];
+
+    for (j = 0; j < NtermRE; j++) {
+      int jre = Nre[j];
+      bre[j][0] += Fv * Ylm->Fr[jre];
+    }
+    for (j = 0; j < NtermIM; j++) {
+      int jim = Nim[j];
+      bim[j][0] += Fv * Ylm->Fi[jim];
+    }
+
+    for (j = 0; j < NtermRE; j++) {
+      int jre = Nre[j];
+      for (k = j; k < NtermRE; k++) {
+        int kre = Nre[k];
+        Are[j][k] += Ylm->Fr[jre] * Ylm->Fr[kre];
+      }
+    }
+
+    for (j = 0; j < NtermIM; j++) {
+      int jim = Nim[j];
+      for (k = j; k < NtermIM; k++) {
+        int kim = Nim[k];
+        Aim[j][k] += Ylm->Fi[jim] * Ylm->Fi[kim];
+      }
+    }
+  }
+
+  for (j = 1; j < NtermRE; j++) {
+    for (k = 0; k < j; k++) {
+      Are[j][k] = Are[k][j];
+    }   
+  }
+  for (j = 1; j < NtermIM; j++) {
+    for (k = 0; k < j; k++) {
+      Aim[j][k] = Aim[k][j];
+    }   
+  }
+
+  if (!dgaussjordan (Are, bre, NtermRE, 1)) {
+    gprint (GP_ERR, "failed to fit data : ill-conditioned matrix\n");
+    return FALSE;
+  }
+  if (!dgaussjordan (Aim, bim, NtermIM, 1)) {
+    gprint (GP_ERR, "failed to fit data : ill-conditioned matrix\n");
+    return FALSE;
+  }
+
+  for (j = 0; j < fit->Nterms; j++) { 
+    fit->Fr[j] = 0.0;
+    fit->Fi[i] = 0.0;
+  }
+  for (j = 0; j < NtermRE; j++) {
+    int jre = Nre[j];
+    // fit->Fr[jre] = bre[j][0] * fit->Nterms / (float) Nfit; // XXX EAM : why is this factor needed?
+    fit->Fr[jre] = bre[j][0];
+  }
+  for (j = 0; j < NtermIM; j++) {
+    int jim = Nim[j];
+    // fit->Fi[jim] = bim[j][0] * fit->Nterms / (float) Nfit; // XXX EAM : why is this factor needed?
+    fit->Fi[jim] = bim[j][0];
+  }
+
+  SHtermsFree (Ylm);
+
+  for (i = 0; i < NtermRE; i++) {
+    free (Are[i]);
+    free (bre[i]);
+  }
+  for (i = 0; i < NtermIM; i++) {
+    free (Aim[i]);
+    free (bim[i]);
+  }
+  free (Are);
+  free (bre);
+  free (Aim);
+  free (bim);
+
   return TRUE;
 }
@@ -542 +705 @@
         double dR = frame->Roff[iD][iR];
         double dD = frame->Doff[iD][iR];
+
+        // do not apply if the fabs(offset) is > 8 arcsec)
+        if ((fabs(dR) > 8.0) || (fabs(dD) > 8.0)) {
+          fprintf (stderr, "skip: %10.6f %10.6f : %7.3f %7.3f\n", average->R, average->D, dR, dD);
+          continue;
+        }
+
         average->R -= dR / 3600.0;
         average->D -= dD / 3600.0;

branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/UpdateObjects.c

@@ -2 +2 @@
 # define PAR_TOOFEW 5
 
-int UpdateObjects_Chips (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats, int i, off_t m, int Nloop);
+int UpdateObjects_Chips (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats, int i, off_t m, int applyGalaxyOffset);
 int UpdateObjects_Stack (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats);
 
@@ -146 +146 @@
 // This function operates on both Measure and MeasureTiny.  In the big stages, this should
 // be called with just MeasureTiny set and Measure == NULL
-int UpdateObjects_Chips (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats, int i, off_t m, int Nloop) {
+int UpdateObjects_Chips (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats, int i, off_t m, int applyGalaxyOffset) {
 
   int setRefColor = areImagesMatched();