Changeset 39585 for trunk/Ohana/src/opihi/cmd.astro/fitpm_irls.c

Timestamp:

May 27, 2016, 5:17:08 PM (10 years ago)

Author:

watersc1

Message:

Added binned initial step method for irls pm/plx fitting.

File:

• trunk/Ohana/src/opihi/cmd.astro/fitpm_irls.c (modified) (8 diffs)

trunk/Ohana/src/opihi/cmd.astro/fitpm_irls.c

@@ -25 +25 @@
   }
 
+  double binning_step = 0.0;
+  if ((N = get_argument (argc, argv, "-binning"))) {
+    remove_argument (N, &argc, argv);
+    binning_step = atof(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+  
+  
   if (argc != 6) {
-    gprint (GP_ERR, "USAGE: fitplx_irls (ra) (dR) (dec) (dD) (mjd) [-mask mask]\n");
+    gprint (GP_ERR, "USAGE: fitplx_irls (ra) (dR) (dec) (dD) (mjd) [-mask mask] [-binning step_size]\n");
     // what about the errors?
     return (FALSE);
@@ -105 +113 @@
 
   PlxFit fit;
-  if (!FitPMonly_IRLS (&fit, X, dX, Y, dY, t, n, VERBOSE)) {
+  if (!FitPMonly_IRLS (&fit, X, dX, Y, dY, t, n, binning_step, VERBOSE)) {
     return FALSE;
   }
@@ -140 +148 @@
 
 /* do we want an init function which does the alloc and a clear function to free? */
-int FitPMonly_IRLS (PlxFit *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts, int VERBOSE) {
+int FitPMonly_IRLS (PlxFit *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts, double binning_step, int VERBOSE) {
 
   int i,j;
@@ -196 +204 @@
   }
   
-  // Solve OLS equation  
-  if (!weighted_LS_PM(T,X,Wx,Y,Wy,Npts,
-                   A,B,VERBOSE)) {
-    // Handle fail case
-    return(FALSE);
+  // Solve OLS equation
+  if (binning_step == 0.0) {
+    if (!weighted_LS_PM(T,X,Wx,Y,Wy,Npts,
+                        A,B,VERBOSE)) {
+      // Handle fail case
+      return(FALSE);
+    }
+  }
+  else {
+    int Nbins;
+    double *Tbin = NULL;
+    double *Xbin = NULL;
+    double *Ybin = NULL;
+    double *WXbin = NULL;
+    double *WYbin = NULL;
+    
+    if (!bin_points(T, X, Wx, Y, Wy, Npts,
+                    &Tbin, &Xbin, &WXbin, &Ybin, &WYbin, &Nbins,  binning_step) ) {
+      return(FALSE);
+    }
+    if (!weighted_LS_PM(Tbin,Xbin,WXbin,Ybin,WYbin, Nbins,
+                        A,B,VERBOSE)) {
+      return(FALSE);
+    }
+    FREE(Tbin);
+    FREE(Xbin);
+    FREE(Ybin);
+    FREE(WXbin);
+    FREE(WYbin);
   }
 
@@ -424 +456 @@
 }
 
+int bin_points (double *T, double *X, double *WX, double *Y, double *WY, int Npts,
+                double **Tbin, double **Xbin, double **WXbin, double **Ybin, double **WYbin, int *Nbins, double binning_step) {
+  double *T_tmp = NULL;
+  int i,j,k,l;
+
+  int Nbins_test;
+  int *bin_test = NULL;
+
+  double *T_int = NULL, *X_int = NULL, *Y_int = NULL, *WX_int = NULL, *WY_int = NULL;
+  double T_min, T_max;
+
+
+  //  printf("In binning code\n");
+
+  // Allocate more bins than we need
+  ALLOCATE(T_tmp,double,Npts);
+  for (i = 0; i < Npts; i++) {
+    T_tmp[i] = T[i];
+  }
+  dsort(T_tmp,Npts);
+
+  Nbins_test = floor((T_tmp[Npts - 1] - T_tmp[0]) / binning_step) + 1;
+  ALLOCATE(bin_test,int,Nbins_test);
+
+  T_min = T_tmp[0];
+  T_max = T_tmp[Npts - 1];
+
+  // printf("Found %d points, putting into %d bins between %f and %f\n",Npts,Nbins_test,T_min,T_max);
+  // Check how many bins are filled.
+  for (j = 0; j < Nbins_test; j++) {
+    bin_test[j] = 0;
+  }
+  for (i = 0; i < Npts; i++) {
+    j = floor((T_tmp[i] - T_tmp[0]) / binning_step);
+    bin_test[j] ++;
+  }
+
+  *Nbins = 0;
+  for (j = 0; j < Nbins_test; j++) {
+    if (bin_test[j] != 0) {
+      *Nbins = *Nbins + 1;
+    }
+  }
+
+  // printf("Using %d actual bins\n",*Nbins);
+
+  double *Tbin_in = NULL;
+  double *Xbin_in = NULL;
+  double *Ybin_in = NULL;
+  double *WXbin_in = NULL;
+  double *WYbin_in = NULL;
+  
+  ALLOCATE(Tbin_in,double,*Nbins);
+  ALLOCATE(Xbin_in,double,*Nbins);
+  ALLOCATE(WXbin_in,double,*Nbins);
+  ALLOCATE(Ybin_in,double,*Nbins);
+  ALLOCATE(WYbin_in,double,*Nbins);
+
+  k = 0;
+  for (j = 0; j < Nbins_test; j++) {
+    if (bin_test[j] == 0) { // No data for this bin
+      continue;
+    }
+    else {
+      // printf("%d bin %d, N = %d\n",k,j,bin_test[j]);
+      // Allocate internal arrays to hold the data that goes into this bin.
+      ALLOCATE(T_int,double, bin_test[j]);
+      ALLOCATE(X_int,double, bin_test[j]);
+      ALLOCATE(Y_int,double, bin_test[j]);
+      ALLOCATE(WX_int,double, bin_test[j]);
+      ALLOCATE(WY_int,double, bin_test[j]);
+
+      // Fill those arrays.
+      l = 0;
+      for (i = 0; i < Npts; i++) {
+        if ((T[i] >= T_min + j * binning_step)&&
+            (T[i] <  T_min + (j + 1) * binning_step)) {
+          T_int[l] = T[i];
+          X_int[l] = X[i];
+          Y_int[l] = Y[i];
+          WX_int[l] = WX[i];
+          WY_int[l] = WY[i];
+          // printf("%d %d %f %f %f %f %f\n",i,l,T[i],X[i],Y[i],WX[i],WY[i]);
+          l++;
+        }
+      } // End loop over input points
+      // printf("%d bin %d, N = %d, l = %d\n",k,j,bin_test[j],l);
+
+      // Make a decision what the binned values should be.
+      if (l == 1) {
+        Tbin_in[k] = T_int[0];
+        Xbin_in[k] = X_int[0];
+        Ybin_in[k] = Y_int[0];
+        WXbin_in[k] = WX_int[0];
+        WYbin_in[k] = WY_int[0];
+      }
+      else {
+        // I think I'm going with medians.
+        dsort(T_int,l);
+        dsort(X_int,l);
+        dsort(Y_int,l);
+        
+        // The median gives the midpoint of all the measurements
+        if ((l % 2) == 0) {
+          Tbin_in[k] = 0.5*(T_int[(int)(0.5*l)] + T_int[(int)(0.5*l) - 1]);
+          Xbin_in[k] = 0.5*(X_int[(int)(0.5*l)] + X_int[(int)(0.5*l) - 1]);
+          Ybin_in[k] = 0.5*(Y_int[(int)(0.5*l)] + Y_int[(int)(0.5*l) - 1]);
+        } else {
+          Tbin_in[k] = T_int[(int)(0.5*l)];
+          Xbin_in[k] = X_int[(int)(0.5*l)];
+          Ybin_in[k] = Y_int[(int)(0.5*l)];
+        }
+        
+        // The scatter between points is probably the most useful measurement for the error.
+        for (i = 0; i < l; i++) {
+          X_int[i] = fabs(X_int[i] - Xbin_in[k]);
+          Y_int[i] = fabs(Y_int[i] - Ybin_in[k]);
+        }
+        dsort(X_int,l);
+        dsort(Y_int,l);
+        
+        if ((l % 2) == 0) {
+          WXbin_in[k] = 1.0 / (0.5*(X_int[(int)(0.5*l)] + X_int[(int)(0.5*l) - 1]));
+          WYbin_in[k] = 1.0 / (0.5*(Y_int[(int)(0.5*l)] + Y_int[(int)(0.5*l) - 1]));
+        } else {
+          WXbin_in[k] = 1.0 / (X_int[(int)(0.5*l)]);
+          WYbin_in[k] = 1.0 / (Y_int[(int)(0.5*l)]);
+        }
+      }
+
+      if (WXbin_in[k] == 0.0) {
+        WXbin_in[k] = WX_int[0];
+      }
+      if (WYbin_in[k] == 0.0) {
+        WYbin_in[k] = WY_int[0];
+      }
+      if (!isfinite(WXbin_in[k])) {
+        WXbin_in[k] = 1.0;
+      }
+      if (!isfinite(WYbin_in[k])) {
+        WYbin_in[k] = 1.0;
+      }
+      
+      //      printf("%f %f %f %f %f\n",Tbin_in[k],Xbin_in[k],Ybin_in[k],WXbin_in[k],WYbin_in[k]);
+      // Increment bin index.
+      k++;
+      
+      // Clean up
+      FREE(T_int);
+      FREE(X_int);
+      FREE(Y_int);
+      FREE(WX_int);
+      FREE(WY_int);
+    } 
+  } // End loop over initial bin set.
+
+  *Tbin = Tbin_in;
+  *Xbin = Xbin_in;
+  *Ybin = Ybin_in;
+  *WXbin = WXbin_in;
+  *WYbin = WYbin_in;
+
+  FREE(T_tmp);
+  FREE(bin_test);
+  
+  return TRUE;
+}
+
+
 double weight_cauchy (double x) {
   double r = x / 2.385;
@@ -450 +651 @@
   dsort(x,N);
 
-  if (N % 2) {
+  if ((N % 2) == 0) {
     median = 0.5*(x[(int)(0.5*N)] + x[(int)(0.5*N) - 1]);
   } else {
@@ -462 +663 @@
   dsort(x,N);
 
-  if (N % 2) {
+  if ((N % 2) == 0) {
     median = 0.5*(x[(int)(0.5*N)] + x[(int)(0.5*N) - 1]);
   } else {
@@ -468 +669 @@
   }
 
+  FREE(x);
   return(median);
 }