Index: /branches/eam_branches/ohana.20170822/src/opihi/dvo/find_dups.c
===================================================================
--- /branches/eam_branches/ohana.20170822/src/opihi/dvo/find_dups.c	(revision 40171)
+++ /branches/eam_branches/ohana.20170822/src/opihi/dvo/find_dups.c	(revision 40172)
@@ -1,12 +1,9 @@
 # include "data.h"
-// NOTE: if there are only a few uniq values, the old algorithm is not bad.  
-// for 10000 uniq values, 30M points take ~20sec in the new algorithm, 
-// 3M points takess 45 sec in the old method.
 
-/* given imageID and detID, make a 64bit joint ID and sort it */ 
-int find_dups (int argc, char **argv) {
+/* given two 32bit IDs, a 64bit joint ID and sort it, along with an index (sequence) vector */ 
+int uniqpair (int argc, char **argv) {
 
-  int Nnew, i, N;
-  Vector *ivec, *ovec;
+  int N;
+  Vector *ID1vec = NULL, *ID2vec = NULL, *SEQvec = NULL, *OUTvec = NULL;
 
   int VERBOSE = FALSE;
@@ -16,8 +13,18 @@
   }
 
-  Vector *cvec = NULL;
+  Vector *SEQdup = NULL;
+  if ((N = get_argument (argc, argv, "-save-dups"))) {
+    remove_argument (N, &argc, argv);
+    if ((SEQdup = SelectVector (argv[N], ANYVECTOR, TRUE)) == NULL) {
+      gprint (GP_ERR, "invalid vector %s\n", "SEQdups");
+      return FALSE;
+    }
+    remove_argument (N, &argc, argv);
+  }
+
+  Vector *CNTvec = NULL;
   if ((N = get_argument (argc, argv, "-c"))) {
     remove_argument (N, &argc, argv);
-    if ((cvec = SelectVector (argv[N], ANYVECTOR, TRUE)) == NULL) {
+    if ((CNTvec = SelectVector (argv[N], ANYVECTOR, TRUE)) == NULL) {
       gprint (GP_ERR, "invalid vector %s\n", argv[N]);
       return FALSE;
@@ -26,56 +33,68 @@
   }
 
-  if (argc != 3) {
-    gprint (GP_ERR, "USAGE: find_dups (ID1) (ID2) (sequence) (out) -c count\n");
+  if (argc != 5) {
+    gprint (GP_ERR, "USAGE: uniqpair (ID1) (ID2) (IDout) -c count\n");
     return (FALSE);
   }
 
   if ((ID1vec = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE);
-  if ((ID2vec = SelectVector (argv[2], ANYVECTOR, TRUE)) == NULL) return (FALSE);
-  if ((SEQvec = SelectVector (argv[3], ANYVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((ID2vec = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((OUTvec = SelectVector (argv[3], ANYVECTOR, TRUE)) == NULL) return (FALSE);
 
-  if ((ID1vec->type == OPIHI_FLT) || (ID2vec->type == OPIHI_FLT) || (ID1vec->type == OPIHI_FLT)) {
-    gprint (GP_ERR, "ERROR: ID and sequence vectors much be int\n");
+  if ((ID1vec->type == OPIHI_FLT) || (ID2vec->type == OPIHI_FLT)) {
+    gprint (GP_ERR, "ERROR: ID vectors much be int\n");
+    return (FALSE);
+  }
+  if (ID1vec->Nelements != ID2vec->Nelements) {
+    gprint (GP_ERR, "ERROR: ID vectors not all the same length\n");
     return (FALSE);
   }
 
-  /* allocate the maximum possible needed */
-  ResetVector (ovec, ivec->type, ivec->Nelements);
-  if (cvec) {
-    ResetVector (cvec, OPIHI_INT, ivec->Nelements);
+  // storage for the output and count vectors
+  int Nout = 0;
+  int NOUT = 10000;
+  ResetVector (OUTvec, OPIHI_INT, NOUT);
+  if (CNTvec) {
+    ResetVector (CNTvec, OPIHI_INT, NOUT);
   }
 
-  Nnew = 0;
-
-  // generate the joined vector
-
+  // storage for the duplicate sequences
+  int NseqDup = 0;
+  int NSEQDUP = 10000;
+  if (SEQdup) {
+    ResetVector (SEQdup, OPIHI_INT, NSEQDUP);
+  }
   
-  // copy the input data to a temporary array to avoid damaging it with sort
+  // incrementing pointers to the input IDs
   opihi_int *ID1 = ID1vec->elements.Int;
   opihi_int *ID2 = ID2vec->elements.Int;
 
+  // generate the joined and sequence vectors
   opihi_int *IDfull = NULL;
+  opihi_int *SEQdata = NULL;
   ALLOCATE (IDfull, opihi_int, ID1vec->Nelements);
-  for (i = 0; i < ID1vec->Nelements; i++) {
+  ALLOCATE (SEQdata, opihi_int, SEQvec[0].Nelements);
+  for (int i = 0; i < ID1vec->Nelements; i++) {
     IDfull[i] = ID1[i] + (ID2[i] << 32);
+    SEQdata[i] = i;
   }
 
-  llsorttwo (IDfull, SEQvec->elements.Int, ivec->Nelements);
+  llsorttwo (IDfull, SEQdata, ID1vec->Nelements);
 
-  Nnew = 0;
-  opihi_int *vtgt = ovec[0].elements.Int;
-
+  opihi_int *vtgt = OUTvec->elements.Int;
   opihi_int *vsrc = IDfull;
-  opihi_int *vidx = SEQvec->elements.Int;
 
   int onePercent = ID1vec->Nelements / 100;
 
-  /* ADD sequence to output
-     do not allocate all memory to start
-   */
-
   struct sigaction *old_sigaction = SetInterrupt();
-  for (i = 0; (i < ID1vec->Nelements) && !interrupt; Nnew++) {
-    vtgt[Nnew] = *vsrc;
+  for (int i = 0; (i < ID1vec->Nelements) && !interrupt; Nout++) {
+    if (Nout >= NOUT) {
+      NOUT += 10000;
+      REALLOCATE (vtgt, opihi_int, NOUT);
+      if (CNTvec) {
+	REALLOCATE (CNTvec->elements.Int, opihi_int, NOUT);
+      }
+    }
+    vtgt[Nout] = *vsrc;
     int Ndup = 0;
     opihi_int lastValue = *vsrc;
@@ -86,15 +105,35 @@
       if (VERBOSE && (i % onePercent == 0)) gprint (GP_ERR, ".");
     }
-    if (cvec) {
-      cvec->elements.Int[Nnew] = Ndup;
+    if (CNTvec) {
+      CNTvec->elements.Int[Nout] = Ndup;
+    }
+    if (SEQdup && (Ndup > 1)) {
+      // we want to save SEQdata values for the duplicates
+      // SEQdata[i-Ndup] .. SEQdata[i-N
+      for (j = i - Ndup; j < i; j++) {
+	SEQdup->elements.Int[NseqDup] = SEQdata[j];
+	NseqDup ++;
+	if (NseqDup >= NSEQDUP) {
+	  NSEQDUP += 10000;
+	  REALLOCATE (SEQdup->elements.Int, opihi_int, NSEQDUP);
+	}
+      }
     }
   }
   ClearInterrupt (old_sigaction);
   if (VERBOSE) gprint (GP_ERR, "\n");
-  free (indata);
 
-  // free up extra memory
-  ResetVector (ovec, ivec->type, Nnew);
-  if (cvec) ResetVector (cvec, OPIHI_INT, Nnew);
+  free (SEQdata);
+  free (IDfull);
+
+  if (SEQdup) {
+    SEQdup->Nelements = NseqDup;
+    REALLOCATE (SEQdup->elements.Int, opihi_int, Nseqdup);
+  }
+
+  // fix references and free up extra memory:
+  OUTvec->elements.Int = vtgt;
+  ResetVector (OUTvec, OPIHI_INT, Nout);
+  if (cvec) ResetVector (cvec, OPIHI_INT, Nout);
 
   return (TRUE);