Index: /branches/eam_branches/ipp-20111122/Ohana/src/opihi/cmd.astro/Makefile
===================================================================
--- /branches/eam_branches/ipp-20111122/Ohana/src/opihi/cmd.astro/Makefile	(revision 32908)
+++ /branches/eam_branches/ipp-20111122/Ohana/src/opihi/cmd.astro/Makefile	(revision 32909)
@@ -58,4 +58,5 @@
 $(SRC)/spec.$(ARCH).o		   \
 $(SRC)/specpairfit.$(ARCH).o	   \
+$(SRC)/spex2dgas.$(ARCH).o	   \
 $(SRC)/mkclusters.$(ARCH).o	   \
 $(SRC)/star.$(ARCH).o		   \
Index: /branches/eam_branches/ipp-20111122/Ohana/src/opihi/cmd.astro/init.c
===================================================================
--- /branches/eam_branches/ipp-20111122/Ohana/src/opihi/cmd.astro/init.c	(revision 32908)
+++ /branches/eam_branches/ipp-20111122/Ohana/src/opihi/cmd.astro/init.c	(revision 32909)
@@ -46,4 +46,5 @@
 int spec                    PROTO((int, char **));
 int specpairfit             PROTO((int, char **));
+int spex2dgas               PROTO((int, char **));
 int mkclusters              PROTO((int, char **));
 int star                    PROTO((int, char **));
@@ -96,4 +97,5 @@
   {1, "spec",        spec,         "extract a spectrum"},
   {1, "specpairfit", specpairfit,  "fit spectrum to another spectrum"},
+  {1, "spex2dgas",   spex2dgas,    "minimize distances in 2D"},
   {1, "mkclusters",  mkclusters,   "group spectra by distance"},
   {1, "star",        star,         "star stats at rough coords"},
Index: /branches/eam_branches/ipp-20111122/Ohana/src/opihi/cmd.astro/spex2dgas.c
===================================================================
--- /branches/eam_branches/ipp-20111122/Ohana/src/opihi/cmd.astro/spex2dgas.c	(revision 32908)
+++ /branches/eam_branches/ipp-20111122/Ohana/src/opihi/cmd.astro/spex2dgas.c	(revision 32909)
@@ -1,11 +1,11 @@
 # include "astro.h"
 
-/* We have N objects with N(N-1)/2 paired distances.
-   We want to find the locus that best describes the observed distances
-   Assume the distances are in a 2D space (though it could be N-D and non-euclidean)
+/* We have N objects with N(N-1)/2 paired distances.  We want to find the 2D distribution that
+   best describes the observed distances.  Assume the distances are in a 2D space, though in
+   principle it could be N-D or even non-euclidean.
 
    Model the points as a gas under pressure.  
 
-   Start with the real positions of the objects at some locations in 2D
+   Start with the a guess for positions of the objects at some locations in 2D
 
    P = (d_now - d_tru)
@@ -17,24 +17,29 @@
    Move based on dP
 
-*/
+   USAGE: spex2dgas incdex1 index2 distance
+
+   outline:
+
+   * load data
+   * generate the unique objects
+   * determine the max distance needed
+   * place each object in the 2D space
+
+   * iterate:
+   ** calculate dP/dX,dPdY for each object
+   ** move each object proportionally to the pressure gradient
+
+   */
 
 typedef struct {
-  int *friend; // indices for all tested relationships
-  float *dist; // distance for this relationship
-  int Ndist;   // number of tested relationships
-  int NDIST;   // number of allocated relationships
-  int group;   // assigned group
-  int Nfriend; // number of friends brought into group
+  int Nindex;   // number of tested relationships
+  int NINDEX;   // number of allocated relationships
+  int *index;  // indices for all other objects
+  float *Dtgt; // target distance for this relationship
+  float *Dcur; // current distance for this relationship
+  float Xo,Yo; // current X,Y position of this object
+  float dPdX;  // pressure in X
+  float dPdY;  // pressure in Y
 } Object;
