Changeset 33037

Timestamp:

Dec 30, 2011, 8:05:04 AM (15 years ago)

Author:

eugene

Message:

avoid double counting the pairs

Location:

branches/eam_branches/ipp-20111122/Ohana/src/opihi/cmd.astro

Files:

• spex1dgas.c (modified) (2 diffs)
• spexseq.c (modified) (13 diffs)

branches/eam_branches/ipp-20111122/Ohana/src/opihi/cmd.astro/spex1dgas.c

@@ -172 +172 @@
     }
 
-    N = index2->elements.Int[i];
-    if (N >= Nobject) abort();
-    m = object[N].Nindex;
-    object[N].index[m] = index1->elements.Int[i];
-    object[N].Dtgt[m] = distance->elements.Flt[i];
-    object[N].Dcur[m] = 0.0;
-    object[N].Nindex ++;
-    if (object[N].Nindex == object[N].NINDEX) {
-      object[N].NINDEX += 100;
-      REALLOCATE (object[N].Dtgt, float, object[N].NINDEX);
-      REALLOCATE (object[N].Dcur, float, object[N].NINDEX);
-      REALLOCATE (object[N].index, int, object[N].NINDEX);
-    }
+    // Above we used the set (N,M,D_N,M) to assigned the distance of obj M to obj N.
+    // Below we use that set to assign the distance of obj N to obj M.  IF we supply all
+    // pairwise distances, then this duplicates all entries, so we need to skip that step
+    // in such a case.
+
+    // XX N = index2->elements.Int[i];
+    // XX if (N >= Nobject) abort();
+    // XX m = object[N].Nindex;
+    // XX object[N].index[m] = index1->elements.Int[i];
+    // XX object[N].Dtgt[m] = distance->elements.Flt[i];
+    // XX object[N].Dcur[m] = 0.0;
+    // XX object[N].Nindex ++;
+    // XX if (object[N].Nindex == object[N].NINDEX) {
+    // XX   object[N].NINDEX += 100;
+    // XX   REALLOCATE (object[N].Dtgt, float, object[N].NINDEX);
+    // XX   REALLOCATE (object[N].Dcur, float, object[N].NINDEX);
+    // XX   REALLOCATE (object[N].index, int, object[N].NINDEX);
+    // XX }
   }
 
@@ -221 +226 @@
   // int idx2 = object[idx1].index[(int)(0.75*object[idx1].Nindex)];
   float A = NAN;
-  for (i = 0; isnan (A) && (i < object[idx1].Nindex); i++) {
-      if (object[idx1].index[i] == idx2) {
-          A = object[idx1].Dtgt[i];
+  for (i = 0; isnan (A) && (i < object[idx2].Nindex); i++) {
+      if (object[idx2].index[i] == idx1) {
+          A = object[idx2].Dtgt[i];
       }
   }

branches/eam_branches/ipp-20111122/Ohana/src/opihi/cmd.astro/spexseq.c

@@ -36 +36 @@
   float Xo, dXo;
   float Yo, dYo;
+  float Srange;
+  int nSeq;
 } Object;
 
@@ -57 +59 @@
 
   // is idx2 closer to pin1 or pin2?
-  N = object[idx2].rindex[pin1];
-  float D1 = object[idx2].Dtgt[N];
-  N = object[idx2].rindex[pin2];
-  float D2 = object[idx2].Dtgt[N];
-
-  int parity = (D1 < D2) ? +1 : -1;
+  // XX N = object[idx2].rindex[pin1];
+  // XX float D1 = object[idx2].Dtgt[N];
+  // XX N = object[idx2].rindex[pin2];
+  // XX float D2 = object[idx2].Dtgt[N];
+  // XX int parity = (D1 < D2) ? +1 : -1;
 
   int Nfriends = 0;
@@ -70 +71 @@
     j = object[idx2].rindex[N];
     assert (object[idx2].index[j] == N);
-    object[idx1].seq[i] = parity*(object[idx2].Dtgt[j] - A);
+    // XX object[idx1].seq[i] = parity*(object[idx2].Dtgt[j] - A);
+    object[idx1].seq[i] = (object[idx2].Dtgt[j] - A);
     Nfriends ++;
   }
+  assert (Nfriends <= object[idx1].Nindex);
+  assert (Nfriends <= Nobject);
   fprintf (stderr, "object %d sequenced by object %d, %d friends sequenced\n", idx1, idx2, Nfriends);
   return;
@@ -144 +148 @@
     }
 
