IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 32909


Ignore:
Timestamp:
Dec 9, 2011, 3:34:31 AM (15 years ago)
Author:
eugene
Message:

adding spex2dgas

Location:
branches/eam_branches/ipp-20111122/Ohana/src/opihi/cmd.astro
Files:
3 edited

Legend:

Unmodified
Added
Removed

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

    r32867 r32909  
    5858$(SRC)/spec.$(ARCH).o              \
    5959$(SRC)/specpairfit.$(ARCH).o       \
     60$(SRC)/spex2dgas.$(ARCH).o         \
    6061$(SRC)/mkclusters.$(ARCH).o        \
    6162$(SRC)/star.$(ARCH).o              \

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

    r32867 r32909  
    4646int spec                    PROTO((int, char **));
    4747int specpairfit             PROTO((int, char **));
     48int spex2dgas               PROTO((int, char **));
    4849int mkclusters              PROTO((int, char **));
    4950int star                    PROTO((int, char **));
    … …  
    9697  {1, "spec",        spec,         "extract a spectrum"},
    9798  {1, "specpairfit", specpairfit,  "fit spectrum to another spectrum"},
     99  {1, "spex2dgas",   spex2dgas,    "minimize distances in 2D"},
    98100  {1, "mkclusters",  mkclusters,   "group spectra by distance"},
    99101  {1, "star",        star,         "star stats at rough coords"},

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

    r32895 r32909  
    11# include "astro.h"
    22
    3 /* We have N objects with N(N-1)/2 paired distances.
    4    We want to find the locus that best describes the observed distances
    5    Assume the distances are in a 2D space (though it could be N-D and non-euclidean)
     3/* We have N objects with N(N-1)/2 paired distances.  We want to find the 2D distribution that
     4   best describes the observed distances.  Assume the distances are in a 2D space, though in
     5   principle it could be N-D or even non-euclidean.
    66
    77   Model the points as a gas under pressure. 
    88
    9    Start with the real positions of the objects at some locations in 2D
     9   Start with the a guess for positions of the objects at some locations in 2D
    1010
    1111   P = (d_now - d_tru)
    … …  
    1717   Move based on dP
    1818
    19 */
     19   USAGE: spex2dgas incdex1 index2 distance
     20
     21   outline:
     22
     23   * load data
     24   * generate the unique objects
     25   * determine the max distance needed
     26   * place each object in the 2D space
     27
     28   * iterate:
     29   ** calculate dP/dX,dPdY for each object
     30   ** move each object proportionally to the pressure gradient
     31
     32   */
    2033
    2134typedef struct {
    22   int *friend; // indices for all tested relationships
    23   float *dist; // distance for this relationship
    24   int Ndist;   // number of tested relationships
    25   int NDIST;   // number of allocated relationships
    26   int group;   // assigned group
    27   int Nfriend; // number of friends brought into group
     35  int Nindex;   // number of tested relationships
     36  int NINDEX;   // number of allocated relationships
     37  int *index;  // indices for all other objects
     38  float *Dtgt; // target distance for this relationship
     39  float *Dcur; // current distance for this relationship
     40  float Xo,Yo; // current X,Y position of this object
     41  float dPdX;  // pressure in X
     42  float dPdY;  // pressure in Y
    2843} Object;
    29 
    30 typedef struct {
    31   int *entry;
    32   int Nentry;
    33   int NENTRY;
    34 } Group;
    35 
    36 static Group  *group = NULL;
    37 static int Ngroup = 0;
    38 static int NGROUP = 0;
    3944
    4045static Object *object = NULL;
    … …  
    4247// static int NOBJECT = 0;
    4348
    44 void add_to_group (int Nobj, float scale, int Nfriends) {
    45 
    46   int j;
    47 
    48   // fprintf (stderr, "add object %d to group %d\n", Nobj, Ngroup);
    49 
    50   object[Nobj].group = Ngroup;
    51 
    52   // add to this group
    53   int N = group[Ngroup].Nentry;
    54   group[Ngroup].entry[N] = Nobj;
    55   group[Ngroup].Nentry ++;
    56   CHECK_REALLOCATE (group[Ngroup].entry, int, group[Ngroup].NENTRY, group[Ngroup].Nentry, 100);
    57 
    58   // add all friends of this object (up to Nfriensd)
    59   // friends are already sorted by distance, so we add closest friends first
    60   for (j = 0; (object[Nobj].Nfriend < Nfriends) && (j < object[Nobj].Ndist); j++) {
    61     if (object[Nobj].dist[j] > scale) continue;
    62 
    63     int Nnew = object[Nobj].friend[j];
    64 
    65     if (object[Nnew].group != -1) continue;
    66 
    67     // found a friend :
    68     object[Nobj].Nfriend ++;
    69     add_to_group (Nnew, scale, Nfriends);
    70   }
     49void sortfriends (float *X, int *IDX1, int N);
     50
     51void get_pressure_gradient (int iObj) {
     52
     53  int i;
     54 
     55  float dPdX = 0.0;
     56  float dPdY = 0.0;
     57
     58  // only use the first N friends
     59  for (i = 0; (i < 10) && (i < object[iObj].Nindex); i++) {
     60    int jObj = object[iObj].index[i];
     61
     62    // only use the 0th object?
     63    // if (jObj != 0) continue;
     64
     65
     66    float Dtgt = object[iObj].Dtgt[i];
     67    float dX = object[jObj].Xo - object[iObj].Xo;
     68    float dY = object[jObj].Yo - object[iObj].Yo;
     69    float Dcur = hypot(dX,dY);
     70
     71    // the force law as a function of (Dcur - Dtgt) : if Dcur is too large, dF is negative
     72    // float dF = (Dcur < 0.01*Dtgt) ? -100.0 : (Dcur - Dtgt) / Dcur; XXX modified spring constant : too crazy
     73    float dF = (Dcur - Dtgt);
     74    dF = MIN (100.0, MAX (-100.0, dF));
     75
     76    float dPdXi = dF * dX / Dcur;
     77    float dPdYi = dF * dY / Dcur;
     78
     79    // if we are too close, then dX/Dcur is too ill-defined, just jump away
     80    if (Dcur < 0.01*Dtgt) {
     81      dPdXi = drand48() - 0.5;
     82      dPdYi = drand48() - 0.5;
     83    }
     84
     85    // fprintf (stderr, "Dcur,Dtgt : %f %f : dX,dY,dP : %f %f : %f : %f %f\n", Dcur, Dtgt, dX, dY, dF, dPdXi, dPdYi);
     86    dPdX += dPdXi;
     87    dPdY += dPdYi;
     88  }
     89  object[iObj].dPdX = dPdX;
     90  object[iObj].dPdY = dPdY;
     91
    7192  return;
    7293}
    7394
    74 void sortvecset (opihi_flt *X, opihi_int *IDX1, opihi_int *IDX2, int N) {
    75 
    76 # define SWAPFUNC(A,B){ opihi_flt tmp; opihi_int itmp; \
    77   tmp  = X[A];    X[A]    = X[B];    X[B]    = tmp; \
    78   itmp = IDX1[A]; IDX1[A] = IDX1[B]; IDX1[B] = itmp; \
    79   itmp = IDX2[A]; IDX2[A] = IDX2[B]; IDX2[B] = itmp; \
     95void move_object (int iObj) {
     96
     97  object[iObj].Xo += 0.25*object[iObj].dPdX;
     98  object[iObj].Yo += 0.25*object[iObj].dPdY;
     99  return;
    80100}
    81101
    82 # define COMPARE(A,B)(X[A] < X[B])
    83 
    84   OHANA_SORT (N, COMPARE, SWAPFUNC);
    85 
    86 # undef SWAPFUNC
    87 # undef COMPARE
    88 
    89 }
    90 
    91 void sortfriends (float *X, int *IDX1, int N) {
    92 
    93 # define SWAPFUNC(A,B){ float tmp; int itmp; \
    94   tmp  = X[A];    X[A]    = X[B];    X[B]    = tmp; \
    95   itmp = IDX1[A]; IDX1[A] = IDX1[B]; IDX1[B] = itmp; \
    96 }
    97 
    98 # define COMPARE(A,B)(X[A] < X[B])
    99 
    100   OHANA_SORT (N, COMPARE, SWAPFUNC);
    101 
    102 # undef SWAPFUNC
    103 # undef COMPARE
    104 
    105 }
    106 
    107 int mkclusters (int argc, char **argv) {
    108  
    109   int i, j;
     102int spex2dgas (int argc, char **argv) {
     103 
     104  int i, iter;
    110105  Vector *index1, *index2, *distance;
    111106
    112   if (argc != 6) goto usage;
     107  // init random numbers
     108  long A, B;
     109  A = time(NULL);
     110  for (B = 0; A == time(NULL); B++);
     111  srand48(B);
     112 
     113  if (argc != 5) goto usage;
    113114
    114115  if ((index1   = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) goto escape;
    115116  if ((index2   = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) goto escape;
    116117  if ((distance = SelectVector (argv[3], OLDVECTOR, TRUE)) == NULL) goto escape;
    117   float scale = atof (argv[4]);
    118   int Nfriends = atoi (argv[5]);
    119 
     118  int Niter = atoi (argv[4]);
    120119  // XXX enforce matching lengths on the three vectors
    121120
    … …  
    123122  CastVector (index2, OPIHI_INT);
    124123
     124  // how many objects do we have?
    125125  Nobject = 0;
    126126  for (i = 0; i < index1->Nelements; i++) {
    … …  
    133133  ALLOCATE (object, Object, Nobject);
    134134  for (i = 0; i < Nobject; i++) {
    135     ALLOCATE (object[i].friend, int, Nobject);
    136     memset (object[i].friend, 0, Nobject*sizeof(int));
    137     ALLOCATE (object[i].dist, float, Nobject);
    138     memset (object[i].dist, 0, Nobject*sizeof(float));
    139     object[i].NDIST = Nobject;
    140     object[i].Ndist = 0;
    141     object[i].group = -1; // unassigned
    142     object[i].Nfriend = 0; // no friends yet
     135    ALLOCATE (object[i].index, int, Nobject);
     136    memset (object[i].index, 0, Nobject*sizeof(int));
     137    ALLOCATE (object[i].Dtgt, float, Nobject);
     138    memset (object[i].Dtgt, 0, Nobject*sizeof(float));
     139    ALLOCATE (object[i].Dcur, float, Nobject);
     140    memset (object[i].Dcur, 0, Nobject*sizeof(float));
     141    object[i].NINDEX = Nobject;
     142    object[i].Nindex = 0;
    143143  }
    144144
    … …  
    147147  for (i = 0; i < index1->Nelements; i++) {
    148148    int N, m;
     149
    149150    N = index1->elements.Int[i];
    150151    if (N >= Nobject) abort();
    151     m = object[N].Ndist;
    152     object[N].friend[m] = index2->elements.Int[i];
    153     object[N].dist[m] = distance->elements.Flt[i];
    154     object[N].Ndist ++;
    155     if (object[N].Ndist == object[N].NDIST) {
    156       object[N].NDIST += 100;
    157       REALLOCATE (object[N].dist, float, object[N].NDIST);
    158       REALLOCATE (object[N].friend, int, object[N].NDIST);
     152    m = object[N].Nindex;
     153    object[N].index[m] = index2->elements.Int[i];
     154    object[N].Dtgt[m] = distance->elements.Flt[i];
     155    object[N].Dcur[m] = 0.0;
     156    object[N].Nindex ++;
     157    if (object[N].Nindex == object[N].NINDEX) {
     158      object[N].NINDEX += 100;
     159      REALLOCATE (object[N].Dtgt, float, object[N].NINDEX);
     160      REALLOCATE (object[N].Dcur, float, object[N].NINDEX);
     161      REALLOCATE (object[N].index, int, object[N].NINDEX);
    159162    }
    160163
    161164    N = index2->elements.Int[i];
    162165    if (N >= Nobject) abort();
    163     m = object[N].Ndist;
    164     object[N].friend[m] = index1->elements.Int[i];
    165     object[N].dist[m] = distance->elements.Flt[i];
    166     object[N].Ndist ++;
    167     if (object[N].Ndist == object[N].NDIST) {
    168       object[N].NDIST += 100;
    169       REALLOCATE (object[N].dist, float, object[N].NDIST);
    170       REALLOCATE (object[N].friend, int, object[N].NDIST);
    171     }
    172   }
    173 
    174   for (i = 0; i < Nobject; i++) {
    175     sortfriends (object[i].dist, object[i].friend, object[i].Ndist);
    176   }
    177 
    178   // generate the set of groups. each group is a list of friends and their friends
    179   Ngroup = 0;
    180   NGROUP = Nobject;
    181   ALLOCATE (group, Group, NGROUP);
    182   for (i = 0; i < NGROUP; i++) {
    183     group[i].Nentry = 0;
    184     group[i].NENTRY = 100;
    185     ALLOCATE (group[i].entry, int, group[i].NENTRY);
    186   }
    187 
    188   sortvecset (distance->elements.Flt, index1->elements.Int, index2->elements.Int, index2->Nelements);
    189  
     166    m = object[N].Nindex;
     167    object[N].index[m] = index1->elements.Int[i];
     168    object[N].Dtgt[m] = distance->elements.Flt[i];
     169    object[N].Dcur[m] = 0.0;
     170    object[N].Nindex ++;
     171    if (object[N].Nindex == object[N].NINDEX) {
     172      object[N].NINDEX += 100;
     173      REALLOCATE (object[N].Dtgt, float, object[N].NINDEX);
     174      REALLOCATE (object[N].Dcur, float, object[N].NINDEX);
     175      REALLOCATE (object[N].index, int, object[N].NINDEX);
     176    }
     177  }
     178
     179  for (i = 0; i < Nobject; i++) {
     180    // sort so closest friends are first
     181    sortfriends (object[i].Dtgt, object[i].index, object[i].Nindex);
     182  }
     183
     184  // find and save the max distance
     185  float Dmax = 0;
    190186  for (i = 0; i < distance->Nelements; i++) {
    191     int Nnew;
    192     int found;
    193 
    194     if (distance->elements.Flt[i] > scale) continue;
    195 
    196     found = FALSE;
    197 
    198     Nnew = index1->elements.Int[i];
    199     if (object[Nnew].group == -1) {
    200       add_to_group (Nnew, scale, Nfriends);
    201       // fprintf (stderr, "new group %d part 1 with %d elements\n", Ngroup, group[Ngroup].Nentry);
    202       found = TRUE;
    203     }
    204 
    205     Nnew = index2->elements.Int[i];
    206     if (object[Nnew].group == -1) {
    207       add_to_group (Nnew, scale, Nfriends);
    208       // fprintf (stderr, "new group %d part 1 with %d elements\n", Ngroup, group[Ngroup].Nentry);
    209       found = TRUE;
    210     }
    211 
    212     if (found) {
    213       fprintf (stderr, "group %d with %d elements\n", Ngroup, group[Ngroup].Nentry);
    214       Ngroup ++;
    215       if (Ngroup >= NGROUP) abort();
    216     }
    217   }
    218 
    219   if (1) {
    220     int Nmiss = 0;
     187    Dmax = MAX (Dmax, distance->elements.Flt[i]);
     188  }
     189
     190  // place the objects at the initial guess locations
     191  // XXX let's try with just a simple grid dividing up the max range
     192  // int Ngrid = sqrt(Nobject);
     193  // float dgrid = 1.5 * Dmax / Ngrid; // XXX remove the fudge factor
     194  for (i = 0; i < Nobject; i++) {
     195    // object[i].Xo = dgrid * (int) (i % Ngrid);
     196    // object[i].Yo = dgrid * (int) (i / Ngrid);
     197    object[i].Xo = Dmax*drand48();
     198    object[i].Yo = Dmax*drand48();
     199    object[i].dPdX = 0.0;
     200    object[i].dPdY = 0.0;
     201  }
     202
     203  for (iter = 0; iter < Niter; iter ++) {
     204    fprintf (stderr, "iter %d\n", iter);
     205
     206    // measure (dP/dX),(dP/dY) for all objects
    221207    for (i = 0; i < Nobject; i++) {
    222       if (object[i].group > -1) continue;
    223       Nmiss ++;
    224       // fprintf (stderr, "object %d not assigned\n", i);
    225       for (j = 0; FALSE && j < object[i].Ndist; j++) {
    226         fprintf (stderr, "friend %d, dist %f\n", object[i].friend[j], object[i].dist[j]);
    227       }
    228     }
    229     fprintf (stderr, "%d objects not assigned\n", Nmiss);
    230   }
     208      get_pressure_gradient (i);
     209      move_object (i);
     210    }
     211
     212    // given (dP/dX),(dP/dY), move each object
     213    // for (i = 0; i < Nobject; i++) {
     214    // }
     215  }
     216
     217  // save the result
     218  FILE *output = fopen ("output.dat", "w");
     219  for (i = 0; i < Nobject; i++) {
     220    fprintf (output, "%f %f : %f %f\n", object[i].Xo, object[i].Yo, object[i].dPdX, object[i].dPdY);
     221  }
     222  fclose (output);
    231223
    232224  return TRUE;
    233225 
    234  escape:
     226escape:
    235227  gprint (GP_ERR, "invalid vector\n");
    236228  return FALSE;
Note: See TracChangeset for help on using the changeset viewer.