Changeset 40406 for trunk

Timestamp:

Apr 25, 2018, 10:29:37 AM (8 years ago)

Author:

eugene

Message:

fix segfault due to missing mask vector

File:

• trunk/Ohana/src/opihi/cmd.astro/fitplx_irls.c (modified) (6 diffs)

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

@@ -12 +12 @@
     if ((mvec = SelectVector (argv[N], ANYVECTOR, TRUE)) == NULL) return (FALSE);
     remove_argument (N, &argc, argv);
-    CastVector (mvec, OPIHI_INT);
   }
 
@@ -81 +80 @@
   double *dD = dDvec->elements.Flt;
 
-  opihi_int *mask = NULL;
-  if (mvec) {
-    mask = mvec->elements.Int;
-  }
-
   // Ntotal : all points supplied by user
   // Nsubset : unmasked points
   int Ntotal = tvec->Nelements; // XXX check other lengths
-  
+
+  if (rvec->Nelements  != Ntotal) ESCAPE ("mis-match in vector lengths (mjd vs ra)  = (%d vs %d)\n", rvec->Nelements, Ntotal);
+  if (dvec->Nelements  != Ntotal) ESCAPE ("mis-match in vector lengths (mjd vs dec) = (%d vs %d)\n", dvec->Nelements, Ntotal);
+  if (dRvec->Nelements != Ntotal) ESCAPE ("mis-match in vector lengths (mjd vs dR)  = (%d vs %d)\n", dRvec->Nelements, Ntotal);
+  if (dDvec->Nelements != Ntotal) ESCAPE ("mis-match in vector lengths (mjd vs dD)  = (%d vs %d)\n", dDvec->Nelements, Ntotal);
+  
+  opihi_int *mask = NULL;
+  if (mvec) {
+    ResetVector (mvec, OPIHI_INT, Ntotal);
+    mask = mvec->elements.Int;
+    for (i = 0; i < Ntotal; i++) { mask[i] = 1; }
+  }
+
+
   double Rmean, Dmean, Tmean;
   PlxSetMeanEpoch (R, D, T, &Rmean, &Dmean, &Tmean, mask, Ntotal);
@@ -135 +142 @@
 
   // update the mask based on the input mask and the outlier limits.
-  double Sum_Wx = 0;
-  double Sum_Wy = 0;
-  mask = mvec->elements.Int;
-  for (i = 0; i < fitdata.Npts; i++) {
-    Sum_Wx += fitdata.Wx[i];
-    Sum_Wy += fitdata.Wy[i];
-  }
-  for (i = 0; i < fitdata.Npts; i++) {
-    // fitdata only includes the previously unmasked points
-    if ((fitdata.Wx[i] < outlier_limit * Sum_Wx / (1.0 * fitdata.Npts))||
-        (fitdata.Wy[i] < outlier_limit * Sum_Wy / (1.0 * fitdata.Npts))) {
-      int n = fitdata.index[i];
-      mask[n] = 0;
+  if (mask) {
+    double Sum_Wx = 0;
+    double Sum_Wy = 0;
+    
+    // calculate the total weight
+    for (i = 0; i < fitdata.Npts; i++) {
+      Sum_Wx += fitdata.Wx[i];
+      Sum_Wy += fitdata.Wy[i];
+    }
+    for (i = 0; i < fitdata.Npts; i++) {
+      // fitdata only includes the previously unmasked points
+      if ((fitdata.Wx[i] < outlier_limit * Sum_Wx / (1.0 * fitdata.Npts))||
+          (fitdata.Wy[i] < outlier_limit * Sum_Wy / (1.0 * fitdata.Npts))) {
+        int n = fitdata.index[i];
+        mask[n] = 0;
       
-      if (VERBOSE == 2) {
+        if (VERBOSE == 2) {
           fprintf (stderr, "%f %f : %f "OPIHI_INT_FMT" : %f %f %f : %f %f %f %f\n", R[n], D[n], T[n], mask[n], fitdata.t[i], fitdata.X[i], fitdata.Y[i], fitdata.Wx[i], fitdata.Wy[i], Sum_Wx, Sum_Wy);
-      }
-    }
-  }
-
-  if (Nresample){
-    // now that the mask has been updated, we need to recalculate mean epoch and positions
-    // XXX make this conditional on actually masking unmasked points above
-    // PlxSetMeanEpoch (R, D, T, &Rmean, &Dmean, &Tmean, mask, Ntotal);
-    // PlxSetEpochPosition (&fitdata, R, D, dR, dD, T, mask, Ntotal, &coords, Tmean);
+        }
+      }
+    }
+  }
+
+  if (Nresample) {
+    // if the mask has been updated, we need to recalculate mean epoch and positions
+    if (mask) {
+      PlxSetMeanEpoch (R, D, T, &Rmean, &Dmean, &Tmean, mask, Ntotal);
+      PlxSetEpochPosition (&fitdata, R, D, dR, dD, T, mask, Ntotal, &coords, Tmean);
+    }
 
     PlxFitData sample;
@@ -218 +229 @@
   // fprintf (stderr, "%f +/- %f | %f %f\n", fit.p, fit.dp, fit.uR, fit.uD);
 
-  Vector *dRresPMP, *dDresPMP, *dRresPLX, *dDresPLX;
+  Vector *dRresPOS, *dDresPOS, *dRresPMP, *dDresPMP, *dRresPLX, *dDresPLX;
 
   // save fit residuals (with only pm removed, and pm and plx removed)
+  if ((dRresPOS = SelectVector ("dRresPOS", ANYVECTOR, TRUE)) == NULL) ESCAPE ("cannot generate vector %s\n", "dRresPOS");
+  if ((dDresPOS = SelectVector ("dDresPOS", ANYVECTOR, TRUE)) == NULL) ESCAPE ("cannot generate vector %s\n", "dDresPOS");
   if ((dRresPMP = SelectVector ("dRresPMP", ANYVECTOR, TRUE)) == NULL) ESCAPE ("cannot generate vector %s\n", "dRresPMP");
   if ((dDresPMP = SelectVector ("dDresPMP", ANYVECTOR, TRUE)) == NULL) ESCAPE ("cannot generate vector %s\n", "dDresPMP");
@@ -226 +239 @@
   if ((dDresPLX = SelectVector ("dDresPLX", ANYVECTOR, TRUE)) == NULL) ESCAPE ("cannot generate vector %s\n", "dDresPLX");
     
+  ResetVector (dRresPOS, OPIHI_FLT, Ntotal);
+  ResetVector (dDresPOS, OPIHI_FLT, Ntotal);
   ResetVector (dRresPMP, OPIHI_FLT, Ntotal);
   ResetVector (dDresPMP, OPIHI_FLT, Ntotal);
@@ -246 +261 @@
     double Yplx = fit.Do + fit.uD*t0;
 
+    dRresPOS->elements.Flt[i] = x0;
+    dDresPOS->elements.Flt[i] = y0;
     dRresPMP->elements.Flt[i] = x0 - Xpmp;
     dDresPMP->elements.Flt[i] = y0 - Ypmp;