-    N = index2->elements.Int[i];
-    if (N >= Nobject) abort();
-    m = object[N].Nindex;
-    object[N].index[m] = index1->elements.Int[i];
-    object[N].Dtgt[m] = distance->elements.Flt[i];
-    object[N].Nindex ++;
-    if (object[N].Nindex == object[N].NINDEX) {
-      object[N].NINDEX += 100;
-      REALLOCATE (object[N].Dtgt, float, object[N].NINDEX);
-      REALLOCATE (object[N].index, int, object[N].NINDEX);
-    }
+    // Above we used the set (N,M,D_N,M) to assigned the distance of obj M to obj N.
+    // Below we use that set to assign the distance of obj N to obj M.  IF we supply all
+    // pairwise distances, then this duplicates all entries, so we need to skip that step
+    // in such a case.
+
+    // XX N = index2->elements.Int[i];
+    // XX if (N >= Nobject) abort();
+    // XX m = object[N].Nindex;
+    // XX object[N].index[m] = index1->elements.Int[i];
+    // XX object[N].Dtgt[m] = distance->elements.Flt[i];
+    // XX object[N].Nindex ++;
+    // XX if (object[N].Nindex == object[N].NINDEX) {
+    // XX   object[N].NINDEX += 100;
+    // XX   REALLOCATE (object[N].Dtgt, float, object[N].NINDEX);
+    // XX   REALLOCATE (object[N].index, int, object[N].NINDEX);
+    // XX }
   }
 
@@ -183 +192 @@
   // find two objects at the opposite extremes
   float Dm1 = object[0].Dtgt[object[0].Nindex - 1];
-  for (i = 0; (object[0].Dtgt[i] < 0.9*Dm1) && (i < object[0].Nindex); i++);
+  for (i = 0; (object[0].Dtgt[i] < 0.75*Dm1) && (i < object[0].Nindex); i++);
   int pin1 = object[0].index[i];
   assert (pin1 >= 0);
@@ -189 +198 @@
 
   float Dm2 = object[pin1].Dtgt[object[pin1].Nindex - 1];
-  for (i = 0; (object[pin1].Dtgt[i] < 0.9*Dm2) && (i < object[pin1].Nindex); i++);
+  for (i = 0; (object[pin1].Dtgt[i] < 0.75*Dm2) && (i < object[pin1].Nindex); i++);
   int pin2 = object[pin1].index[i];
   assert (pin2 >= 0);
@@ -223 +232 @@
       float dS2m = 0.0;
       int nS = 0;
+
+      float Smin = +1000.0;
+      float Smax = -1000.0;
 
       for (k = 0; k < object[Nj].Nindex; k++) {
@@ -231 +243 @@
         if (isnan(object[i].seq[Nsi])) continue;
         
+        Smin = MIN(object[Nj].seq[k], Smin);
+        Smax = MAX(object[Nj].seq[k], Smax);
+
         float dSp = object[Nj].seq[k] - object[i].seq[Nsi];
         float dSm = -object[Nj].seq[k] - object[i].seq[Nsi];
@@ -243 +258 @@
       float dSp = sqrt(dS2p / nS - Sp*Sp);
       float dSm = sqrt(dS2m / nS - Sm*Sm);
-      fprintf (stderr, "%d: %f +/- %f or %f +/- %f\n", Nj, Sp, dSp, Sm, dSm);
-
-      // if (dSm > dSp) {
-      if (1) {
-        // parity is always 1 (based on pins)
+      fprintf (stderr, "%d: %f +/- %f or %f +/- %f : %f range, %d matches\n", Nj, Sp, dSp, Sm, dSm, Smax - Smin, nS);
+
+      if (dSm > dSp) {
+        // parity is always 1 if we use the pins (but it does not seem to work)
         for (k = 0; k < object[Nj].Nindex; k++) {
           if (isnan(object[Nj].seq[k])) continue;
@@ -263 +277 @@
         object[Nj].dSo = dSm;
       }
+      object[Nj].Srange = Smax - Smin;
+      object[Nj].nSeq = nS;
       object[Nj].shifted = TRUE;
       Nfixed ++;
     }
-    fprintf (stderr, "%d used to reconcile %d friends\n", i, Nfixed);
+    // fprintf (stderr, "%d used to reconcile %d friends\n", i, Nfixed);
   }
 
@@ -318 +334 @@
     Vector *outdXo   = SelectVector ("sp1d_dXo", ANYVECTOR, TRUE); if (!outXo) goto escape;
 
+    Vector *outSr   = SelectVector ("sp1d_Sr", ANYVECTOR, TRUE); if (!outSr) goto escape;
+    Vector *outSn   = SelectVector ("sp1d_nS", ANYVECTOR, TRUE); if (!outSn) goto escape;
+
     ResetVector (outindex, OPIHI_INT, Nobject);
     ResetVector (outXo, OPIHI_FLT, Nobject);
@@ -325 +344 @@
     ResetVector (outdSo, OPIHI_FLT, Nobject);
     ResetVector (outdXo, OPIHI_FLT, Nobject);
+
+    ResetVector (outSr, OPIHI_FLT, Nobject);
+    ResetVector (outSn, OPIHI_INT, Nobject);
 
     for (i = 0; i < Nobject; i++) {
@@ -333 +355 @@
       outdSo->elements.Flt[i] = object[i].dSo;
       outdXo->elements.Flt[i] = object[i].dXo;
+      outSr->elements.Flt[i] = object[i].Srange;
+      outSn->elements.Int[i] = object[i].nSeq;
     }
   }