-
-typedef struct {
-  int *entry;
-  int Nentry;
-  int NENTRY;
-} Group;
-
-static Group  *group = NULL;
-static int Ngroup = 0;
-static int NGROUP = 0;
 
 static Object *object = NULL;
@@ -42,80 +47,74 @@
 // static int NOBJECT = 0;
 
-void add_to_group (int Nobj, float scale, int Nfriends) {
-
-  int j;
-
-  // fprintf (stderr, "add object %d to group %d\n", Nobj, Ngroup);
-
-  object[Nobj].group = Ngroup;
-
-  // add to this group
-  int N = group[Ngroup].Nentry;
-  group[Ngroup].entry[N] = Nobj;
-  group[Ngroup].Nentry ++;
-  CHECK_REALLOCATE (group[Ngroup].entry, int, group[Ngroup].NENTRY, group[Ngroup].Nentry, 100);
-
-  // add all friends of this object (up to Nfriensd)
-  // friends are already sorted by distance, so we add closest friends first
-  for (j = 0; (object[Nobj].Nfriend < Nfriends) && (j < object[Nobj].Ndist); j++) {
-    if (object[Nobj].dist[j] > scale) continue;
-
-    int Nnew = object[Nobj].friend[j];
-
-    if (object[Nnew].group != -1) continue;
-
-    // found a friend : 
-    object[Nobj].Nfriend ++;
-    add_to_group (Nnew, scale, Nfriends);
-  }
+void sortfriends (float *X, int *IDX1, int N);
+
+void get_pressure_gradient (int iObj) {
+
+  int i;
+  
+  float dPdX = 0.0;
+  float dPdY = 0.0;
+
+  // only use the first N friends
+  for (i = 0; (i < 10) && (i < object[iObj].Nindex); i++) {
+    int jObj = object[iObj].index[i];
+
+    // only use the 0th object?
+    // if (jObj != 0) continue;
+
+
+    float Dtgt = object[iObj].Dtgt[i];
+    float dX = object[jObj].Xo - object[iObj].Xo;
+    float dY = object[jObj].Yo - object[iObj].Yo;
+    float Dcur = hypot(dX,dY);
+
+    // the force law as a function of (Dcur - Dtgt) : if Dcur is too large, dF is negative
+    // float dF = (Dcur < 0.01*Dtgt) ? -100.0 : (Dcur - Dtgt) / Dcur; XXX modified spring constant : too crazy
+    float dF = (Dcur - Dtgt);
+    dF = MIN (100.0, MAX (-100.0, dF));
+
+    float dPdXi = dF * dX / Dcur;
+    float dPdYi = dF * dY / Dcur;
+
+    // if we are too close, then dX/Dcur is too ill-defined, just jump away
+    if (Dcur < 0.01*Dtgt) {
+      dPdXi = drand48() - 0.5;
+      dPdYi = drand48() - 0.5;
+    }
+
+    // fprintf (stderr, "Dcur,Dtgt : %f %f : dX,dY,dP : %f %f : %f : %f %f\n", Dcur, Dtgt, dX, dY, dF, dPdXi, dPdYi);
+    dPdX += dPdXi;
+    dPdY += dPdYi;
+  }
+  object[iObj].dPdX = dPdX;
+  object[iObj].dPdY = dPdY;
+
   return;
 }
 
-void sortvecset (opihi_flt *X, opihi_int *IDX1, opihi_int *IDX2, int N) {
-
-# define SWAPFUNC(A,B){ opihi_flt tmp; opihi_int itmp; \
-  tmp  = X[A];    X[A]    = X[B];    X[B]    = tmp; \
-  itmp = IDX1[A]; IDX1[A] = IDX1[B]; IDX1[B] = itmp; \
-  itmp = IDX2[A]; IDX2[A] = IDX2[B]; IDX2[B] = itmp; \
+void move_object (int iObj) {
+
+  object[iObj].Xo += 0.25*object[iObj].dPdX;
+  object[iObj].Yo += 0.25*object[iObj].dPdY;
+  return;
 }
 
-# define COMPARE(A,B)(X[A] < X[B])
-
-  OHANA_SORT (N, COMPARE, SWAPFUNC);
-
-# undef SWAPFUNC
-# undef COMPARE
-
-}
-
-void sortfriends (float *X, int *IDX1, int N) {
-
-# define SWAPFUNC(A,B){ float tmp; int itmp; \
-  tmp  = X[A];    X[A]    = X[B];    X[B]    = tmp; \
-  itmp = IDX1[A]; IDX1[A] = IDX1[B]; IDX1[B] = itmp; \
-}
-
-# define COMPARE(A,B)(X[A] < X[B])
-
-  OHANA_SORT (N, COMPARE, SWAPFUNC);
-
-# undef SWAPFUNC
-# undef COMPARE
-
-}
-
-int mkclusters (int argc, char **argv) {
-  
-  int i, j;
+int spex2dgas (int argc, char **argv) {
+  
+  int i, iter;
   Vector *index1, *index2, *distance;
 
-  if (argc != 6) goto usage;
+  // init random numbers
+  long A, B;
+  A = time(NULL);
+  for (B = 0; A == time(NULL); B++);
+  srand48(B);
+ 
+  if (argc != 5) goto usage;
 
   if ((index1   = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) goto escape;
   if ((index2   = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) goto escape;
   if ((distance = SelectVector (argv[3], OLDVECTOR, TRUE)) == NULL) goto escape;
-  float scale = atof (argv[4]);
-  int Nfriends = atoi (argv[5]);
-
+  int Niter = atoi (argv[4]);
   // XXX enforce matching lengths on the three vectors
 
@@ -123,4 +122,5 @@
   CastVector (index2, OPIHI_INT);
 
+  // how many objects do we have?
   Nobject = 0;
   for (i = 0; i < index1->Nelements; i++) {
@@ -133,12 +133,12 @@
   ALLOCATE (object, Object, Nobject);
   for (i = 0; i < Nobject; i++) {
-    ALLOCATE (object[i].friend, int, Nobject);
-    memset (object[i].friend, 0, Nobject*sizeof(int));
-    ALLOCATE (object[i].dist, float, Nobject);
-    memset (object[i].dist, 0, Nobject*sizeof(float));
-    object[i].NDIST = Nobject;
-    object[i].Ndist = 0;
-    object[i].group = -1; // unassigned
-    object[i].Nfriend = 0; // no friends yet
+    ALLOCATE (object[i].index, int, Nobject);
+    memset (object[i].index, 0, Nobject*sizeof(int));
+    ALLOCATE (object[i].Dtgt, float, Nobject);
+    memset (object[i].Dtgt, 0, Nobject*sizeof(float));
+    ALLOCATE (object[i].Dcur, float, Nobject);
+    memset (object[i].Dcur, 0, Nobject*sizeof(float));
+    object[i].NINDEX = Nobject;
+    object[i].Nindex = 0;
   }
 
@@ -147,90 +147,82 @@
   for (i = 0; i < index1->Nelements; i++) {
     int N, m;
+
     N = index1->elements.Int[i];
     if (N >= Nobject) abort();
-    m = object[N].Ndist;
-    object[N].friend[m] = index2->elements.Int[i];
-    object[N].dist[m] = distance->elements.Flt[i];
-    object[N].Ndist ++;
-    if (object[N].Ndist == object[N].NDIST) {
-      object[N].NDIST += 100;
-      REALLOCATE (object[N].dist, float, object[N].NDIST);
-      REALLOCATE (object[N].friend, int, object[N].NDIST);
+    m = object[N].Nindex;
+    object[N].index[m] = index2->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);
     }
 
     N = index2->elements.Int[i];
     if (N >= Nobject) abort();
-    m = object[N].Ndist;
-    object[N].friend[m] = index1->elements.Int[i];
-    object[N].dist[m] = distance->elements.Flt[i];
-    object[N].Ndist ++;
-    if (object[N].Ndist == object[N].NDIST) {
-      object[N].NDIST += 100;
-      REALLOCATE (object[N].dist, float, object[N].NDIST);
-      REALLOCATE (object[N].friend, int, object[N].NDIST);
-    }
-  }
-
-  for (i = 0; i < Nobject; i++) {
-    sortfriends (object[i].dist, object[i].friend, object[i].Ndist);
-  }
-
-  // generate the set of groups. each group is a list of friends and their friends
-  Ngroup = 0;
-  NGROUP = Nobject;
-  ALLOCATE (group, Group, NGROUP);
-  for (i = 0; i < NGROUP; i++) {
-    group[i].Nentry = 0;
-    group[i].NENTRY = 100;
-    ALLOCATE (group[i].entry, int, group[i].NENTRY);
-  }
-
-  sortvecset (distance->elements.Flt, index1->elements.Int, index2->elements.Int, index2->Nelements);
-  
+    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);
+    }
+  }
+
+  for (i = 0; i < Nobject; i++) {
+    // sort so closest friends are first
+    sortfriends (object[i].Dtgt, object[i].index, object[i].Nindex);
+  }
+
+  // find and save the max distance
+  float Dmax = 0;
   for (i = 0; i < distance->Nelements; i++) {
-    int Nnew;
-    int found;
-
-    if (distance->elements.Flt[i] > scale) continue;
-
-    found = FALSE;
-
-    Nnew = index1->elements.Int[i];
-    if (object[Nnew].group == -1) {
-      add_to_group (Nnew, scale, Nfriends);
-      // fprintf (stderr, "new group %d part 1 with %d elements\n", Ngroup, group[Ngroup].Nentry);
-      found = TRUE;
-    }
-
-    Nnew = index2->elements.Int[i];
-    if (object[Nnew].group == -1) {
-      add_to_group (Nnew, scale, Nfriends);
-      // fprintf (stderr, "new group %d part 1 with %d elements\n", Ngroup, group[Ngroup].Nentry);
-      found = TRUE;
-    }
-
-    if (found) {
-      fprintf (stderr, "group %d with %d elements\n", Ngroup, group[Ngroup].Nentry);
-      Ngroup ++;
-      if (Ngroup >= NGROUP) abort();
-    }
-  }
-
-  if (1) {
-    int Nmiss = 0;
+    Dmax = MAX (Dmax, distance->elements.Flt[i]);
+  }
+
+  // place the objects at the initial guess locations
+  // XXX let's try with just a simple grid dividing up the max range
+  // int Ngrid = sqrt(Nobject);
+  // float dgrid = 1.5 * Dmax / Ngrid; // XXX remove the fudge factor
+  for (i = 0; i < Nobject; i++) {
+    // object[i].Xo = dgrid * (int) (i % Ngrid);
+    // object[i].Yo = dgrid * (int) (i / Ngrid);
+    object[i].Xo = Dmax*drand48();
+    object[i].Yo = Dmax*drand48();
+    object[i].dPdX = 0.0;
+    object[i].dPdY = 0.0;
+  }
+
+  for (iter = 0; iter < Niter; iter ++) {
+    fprintf (stderr, "iter %d\n", iter);
+
+    // measure (dP/dX),(dP/dY) for all objects
     for (i = 0; i < Nobject; i++) {
-      if (object[i].group > -1) continue;
-      Nmiss ++;
-      // fprintf (stderr, "object %d not assigned\n", i);
-      for (j = 0; FALSE && j < object[i].Ndist; j++) {
-	fprintf (stderr, "friend %d, dist %f\n", object[i].friend[j], object[i].dist[j]);
-      }
-    }
-    fprintf (stderr, "%d objects not assigned\n", Nmiss);
-  } 
+      get_pressure_gradient (i);
+      move_object (i);
+    }
+
+    // given (dP/dX),(dP/dY), move each object 
+    // for (i = 0; i < Nobject; i++) {
+    // }
+  }
+
+  // save the result
+  FILE *output = fopen ("output.dat", "w");
+  for (i = 0; i < Nobject; i++) {
+    fprintf (output, "%f %f : %f %f\n", object[i].Xo, object[i].Yo, object[i].dPdX, object[i].dPdY);
+  }
+  fclose (output);
 
   return TRUE;
   
- escape: 
+escape: 
   gprint (GP_ERR, "invalid vector\n");
   return FALSE;