Index: /trunk/Ohana/src/opihi/Makefile =================================================================== --- /trunk/Ohana/src/opihi/Makefile (revision 33661) +++ /trunk/Ohana/src/opihi/Makefile (revision 33662) @@ -9,4 +9,5 @@ dimm: lib.data lib.shell cmd.basic cmd.data cmd.astro dvo: lib.data lib.shell cmd.basic cmd.data cmd.astro +dvo_client: lib.data lib.shell cmd.basic cmd.data cmd.astro pantasks: lib.data lib.shell cmd.basic cmd.data cmd.astro pantasks_server: lib.data lib.shell cmd.basic cmd.data cmd.astro @@ -18,5 +19,5 @@ PROGRAM = mana dvo pantasks pclient pcontrol -SPECIAL = pantasks_client pantasks_server +SPECIAL = pantasks_client pantasks_server dvo_client EXTRAS = dimm @@ -24,6 +25,4 @@ all: for i in $(PROGRAM) $(SPECIAL); do $(MAKE) $$i || exit; done - $(MAKE) pantasks_client - $(MAKE) pantasks_server libs: @@ -35,6 +34,4 @@ install: for i in $(PROGRAM) $(SPECIAL); do $(MAKE) $$i.install || exit; done - $(MAKE) pantasks_client.install - $(MAKE) pantasks_server.install extras-install: @@ -63,4 +60,10 @@ if [ -d pantasks ]; then (cd pantasks && $(MAKE) $@); fi +dvo_client: + if [ -d dvo ]; then (cd dvo && $(MAKE) $@); fi + +dvo_client.install: dvo_client + if [ -d dvo ]; then (cd dvo && $(MAKE) $@); fi + $(PROGRAM) $(LIBS) $(EXTRAS): if [ -d "$@" ]; then (cd $@ && $(MAKE)); fi Index: /trunk/Ohana/src/opihi/cmd.astro/Makefile =================================================================== --- /trunk/Ohana/src/opihi/cmd.astro/Makefile (revision 33661) +++ /trunk/Ohana/src/opihi/cmd.astro/Makefile (revision 33662) @@ -32,7 +32,8 @@ $(SRC)/fixwrap.$(ARCH).o \ $(SRC)/fixcols.$(ARCH).o \ +$(SRC)/fiximage.$(ARCH).o \ $(SRC)/gauss.$(ARCH).o \ -$(SRC)/getvel.$(ARCH).o \ -$(SRC)/getlst.$(ARCH).o \ +$(SRC)/getvel.$(ARCH).o \ +$(SRC)/getlst.$(ARCH).o \ $(SRC)/medianmap.$(ARCH).o \ $(SRC)/mkgauss.$(ARCH).o \ @@ -51,12 +52,17 @@ $(SRC)/precess.$(ARCH).o \ $(SRC)/profile.$(ARCH).o \ -$(SRC)/radec.$(ARCH).o \ -$(SRC)/region.$(ARCH).o \ +$(SRC)/radec.$(ARCH).o \ +$(SRC)/region.$(ARCH).o \ $(SRC)/rotcurve.$(ARCH).o \ $(SRC)/scale.$(ARCH).o \ $(SRC)/sexigesimal.$(ARCH).o \ $(SRC)/spec.$(ARCH).o \ +$(SRC)/specpairfit.$(ARCH).o \ +$(SRC)/spexseq.$(ARCH).o \ +$(SRC)/spex1dgas.$(ARCH).o \ +$(SRC)/spex2dgas.$(ARCH).o \ +$(SRC)/mkclusters.$(ARCH).o \ $(SRC)/star.$(ARCH).o \ -$(SRC)/transform.$(ARCH).o \ +$(SRC)/transform.$(ARCH).o \ $(SRC)/imsub.$(ARCH).o \ $(SRC)/imfit.$(ARCH).o \ Index: /trunk/Ohana/src/opihi/cmd.astro/coords.c =================================================================== --- /trunk/Ohana/src/opihi/cmd.astro/coords.c (revision 33661) +++ /trunk/Ohana/src/opihi/cmd.astro/coords.c (revision 33662) @@ -83,4 +83,5 @@ } if (mode == PIXEL) { + Xin = ohana_normalize_angle_to_midpoint (Xin, coords.crval1); RD_to_XY (&Xout, &Yout, Xin, Yin, &coords); if (!Quiet) gprint (GP_LOG, "%7.2f %7.2f\n", Xout, Yout); @@ -92,5 +93,9 @@ if (mode == SKY) { for (i = 0; i < xvec[0].Nelements; i++) { - XY_to_RD (&xvec[0].elements.Flt[i], &yvec[0].elements.Flt[i], xvec[0].elements.Flt[i], yvec[0].elements.Flt[i], &coords); + double Xin = xvec[0].elements.Flt[i]; + double Yin = yvec[0].elements.Flt[i]; + xvec[0].elements.Flt[i] = NAN; + yvec[0].elements.Flt[i] = NAN; + XY_to_RD (&xvec[0].elements.Flt[i], &yvec[0].elements.Flt[i], Xin, Yin, &coords); } return (TRUE); @@ -98,5 +103,10 @@ if (mode == PIXEL) { for (i = 0; i < xvec[0].Nelements; i++) { - RD_to_XY (&xvec[0].elements.Flt[i], &yvec[0].elements.Flt[i], xvec[0].elements.Flt[i], yvec[0].elements.Flt[i], &coords); + double Xin = xvec[0].elements.Flt[i]; + double Yin = yvec[0].elements.Flt[i]; + xvec[0].elements.Flt[i] = NAN; + yvec[0].elements.Flt[i] = NAN; + Xin = ohana_normalize_angle_to_midpoint (Xin, coords.crval1); + RD_to_XY (&xvec[0].elements.Flt[i], &yvec[0].elements.Flt[i], Xin, Yin, &coords); } return (TRUE); Index: /trunk/Ohana/src/opihi/cmd.astro/fiximage.c =================================================================== --- /trunk/Ohana/src/opihi/cmd.astro/fiximage.c (revision 33662) +++ /trunk/Ohana/src/opihi/cmd.astro/fiximage.c (revision 33662) @@ -0,0 +1,91 @@ +# include "astro.h" + +int fiximage (int argc, char **argv) { + + int ix, iy, N; + Buffer *in, *ct, *mask; + + int VERBOSE = FALSE; + if ((N = get_argument (argc, argv, "-v"))) { + VERBOSE = TRUE; + remove_argument (N, &argc, argv); + } + + if (argc != 4) { + gprint (GP_ERR, "USAGE: fiximage (mask)\n"); + return (FALSE); + } + + if ((in = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE); + if ((ct = SelectBuffer (argv[2], OLDBUFFER, TRUE)) == NULL) return (FALSE); + if ((mask = SelectBuffer (argv[3], ANYBUFFER, TRUE)) == NULL) return (FALSE); + + int Nx = in[0].matrix.Naxis[0]; + int Ny = in[0].matrix.Naxis[1]; + + gfits_free_matrix (&mask[0].matrix); + gfits_free_header (&mask[0].header); + CreateBuffer (mask, Nx, Ny, -32, 0.0, 1.0); + strcpy (mask[0].file, "(empty)"); + memset (mask[0].matrix.buffer, 0, Nx*Ny*sizeof(float)); + // 0 -- init, cannot be repaired + // 1 -- does not need to be repaired + // 2 -- repair x only + // 4 -- repair y only + // 6 -- repair x & y only + + float *inB = (float *) in[0].matrix.buffer; + float *ctB = (float *) ct[0].matrix.buffer; + float *mkB = (float *) mask[0].matrix.buffer; + + // loop over input ct and find pixels that can be repaired. + for (iy = 1; iy < Ny - 1; iy++) { + for (ix = 1; ix < Nx - 1; ix++) { + + int Npix = ix + Nx*iy; + + // does not need to be repaired + if (ctB[Npix]) { mkB[Npix] = 1.0; continue; } + + // if we have all 4 valid neighbors, we can fix + if (ctB[Npix - 1 ] && ctB[Npix + 1 ]) mkB[Npix] += 2.0; + if (ctB[Npix - Nx] && ctB[Npix + Nx]) mkB[Npix] += 4.0; + } + } + + // loop over input ct and repair the reparable pixels + for (iy = 1; iy < Ny - 1; iy++) { + for (ix = 1; ix < Nx - 1; ix++) { + + int Npix = ix + Nx*iy; + + if (mkB[Npix] == 0.0) continue; // cannot be repaired + if (mkB[Npix] == 1.0) continue; // does not need to be repaired + + float Vxm, Vxp, Vym, Vyp, Value; + if ((mkB[Npix] == 2) || (mkB[Npix] == 6)) { + Vxm = inB[Npix - 1 ] / ctB[Npix - 1 ]; + Vxp = inB[Npix + 1 ] / ctB[Npix + 1 ]; + } + if ((mkB[Npix] == 4) || (mkB[Npix] == 6)) { + Vym = inB[Npix - Nx] / ctB[Npix - Nx]; + Vyp = inB[Npix + Nx] / ctB[Npix + Nx]; + } + + if (mkB[Npix] == 2) { + Value = 0.5*(Vxm + Vxp); + } + if (mkB[Npix] == 4) { + Value = 0.5*(Vym + Vyp); + } + if (mkB[Npix] == 6) { + Value = 0.25*(Vxm + Vxp + Vym + Vyp); + } + + ctB[Npix] = 1; + inB[Npix] = Value; + } + } + + return TRUE; +} Index: /trunk/Ohana/src/opihi/cmd.astro/init.c =================================================================== --- /trunk/Ohana/src/opihi/cmd.astro/init.c (revision 33661) +++ /trunk/Ohana/src/opihi/cmd.astro/init.c (revision 33662) @@ -15,4 +15,5 @@ int flux PROTO((int, char **)); int fixwrap PROTO((int, char **)); +int fiximage PROTO((int, char **)); int fixcols PROTO((int, char **)); int fixrows PROTO((int, char **)); @@ -45,4 +46,9 @@ int sexigesimal PROTO((int, char **)); int spec PROTO((int, char **)); +int specpairfit PROTO((int, char **)); +int spexseq PROTO((int, char **)); +int spex1dgas PROTO((int, char **)); +int spex2dgas PROTO((int, char **)); +int mkclusters PROTO((int, char **)); int star PROTO((int, char **)); int times PROTO((int, char **)); @@ -64,4 +70,5 @@ {1, "flux", flux, "flux in a convex contour"}, {1, "fixwrap", fixwrap, "fix megacam over-wrapped pixels"}, + {1, "fiximage", fiximage, "fix pixels in an image by interpolation"}, {1, "fixcols", fixcols, "fix bad columns by comparing with others"}, {1, "fixrows", fixrows, "fix bad rows by comparing with others"}, @@ -93,4 +100,9 @@ {1, "sexigesimal", sexigesimal, "convert to/from sexigesimal/decimal"}, {1, "spec", spec, "extract a spectrum"}, + {1, "specpairfit", specpairfit, "fit spectrum to another spectrum"}, + {1, "spexseq", spexseq, "generate the spectral sequence"}, + {1, "spex1dgas", spex1dgas, "minimize distances in 1D"}, + {1, "spex2dgas", spex2dgas, "minimize distances in 2D"}, + {1, "mkclusters", mkclusters, "group spectra by distance"}, {1, "star", star, "star stats at rough coords"}, {1, "transform", transform, "geometric transformation of image"}, Index: /trunk/Ohana/src/opihi/cmd.astro/mkclusters.c =================================================================== --- /trunk/Ohana/src/opihi/cmd.astro/mkclusters.c (revision 33662) +++ /trunk/Ohana/src/opihi/cmd.astro/mkclusters.c (revision 33662) @@ -0,0 +1,233 @@ +# include "astro.h" + +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 +} 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; +static int Nobject = 0; +// 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); + } + 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; \ +} + +# 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; + Vector *index1, *index2, *distance; + + if (argc != 6) 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]); + + // XXX enforce matching lengths on the three vectors + + CastVector (index1, OPIHI_INT); + CastVector (index2, OPIHI_INT); + + Nobject = 0; + for (i = 0; i < index1->Nelements; i++) { + Nobject = MAX (Nobject, index1->elements.Int[i]); + Nobject = MAX (Nobject, index2->elements.Int[i]); + } + Nobject ++; // after the loop, Nobject has the value of the highest index, not the count + + // allocate the list of Object -- these list the possible friends of the given object[i] + 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 + } + + // generate the set of vectors of all distances for each entry + // the object seq number (object[i]) will be used to match the index1 and index2 values + 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); + } + + 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); + + 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; + 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); + } + + return TRUE; + + escape: + gprint (GP_ERR, "invalid vector\n"); + return FALSE; + +usage: + gprint (GP_ERR, "USAGE: mkclusters (index1) (index2) (distance) scale\n"); + return FALSE; +} + +/* this function takes 3 vectors: a set of distances and a pair of indicies identifying the end points. + + the goal is to generate groups of the end points that are relatively closer. + + note that the distances may not be Euclidean: nothing can be assumed about relationships between d(1,2), d(2,3), and d(1,3) + + the indicies should be sequential. the algorithm does not require it, but storage is much cleaner if it is + +*/ Index: /trunk/Ohana/src/opihi/cmd.astro/specpairfit.c =================================================================== --- /trunk/Ohana/src/opihi/cmd.astro/specpairfit.c (revision 33662) +++ /trunk/Ohana/src/opihi/cmd.astro/specpairfit.c (revision 33662) @@ -0,0 +1,61 @@ +# include "astro.h" + +int specpairfit (int argc, char **argv) { + + int i; + Vector *flux1, *flux2, *dflux1, *dflux2, *window; + + if (argc != 7) goto usage; + + if ((flux1 = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) goto escape; + if ((dflux1 = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) goto escape; + if ((flux2 = SelectVector (argv[3], OLDVECTOR, TRUE)) == NULL) goto escape; + if ((dflux2 = SelectVector (argv[4], OLDVECTOR, TRUE)) == NULL) goto escape; + if ((window = SelectVector (argv[5], OLDVECTOR, TRUE)) == NULL) goto escape; + + // XXX enforce matching lengths on the 6 vectors + + int mask = atoi (argv[6]); + CastVector (window, OPIHI_INT); + + // minimize (flux1 - flux2*alpha) in window defined by mask + // note that the mask is a SELECTION mask not an EXCLUSION mask + + double F12 = 0.0; + double F22 = 0.0; + for (i = 0; i < flux1->Nelements; i++) { + if ((mask & window->elements.Int[i]) == 0) continue; + double weight = 1.0 / (SQ(dflux1->elements.Flt[i]) + SQ(dflux2->elements.Flt[i])); + F12 += flux1->elements.Flt[i] * flux2->elements.Flt[i] * weight; + F22 += flux2->elements.Flt[i] * flux2->elements.Flt[i] * weight; + } + + double Ao = F12 / F22; + double dA = sqrt(1.0 / F22); + + int Ndof = -1; // 1 parameter fit + double chisq = 0.0; + for (i = 0; i < flux1->Nelements; i++) { + if ((mask & window->elements.Int[i]) == 0) continue; + double weight = 1.0 / (SQ(dflux1->elements.Flt[i]) + SQ(dflux2->elements.Flt[i])); + chisq += SQ(flux1->elements.Flt[i] - Ao * flux2->elements.Flt[i]) * weight; + Ndof ++; + } + + double chisqNu = chisq / Ndof; + + // fprintf (stderr, "Ao: %f +/- %f, chisq: %f, chisq_nu : %f for %d dof\n", Ao, dA, chisq, chisqNu, Ndof); + set_variable ("Ao", Ao); + set_variable ("dA", dA); + set_variable ("Xv", chisqNu); + set_variable ("Nd", Ndof); + return (TRUE); + + escape: + gprint (GP_ERR, "invalid vector\n"); + return (FALSE); + +usage: + gprint (GP_ERR, "USAGE: specpairfit (flux1) (dflux1) (flux2) (dflux2) (window) (mask) [options]\n"); + return FALSE; +} Index: /trunk/Ohana/src/opihi/cmd.astro/spex1dgas.c =================================================================== --- /trunk/Ohana/src/opihi/cmd.astro/spex1dgas.c (revision 33662) +++ /trunk/Ohana/src/opihi/cmd.astro/spex1dgas.c (revision 33662) @@ -0,0 +1,438 @@ +# include "astro.h" + +/* We have N objects with N(N-1)/2 paired distances. We want to find the 1D distribution that + best describes the observed distances. Assume the distances are in a 1D space, though in + principle it could be N-D or even non-euclidean. + + Model the points as a gas under pressure. + + Start with the a guess for positions of the objects at some locations in 2D + + P = (d_now - d_tru) + + Iterate over the points and move in X,Y + + Find (dP/dX,dP/dY) for modest moves in X,Y + + Move based on dP + + USAGE: spex1dgas 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 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; + +static Object *object = NULL; +static int Nobject = 0; +// static int NOBJECT = 0; + +void sortfriends (float *X, int *IDX1, int N); + +static void get_pressure_gradient (int iObj, int nCloseMax, float farFrac, int nearNeighbors, float maxPressure) { + + int i; + + float dPdX = 0.0; + + // only use the first N friends + for (i = 0; (i < object[iObj].Nindex); i++) { + int jObj = object[iObj].index[i]; + + // some options: + // if the iterations are small, we should only worry about getting the near neighbors right + // if the iterations are large, we should add in more distant objects + if (nearNeighbors && (i >= nCloseMax)) break; + if (!nearNeighbors && (i >= nCloseMax)) { + if (drand48() > farFrac) continue; + } + + float Dtgt = object[iObj].Dtgt[i]; + float dX = object[jObj].Xo - object[iObj].Xo; + float Dcur = dX; + + // 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 (maxPressure, MAX (-maxPressure, dF)); + + float dPdXi; + if (fabs(Dcur) < 1e-6) { + dPdXi = 0.0; + } else { + dPdXi = dF; + } + + if (isnan(Dtgt) || isnan(dX) || isnan(Dcur) || isnan(dF) || isnan(dPdXi)) abort(); + dPdX += dPdXi; + } + object[iObj].dPdX = dPdX; + object[iObj].dPdY = 0.0; + + return; +} + +static void move_object (int iObj) { + + object[iObj].Xo += 0.25*object[iObj].dPdX; + return; +} + +int spex1dgas (int argc, char **argv) { + + int i, j, iter; + Vector *index1, *index2, *distance; + float *XoList; + int *IDList, *MidObj; + + MidObj = NULL; + IDList = NULL; + XoList = NULL; + { + // init random numbers + long A, B; + A = time(NULL); + for (B = 0; A == time(NULL); B++); + srand48(B); + } + + if (argc != 11) goto usage; + + // XXX enforce matching lengths on the three vectors + 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; + + int Niter = atoi (argv[4]); + int nCloseMax = atoi (argv[5]); + int nCloseIter = atoi (argv[6]); + float farFrac = atof (argv[7]); + float maxPressure = atof (argv[8]); + + int idx1 = atoi(argv[9]); + int idx2 = atoi(argv[10]); + + CastVector (index1, OPIHI_INT); + CastVector (index2, OPIHI_INT); + + // how many objects do we have? + Nobject = 0; + for (i = 0; i < index1->Nelements; i++) { + Nobject = MAX (Nobject, index1->elements.Int[i]); + Nobject = MAX (Nobject, index2->elements.Int[i]); + } + Nobject ++; // after the loop, Nobject has the value of the highest index, not the count + + // allocate the list of Object -- these list the possible friends of the given object[i] + ALLOCATE (object, Object, Nobject); + for (i = 0; i < Nobject; i++) { + 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; + } + + // generate the set of vectors of all distances for each entry + // the object seq number (object[i]) will be used to match the index1 and index2 values + for (i = 0; i < index1->Nelements; i++) { + int N, m; + + N = index1->elements.Int[i]; + if (N >= Nobject) abort(); + 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); + } + + // 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 } + } + + 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++) { + 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 + // XXX for (i = 0; i < Nobject; i++) { + // XXX // object[i].Xo = dgrid * (int) (i % Ngrid); + // XXX // object[i].Yo = dgrid * (int) (i / Ngrid); + // XXX object[i].Xo = Dmax*drand48(); + // XXX object[i].Yo = Dmax*drand48(); + // XXX object[i].dPdX = 0.0; + // XXX object[i].dPdY = 0.0; + // XXX } + + // place the objects at the initial guess locations + // use object 0 and its most distant friend to constrain: + + // we are going to use the projected X distance as the location in X, and ignore the Y + // coordinate + + // need to get pretty close to a good starting point. start with a guess using the first entry & its extremum + // int idx1 = object[0].index[(int)(0.75*object[0].Nindex)]; + // int idx2 = object[idx1].index[(int)(0.75*object[idx1].Nindex)]; + float A = NAN; + for (i = 0; isnan (A) && (i < object[idx2].Nindex); i++) { + if (object[idx2].index[i] == idx1) { + A = object[idx2].Dtgt[i]; + } + } + fprintf (stderr, "idx1: %d, idx2: %d, A: %f\n", idx1, idx2, A); + set_variable ("SP1D_A", A); + object[idx1].Xo = 0.0; + object[idx1].Yo = 0.0; + object[idx1].dPdX = 0.0; + object[idx1].dPdY = 0.0; + + object[idx2].Xo = A; + object[idx2].Yo = 0.0; + object[idx2].dPdX = 0.0; + object[idx2].dPdY = 0.0; + + Vector *tmpB = SelectVector ("sp1d_B", ANYVECTOR, TRUE); + Vector *tmpC = SelectVector ("sp1d_C", ANYVECTOR, TRUE); + ResetVector (tmpB, OPIHI_FLT, Nobject); + ResetVector (tmpC, OPIHI_FLT, Nobject); + + for (i = 0; i < Nobject; i++) { + if (i == idx1) continue; + if (i == idx2) continue; + float B = NAN; + float C = NAN; + for (j = 0; (isnan(B) || isnan(C)) && (j < object[i].Nindex); j++) { + if (object[i].index[j] == idx1) { B = object[i].Dtgt[j]; } + if (object[i].index[j] == idx2) { C = object[i].Dtgt[j]; } + } + if (isnan(B) || isnan(C)) abort(); + + // C^2 = A^2 + B^2 - 2AB cos(t) + // Xo = B cos(t) + // Xo = A^2 + B^2 - C^2 / 2 A + float Xo = (SQ(A) + SQ(B) - SQ(C)) / (2*A); + + object[i].Xo = Xo; + object[i].dPdX = 0.0; + object[i].Yo = 0.0; + object[i].dPdY = 0.0; + tmpB->elements.Flt[i] = B; + tmpC->elements.Flt[i] = C; + } + tmpB->elements.Flt[idx1] = 0; + tmpB->elements.Flt[idx2] = A; + tmpC->elements.Flt[idx1] = A; + tmpC->elements.Flt[idx2] = 0; + + if (Niter >= -1) { + // now choose several in the middle of the range, and find the mean distance for each relative to that group + ALLOCATE (XoList, float, Nobject); + ALLOCATE (IDList, int, Nobject); + for (i = 0; i < Nobject; i++) { + XoList[i] = object[i].Xo; + IDList[i] = i; + } + sortfriends (XoList, IDList, Nobject); + float XoMid = XoList[(int)(Nobject/2)]; + +# define NMID 5 + ALLOCATE (MidObj, int, NMID); + MidObj[0] = IDList[(int)(Nobject/2) + 0]; + MidObj[1] = IDList[(int)(Nobject/2) + 1]; + MidObj[2] = IDList[(int)(Nobject/2) - 1]; + MidObj[3] = IDList[(int)(Nobject/2) + 2]; + MidObj[4] = IDList[(int)(Nobject/2) - 2]; + + for (i = 0; i < Nobject; i++) { + for (j = 0; j < NMID; j++) { + if (MidObj[j] == i) goto skip_object; + } + float Dmean = 0.0; + for (j = 0; j < NMID; j++) { + Dmean += object[i].Dtgt[MidObj[j]]; // desired distance between object (i) and object in middle section + } + Dmean = Dmean / NMID; + if (object[i].Xo < XoMid) { + object[i].Xo = XoMid - Dmean; + } else { + object[i].Xo = XoMid + Dmean; + } + skip_object: + continue; + } + } + + if (Niter >= 0) { + float XoMean = 0.0; + // now choose one end of the range, and find the mean distance for each relative to THAT group + for (i = 0; i < Nobject; i++) { + XoList[i] = object[i].Xo; + IDList[i] = i; + XoMean += object[i].Xo; + } + XoMean = XoMean / Nobject; + sortfriends (XoList, IDList, Nobject); + float XoMin = XoList[0]; + + MidObj[0] = IDList[0]; + MidObj[1] = IDList[1]; + MidObj[2] = IDList[2]; + MidObj[3] = IDList[3]; + MidObj[4] = IDList[4]; + + for (i = 0; i < Nobject; i++) { + for (j = 0; j < NMID; j++) { + if (MidObj[j] == i) goto skip_object_p2; + } + float Dmean = 0.0; + for (j = 0; j < NMID; j++) { + Dmean += object[i].Dtgt[MidObj[j]]; // desired distance between object (i) and object in middle section + } + Dmean = Dmean / NMID; + if (object[i].Xo < XoMin) { + object[i].Xo = XoMin - Dmean; + } else { + object[i].Xo = XoMin + Dmean; + } + skip_object_p2: + continue; + } + } + + for (iter = 0; iter < Niter; iter ++) { + fprintf (stderr, "iter %d\n", iter); + + int nearNeighbors = (iter < nCloseIter); + + // save the result + if (0) { + char name[64]; + snprintf (name, 64, "output.%02d.dat", iter); + FILE *output = fopen (name, "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); + } + + // measure (dP/dX) & move object + for (i = 0; i < Nobject; i++) { + get_pressure_gradient (i, nCloseMax, farFrac, nearNeighbors, maxPressure); + move_object (i); + } + } + + // find the mean Xo, Yo positions for all objects + float XoSum = 0.0; + for (i = 0; i < Nobject; i++) { + float dY2sum = 0.0; + int Npts = 0; + for (j = 0; (j < nCloseMax) && (j < object[i].Nindex); j++) { + int idx = object[i].index[j]; + float dR2 = SQ(object[i].Dtgt[j]); + float dX2 = SQ(object[i].Xo - object[idx].Xo); + float dY2 = dR2 - dX2; + dY2sum += dY2; + Npts ++; + } + XoSum += object[i].Xo; + object[i].Yo = (dY2sum > 0.0) ? sqrt(fabs(dY2sum / Npts)) : -sqrt(fabs(dY2sum / Npts)); + } + float XoMean = XoSum / Nobject; + + // remove the mean Xo (uninteresting) + for (i = 0; i < Nobject; i++) { + object[i].Xo -= XoMean; + } + + // save the result + { + Vector *outindex = SelectVector ("sp1d_idx", ANYVECTOR, TRUE); if (!outindex) goto escape; + Vector *outXo = SelectVector ("sp1d_Xo", ANYVECTOR, TRUE); if (!outXo) goto escape; + Vector *outYo = SelectVector ("sp1d_Yo", ANYVECTOR, TRUE); if (!outYo) goto escape; + + ResetVector (outindex, OPIHI_INT, Nobject); + ResetVector (outXo, OPIHI_FLT, Nobject); + ResetVector (outYo, OPIHI_FLT, Nobject); + + for (i = 0; i < Nobject; i++) { + outindex->elements.Int[i] = i; + outXo->elements.Flt[i] = object[i].Xo; + outYo->elements.Flt[i] = object[i].Yo; + } + } + + return TRUE; + +escape: + gprint (GP_ERR, "invalid vector\n"); + return FALSE; + +usage: + gprint (GP_ERR, "USAGE: spex2dgas (index1) (index2) (distance) (Niter) (nCloseMax) (nCloseIter) (farFrac) (maxPressure) (idx1) (idx2)\n"); + return FALSE; +} + +/* this function takes 3 vectors: a set of distances and a pair of indicies identifying the end points. + + the goal is to generate groups of the end points that are relatively closer. + + note that the distances may not be Euclidean: nothing can be assumed about relationships between d(1,2), d(2,3), and d(1,3) + + the indicies should be sequential. the algorithm does not require it, but storage is much cleaner if it is + +*/ Index: /trunk/Ohana/src/opihi/cmd.astro/spex2dgas.c =================================================================== --- /trunk/Ohana/src/opihi/cmd.astro/spex2dgas.c (revision 33662) +++ /trunk/Ohana/src/opihi/cmd.astro/spex2dgas.c (revision 33662) @@ -0,0 +1,344 @@ +# include "astro.h" + +/* 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 a guess for positions of the objects at some locations in 2D + + P = (d_now - d_tru) + + Iterate over the points and move in X,Y + + Find (dP/dX,dP/dY) for modest moves in X,Y + + 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 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; + +static Object *object = NULL; +static int Nobject = 0; +// static int NOBJECT = 0; + +void sortfriends (float *X, int *IDX1, int N); + +static void get_pressure_gradient (int iObj, int nCloseMax, float farFrac, int nearNeighbors, float maxPressure) { + + int i; + + float dPdX = 0.0; + float dPdY = 0.0; + + // only use the first N friends + for (i = 0; (i < object[iObj].Nindex); i++) { + int jObj = object[iObj].index[i]; + + // some options: + // if the iterations are small, we should only worry about getting the near neighbors right + // if the iterations are large, we should add in more distant objects + if (nearNeighbors && (i >= nCloseMax)) break; + if (!nearNeighbors && (i >= nCloseMax)) { + if (drand48() > farFrac) 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 (maxPressure, MAX (-maxPressure, dF)); + + float dPdXi, dPdYi; + if (fabs(Dcur) < 1e-6) { + dPdXi = 0.0; + dPdYi = 0.0; + } else { + dPdXi = dF * dX / Dcur; + dPdYi = dF * dY / Dcur; + } + + if (isnan(Dtgt) || isnan(dX) || isnan(dY) || isnan(Dcur) || isnan(dF) || isnan(dPdXi) || isnan(dPdYi)) abort(); + + // if we are too close, then dX/Dcur is too ill-defined, just jump away by 5% of Dtgt + if (Dcur < 0.01*Dtgt) { + dPdXi = (drand48() - 0.5)*0.1*Dtgt; + dPdYi = (drand48() - 0.5)*0.1*Dtgt; + } + + if (i >= nCloseMax) { + 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; +} + +static void move_object (int iObj) { + + object[iObj].Xo += 0.25*object[iObj].dPdX; + object[iObj].Yo += 0.25*object[iObj].dPdY; + return; +} + +int spex2dgas (int argc, char **argv) { + + int i, j, iter, IoMax; + Vector *index1, *index2, *distance; + float XoMax, YoMax; + + { + // init random numbers + long A, B; + A = time(NULL); + for (B = 0; A == time(NULL); B++); + srand48(B); + } + + if (argc != 9) 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; + int Niter = atoi (argv[4]); + int nCloseMax = atoi (argv[5]); + int nCloseIter = atoi (argv[6]); + float farFrac = atof (argv[7]); + float maxPressure = atof (argv[8]); + // XXX enforce matching lengths on the three vectors + + CastVector (index1, OPIHI_INT); + CastVector (index2, OPIHI_INT); + + // how many objects do we have? + Nobject = 0; + for (i = 0; i < index1->Nelements; i++) { + Nobject = MAX (Nobject, index1->elements.Int[i]); + Nobject = MAX (Nobject, index2->elements.Int[i]); + } + Nobject ++; // after the loop, Nobject has the value of the highest index, not the count + + // allocate the list of Object -- these list the possible friends of the given object[i] + ALLOCATE (object, Object, Nobject); + for (i = 0; i < Nobject; i++) { + 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; + } + + // generate the set of vectors of all distances for each entry + // the object seq number (object[i]) will be used to match the index1 and index2 values + for (i = 0; i < index1->Nelements; i++) { + int N, m; + + N = index1->elements.Int[i]; + if (N >= Nobject) abort(); + 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].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++) { + 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 + // XXX for (i = 0; i < Nobject; i++) { + // XXX // object[i].Xo = dgrid * (int) (i % Ngrid); + // XXX // object[i].Yo = dgrid * (int) (i / Ngrid); + // XXX object[i].Xo = Dmax*drand48(); + // XXX object[i].Yo = Dmax*drand48(); + // XXX object[i].dPdX = 0.0; + // XXX object[i].dPdY = 0.0; + // XXX } + + // place the objects at the initial guess locations + // use object 0 and its most distant friend to constrain: + + int idx1 = 0; + int idx2 = object[idx1].index[object[idx1].Nindex-1]; + float A = object[idx1].Dtgt[object[idx1].Nindex-1]; + object[idx1].Xo = 0.0; + object[idx1].Yo = 0.0; + object[idx1].dPdX = 0.0; + object[idx1].dPdY = 0.0; + + object[idx2].Xo = A; + object[idx2].Yo = 0.0; + object[idx2].dPdX = 0.0; + object[idx2].dPdY = 0.0; + + // choose the correct Y side by comparing to the distance from the most deviant + YoMax = 0; + XoMax = 0; + IoMax = 0; + + for (i = 0; i < Nobject; i++) { + if (i == idx1) continue; + if (i == idx2) continue; + float B = NAN; + float C = NAN; + for (j = 0; (isnan(B) || isnan(C)) && (j < object[i].Nindex); j++) { + if (object[i].index[j] == idx1) { B = object[i].Dtgt[j]; } + if (object[i].index[j] == idx2) { C = object[i].Dtgt[j]; } + } + if (isnan(B) || isnan(C)) abort(); + + float Xo = (SQ(A) + SQ(B) - SQ(C)) / (2*A); + float Y2 = SQ(B) - SQ(Xo); + + float Yo = (Y2 < 0) ? 0.0 : sqrt(Y2); + if (isnan(Yo)) abort(); + + object[i].Xo = Xo; + object[i].Yo = Yo; + object[i].dPdX = 0.0; + object[i].dPdY = 0.0; + if (Yo > YoMax) { + YoMax = Yo; + XoMax = Xo; + IoMax = i; + } + } + + for (i = 0; i < Nobject; i++) { + if (i == idx1) continue; + if (i == idx2) continue; + if (i == IoMax) continue; + for (j = 0; j < object[i].Nindex; j++) { + if (object[i].index[j] == IoMax) { + float d02 = SQ(object[i].Dtgt[j]); + float dX = XoMax - object[i].Xo; + float dY1 = YoMax - object[i].Yo; + float dY2 = YoMax + object[i].Yo; + float d12 = SQ(dX) + SQ(dY1); + float d22 = SQ(dX) + SQ(dY2); + if (fabs(d12 - d02) > fabs(d22 - d02)) { + object[i].Yo = -object[i].Yo; + break; + } + } + } + } + + for (iter = 0; iter < Niter; iter ++) { + fprintf (stderr, "iter %d\n", iter); + + int nearNeighbors = (iter < nCloseIter); + + // save the result + char name[64]; + snprintf (name, 64, "output.%02d.dat", iter); + FILE *output = fopen (name, "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); + + // measure (dP/dX),(dP/dY) for all objects + for (i = 0; i < Nobject; i++) { + get_pressure_gradient (i, nCloseMax, farFrac, nearNeighbors, maxPressure); + 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: + gprint (GP_ERR, "invalid vector\n"); + return FALSE; + +usage: + gprint (GP_ERR, "USAGE: spex2dgas (index1) (index2) (distance) (Niter) (nCloseMax) (nCloseIter) (farFrac) (maxPressure)\n"); + return FALSE; +} + +/* this function takes 3 vectors: a set of distances and a pair of indicies identifying the end points. + + the goal is to generate groups of the end points that are relatively closer. + + note that the distances may not be Euclidean: nothing can be assumed about relationships between d(1,2), d(2,3), and d(1,3) + + the indicies should be sequential. the algorithm does not require it, but storage is much cleaner if it is + +*/ Index: /trunk/Ohana/src/opihi/cmd.astro/spexseq.c =================================================================== --- /trunk/Ohana/src/opihi/cmd.astro/spexseq.c (revision 33662) +++ /trunk/Ohana/src/opihi/cmd.astro/spexseq.c (revision 33662) @@ -0,0 +1,381 @@ +# include "astro.h" + +/* We have N objects with N(N-1)/2 paired chisq. We want to find the 1D sequence / index that + best describes the observed distances. + + A sequence is best described by the distances from objects which are only modestly close. + If we get too close to the reference, we (a) hit a non-euclidean (or N-D) minimum (A-B ~ A-C + ~ B-C) and (b) the positive and negative sides of the sequence are ambiguous. If we get too + far away, the distances are all either equivalent or less consistent + + 1) find the dynamic range of distances + + 2) choose object (1) and select a friend (2) > a fixed fraction of the full dynamic range (say, + 10-20%) + + 3) all nearby friends of (1) (D < f MaxDist) have sequences assigned based on the distance + to (2). center based on the distance (1)-(2) ie, (1) gets a sequence value of 0.0 + + 4) choose another object (1') and repeat [need to choose (2') on the same side as (2)?] + + ** reconciliation is a bit tricky + ** edge effects imply that we prefer to start in the center + ** is this well-defined? + + */ + +typedef struct { + int Nindex; // number of tested relationships + int NINDEX; // number of allocated relationships + int *index; // index for distances to other objects: object[i].Dtgt[j] is distance from object[i] to object[idx] where idx = object[i].index[j] (sorted by Dtgt) + int *rindex; // reverse index for Dtgt: distance from object[i] to object[idx] is object[i].Dtgt[j] where object[i].rindex[idx] = j + int shifted; // has this object already had the sequence adjusted? + float *Dtgt; // target distance for this relationship + float *seq; // sequence value for this friend + float So, dSo; + float Xo, dXo; + float Yo, dYo; + float Srange; + int nSeq; +} Object; + +void sortfriends (float *X, int *IDX1, int N); + +void get_sequence (Object *object, int Nobject, int idx1, float Dmax, float f1, float f2, int pin1, int pin2) { + + int i, j, N; + + // object (1) is provided above; choose object (2) + int idx2 = -1; + float A = NAN; + for (i = 0; (idx2 == -1) && (i < object[idx1].Nindex); i++) { + if (object[idx1].Dtgt[i] > f1*Dmax) { + idx2 = object[idx1].index[i]; + } + A = object[idx1].Dtgt[i]; + } + // fprintf (stderr, "idx1: %d, idx2: %d, A: %f\n", idx1, idx2, A); + // set_variable ("SP1D_A", A); + + // is idx2 closer to pin1 or pin2? + // 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; + for (i = 0; (object[idx1].Dtgt[i] < f2*Dmax) && (i < object[idx1].Nindex); i++) { + // what is the distance from this friend to object (2)? + N = object[idx1].index[i]; + j = object[idx2].rindex[N]; + assert (object[idx2].index[j] == N); + // 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; +} + +int spexseq (int argc, char **argv) { + + int i, j, k; + Vector *index1, *index2, *distance; + Object *object = NULL; + int Nobject = 0; + + { + // init random numbers + long A, B; + A = time(NULL); + for (B = 0; A == time(NULL); B++); + srand48(B); + } + + if (argc != 6) goto usage; + + // XXX enforce matching lengths on the three vectors + 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 f1 = atof (argv[4]); + float f2 = atof (argv[5]); + + // int idx1 = atoi(argv[9]); + // int idx2 = atoi(argv[10]); + + CastVector (index1, OPIHI_INT); + CastVector (index2, OPIHI_INT); + + // how many objects do we have? + Nobject = 0; + for (i = 0; i < index1->Nelements; i++) { + Nobject = MAX (Nobject, index1->elements.Int[i]); + Nobject = MAX (Nobject, index2->elements.Int[i]); + } + Nobject ++; // after the loop, Nobject has the value of the highest index, not the count + + // allocate the list of Object -- these list the possible friends of the given object[i] + ALLOCATE (object, Object, Nobject); + for (i = 0; i < Nobject; i++) { + 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)); + object[i].NINDEX = Nobject; + object[i].Nindex = 0; + } + + // generate the set of vectors of all distances for each entry + // the object seq number (object[i]) will be used to match the index1 and index2 values + for (i = 0; i < index1->Nelements; i++) { + int N, m; + + N = index1->elements.Int[i]; + if (N >= Nobject) abort(); + m = object[N].Nindex; + object[N].index[m] = index2->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 } + } + + for (i = 0; i < Nobject; i++) { + // sort so closest friends are first + sortfriends (object[i].Dtgt, object[i].index, object[i].Nindex); + object[i].shifted = FALSE; + object[i].So = NAN; + object[i].dSo = NAN; + object[i].Xo = NAN; + object[i].dXo = NAN; + // generate the reverse index + REALLOCATE (object[i].rindex, int, object[i].Nindex); + REALLOCATE (object[i].seq, float, object[i].Nindex); + for (j = 0; j < object[i].Nindex; j++) { + int N = object[i].index[j]; + object[i].rindex[N] = j; + object[i].seq[j] = NAN; + } + } + + // find and save the max distance + float Dmax = 0; + for (i = 0; i < distance->Nelements; i++) { + Dmax = MAX (Dmax, distance->elements.Flt[i]); + } + + // find two objects at the opposite extremes + float Dm1 = object[0].Dtgt[object[0].Nindex - 1]; + for (i = 0; (object[0].Dtgt[i] < 0.75*Dm1) && (i < object[0].Nindex); i++); + int pin1 = object[0].index[i]; + assert (pin1 >= 0); + assert (pin1 < Nobject); + + float Dm2 = object[pin1].Dtgt[object[pin1].Nindex - 1]; + for (i = 0; (object[pin1].Dtgt[i] < 0.75*Dm2) && (i < object[pin1].Nindex); i++); + int pin2 = object[pin1].index[i]; + assert (pin2 >= 0); + assert (pin2 < Nobject); + + // generate sequences for the near friends of all objects based on a somewhat distant friend + for (i = 0; i < Nobject; i++) { + get_sequence (object, Nobject, i, Dmax, f1, f2, pin1, pin2); + } + + // now we need to reconcile these sequences. + // start with object 0, find shift of all friends relative to 0 + object[0].shifted = TRUE; + for (i = 0; i < Nobject; i++) { + + // don't use unreconciled objects as a reference (we'll catch them eventually) + if (!object[i].shifted) continue; + + int Nfixed = 0; + + // shift all of the close friends of object (i) wrt object (i) + for (j = 0; (j < object[i].Nindex) && (object[i].Dtgt[j] < f2*Dmax); j++) { + + // loop over the sequence for object Nj + int Nj = object[i].index[j]; + + // only reconcile once + if (object[Nj].shifted == TRUE) continue; + + float dS1p = 0.0; + float dS2p = 0.0; + float dS1m = 0.0; + float dS2m = 0.0; + int nS = 0; + + float Smin = +1000.0; + float Smax = -1000.0; + + for (k = 0; k < object[Nj].Nindex; k++) { + int Nsj = object[Nj].index[k]; + int Nsi = object[i].rindex[Nsj]; + + if (isnan(object[Nj].seq[k])) continue; + 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]; + dS1p += dSp; + dS2p += dSp*dSp; + dS1m += dSm; + dS2m += dSm*dSm; + nS ++; + } + float Sp = dS1p / nS; + float Sm = dS1m / nS; + float dSp = sqrt(dS2p / nS - Sp*Sp); + float dSm = sqrt(dS2m / nS - Sm*Sm); + 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; + object[Nj].seq[k] -= Sp; + } + object[Nj].So = Sp; + object[Nj].dSo = dSp; + } else { + // flip the sequence and offset + for (k = 0; k < object[Nj].Nindex; k++) { + if (isnan(object[Nj].seq[k])) continue; + object[Nj].seq[k] = -object[Nj].seq[k] - Sm; + } + object[Nj].So = Sm; + 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); + } + + int Nfixed = 0; + int Nunfix = 0; + for (i = 0; i < Nobject; i++) { + if (object[i].shifted) { + Nfixed ++; + } else { + Nunfix ++; + } + } + fprintf (stderr, "%d fixed, %d unfixed\n", Nfixed, Nunfix); + + + // calculate the mean Xo for each object using the values of the sequences calculated + for (i = 0; i < Nobject; i++) { + + if (!object[i].shifted) { + fprintf (stderr, "this object was not reconciled\n"); + } + float XoSum = 0.0; + float XoS2 = 0.0; + int XoNum = 0; + + // find all the measurements of this object's sequence (from its friends values) + for (j = 0; j < Nobject; j++) { + int Ns = object[j].rindex[i]; + if (isnan(object[j].seq[Ns])) continue; + XoSum += object[j].seq[Ns]; + XoS2 += SQ(object[j].seq[Ns]); + XoNum ++; + } + float Xo = XoSum / XoNum; + float dXo = sqrt(XoS2 / XoNum - Xo*Xo); + object[i].Xo = Xo; + object[i].dXo = dXo; + } + + // XXX seq is effectively Xo, calculate Yo for the objects based on the local + // objects + + // save the result + { + Vector *outindex = SelectVector ("sp1d_idx", ANYVECTOR, TRUE); if (!outindex) goto escape; + Vector *outXo = SelectVector ("sp1d_Xo", ANYVECTOR, TRUE); if (!outXo) goto escape; + Vector *outYo = SelectVector ("sp1d_Yo", ANYVECTOR, TRUE); if (!outYo) goto escape; + + Vector *outSo = SelectVector ("sp1d_So", ANYVECTOR, TRUE); if (!outXo) goto escape; + Vector *outdSo = SelectVector ("sp1d_dSo", ANYVECTOR, TRUE); if (!outYo) goto escape; + 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); + ResetVector (outYo, OPIHI_FLT, Nobject); + + ResetVector (outSo, OPIHI_FLT, Nobject); + 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++) { + outindex->elements.Int[i] = i; + outXo->elements.Flt[i] = object[i].Xo; + outYo->elements.Flt[i] = object[i].Yo; + outSo->elements.Flt[i] = object[i].So; + 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; + } + } + + return TRUE; + +escape: + gprint (GP_ERR, "invalid vector\n"); + return FALSE; + +usage: + gprint (GP_ERR, "USAGE: spexseq (index1) (index2) (distance) (Niter) (nCloseMax) (nCloseIter) (farFrac) (maxPressure) (idx1) (idx2)\n"); + return FALSE; +} + +/* this function takes 3 vectors: a set of distances and a pair of indicies identifying the end points. + + the goal is to generate groups of the end points that are relatively closer. + + note that the distances may not be Euclidean: nothing can be assumed about relationships between d(1,2), d(2,3), and d(1,3) + + the indicies should be sequential. the algorithm does not require it, but storage is much cleaner if it is + +*/ Index: /trunk/Ohana/src/opihi/cmd.basic/list.c =================================================================== --- /trunk/Ohana/src/opihi/cmd.basic/list.c (revision 33661) +++ /trunk/Ohana/src/opihi/cmd.basic/list.c (revision 33662) @@ -19,4 +19,24 @@ remove_argument (N, &argc, argv); RunCommand = TRUE; + } + + if ((N = get_argument (argc, argv, "-vectors"))) { + remove_argument (N, &argc, argv); + if (argc != 2) { + gprint (GP_ERR, "USAGE: list (root) -vectors\n"); + return (FALSE); + } + ListVectorsToList (argv[1]); + return TRUE; + } + + if ((N = get_argument (argc, argv, "-buffers"))) { + remove_argument (N, &argc, argv); + if (argc != 2) { + gprint (GP_ERR, "USAGE: list (root) -buffers\n"); + return (FALSE); + } + ListBuffersToList (argv[1]); + return TRUE; } @@ -85,6 +105,14 @@ } + // old list must exist, or give an error sprintf (line, "%s:n", argv[3]); N = get_int_variable (line, &found); + if (!found) { + gprint (GP_ERR, "USAGE: list (newlist) -copy (oldlist) ...\n"); + gprint (GP_ERR, "ERROR: missing input list\n"); + return (FALSE); + } + + sprintf (line, "%s:n", argv[1]); set_int_variable (line, N); @@ -118,5 +146,4 @@ // remove the single named entry from the list (finds entry with given name, reduces list length by one) - // return an error if -add is given with no other args if ((argc > 2) && (!strcmp (argv[2], "-del"))) { if (argc != 4) { @@ -155,8 +182,13 @@ if (argc != 2) { - gprint (GP_ERR, "USAGE: list (root) : supply list data, terminate with 'END'\n"); - gprint (GP_ERR, "USAGE: list (root) -x (command) : create list from shell output\n"); - gprint (GP_ERR, "USAGE: list (root) -split (words) : create list from words\n"); - gprint (GP_ERR, "USAGE: list (root) -add (words) : extend a list\n"); + gprint (GP_ERR, "USAGE: list (root) : supply list data, terminate with 'END'\n"); + gprint (GP_ERR, "USAGE: list (root) -x (command) : create list from shell output\n"); + gprint (GP_ERR, "USAGE: list (root) -vectors : create list from vector names\n"); + gprint (GP_ERR, "USAGE: list (root) -buffers : create list from buffer names\n"); + gprint (GP_ERR, "USAGE: list (root) -split (words) : create list from words\n"); + gprint (GP_ERR, "USAGE: list (root) -splitbychar (char) (word) [(words)..] : create list from words\n"); + gprint (GP_ERR, "USAGE: list (root) -add (words) : extend a list\n"); + gprint (GP_ERR, "USAGE: list (root) -copy (list) : copy a list to a new name\n"); + gprint (GP_ERR, "USAGE: list (root) -del (word) : delete the entry by value\n"); return (FALSE); } Index: /trunk/Ohana/src/opihi/cmd.basic/shell.c =================================================================== --- /trunk/Ohana/src/opihi/cmd.basic/shell.c (revision 33661) +++ /trunk/Ohana/src/opihi/cmd.basic/shell.c (revision 33662) @@ -13,5 +13,5 @@ int exit_status; int wait_status; - int result, length; + int result; char **args, *shell; struct timeval start, now; @@ -33,16 +33,5 @@ args[0] = shell; args[1] = cmdflag; - - length = 0; - for (i = 1; i < argc; i++) { - length += strlen(argv[i]) + 1; - } - - ALLOCATE (args[2], char, length); - args[2][0] = 0; - for (i = 1; i < argc; i++) { - strcat (args[2], argv[i]); - if (i < argc - 1) strcat (args[2], " "); - } + args[2] = paste_args (argc - 1, &argv[1]); args[3] = NULL; Index: /trunk/Ohana/src/opihi/cmd.data/Makefile =================================================================== --- /trunk/Ohana/src/opihi/cmd.data/Makefile (revision 33661) +++ /trunk/Ohana/src/opihi/cmd.data/Makefile (revision 33662) @@ -123,4 +123,5 @@ $(SRC)/swapbytes.$(ARCH).o \ $(SRC)/textline.$(ARCH).o \ +$(SRC)/threshold.$(ARCH).o \ $(SRC)/tv.$(ARCH).o \ $(SRC)/tvchannel.$(ARCH).o \ Index: /trunk/Ohana/src/opihi/cmd.data/cumulative.c =================================================================== --- /trunk/Ohana/src/opihi/cmd.data/cumulative.c (revision 33661) +++ /trunk/Ohana/src/opihi/cmd.data/cumulative.c (revision 33662) @@ -23,4 +23,6 @@ opihi_flt *Vi = ivec[0].elements.Flt; *Vo = *Vi; + Vi++; + Vo++; for (i = 1; i < ivec[0].Nelements; i++, Vi++, Vo++) { *Vo = Vo[-1] + *Vi; @@ -29,4 +31,6 @@ opihi_int *Vi = ivec[0].elements.Int; *Vo = *Vi; + Vi++; + Vo++; for (i = 1; i < ivec[0].Nelements; i++, Vi++, Vo++) { *Vo = Vo[-1] + *Vi; Index: /trunk/Ohana/src/opihi/cmd.data/cursor.c =================================================================== --- /trunk/Ohana/src/opihi/cmd.data/cursor.c (revision 33661) +++ /trunk/Ohana/src/opihi/cmd.data/cursor.c (revision 33662) @@ -4,5 +4,5 @@ char string[20], key[20], *name; - int i, N, kapa; + int i, N, kapa, VERBOSE; double X, Y, R, D, Z; void *oldsignal; @@ -13,4 +13,9 @@ // if ((N = get_argument (argc, argv, "-g"))) { // if ((N = get_argument (argc, argv, "-i"))) { + + VERBOSE = TRUE; + if ((N = get_argument (argc, argv, "-a"))) { + VERBOSE = FALSE; + } name = NULL; @@ -53,5 +58,5 @@ set_str_variable ("KEY", key); - gprint (GP_LOG, "%s %f %f %f %f %f\n", key, X, Y, Z, R, D); + if (VERBOSE) gprint (GP_LOG, "%s %f %f %f %f %f\n", key, X, Y, Z, R, D); if (!strcasecmp (key, "Q")) break; Index: /trunk/Ohana/src/opihi/cmd.data/cut.c =================================================================== --- /trunk/Ohana/src/opihi/cmd.data/cut.c (revision 33661) +++ /trunk/Ohana/src/opihi/cmd.data/cut.c (revision 33662) @@ -48,6 +48,6 @@ case 'X': /* create output vectors */ - ResetVector (xvec, OPIHI_FLT, MAX (nx, 1)); - ResetVector (yvec, OPIHI_FLT, MAX (nx, 1)); + ResetVector (xvec, OPIHI_FLT, nx); + ResetVector (yvec, OPIHI_FLT, nx); bzero (yvec[0].elements.Flt, nx*sizeof(opihi_flt)); for (i = 0; i < nx; i++) { @@ -80,6 +80,6 @@ case 'y': case 'Y': - ResetVector (xvec, OPIHI_FLT, MAX (ny, 1)); - ResetVector (yvec, OPIHI_FLT, MAX (ny, 1)); + ResetVector (xvec, OPIHI_FLT, ny); + ResetVector (yvec, OPIHI_FLT, ny); bzero (yvec[0].elements.Flt, ny*sizeof(opihi_flt)); for (i = 0; i < ny; i++) { Index: /trunk/Ohana/src/opihi/cmd.data/imcut.c =================================================================== --- /trunk/Ohana/src/opihi/cmd.data/imcut.c (revision 33661) +++ /trunk/Ohana/src/opihi/cmd.data/imcut.c (revision 33662) @@ -37,6 +37,6 @@ dY = dY / L; - ResetVector (xvec, OPIHI_FLT, MAX (L, 1)); - ResetVector (yvec, OPIHI_FLT, MAX (L, 1)); + ResetVector (xvec, OPIHI_FLT, L); + ResetVector (yvec, OPIHI_FLT, L); V = (float *)buf[0].matrix.buffer; Index: /trunk/Ohana/src/opihi/cmd.data/init.c =================================================================== --- /trunk/Ohana/src/opihi/cmd.data/init.c (revision 33661) +++ /trunk/Ohana/src/opihi/cmd.data/init.c (revision 33662) @@ -61,4 +61,5 @@ int header PROTO((int, char **)); int list_vectors PROTO((int, char **)); +int vtype PROTO((int, char **)); int load PROTO((int, char **)); int lookup PROTO((int, char **)); @@ -111,4 +112,5 @@ int swapbytes PROTO((int, char **)); int textline PROTO((int, char **)); +int threshold PROTO((int, char **)); int tv PROTO((int, char **)); int tvchannel PROTO((int, char **)); @@ -138,4 +140,8 @@ int zplot PROTO((int, char **)); int zcplot PROTO((int, char **)); + +// ??? +// int mtype PROTO((int, char **)); +// {1, "mtype", mtype, "return the type of the defined buffer"}, static Command cmds[] = { @@ -257,4 +263,5 @@ {1, "swapbytes", swapbytes, "byte swap thing"}, {1, "textline", textline, "write text line on graph"}, + {1, "threshold", threshold, "find (interpolate) location of transition"}, {1, "tv", tv, "display an image on the Kii window"}, {1, "tvchannel", tvchannel, "set the current tv channel"}, @@ -269,4 +276,5 @@ {1, "vclip", vclip, "clip values in a vector to be within a range"}, {1, "vectors", list_vectors, "list vectors"}, + {1, "vtype", vtype, "return the vector type (FLT or INT)"}, {1, "vgauss", vgauss, "fit a Gaussian to a vector"}, {1, "vellipse", vellipse, "fit a Ellipse to a vector pair"}, Index: /trunk/Ohana/src/opihi/cmd.data/line.c =================================================================== --- /trunk/Ohana/src/opihi/cmd.data/line.c (revision 33661) +++ /trunk/Ohana/src/opihi/cmd.data/line.c (revision 33662) @@ -3,7 +3,14 @@ int line (int argc, char **argv) { - int kapa; + int kapa, N; Graphdata graphmode; float x[2], y[2]; + + /* FracPositions uses coordinates of 0-1 relative to axis range */ + int FracPositions = FALSE; + if ((N = get_argument (argc, argv, "-frac"))) { + remove_argument (N, &argc, argv); + FracPositions = TRUE; + } if (!style_args (&graphmode, &argc, argv, &kapa)) return FALSE; @@ -18,4 +25,11 @@ y[1] = atof(argv[5]); + if (FracPositions) { + x[0] = x[0] * (graphmode.xmax - graphmode.xmin) + graphmode.xmin; + y[0] = y[0] * (graphmode.ymax - graphmode.ymin) + graphmode.ymin; + x[1] = x[1] * (graphmode.xmax - graphmode.xmin) + graphmode.xmin; + y[1] = y[1] * (graphmode.ymax - graphmode.ymin) + graphmode.ymin; + } + /* set point style and errorbar mode (these are NOT sticky) */ graphmode.style = 0; Index: /trunk/Ohana/src/opihi/cmd.data/list_buffers.c =================================================================== --- /trunk/Ohana/src/opihi/cmd.data/list_buffers.c (revision 33661) +++ /trunk/Ohana/src/opihi/cmd.data/list_buffers.c (revision 33662) @@ -9,2 +9,43 @@ } + +# if (0) +int mtype (int argc, char **argv) { + + Variable = (char *) NULL; + if ((N = get_argument (argc, argv, "-var"))) { + remove_argument (N, &argc, argv); + Variable = strcreate (argv[N]); + remove_argument (N, &argc, argv); + } + + if (argc != 2) { + gprint (GP_ERR, "USAGE: mtype (buffer) [-var out]\n"); + return (FALSE); + } + + Buffer *buf = SelectBuffer (argv[1], OLDBUFFER, FALSE); + if (!buf) { + gprint (GP_ERR, "unknown buffer %s\n", argv[1]); + free (Variable); + return FALSE; + } + + if (buf->type == OPIHI_FLT) { + if (Variable) { + set_str_variable (Variable, "FLT"); + } else { + gprint (GP_LOG, "%s : FLT\n", argv[1]); + } + } else { + if (Variable) { + set_str_variable (Variable, "INT"); + } else { + gprint (GP_LOG, "%s : INT\n", argv[1]); + } + } + if (Variable) free (Variable); + + return (TRUE); +} +# endif Index: /trunk/Ohana/src/opihi/cmd.data/list_vectors.c =================================================================== --- /trunk/Ohana/src/opihi/cmd.data/list_vectors.c (revision 33661) +++ /trunk/Ohana/src/opihi/cmd.data/list_vectors.c (revision 33662) @@ -7,2 +7,43 @@ } +int vtype (int argc, char **argv) { + + int N; + + char *Variable = (char *) NULL; + if ((N = get_argument (argc, argv, "-var"))) { + remove_argument (N, &argc, argv); + Variable = strcreate (argv[N]); + remove_argument (N, &argc, argv); + } + + if (argc != 2) { + gprint (GP_ERR, "USAGE: vtype (vector) [-var out]\n"); + return (FALSE); + } + + Vector *vec = SelectVector (argv[1], OLDVECTOR, FALSE); + if (!vec) { + gprint (GP_ERR, "unknown vector %s\n", argv[1]); + free (Variable); + return FALSE; + } + + if (vec->type == OPIHI_FLT) { + if (Variable) { + set_str_variable (Variable, "FLT"); + } else { + gprint (GP_LOG, "%s : FLT\n", argv[1]); + } + } else { + if (Variable) { + set_str_variable (Variable, "INT"); + } else { + gprint (GP_LOG, "%s : INT\n", argv[1]); + } + } + if (Variable) free (Variable); + + return (TRUE); +} + Index: /trunk/Ohana/src/opihi/cmd.data/read_vectors.c =================================================================== --- /trunk/Ohana/src/opihi/cmd.data/read_vectors.c (revision 33661) +++ /trunk/Ohana/src/opihi/cmd.data/read_vectors.c (revision 33662) @@ -21,19 +21,23 @@ } +// vector types +enum {COLTYPE_NONE, COLTYPE_FLT, COLTYPE_INT, COLTYPE_TIME}; + int read_vectors (int argc, char **argv) { - int i, j, Nskip, Nvec, *col, done, status, IsCSV, VERBOSE; - int Nbytes, nbytes, Nstart, NELEM, N, nread; + int TimeFormat; + time_t TimeReference; + int i, j, Nskip, Narg, Nvec, *col, IsCSV, VERBOSE; + int Nbytes, Nstart, NELEM, Nelem, nread, *coltype; char *colstr, *c0, *c1, *buffer, *extname; - double value; Vector **vec; /* auto-sense table type */ - if ((N = get_argument (argc, argv, "-fits"))) { - remove_argument (N, &argc, argv); - extname = strcreate (argv[N]); + if ((Narg = get_argument (argc, argv, "-fits"))) { + remove_argument (Narg, &argc, argv); + extname = strcreate (argv[Narg]); if (extname == (char *) NULL) return (FALSE); - remove_argument (N, &argc, argv); - status = read_table_vectors (argc, argv, extname); + remove_argument (Narg, &argc, argv); + int status = read_table_vectors (argc, argv, extname); free (extname); return (status); @@ -41,19 +45,19 @@ Nskip = 0; - if ((N = get_argument (argc, argv, "-skip"))) { - remove_argument (N, &argc, argv); - Nskip = atof (argv[N]); - remove_argument (N, &argc, argv); + if ((Narg = get_argument (argc, argv, "-skip"))) { + remove_argument (Narg, &argc, argv); + Nskip = atof (argv[Narg]); + remove_argument (Narg, &argc, argv); } IsCSV = FALSE; - if ((N = get_argument (argc, argv, "-csv"))) { - remove_argument (N, &argc, argv); + if ((Narg = get_argument (argc, argv, "-csv"))) { + remove_argument (Narg, &argc, argv); IsCSV = TRUE; } VERBOSE = FALSE; - if ((N = get_argument (argc, argv, "-v"))) { - remove_argument (N, &argc, argv); + if ((Narg = get_argument (argc, argv, "-v"))) { + remove_argument (Narg, &argc, argv); VERBOSE = TRUE; } @@ -65,4 +69,7 @@ /* read name N name N */ + // do this only optionally? + GetTimeFormat (&TimeReference, &TimeFormat); + if (f == (FILE *) NULL) { gprint (GP_ERR, "no open file for read\n"); @@ -71,14 +78,31 @@ fseeko (f, 0LL, SEEK_SET); -// if (IsCSV) { -// status = read_vectors_csv (argc, argv, Nskip, f); -// return status; -// } - Nvec = (argc - 1) / 2; ALLOCATE (vec, Vector *, Nvec); ALLOCATE (col, int, Nvec); + ALLOCATE (coltype, int, Nvec); for (i = 0; i < Nvec; i++) { + + // interpret the column names including type flags + // XXX review the grammar before releasing this : is foo:type best, or is something else needed? + // Note the conflict wrt list entries + // the name may be of the form foo:type, where type may be one of : int, float, time + + coltype[i] = COLTYPE_FLT; + char *colname = argv[2*i + 1]; + char *ptr = strchr (colname, ':'); + if (ptr) { + // split out colname and type + *ptr = 0; + ptr ++; + if (!ptr) goto bad_colname; + coltype[i] = COLTYPE_NONE; + if (!strcasecmp(ptr, "float")) { coltype[i] = COLTYPE_FLT; } + if (!strcasecmp(ptr, "int")) { coltype[i] = COLTYPE_INT; } + if (!strcasecmp(ptr, "time")) { coltype[i] = COLTYPE_TIME; } + if (!coltype[i]) goto bad_colname; + } + if ((vec[i] = SelectVector (argv[2*i + 1], ANYVECTOR, TRUE)) == NULL) { gprint (GP_ERR, "USAGE: read name N name N ...\n"); @@ -87,5 +111,5 @@ return (FALSE); } - // XXX we could allow flags (eg, N.i) to specify INT vs FLT types vectors... + colstr = argv[2*i+2]; for (j = 0; j < strlen (colstr); j++) { @@ -119,5 +143,10 @@ NELEM = 1000; for (i = 0; i < Nvec; i++) { - ResetVector (vec[i], OPIHI_FLT, NELEM); + if (coltype[i] == COLTYPE_INT) { + ResetVector (vec[i], OPIHI_INT, NELEM); + } else { + // note that COLTYPE_TIME is a type of float + ResetVector (vec[i], OPIHI_FLT, NELEM); + } } @@ -128,8 +157,11 @@ } - Nstart = 0; - N = 0; - done = FALSE; - while (!done) { + // we have a working buffer read from the file. we parse the lines in the working buffer + // until we reach the last chunk without an EOL char. at that point, we shift the start + // of the last (partial) line to the start of the buffer and re-fill. + + Nstart = 0; // location of the last valid byte in the buffer (start filling here) + Nelem = 0; // number of valid rows read (vector elements) + while (TRUE) { Nbytes = 0x10000 - Nstart; bzero (&buffer[Nstart], Nbytes); @@ -139,61 +171,83 @@ break; } - if (nread == 0) break; - nbytes = nread + Nstart; + if (nread == 0) break; // end of the file + // nbytes = nread + Nstart; - status = TRUE; - c0 = buffer; - while (status) { - c1 = strchr (c0, '\n'); + int bufferStatus = TRUE; + c0 = buffer; // c0 always marks the start of a line + while (bufferStatus) { + c1 = strchr (c0, '\n'); // find the end of this current line if (c1 == (char *) NULL) { Nstart = strlen (c0); memmove (buffer, c0, Nstart); - status = FALSE; - } else { - *c1 = 0; - } - if ((*c0 != '#') && (*c0 != '!')) { - for (i = 0; (i < Nvec) && status; i++) { + bufferStatus = FALSE; + continue; + } + *c1 = 0; // mark the end of the line + + if (*c0 == '#') { c0 = c1 + 1; continue; } + if (*c0 == '!') { c0 = c1 + 1; continue; } + + // parse the vectors in this line. this code is a bit inefficient: each column + // requires a separate pass through the line. + + int lineStatus = TRUE; + for (i = 0; i < Nvec; i++) { + int ivalue; + double dvalue; + time_t tvalue; + int readStatus = FALSE; + // need to make the if cases for coltype[i] + switch (coltype[i]) { + case COLTYPE_INT: + readStatus = IsCSV ? iparse_csv (&ivalue, col[i], c0) : iparse (&ivalue, col[i], c0); + vec[i][0].elements.Int[Nelem] = readStatus ? ivalue : 0; + break; + case COLTYPE_FLT: + readStatus = IsCSV ? dparse_csv (&dvalue, col[i], c0) : dparse (&dvalue, col[i], c0); + vec[i][0].elements.Flt[Nelem] = readStatus ? dvalue : NAN; + break; + case COLTYPE_TIME: + readStatus = IsCSV ? tparse_csv (&tvalue, col[i], c0) : tparse (&tvalue, col[i], c0); + dvalue = TimeValue (tvalue, TimeReference, TimeFormat); + vec[i][0].elements.Flt[Nelem] = readStatus ? dvalue : NAN; + break; + } + if (!readStatus && VERBOSE) { if (IsCSV) { - status = dparse_csv (&value, col[i], c0); + gprint (GP_ERR, "suspect field: %d (%s) in %s\n", col[i], argv[2*i+2], c0); } else { - status = dparse (&value, col[i], c0); + gprint (GP_ERR, "suspect field: %d in %s\n", col[i], c0); } - if (status) { - vec[i][0].elements.Flt[N] = value; - } else { - vec[i][0].elements.Flt[N] = NAN; - if (VERBOSE) { - if (IsCSV) { - gprint (GP_ERR, "suspect field: %d (%s) in %s\n", col[i], argv[2*i+2], c0); - } else { - gprint (GP_ERR, "suspect field: %d in %s\n", col[i], c0); - } - } - } - } - if (status) { - N++; - } else { - if (VERBOSE && FALSE) { - char temp[32]; - strncpy (temp, c0, 32); - temp[31] = 0; - gprint (GP_ERR, "skip line %s\n\n", temp); - } - } - } - c0 = c1 + 1; - if (N == NELEM) { + } + lineStatus &= readStatus; + } + if (!lineStatus && VERBOSE) { + char temp[32]; + strncpy (temp, c0, 32); + temp[31] = 0; + gprint (GP_ERR, "skip line %s\n\n", temp); + } + Nelem ++; + if (Nelem == NELEM) { NELEM += 1000; for (i = 0; i < Nvec; i++) { - REALLOCATE (vec[i][0].elements.Flt, opihi_flt, NELEM); - } - } + if (coltype[i] == COLTYPE_INT) { + REALLOCATE (vec[i][0].elements.Int, opihi_int, NELEM); + } else { + REALLOCATE (vec[i][0].elements.Flt, opihi_flt, NELEM); + } + } + } + c0 = c1 + 1; } } for (i = 0; i < Nvec; i++) { - REALLOCATE (vec[i][0].elements.Flt, opihi_flt, MAX (N,1)); - vec[i][0].Nelements = N; + if (coltype[i] == COLTYPE_INT) { + REALLOCATE (vec[i][0].elements.Int, opihi_int, MAX (Nelem,1)); + } else { + REALLOCATE (vec[i][0].elements.Flt, opihi_flt, MAX (Nelem,1)); + } + vec[i][0].Nelements = Nelem; } @@ -242,4 +296,12 @@ Nextend = atoi (extname); } + + // XXX ReadAll needs: deal with Extnum vs Extname, save vectors, etc + // ReadAll = FALSE; + // if ((N = get_argument (argc, argv, "-all"))) { + // remove_argument (N, &argc, argv); + // ReadAll = atoi (extname); + // if (argc != 1) ESCAPE ("-all option cannot be mixed with selected field"); + // } if (argc < 2) ESCAPE ("USAGE: read -fits extension [-extnum] [-keyword key] name name ..."); @@ -314,9 +376,4 @@ for (i = 1; i < argc; i++) { void *data; - char *Pc; - short *Ps; - int *Pi; - float *Pf; - double *Pd; int Nval; char name[80]; @@ -345,44 +402,52 @@ sprintf (name, "%s:%d", argv[i], j); if ((vec[j] = SelectVector (name, ANYVECTOR, TRUE)) == NULL) ESCAPE ("bad vector name"); - ResetVector (vec[j], vecType, MAX (Ny,1)); - } - + ResetVector (vec[j], vecType, Ny); + } + + if (!strcmp (type, "char")) { + char *Ptr = data; + for (j = 0; j < Ny; j++) { + for (k = 0; k < Nval; k++, Ptr++) { + vec[k][0].elements.Int[j] = *Ptr; + } + } + } + if (!strcmp (type, "short")) { + short *Ptr = data; + for (j = 0; j < Ny; j++) { + for (k = 0; k < Nval; k++, Ptr++) { + vec[k][0].elements.Int[j] = *Ptr; + } + } + } + if (!strcmp (type, "int")) { + int *Ptr = data; + for (j = 0; j < Ny; j++) { + for (k = 0; k < Nval; k++, Ptr++) { + vec[k][0].elements.Int[j] = *Ptr; + } + } + } + if (!strcmp (type, "int64_t")) { + int64_t *Ptr = data; + for (j = 0; j < Ny; j++) { + for (k = 0; k < Nval; k++, Ptr++) { + vec[k][0].elements.Int[j] = *Ptr; + } + } + } + if (!strcmp (type, "float")) { + float *Ptr = data; + for (j = 0; j < Ny; j++) { + for (k = 0; k < Nval; k++, Ptr++) { + vec[k][0].elements.Flt[j] = *Ptr; + } + } + } if (!strcmp (type, "double")) { - Pd = (double *) data; - for (j = 0; j < Ny; j++) { - for (k = 0; k < Nval; k++, Pd++) { - vec[k][0].elements.Flt[j] = *Pd; - } - } - } - if (!strcmp (type, "float")) { - Pf = (float *) data; - for (j = 0; j < Ny; j++) { - for (k = 0; k < Nval; k++, Pf++) { - vec[k][0].elements.Flt[j] = *Pf; - } - } - } - if (!strcmp (type, "int")) { - Pi = (int *) data; - for (j = 0; j < Ny; j++) { - for (k = 0; k < Nval; k++, Pi++) { - vec[k][0].elements.Int[j] = *Pi; - } - } - } - if (!strcmp (type, "short")) { - Ps = (short *) data; - for (j = 0; j < Ny; j++) { - for (k = 0; k < Nval; k++, Ps++) { - vec[k][0].elements.Int[j] = *Ps; - } - } - } - if (!strcmp (type, "char")) { - Pc = (char *) data; - for (j = 0; j < Ny; j++) { - for (k = 0; k < Nval; k++, Pc++) { - vec[k][0].elements.Int[j] = *Pc; + double *Ptr = data; + for (j = 0; j < Ny; j++) { + for (k = 0; k < Nval; k++, Ptr++) { + vec[k][0].elements.Flt[j] = *Ptr; } } @@ -396,104 +461,2 @@ return (TRUE); } - -# if (0) -int read_vectors_csv (int argc, char **argv, int Nskip, FILE *f) { - - int i, j, Nvec, *col, done, status; - int Nbytes, nbytes, Nstart, NELEM, N, nread; - char *colstr, *c0, *c1, *buffer, *extname; - double value; - Vector **vec; - - Nvec = (argc - 1) / 2; - ALLOCATE (vec, Vector *, Nvec); - ALLOCATE (col, int, Nvec); - - for (i = 0; i < Nvec; i++) { - if ((vec[i] = SelectVector (argv[2*i + 1], ANYVECTOR, TRUE)) == NULL) { - gprint (GP_ERR, "USAGE: read name N name N ...\n"); - free (vec); - free (col); - return (FALSE); - } - // XXX we could allow flags (eg, N.i) to specify INT vs FLT types vectors... - colstr = argv[2*i+2]; - for (j = 0; j < strlen (colstr); j++) { - if (!isdigit(colstr[j])) { - gprint (GP_ERR, "USAGE: read name N name N ...\n"); - free (vec); - free (col); - return (FALSE); - } - } - col[i] = atof (colstr); - } - - // currently, all read vectors are forced to be type FLT - NELEM = 1000; - for (i = 0; i < Nvec; i++) { - ResetVector (vec[i], OPIHI_FLT, NELEM); - } - - ALLOCATE (buffer, char, 0x10001); - bzero (buffer, 0x10001); - for (i = 0; i < Nskip; i++) { - scan_line (f, buffer); - } - - Nstart = 0; - N = 0; - done = FALSE; - while (!done) { - Nbytes = 0x10000 - Nstart; - bzero (&buffer[Nstart], Nbytes); - nread = fread (&buffer[Nstart], 1, Nbytes, f); - if (ferror (f)) { - perror ("error reading data file"); - break; - } - if (nread == 0) break; - nbytes = nread + Nstart; - - status = TRUE; - c0 = buffer; - while (status) { - c1 = strchr (c0, '\n'); - if (c1 == (char *) NULL) { - Nstart = strlen (c0); - memmove (buffer, c0, Nstart); - status = FALSE; - } else { - *c1 = 0; - } - if ((*c0 != '#') && (*c0 != '!')) { - for (i = 0; (i < Nvec) && status; i++) { - status = dparse_csv (&value, col[i], c0); - vec[i][0].elements.Flt[N] = value; - if (!status) vec[i][0].elements.Flt[N] = NAN; - } - if (status) N++; - } - c0 = c1 + 1; - if (N == NELEM) { - NELEM += 1000; - for (i = 0; i < Nvec; i++) { - REALLOCATE (vec[i][0].elements.Flt, opihi_flt, NELEM); - } - } - - } - } - for (i = 0; i < Nvec; i++) { - REALLOCATE (vec[i][0].elements.Flt, opihi_flt, MAX (N,1)); - vec[i][0].Nelements = N; - } - - free (vec); - free (col); - free (buffer); - return (TRUE); - -} - -# endif Index: /trunk/Ohana/src/opihi/cmd.data/reindex.c =================================================================== --- /trunk/Ohana/src/opihi/cmd.data/reindex.c (revision 33661) +++ /trunk/Ohana/src/opihi/cmd.data/reindex.c (revision 33662) @@ -22,5 +22,5 @@ // ovec matches ivec in type and xvec in size (xvec need not have all ivec elements, and may have duplicates - ResetVector (ovec, ivec->type, MAX (xvec[0].Nelements, 1)); + ResetVector (ovec, ivec->type, xvec[0].Nelements); Nmax = ivec[0].Nelements - 1; @@ -51,5 +51,5 @@ // free up unused memory - ResetVector (ovec, ivec->type, MAX (Npts, 1)); + ResetVector (ovec, ivec->type, Npts); return (TRUE); Index: /trunk/Ohana/src/opihi/cmd.data/select.c =================================================================== --- /trunk/Ohana/src/opihi/cmd.data/select.c (revision 33661) +++ /trunk/Ohana/src/opihi/cmd.data/select.c (revision 33662) @@ -28,7 +28,7 @@ if ((in1->type == OPIHI_INT) && (in2->type == OPIHI_INT)) { - ResetVector (ovec, OPIHI_INT, MAX (tvec[0].Nelements, 1)); + ResetVector (ovec, OPIHI_INT, tvec[0].Nelements); } else { - ResetVector (ovec, OPIHI_FLT, MAX (tvec[0].Nelements, 1)); + ResetVector (ovec, OPIHI_FLT, tvec[0].Nelements); } Index: /trunk/Ohana/src/opihi/cmd.data/test/read.sh =================================================================== --- /trunk/Ohana/src/opihi/cmd.data/test/read.sh (revision 33662) +++ /trunk/Ohana/src/opihi/cmd.data/test/read.sh (revision 33662) @@ -0,0 +1,145 @@ + +list tests + #test1 + #test2 + #test3 + #test4 + test5 + test6 +end + +# basic read +macro test1 + + $PASS = 1 + + data read.t1.dat + read t1 1 t2 2 t3 3 + + if (t1[0] != 1) set PASS = 0 + if (t1[1] != 2) set PASS = 0 + if (t1[2] != 4) set PASS = 0 + + if (t2[0] != 5) set PASS = 0 + if (t2[1] != 3) set PASS = 0 + if (t2[2] != 5) set PASS = 0 + + if (t3[0] != 8) set PASS = 0 + if (t3[1] != 9) set PASS = 0 + if (t3[2] != 7) set PASS = 0 +end + +# basic read with int types +macro test2 + + $PASS = 1 + + data read.t1.dat + read t1:int 1 t2:int 2 t3:int 3 + + # note: no way to verify the mode of a vector + vectors + + if (t1[0] != 1) set PASS = 0 + if (t1[1] != 2) set PASS = 0 + if (t1[2] != 4) set PASS = 0 + + if (t2[0] != 5) set PASS = 0 + if (t2[1] != 3) set PASS = 0 + if (t2[2] != 5) set PASS = 0 + + if (t3[0] != 8) set PASS = 0 + if (t3[1] != 9) set PASS = 0 + if (t3[2] != 7) set PASS = 0 +end + +# basic read with comments +macro test3 + + $PASS = 1 + + data read.t2.dat + read t1 1 t2 2 t3 3 + + if (t1[0] != 1) set PASS = 0 + if (t1[1] != 2) set PASS = 0 + if (t1[2] != 4) set PASS = 0 + + if (t2[0] != 5) set PASS = 0 + if (t2[1] != 3) set PASS = 0 + if (t2[2] != 5) set PASS = 0 + + if (t3[0] != 8) set PASS = 0 + if (t3[1] != 9) set PASS = 0 + if (t3[2] != 7) set PASS = 0 +end + +# basic read with comments & int types +macro test4 + + $PASS = 1 + + data read.t2.dat + read t1:int 1 t2:int 2 t3:int 3 + + # note: no way to verify the mode of a vector + vectors + + if (t1[0] != 1) set PASS = 0 + if (t1[1] != 2) set PASS = 0 + if (t1[2] != 4) set PASS = 0 + + if (t2[0] != 5) set PASS = 0 + if (t2[1] != 3) set PASS = 0 + if (t2[2] != 5) set PASS = 0 + + if (t3[0] != 8) set PASS = 0 + if (t3[1] != 9) set PASS = 0 + if (t3[2] != 7) set PASS = 0 +end + +# read with invalid entries +macro test5 + + $PASS = 1 + + data read.t3.dat + read t1 1 t2 2 t3 3 + + if (t1[0] != 1) set PASS = 0 + if (t1[1] != 2) set PASS = 0 + if (t1[2] != 4) set PASS = 0 + + if (t2[0] != 5) set PASS = 0 + if (not(isnan(t2[1]))) set PASS = 0 + if (t2[2] != 5) set PASS = 0 + + if (t3[0] != 8) set PASS = 0 + if (t3[1] != 9) set PASS = 0 + if (not(isnan(t3[2]))) set PASS = 0 +end + +# read with invalid entries & int types +macro test6 + + $PASS = 1 + + data read.t3.dat + read t1:int 1 t2:int 2 t3:int 3 + + # note: no way to verify the mode of a vector + vectors + + if (t1[0] != 1) set PASS = 0 + if (t1[1] != 2) set PASS = 0 + if (t1[2] != 4) set PASS = 0 + + if (t2[0] != 5) set PASS = 0 + if (t2[1] != 0) set PASS = 0 + if (t2[2] != 5) set PASS = 0 + + if (t3[0] != 8) set PASS = 0 + if (t3[1] != 9) set PASS = 0 + if (t3[2] != 0) set PASS = 0 +end + Index: /trunk/Ohana/src/opihi/cmd.data/test/read.t1.dat =================================================================== --- /trunk/Ohana/src/opihi/cmd.data/test/read.t1.dat (revision 33662) +++ /trunk/Ohana/src/opihi/cmd.data/test/read.t1.dat (revision 33662) @@ -0,0 +1,3 @@ +1 5 8 +2 3 9 +4 5 7 Index: /trunk/Ohana/src/opihi/cmd.data/test/read.t2.dat =================================================================== --- /trunk/Ohana/src/opihi/cmd.data/test/read.t2.dat (revision 33662) +++ /trunk/Ohana/src/opihi/cmd.data/test/read.t2.dat (revision 33662) @@ -0,0 +1,5 @@ +# comment line +1 5 8 +2 3 9 +! comment line 2 +4 5 7 Index: /trunk/Ohana/src/opihi/cmd.data/test/read.t3.dat =================================================================== --- /trunk/Ohana/src/opihi/cmd.data/test/read.t3.dat (revision 33662) +++ /trunk/Ohana/src/opihi/cmd.data/test/read.t3.dat (revision 33662) @@ -0,0 +1,3 @@ +1 5 8 +2 X 9 +4 5 Y Index: /trunk/Ohana/src/opihi/cmd.data/test/read.t4.dat =================================================================== --- /trunk/Ohana/src/opihi/cmd.data/test/read.t4.dat (revision 33662) +++ /trunk/Ohana/src/opihi/cmd.data/test/read.t4.dat (revision 33662) @@ -0,0 +1,7 @@ +1 5 8 2000/01/01,00:00:00 +2 3 9 2000/01/01,01:00:00 +4 5 7 2000/01/01,02:00:00 +4 5 7 now +4 5 7 today +4 5 7 2451544.750000j +4 5 7 51544.500000J Index: /trunk/Ohana/src/opihi/cmd.data/threshold.c =================================================================== --- /trunk/Ohana/src/opihi/cmd.data/threshold.c (revision 33662) +++ /trunk/Ohana/src/opihi/cmd.data/threshold.c (revision 33662) @@ -0,0 +1,88 @@ +# include "data.h" + +int threshold (int argc, char **argv) { + + int N, QUIET; + double value; + Vector *vecx, *vecy; + + QUIET = FALSE; + if ((N = get_argument (argc, argv, "-q"))) { + QUIET = TRUE; + remove_argument (N, &argc, argv); + } + + if (argc != 4) { + gprint (GP_ERR, "USAGE: threshold (value)\n"); + gprint (GP_ERR, " find the x coordinate at which we pass the specified value\n"); + gprint (GP_ERR, " by default, y must be monotonically increasing\n"); + return (FALSE); + } + + if ((vecx = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE); + if ((vecy = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) return (FALSE); + + value = atof (argv[3]); + + // if (argc == 6) { + // start = atof (argv[4]); + // end = atof (argv[5]); + // } else { + // start = (vecx[0].type == OPIHI_FLT) ? vecx[0].elements.Flt[0] : vecx[0].elements.Int[0]; + // end = (vecx[0].type == OPIHI_FLT) ? vecx[0].elements.Flt[vecx[0].Nelements - 1] : vecx[0].elements.Int[vecx[0].Nelements - 1]; + // } + + int isFltX = (vecx[0].type == OPIHI_FLT); + int isFltY = (vecy[0].type == OPIHI_FLT); + + // use bisection to find the value + int Nlo, Nhi; + int Nelements = vecx[0].Nelements; + + // find the last entry before start + Nlo = 0; + Nhi = Nelements - 1; + while (Nhi - Nlo > 10) { + N = 0.5*(Nlo + Nhi); + double testval = isFltY ? vecy[0].elements.Flt[N] : vecy[0].elements.Int[N]; + if (testval < value) { + Nlo = MAX(N, 0); + } else { + Nhi = MIN(N, Nelements - 1); + } + } + // v[Nlo] < value <= v[Nhi] + for (N = Nlo; N <= Nhi; N++) { + double testval = isFltY ? vecy[0].elements.Flt[N] : vecy[0].elements.Int[N]; + if (testval > value) { + Nhi = N; + break; + } + } + // v[Nhi] is transition bin + + double x0, x1, y0, y1, Xvalue; + if (Nhi == 0) { + // interpolate to value: + y0 = isFltY ? vecy[0].elements.Flt[Nhi] : vecy[0].elements.Int[Nhi]; + y1 = isFltY ? vecy[0].elements.Flt[Nhi+1] : vecy[0].elements.Int[Nhi+1]; + x0 = isFltX ? vecx[0].elements.Flt[Nhi] : vecy[0].elements.Int[Nhi]; + x1 = isFltX ? vecx[0].elements.Flt[Nhi+1] : vecy[0].elements.Int[Nhi+1]; + Xvalue = (value - y0) * (x1 - x0) / (y1 - y0) + x0; + } else { + // interpolate to value: + y0 = isFltY ? vecy[0].elements.Flt[Nhi-1] : vecy[0].elements.Int[Nhi-1]; + y1 = isFltY ? vecy[0].elements.Flt[Nhi] : vecy[0].elements.Int[Nhi]; + x0 = isFltX ? vecx[0].elements.Flt[Nhi-1] : vecy[0].elements.Int[Nhi-1]; + x1 = isFltX ? vecx[0].elements.Flt[Nhi] : vecy[0].elements.Int[Nhi]; + Xvalue = (value - y0) * (x1 - x0) / (y1 - y0) + x0; + } + + set_variable ("threshval", y1); + set_int_variable ("threshbin", Nhi); + set_variable ("threshold", Xvalue); + + if (!QUIET) gprint (GP_LOG, "theshold %f (bin %d is %f)\n", Xvalue, Nhi, y1); + + return (TRUE); +} Index: /trunk/Ohana/src/opihi/cmd.data/write_vectors.c =================================================================== --- /trunk/Ohana/src/opihi/cmd.data/write_vectors.c (revision 33661) +++ /trunk/Ohana/src/opihi/cmd.data/write_vectors.c (revision 33662) @@ -7,5 +7,5 @@ FILE *f; char **fmtlist, *fmttype; - char *p0, *p1, *p2, *format; + char *p0, *p1, *p2, *format, *FITS; Vector **vec; @@ -15,4 +15,12 @@ remove_argument (N, &argc, argv); format = strcreate (argv[N]); + remove_argument (N, &argc, argv); + } + + /* option generate a FITS output table */ + FITS = NULL; + if ((N = get_argument (argc, argv, "-fits"))) { + remove_argument (N, &argc, argv); + FITS = strcreate (argv[N]); remove_argument (N, &argc, argv); } @@ -29,22 +37,9 @@ } - /* open file for outuput */ - if (append) { - f = fopen (argv[1], "a"); - } else { - f = fopen (argv[1], "w"); - } - if (f == (FILE *) NULL) { - gprint (GP_ERR, "can't open file for write\n"); - return (FALSE); - } - /* find number of output vectors */ Nvec = (argc - 2); if (Nvec < 1) { - gprint (GP_ERR, "USAGE: write (file) vector vector ...\n"); - fclose (f); - fflush (f); - return (FALSE); + gprint (GP_ERR, "USAGE: write (file) vector vector ...\n"); + return (FALSE); } ALLOCATE (vec, Vector *, Nvec); @@ -56,6 +51,4 @@ gprint (GP_ERR, "USAGE: write (file) vector vector ...\n"); free (vec); - fclose (f); - fflush (f); return (FALSE); } @@ -68,8 +61,23 @@ gprint (GP_ERR, "error: vectors must all be the same size\n"); free (vec); - fclose (f); - fflush (f); return (FALSE); } + } + + if (FITS) { + int status = WriteVectorTableFITS (argv[1], FITS, vec, Nvec, append, format); + free (vec); + return status; + } + + /* open file for outuput */ + if (append) { + f = fopen (argv[1], "a"); + } else { + f = fopen (argv[1], "w"); + } + if (f == (FILE *) NULL) { + gprint (GP_ERR, "can't open file for write\n"); + return (FALSE); } @@ -88,5 +96,5 @@ fclose (f); free (vec); - fflush (f); + fflush (f); return (TRUE); } @@ -111,5 +119,5 @@ free (format); fclose (f); - fflush (f); + fflush (f); return (FALSE); } Index: /trunk/Ohana/src/opihi/doc/dvo.parallel.txt =================================================================== --- /trunk/Ohana/src/opihi/doc/dvo.parallel.txt (revision 33662) +++ /trunk/Ohana/src/opihi/doc/dvo.parallel.txt (revision 33662) @@ -0,0 +1,12 @@ + +2012.02.17 + + To parallelize dvo queries, like mextract, I need a version of dvo + that takes acts as the remote client. It should take a single dvo + command on the command line and operate on it on the remote host. + eg: + +dvo_client mextract (reset of line) + +it would dump out a result table somewhere (what is loaded?) + Index: /trunk/Ohana/src/opihi/dvo/Makefile =================================================================== --- /trunk/Ohana/src/opihi/dvo/Makefile (revision 33661) +++ /trunk/Ohana/src/opihi/dvo/Makefile (revision 33662) @@ -23,8 +23,8 @@ $(SRC)/compare.$(ARCH).o \ $(SRC)/dvomisc.$(ARCH).o \ +$(SRC)/dvo_host_utils.$(ARCH).o \ $(SRC)/region_list.$(ARCH).o \ $(SRC)/find_matches.$(ARCH).o \ -$(SRC)/photometry.$(ARCH).o \ -$(SRC)/dvo.$(ARCH).o +$(SRC)/photometry.$(ARCH).o cmds = \ @@ -50,6 +50,8 @@ $(SRC)/fitsed.$(ARCH).o \ $(SRC)/gcat.$(ARCH).o \ +$(SRC)/catlist.$(ARCH).o \ $(SRC)/gimages.$(ARCH).o \ $(SRC)/gstar.$(ARCH).o \ +$(SRC)/hosts.$(ARCH).o \ $(SRC)/images.$(ARCH).o \ $(SRC)/imbox.$(ARCH).o \ @@ -71,4 +73,5 @@ $(SRC)/paverage.$(ARCH).o \ $(SRC)/procks.$(ARCH).o \ +$(SRC)/remote.$(ARCH).o \ $(SRC)/skycat.$(ARCH).o \ $(SRC)/skycoverage.$(ARCH).o \ @@ -88,5 +91,12 @@ dvo: $(BIN)/dvo.$(ARCH) $(SRC)/dvo.$(ARCH).o : $(libs) -$(BIN)/dvo.$(ARCH) : $(funcs) $(cmds) +$(BIN)/dvo.$(ARCH) : $(funcs) $(cmds) $(SRC)/dvo.$(ARCH).o + +dvo_client: $(BIN)/dvo_client.$(ARCH) +$(SRC)/dvo_client.$(ARCH).o : $(libs) +$(BIN)/dvo_client.$(ARCH) : $(funcs) $(cmds) $(SRC)/dvo_client.$(ARCH).o + +dvo.install: $(DESTBIN)/dvo +dvo_client.install: $(DESTBIN)/dvo_client install: $(DESTBIN)/dvo help modules @@ -98,5 +108,5 @@ modules: dvo.modules -.PHONY: dvo +.PHONY: dvo dvo_client # deprecated functions: verify & delete Index: /trunk/Ohana/src/opihi/dvo/avextract.c =================================================================== --- /trunk/Ohana/src/opihi/dvo/avextract.c (revision 33661) +++ /trunk/Ohana/src/opihi/dvo/avextract.c (revision 33662) @@ -37,4 +37,34 @@ } + int PARALLEL = FALSE; + if ((N = get_argument (argc, argv, "-parallel"))) { + remove_argument (N, &argc, argv); + PARALLEL = TRUE; + } + + // use remote tables, but not dvo_client.. + int PARALLEL_LOCAL = FALSE; + HostTable *table = NULL; + if ((N = get_argument (argc, argv, "-parallel-local"))) { + remove_argument (N, &argc, argv); + PARALLEL_LOCAL = TRUE; + + char *CATDIR = GetCATDIR(); + if (!CATDIR) { + gprint (GP_ERR, "CATDIR is not set\n"); + return FALSE; + } + SkyTable *sky = GetSkyTable(); + if (!sky) { + gprint (GP_ERR, "failed to load sky table for database\n"); + return FALSE; + } + table = HostTableLoad (CATDIR, sky->hosts); + if (!table) { + gprint (GP_ERR, "ERROR: failure reading Host Table %s for database %s\n", sky->hosts, CATDIR); + return FALSE; + } + } + dvo_catalog_init (&catalog, TRUE); @@ -43,5 +73,7 @@ Nsecfilt = GetPhotcodeNsecfilt (); - // parse skyregion options + // parse skyregion options. NOTE: this is stripped off in parallel operation and always + // defined for the client via the -skyregion option. The dvo_client parses this + // argument in the main program, before it is passed to the command (like mextract) if ((selection = SetRegionSelection (&argc, argv)) == NULL) { print_error(); @@ -85,4 +117,17 @@ /* load region corresponding to selection above */ if ((skylist = SelectRegions (selection)) == NULL) goto escape; + + // this does all the work of re-packaging the command, calling it on the remote machines, then loading in the results + if (PARALLEL && !HOST_ID) { + int status = HostTableParallelOps (argc, argv, RESULT_FILE, 0, VERBOSE); + + dbFreeFields (fields, Nfields); + dbFreeStack (stack, Nstack); + free (stack); + FreeSkyRegionSelection (selection); + dvo_catalog_free (&catalog); + + return status; + } /* create output storage vectors */ @@ -114,6 +159,17 @@ interrupt = FALSE; for (i = 0; (i < skylist[0].Nregions) && !interrupt; i++) { + + // does this host ID match the desired location for the table? + if (!HostTableTestHost(skylist[0].regions[i], HOST_ID)) continue; + /* lock, load, unlock catalog */ - catalog.filename = skylist[0].filename[i]; + char hostfile[1024]; + if (PARALLEL_LOCAL) { + int hostID = (skylist[0].regions[i]->hostFlags & DATA_USE_BCK) ? skylist[0].regions[i]->backupID : skylist[0].regions[i]->hostID; + int seq = table->index[hostID]; + HOSTDIR = table->hosts[seq].pathname; + } + snprintf (hostfile, 1024, "%s/%s.cpt", HOSTDIR, skylist[0].regions[i]->name); + catalog.filename = (HOST_ID || PARALLEL_LOCAL) ? hostfile : skylist[0].filename[i]; catalog.catflags = LOAD_AVES | LOAD_SECF; if (needMeasures) { @@ -169,12 +225,14 @@ interrupt = FALSE; for (n = 0; n < Nreturn; n++) { - vec[n][0].Nelements = Npts; - if (vec[n][0].type == OPIHI_FLT) { - REALLOCATE (vec[n][0].elements.Flt, opihi_flt, MAX(1,Npts)); - } else { - REALLOCATE (vec[n][0].elements.Int, opihi_int, MAX(1,Npts)); - } - } - + ResetVector (vec[n], fields[n].type, Npts); + } + + // write vectors to a table (this is used by parallel dvo operations, but can be used elsewhere) + if (RESULT_FILE) { + int status = WriteVectorTableFITS (RESULT_FILE, "RESULT", vec, Nreturn, FALSE, NULL); + if (!status) goto escape; + } + + if (table) free (table); if (vec) free (vec); if (values) free (values); @@ -187,4 +245,5 @@ escape: + if (table) free (table); if (vec) free (vec); if (values) free (values); Index: /trunk/Ohana/src/opihi/dvo/avmatch.c =================================================================== --- /trunk/Ohana/src/opihi/dvo/avmatch.c (revision 33661) +++ /trunk/Ohana/src/opihi/dvo/avmatch.c (revision 33662) @@ -14,7 +14,7 @@ off_t i, j, n, m, *index; - int N, Ncat, Npts, NPTS, last, Nfields, Nsecfilt; + int N, Ncat, Npts, NPTS, last, Nfields, Nsecfilt, Ninvec; int VERBOSE; - char name[1024]; + char name[1024], *found; void *Signal; float RADIUS; @@ -22,5 +22,5 @@ Catalog catalog; - Vector **vec, *RAvec, *DECvec; + Vector **vec, **invec, *RAvec, *DECvec; dbField *fields; dbValue *values; @@ -29,11 +29,12 @@ /* defaults */ vec = NULL; + invec = NULL; fields = NULL; values = NULL; skylist = NULL; + Ninvec = 0; if ((N = get_argument (argc, argv, "-h"))) goto help; if ((N = get_argument (argc, argv, "--help"))) goto help; - if (argc < 5) goto help; VERBOSE = FALSE; @@ -43,4 +44,21 @@ } + int PARALLEL = FALSE; + if ((N = get_argument (argc, argv, "-parallel"))) { + remove_argument (N, &argc, argv); + PARALLEL = TRUE; + } + + // dump results directly to fits file (esp for parallel dvo) + char *CoordsFile = NULL; + if ((N = get_argument (argc, argv, "-coords"))) { + remove_argument (N, &argc, argv); + CoordsFile = strcreate(argv[N]); + remove_argument (N, &argc, argv); + } + + if (!CoordsFile && (argc < 5)) goto help; + if ( CoordsFile && (argc < 3)) goto help; + dvo_catalog_init (&catalog, TRUE); @@ -49,12 +67,55 @@ Nsecfilt = GetPhotcodeNsecfilt (); + // this does all the work of re-packaging the command, calling it on the remote machines, then loading in the results + RAvec = NULL; + DECvec = NULL; + if (PARALLEL && !HOST_ID) { + if (!CoordsFile) { + // get vectors corresponding to coordinates of interest + if ((RAvec = SelectVector (argv[1], ANYVECTOR, TRUE)) == NULL) goto help; + if ((DECvec = SelectVector (argv[2], ANYVECTOR, TRUE)) == NULL) goto help; + + ALLOCATE (vec, Vector *, 2); + vec[0] = RAvec; + vec[1] = DECvec; + + CoordsFile = abspath("coords.fits", 1024); + int status = WriteVectorTableFITS (CoordsFile, "COORDS", vec, 2, FALSE, NULL); + if (!status) goto escape; + } + + char *targv1 = argv[1]; + char *targv2 = argv[2]; + argv[1] = strcreate ("-coords"); + argv[2] = strcreate (CoordsFile); + free (CoordsFile); + + // I need to pass the RA & DEC vectors to the remote clients... + int status = HostTableParallelOps (argc, argv, RESULT_FILE, RAvec->Nelements, VERBOSE); + if (vec) free (vec); + + free (argv[1]); + free (argv[2]); + argv[1] = targv1; + argv[2] = targv2; + + return status; + } + // get vectors corresponding to coordinates of interest - if ((RAvec = SelectVector (argv[1], ANYVECTOR, TRUE)) == NULL) goto help; - if ((DECvec = SelectVector (argv[2], ANYVECTOR, TRUE)) == NULL) goto help; - RADIUS = atof (argv[3]); - - // strip off RA, DEC, RADIUS arguments - remove_argument (1, &argc, argv); - remove_argument (1, &argc, argv); + if (CoordsFile) { + // read RAvec, DECvec from coords file (1st 2 fields?) + Ninvec = 0; + invec = ReadVectorTableFITS (CoordsFile, "COORDS", &Ninvec); + RAvec = invec[0]; + DECvec = invec[1]; + } else { + if ((RAvec = SelectVector (argv[1], ANYVECTOR, TRUE)) == NULL) goto help; + if ((DECvec = SelectVector (argv[2], ANYVECTOR, TRUE)) == NULL) goto help; + // strip off RA & DEC + remove_argument (1, &argc, argv); + remove_argument (1, &argc, argv); + } + RADIUS = atof (argv[1]); remove_argument (1, &argc, argv); @@ -92,4 +153,6 @@ } ALLOCATE (index, off_t, NPTS); + ALLOCATE (found, char, NPTS); + memset (found, 0, NPTS*sizeof(char)); // grab data from all selected sky regions @@ -97,6 +160,12 @@ interrupt = FALSE; for (i = 0; (i < skylist[0].Nregions) && !interrupt; i++) { + + // does this host ID match the desired location for the table? + if (!HostTableTestHost(skylist[0].regions[i], HOST_ID)) continue; + /* lock, load, unlock catalog */ - catalog.filename = skylist[0].filename[i]; + char hostfile[1024]; + snprintf (hostfile, 1024, "%s/%s.cpt", HOSTDIR, skylist[0].regions[i]->name); + catalog.filename = HOST_ID ? hostfile : skylist[0].filename[i]; catalog.catflags = LOAD_AVES | LOAD_MEAS | LOAD_SECF; catalog.Nsecfilt = 0; @@ -124,5 +193,4 @@ if (Ncat == -2) continue; - // XXX set a 'found' vector to double check we catch everything? m = catalog.average[Ncat].measureOffset; @@ -141,4 +209,5 @@ } } + found[Npts] = TRUE; } dvo_catalog_free (&catalog); @@ -147,6 +216,61 @@ interrupt = FALSE; + // write vectors to a table (this is used by parallel dvo operations, but can be used elsewhere) + // only write the fields which were in a valid catalog + if (RESULT_FILE) { + // extend the array by one to hold index array + REALLOCATE (vec, Vector *, Nfields + 1); + vec[Nfields] = InitVector(); + strcpy (vec[Nfields]->name, "index"); + ResetVector (vec[Nfields], OPIHI_INT, NPTS); + Vector *idxVec = vec[Nfields]; + + // only write out the rows which were found + Npts = 0; + for (i = 0; i < NPTS; i++) { + if (!found[i]) continue; + idxVec->elements.Int[Npts] = i; + Npts ++; + } + int Nfound = Npts; + idxVec->Nelements = Nfound; + + fprintf (stderr, "found %d of %d pts\n", Nfound, NPTS); + + for (i = 0; i < Nfields; i++) { + if (vec[i][0].type == OPIHI_FLT) { + opihi_flt *tmp = NULL; + ALLOCATE (tmp, opihi_flt, Nfound); + Npts = 0; + for (j = 0; j < NPTS; j++) { + if (!found[j]) continue; + tmp[Npts] = vec[i][0].elements.Flt[j]; + Npts++; + } + free (vec[i][0].elements.Flt); + vec[i][0].elements.Flt = tmp; + } else { + opihi_int *tmp = NULL; + ALLOCATE (tmp, opihi_int, Nfound); + Npts = 0; + for (j = 0; j < NPTS; j++) { + if (!found[j]) continue; + tmp[Npts] = vec[i][0].elements.Int[j]; + Npts++; + } + free (vec[i][0].elements.Int); + vec[i][0].elements.Int = tmp; + } + vec[i][0].Nelements = Nfound; + } + int status = WriteVectorTableFITS (RESULT_FILE, "RESULT", vec, Nfields + 1, FALSE, NULL); + free (vec[Nfields]->elements.Int); + free (vec[Nfields]); + if (!status) goto escape; + } + if (vec) free (vec); if (values) free (values); + if (invec) FreeVectorArray (invec, Ninvec); dbFreeFields (fields, Nfields); SkyListFree (skylist); @@ -156,4 +280,5 @@ if (vec) free (vec); if (values) free (values); + if (invec) FreeVectorArray (invec, Ninvec); dbFreeFields (fields, Nfields); SkyListFree (skylist); @@ -162,4 +287,5 @@ help: gprint (GP_ERR, "USAGE: avmatch (RA) (DEC) (RADIUS) field[,field,field...]\n"); + gprint (GP_ERR, " OR: avmatch -coords (filename.fits) (RADIUS) field[,field,field...]\n"); if ((argc > N + 1) && !strcasecmp (argv[N+1], "fields")) { Index: /trunk/Ohana/src/opihi/dvo/catlist.c =================================================================== --- /trunk/Ohana/src/opihi/dvo/catlist.c (revision 33662) +++ /trunk/Ohana/src/opihi/dvo/catlist.c (revision 33662) @@ -0,0 +1,105 @@ +# include "dvoshell.h" + +// find the catalog files which are in a given region, possibly only by host +int catlist (int argc, char **argv) { + + int i, N; + struct stat filestat; + + int VERBOSE = FALSE; + if ((N = get_argument (argc, argv, "-v"))) { + remove_argument (N, &argc, argv); + VERBOSE = TRUE; + } + int ShowAll = FALSE; + if ((N = get_argument (argc, argv, "-all"))) { + remove_argument (N, &argc, argv); + ShowAll = TRUE; + } + int Depth = -1; + if ((N = get_argument (argc, argv, "-depth"))) { + remove_argument (N, &argc, argv); + Depth = atoi (argv[N]); + remove_argument (N, &argc, argv); + } + int ThisHost = FALSE; + if ((N = get_argument (argc, argv, "-this-host"))) { + remove_argument (N, &argc, argv); + ThisHost = TRUE; + } + + // use remote tables, but not dvo_client.. + int PARALLEL_LOCAL = FALSE; + HostTable *table = NULL; + if ((N = get_argument (argc, argv, "-parallel-local"))) { + remove_argument (N, &argc, argv); + PARALLEL_LOCAL = TRUE; + } + + if (argc != 5) { + gprint (GP_ERR, "USAGE: catlist Rmin Rmax Dmin Dmax\n"); + return (FALSE); + } + + float Rmin = atof(argv[1]); + float Rmax = atof(argv[2]); + float Dmin = atof(argv[3]); + float Dmax = atof(argv[4]); + + SkyTable *sky = GetSkyTable (); + SkyList *skylist = SkyListByBounds (sky, Depth, Rmin, Rmax, Dmin, Dmax); + + SkyRegion **regions = skylist[0].regions; + + + if (PARALLEL_LOCAL) { + char *CATDIR = GetCATDIR(); + if (!CATDIR) { + gprint (GP_ERR, "CATDIR is not set\n"); + return FALSE; + } + table = HostTableLoad (CATDIR, sky->hosts); + if (!table) { + gprint (GP_ERR, "ERROR: failure reading Host Table %s for database %s\n", sky->hosts, CATDIR); + return FALSE; + } + } + + int Nregion = 0; + for (i = 0; i < skylist[0].Nregions; i++) { + + // skip tables that are not on this host (if -this-host supplied) + if (ThisHost && !HostTableTestHost (regions[i], HOST_ID)) continue; + + if (PARALLEL_LOCAL) { + int hostID = (skylist[0].regions[i]->hostFlags & DATA_USE_BCK) ? skylist[0].regions[i]->backupID : skylist[0].regions[i]->hostID; + int seq = table->index[hostID]; + HOSTDIR = table->hosts[seq].pathname; + } + char hostfile[1024]; + snprintf (hostfile, 1024, "%s/%s.cpt", HOSTDIR, skylist[0].regions[i]->name); + char *filename = (HOST_ID || PARALLEL_LOCAL) ? hostfile : skylist[0].filename[i]; + + if (ShowAll || (stat (filename, &filestat) != -1)) { + if (VERBOSE) gprint (GP_ERR, "%3d %s %6.2f - %6.2f, %6.2f - %6.2f\n", i, regions[i][0].name, + regions[i][0].Rmin, regions[i][0].Rmax, regions[i][0].Dmin, regions[i][0].Dmax); + + char name[64], nameRmin[64], nameRmax[64], nameDmin[64], nameDmax[64]; + snprintf (name, 64, "region:%d", Nregion); + snprintf (nameRmin, 64, "region_Rmin:%d", Nregion); + snprintf (nameRmax, 64, "region_Rmax:%d", Nregion); + snprintf (nameDmin, 64, "region_Dmin:%d", Nregion); + snprintf (nameDmax, 64, "region_Dmax:%d", Nregion); + set_str_variable (name, regions[i][0].name); + set_variable (nameRmin, regions[i][0].Rmin); + set_variable (nameRmax, regions[i][0].Rmax); + set_variable (nameDmin, regions[i][0].Dmin); + set_variable (nameDmax, regions[i][0].Dmax); + Nregion ++; + } + } + set_int_variable ("region:n", Nregion); + + if (table) free (table); + return (TRUE); +} Index: /trunk/Ohana/src/opihi/dvo/catstats.c =================================================================== --- /trunk/Ohana/src/opihi/dvo/catstats.c (revision 33662) +++ /trunk/Ohana/src/opihi/dvo/catstats.c (revision 33662) @@ -0,0 +1,91 @@ +# include "dvoshell.h" + +// return some stats on the specified catalog files +int catstats (int argc, char **argv) { + + int i, N; + struct stat filestat; + + int VERBOSE = FALSE; + if ((N = get_argument (argc, argv, "-v"))) { + remove_argument (N, &argc, argv); + VERBOSE = TRUE; + } + int ShowAll = FALSE; + if ((N = get_argument (argc, argv, "-all"))) { + remove_argument (N, &argc, argv); + ShowAll = TRUE; + } + int Depth = -1; + if ((N = get_argument (argc, argv, "-depth"))) { + remove_argument (N, &argc, argv); + Depth = atoi (argv[N]); + remove_argument (N, &argc, argv); + } + int ThisHost = FALSE; + if ((N = get_argument (argc, argv, "-this-host"))) { + remove_argument (N, &argc, argv); + ThisHost = TRUE; + } + + // use remote tables, but not dvo_client.. + int PARALLEL_LOCAL = FALSE; + HostTable *table = NULL; + if ((N = get_argument (argc, argv, "-parallel-local"))) { + remove_argument (N, &argc, argv); + PARALLEL_LOCAL = TRUE; + } + + if (argc != 5) { + gprint (GP_ERR, "USAGE: catstats (name)\n"); + return (FALSE); + } + + SkyTable *sky = GetSkyTable (); + SkyList *skylist = SkyListByBounds (sky, Depth, Rmin, Rmax, Dmin, Dmax); + SkyRegion **regions = skylist[0].regions; + + + if (PARALLEL_LOCAL) { + char *CATDIR = GetCATDIR(); + if (!CATDIR) { + gprint (GP_ERR, "CATDIR is not set\n"); + return FALSE; + } + table = HostTableLoad (CATDIR, sky->hosts); + if (!table) { + gprint (GP_ERR, "ERROR: failure reading Host Table %s for database %s\n", sky->hosts, CATDIR); + return FALSE; + } + } + + int Nregion = 0; + for (i = 0; i < skylist[0].Nregions; i++) { + + // skip tables that are not on this host (if -this-host supplied) + if (ThisHost && !HostTableTestHost (regions[i], HOST_ID)) continue; + + if (PARALLEL_LOCAL) { + int hostID = (skylist[0].regions[i]->hostFlags & DATA_USE_BCK) ? skylist[0].regions[i]->backupID : skylist[0].regions[i]->hostID; + int seq = table->index[hostID]; + HOSTDIR = table->hosts[seq].pathname; + } + char hostfile[1024]; + snprintf (hostfile, 1024, "%s/%s.cpt", HOSTDIR, skylist[0].regions[i]->name); + char *filename = (HOST_ID || PARALLEL_LOCAL) ? hostfile : skylist[0].filename[i]; + + if (ShowAll || (stat (filename, &filestat) != -1)) { + if (VERBOSE) gprint (GP_ERR, "%3d %s %6.2f - %6.2f, %6.2f - %6.2f\n", i, regions[i][0].name, + regions[i][0].Rmin, regions[i][0].Rmax, regions[i][0].Dmin, regions[i][0].Dmax); + + char name[64]; + snprintf (name, 64, "region:%d", Nregion); + set_str_variable (name, regions[i][0].name); + Nregion ++; + } + } + set_int_variable ("region:n", Nregion); + + if (table) free (table); + return (TRUE); +} Index: /trunk/Ohana/src/opihi/dvo/dvo.c.in =================================================================== --- /trunk/Ohana/src/opihi/dvo/dvo.c.in (revision 33661) +++ /trunk/Ohana/src/opihi/dvo/dvo.c.in (revision 33662) @@ -29,4 +29,5 @@ set_str_variable ("PROMPT", opihi_prompt); set_str_variable ("RCFILE", opihi_rcfile); + set_int_variable ("DVO_CLIENT", FALSE); { @@ -44,4 +45,9 @@ free (modules); } + + // dvo_client uses the following; they should be NULL for dvo + HOST_ID = 0; + HOSTDIR = NULL; + RESULT_FILE = NULL; return; Index: /trunk/Ohana/src/opihi/dvo/dvo_client.c =================================================================== --- /trunk/Ohana/src/opihi/dvo/dvo_client.c (revision 33662) +++ /trunk/Ohana/src/opihi/dvo/dvo_client.c (revision 33662) @@ -0,0 +1,134 @@ +# include "dvoshell.h" + +int input PROTO((int, char **)); + +/* program-dependent initialization */ +void program_init (int *argc, char **argv) { + + int N; + + /* load the commands used by this implementation */ + InitBasic (); + InitData (); + InitAstro (); + InitDVO (); + + { + char *helpdir; + char *modules; + static char *datadir = "@DATADIR@"; + ALLOCATE (helpdir, char, strlen(datadir) + strlen("/help") + 2); + sprintf (helpdir, "%s/help", datadir); + set_str_variable ("HELPDIR", helpdir); + free (helpdir); + ALLOCATE (modules, char, strlen(datadir) + strlen("/modules") + 2); + sprintf (modules, "%s/modules", datadir); + set_str_variable ("MODULES:0", modules); + set_int_variable ("MODULES:n", 1); + free (modules); + } + + set_int_variable ("DVO_CLIENT", TRUE); + set_int_variable ("SCRIPT", FALSE); + + // dvo_client should have 2 standard arguments: -hostID and -hostdir + HOST_ID = 0; + if ((N = get_argument (*argc, argv, "-hostID"))) { + remove_argument (N, argc, argv); + HOST_ID = atoi (argv[N]);; + remove_argument (N, argc, argv); + } + + HOSTDIR = NULL; + if ((N = get_argument (*argc, argv, "-hostdir"))) { + remove_argument (N, argc, argv); + HOSTDIR = strcreate (argv[N]);; + remove_argument (N, argc, argv); + } + + RESULT_FILE = NULL; + if ((N = get_argument (*argc, argv, "-result"))) { + remove_argument (N, argc, argv); + RESULT_FILE = strcreate (argv[N]);; + remove_argument (N, argc, argv); + } + + // parse -skyregion option (used by most commands) + set_skyregion (0.0, 360.0, -90.0, +90.0); + if ((N = get_argument (*argc, argv, "-skyregion"))) { + if (N + 4 >= *argc) { + gprint (GP_ERR, "USAGE: -skyregion (RA) (RA) (DEC) (DEC)\n"); + exit (1); + } + remove_argument (N, argc, argv); + set_skyregion (atof(argv[N]), atof(argv[N+1]), atof(argv[N+2]), atof(argv[N+3])); + remove_argument (N, argc, argv); + remove_argument (N, argc, argv); + remove_argument (N, argc, argv); + remove_argument (N, argc, argv); + } + + // parse -time option + if ((N = get_argument (*argc, argv, "-time"))) { + if (N + 2 >= *argc) { + gprint (GP_ERR, "USAGE: -time (TimeRef) (TimeFormat)\n"); + exit (1); + } + remove_argument (N, argc, argv); + set_str_variable ("TIMEREF", argv[N]); + remove_argument (N, argc, argv); + set_str_variable ("TIMEFORMAT", argv[N]); + remove_argument (N, argc, argv); + } + + return; +} + +// dvo_client should be called like a dvo command; +int main (int argc, char **argv) { + + // parse out whatever might be needed up front + general_init (&argc, argv); + program_init (&argc, argv); + + if (argc < 2) { + fprintf (stderr, "USAGE: dvo_client (command) (options)\n"); + exit (3); + } + + int status = input (2, &argv[0]); + if (!status) exit (2); + exit (0); + +# if (0) + // identiy the comm + Command *cmd = MatchCommand (argv[1], TRUE, TRUE); + if (cmd == NULL) { + fprintf (stderr, "unknown command %s\n", argv[1]); + exit (1); + } + + // argv[0] is dvo_client; we want to pass in argv[1].. as argv[0] + int status = (*cmd[0].func) (argc - 1, argv + 1); + if (!status) exit (2); + exit (0); +# endif +} + +/* standard welcome message */ +void welcome () { + gprint (GP_ERR, "starting dvo_client...\n"); +} + +/* add program-dependent exit functions here */ +void cleanup () { + QuitKapa (); + ConfigFree (); + + FreeBasic (); + FreeData (); + FreeAstro (); + FreeDVO (); + return; +} + Index: /trunk/Ohana/src/opihi/dvo/dvo_host_utils.c =================================================================== --- /trunk/Ohana/src/opihi/dvo/dvo_host_utils.c (revision 33662) +++ /trunk/Ohana/src/opihi/dvo/dvo_host_utils.c (revision 33662) @@ -0,0 +1,185 @@ +# include "dvoshell.h" + +# define DEBUG 0 +# define PARALLEL_MANUAL 0 +# define PARALLEL_SERIAL 0 +# define MAX_PATH_LENGTH 1024 + +int HostTableLaunchJobs (HostTable *table, char *basecmd, char *options) { + + char uniquer[12]; + int TIME = time(NULL); + int PID = getpid(); + snprintf (uniquer, 12, "%05d.%05d", PID, TIME % 100000); + + int i; + for (i = 0; i < table->Nhosts; i++) { + + // ensure that the paths are absolute path names + char *tmppath = abspath (table->hosts[i].pathname, MAX_PATH_LENGTH); + free (table->hosts[i].pathname); + table->hosts[i].pathname = tmppath; + + // need to save the results filename with the uniquer + // XXX a bit of a waste (but only 1024 * 60 bytes or so + ALLOCATE (table->hosts[i].results, char, MAX_PATH_LENGTH); + snprintf (table->hosts[i].results, MAX_PATH_LENGTH, "%s/dvo.results.%s.fits", table->hosts[i].pathname, uniquer); + + char commandBase[MAX_PATH_LENGTH]; + snprintf (commandBase, MAX_PATH_LENGTH, "dvo.command.%s.txt", uniquer); + char *commandFile = abspath(commandBase, MAX_PATH_LENGTH); + + FILE *f = fopen (commandFile, "w"); + fprintf (f, "%s\n", basecmd); + fclose (f); + + char command[1024]; + snprintf (command, 1024, "dvo_client %s -result %s %s -hostID %d -hostdir %s", commandFile, table->hosts[i].results, options, table->hosts[i].hostID, table->hosts[i].pathname); + free (commandFile); + + if (DEBUG || PARALLEL_MANUAL) fprintf (stderr, "command: %s\n", command); + + if (PARALLEL_MANUAL) { + continue; + } + + if (PARALLEL_SERIAL) { + int status = system (command); + if (status) { + fprintf (stderr, "ERROR running relphot_client\n"); + exit (2); + } + } else { + // launch the job on the remote machine (no handshake) + int errorInfo = 0; + int pid = rconnect ("ssh", table->hosts[i].hostname, command, table->hosts[i].stdio, &errorInfo, FALSE); + if (!pid) { + if (DEBUG) fprintf (stderr, "failure to start %s (error %d)\n", table->hosts[i].hostname, errorInfo); + exit (1); + } + table->hosts[i].pid = pid; // save for future reference + } + } + return TRUE; +} + +// bundle the arguments into a command and pass to dvo_client. this implementation +// expects there to be a result file from the clients, and to load this into vectors in +// the main shell. 'Nelements' is a temp hack : for most commands, the result vectors are +// concatenated, but for avmatch, the vectors are merged by index into a pre-known +// length. this is probably not a solution to a general problem.. +int HostTableParallelOps (int argc, char **argv, char *ResultFile, int Nelements, int VERBOSE) { + + int i; + + // load the list of hosts + SkyTable *sky = GetSkyTable(); + if (!sky) { + gprint (GP_ERR, "failed to load sky table for database\n"); + return FALSE; + } + + char *tmppath = GetCATDIR (); + if (!tmppath) { + gprint (GP_ERR, "failed to get CATDIR for database\n"); + return FALSE; + } + + char *CATDIR = abspath (tmppath, MAX_PATH_LENGTH); + if (!CATDIR) { + gprint (GP_ERR, "failed to make an absolute path from %s (too long)\n", tmppath); + return FALSE; + } + + HostTable *table = HostTableLoad (CATDIR, sky->hosts); + if (!table) { + gprint (GP_ERR, "ERROR: failure reading Host Table %s for database %s\n", sky->hosts, CATDIR); + return FALSE; + } + + // other things I need to append? + char *basecmd = paste_args (argc, argv); + + // determine the sky region + double Rmin, Rmax, Dmin, Dmax; + get_skyregion (&Rmin, &Rmax, &Dmin, &Dmax); + + // determine time reference and format + char *TimeRef = get_variable ("TIMEREF"); + if (!TimeRef) { + gprint (GP_ERR, "failed to find TIMEREF variable\n"); + return FALSE; + } + + char *TimeFormat = get_variable ("TIMEFORMAT"); + if (!TimeFormat) { + gprint (GP_ERR, "failed to find TIMEFORMAT variable\n"); + return FALSE; + } + + char tmp; + char *options = NULL; + int length = snprintf (&tmp, 0, "-D CATDIR %s -time %s %s -skyregion %f %f %f %f", CATDIR, TimeRef, TimeFormat, Rmin, Rmax, Dmin, Dmax); + + ALLOCATE (options, char, length); + snprintf (options, length, "-D CATDIR %s -time %s %s -skyregion %f %f %f %f", CATDIR, TimeRef, TimeFormat, Rmin, Rmax, Dmin, Dmax); + + // launch this command remotely + HostTableLaunchJobs (table, basecmd, options); + free (options); + free (basecmd); + + if (PARALLEL_MANUAL) { + gprint (GP_ERR, "run the relphot_client commands above. when these are done, hit return\n"); + getchar(); + } + if (!PARALLEL_MANUAL && !PARALLEL_SERIAL) { + HostTableWaitJobsGetIO (table, __FILE__, __LINE__, VERBOSE); + } + + // load fields from file + int Nvec = 0; + Vector **vec = NULL; + for (i = 0; i < table->Nhosts; i++) { + if (table->hosts[i].status) continue; + + int Ninvec = 0; + Vector **invec = ReadVectorTableFITS (table->hosts[i].results, "RESULT", &Ninvec); + if (!invec) { + // failed to read the file, now what? + gprint (GP_ERR, "failed to read remote result file : %s\n", table->hosts[i].results); + free (table->hosts[i].results); + table->hosts[i].results = NULL; + continue; + } + free (table->hosts[i].results); + table->hosts[i].results = NULL; + + if (Nelements == 0) { + vec = MergeVectors (vec, &Nvec, invec, Ninvec); + if (vec != invec) { + FreeVectorArray (invec, Ninvec); + } + } else { + vec = MergeVectorsByIndex (vec, &Nvec, invec, Ninvec, Nelements); + FreeVectorArray (invec, Ninvec); + } + } + + // write vectors to a table (this is used by parallel dvo operations, but can be used elsewhere) + if (ResultFile) { + int status = WriteVectorTableFITS (ResultFile, "RESULT", vec, Nvec, FALSE, NULL); + if (!status) { + gprint (GP_ERR, "failed to write result file %s\n", ResultFile); + return FALSE; + } + } + + for (i = 0; i < Nvec; i++) { + AssignVector (vec[i], vec[i]->name, ANYVECTOR, TRUE); + } + free (vec); + + free (table); + return TRUE; +} Index: /trunk/Ohana/src/opihi/dvo/gcat.c =================================================================== --- /trunk/Ohana/src/opihi/dvo/gcat.c (revision 33661) +++ /trunk/Ohana/src/opihi/dvo/gcat.c (revision 33662) @@ -3,5 +3,5 @@ int gcat (int argc, char **argv) { - int i; + int i, N; struct stat filestat; double Ra, Dec, Radius; @@ -10,6 +10,22 @@ void *Signal; + int ShowHost = FALSE; + if ((N = get_argument (argc, argv, "-host"))) { + remove_argument (N, &argc, argv); + ShowHost = TRUE; + } + int ShowBackup = FALSE; + if ((N = get_argument (argc, argv, "-backup"))) { + remove_argument (N, &argc, argv); + ShowBackup = TRUE; + } + int ShowFlags = FALSE; + if ((N = get_argument (argc, argv, "-flags"))) { + remove_argument (N, &argc, argv); + ShowFlags = TRUE; + } + if ((argc != 3) && (argc != 4)) { - gprint (GP_ERR, "USAGE: gcat RA DEC [Radius]\n"); + gprint (GP_ERR, "USAGE: gcat RA DEC [Radius] [-host] [-backup] [-flags]\n"); return (FALSE); } @@ -31,9 +47,26 @@ for (i = 0; (i < skylist[0].Nregions) && !interrupt; i++) { + gprint (GP_ERR, "%3d %s", i, skylist[0].regions[i][0].name); if (stat (skylist[0].filename[i], &filestat) != -1) { - gprint (GP_ERR, "%3d %s *\n", i, skylist[0].regions[i][0].name); + gprint (GP_ERR, " *"); } else { - gprint (GP_ERR, "%3d %s\n", i, skylist[0].regions[i][0].name); + gprint (GP_ERR, " "); } + if (ShowHost) { + gprint (GP_ERR, " %3d", skylist[0].regions[i][0].hostID); + } else { + gprint (GP_ERR, " "); + } + if (ShowBackup) { + gprint (GP_ERR, " %3d", skylist[0].regions[i][0].backupID); + } else { + gprint (GP_ERR, " "); + } + if (ShowFlags) { + gprint (GP_ERR, " 0x%04x", skylist[0].regions[i][0].hostFlags); + } else { + gprint (GP_ERR, " "); + } + gprint (GP_ERR, "\n"); set_str_variable ("CATNAME", skylist[0].filename[i]); } Index: /trunk/Ohana/src/opihi/dvo/gstar.c =================================================================== --- /trunk/Ohana/src/opihi/dvo/gstar.c (revision 33661) +++ /trunk/Ohana/src/opihi/dvo/gstar.c (revision 33662) @@ -68,4 +68,10 @@ } + int PARALLEL = FALSE; + if ((N = get_argument (argc, argv, "-parallel"))) { + remove_argument (N, &argc, argv); + PARALLEL = TRUE; + } + if (argc != 4) { gprint (GP_ERR, "USAGE: gstar RA DEC Radius [-m]\n"); @@ -89,6 +95,25 @@ } + HostTable *table = NULL; + if (PARALLEL) { + char *CATDIR = GetCATDIR(); + table = HostTableLoad (CATDIR, sky->hosts); + if (!table) { + gprint (GP_ERR, "ERROR: failure reading Host Table %s for database %s\n", sky->hosts, CATDIR); + return FALSE; + } + + SkyRegion *region = skylist[0].regions[0]; + int hostID = (region->hostFlags & DATA_USE_BCK) ? region->backupID : region->hostID; + int index = table->index[hostID]; + + char hostfile[1024]; + snprintf (hostfile, 1024, "%s/%s.cpt", table->hosts[index].pathname, region->name); + catalog.filename = hostfile; + } else { + catalog.filename = skylist[0].filename[0]; + } + /* lock, load, unlock catalog */ - catalog.filename = skylist[0].filename[0]; catalog.catflags = GetMeasures ? LOAD_AVES | LOAD_MEAS | LOAD_SECF : LOAD_AVES | LOAD_SECF; catalog.Nsecfilt = 0; @@ -177,5 +202,13 @@ gprint (GP_LOG, "%f ", catalog.average[k].duD); gprint (GP_LOG, "%f ", catalog.average[k].P); - gprint (GP_LOG, "%f ", catalog.average[k].dP); + gprint (GP_LOG, "%f ", catalog.average[k].dP); + + gprint (GP_LOG, "%f ", catalog.average[k].ChiSqAve); + gprint (GP_LOG, "%f ", catalog.average[k].ChiSqPM); + gprint (GP_LOG, "%f ", catalog.average[k].ChiSqPar); + + date = ohana_sec_to_date (catalog.average[k].Tmean); + gprint (GP_LOG, "%20s ", date); + gprint (GP_LOG, "%f ", catalog.average[k].Trange / 86400.0); } @@ -197,6 +230,24 @@ for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, 2); gprint (GP_LOG, "\n"); - } - + + if (FULL_OUTPUT) { + /* Map */ + for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, 5); + gprint (GP_LOG, "\n"); + + /* M_20 */ + for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, 6); + gprint (GP_LOG, "\n"); + + /* M_80 */ + for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, 7); + gprint (GP_LOG, "\n"); + + /* secfilt flags */ + for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, 4); + gprint (GP_LOG, "\n"); + + } + } if (GetMeasures || SaveVectors) { m = catalog.average[k].measureOffset; @@ -218,6 +269,6 @@ gprint (GP_LOG, "%7.4f ", catalog.measure[m].dR); gprint (GP_LOG, "%7.4f ", catalog.measure[m].dD); - gprint (GP_LOG, "%4x ", catalog.measure[m].photFlags); - gprint (GP_LOG, "%3x ", catalog.measure[m].dbFlags); + gprint (GP_LOG, "0x%08x ", catalog.measure[m].photFlags); + gprint (GP_LOG, "0x%08x ", catalog.measure[m].dbFlags); gprint (GP_LOG, "%5d ", catalog.measure[m].photcode); gprint (GP_LOG, "%-20s ", GetPhotcodeNamebyCode (catalog.measure[m].photcode)); @@ -226,23 +277,23 @@ if (FULL_OUTPUT) { - gprint (GP_LOG, "%f ", catalog.measure[m].Mcal); - gprint (GP_LOG, "%f ", catalog.measure[m].Map); - gprint (GP_LOG, "%f ", pow(10.0, 0.4*catalog.measure[m].dt)); - gprint (GP_LOG, "%f ", 1.0 + catalog.measure[m].airmass); - gprint (GP_LOG, "%f ", catalog.measure[m].az); - gprint (GP_LOG, "%f ", catalog.measure[m].Xccd); - gprint (GP_LOG, "%f ", catalog.measure[m].Yccd); - gprint (GP_LOG, "%f ", FromShortPixels(catalog.measure[m].dXccd)); - gprint (GP_LOG, "%f ", FromShortPixels(catalog.measure[m].dYccd)); - gprint (GP_LOG, "%f ", catalog.measure[m].Sky); - gprint (GP_LOG, "%f ", catalog.measure[m].dSky); - gprint (GP_LOG, "%d ", catalog.measure[m].averef); - gprint (GP_LOG, "%d ", catalog.measure[m].detID); - gprint (GP_LOG, "%d ", catalog.measure[m].imageID); - gprint (GP_LOG, "%f ", catalog.measure[m].psfQual); - gprint (GP_LOG, "%f ", catalog.measure[m].psfChisq); - gprint (GP_LOG, "%f ", catalog.measure[m].crNsigma); - gprint (GP_LOG, "%f ", catalog.measure[m].extNsigma); - gprint (GP_LOG, "%f ", FromShortDegrees(catalog.measure[m].theta)); + gprint (GP_LOG, "%6.3f ", catalog.measure[m].Mcal); + gprint (GP_LOG, "%6.3f ", catalog.measure[m].Map); + gprint (GP_LOG, "%5.1f ", pow(10.0, 0.4*catalog.measure[m].dt)); + gprint (GP_LOG, "%5.3f ", catalog.measure[m].airmass); + gprint (GP_LOG, "%6.1f ", catalog.measure[m].az); + gprint (GP_LOG, "%6.1f ", catalog.measure[m].Xccd); + gprint (GP_LOG, "%6.1f ", catalog.measure[m].Yccd); + gprint (GP_LOG, "%3.1f ", FromShortPixels(catalog.measure[m].dXccd)); + gprint (GP_LOG, "%3.1f ", FromShortPixels(catalog.measure[m].dYccd)); + gprint (GP_LOG, "%6.1f ", catalog.measure[m].Sky); + gprint (GP_LOG, "%5.1f ", catalog.measure[m].dSky); + gprint (GP_LOG, "%8d ", catalog.measure[m].averef); + gprint (GP_LOG, "0x%08x ", catalog.measure[m].detID); + gprint (GP_LOG, "0x%08x ", catalog.measure[m].imageID); + gprint (GP_LOG, "%.3f ", catalog.measure[m].psfQual); + gprint (GP_LOG, "%.1f ", catalog.measure[m].psfChisq); + gprint (GP_LOG, "%.1f ", catalog.measure[m].crNsigma); + gprint (GP_LOG, "%.1f ", catalog.measure[m].extNsigma); + gprint (GP_LOG, "%.1f ", FromShortDegrees(catalog.measure[m].theta)); } gprint (GP_LOG, "\n"); @@ -295,14 +346,14 @@ void print_double (double value) { if (isnan(value)) - gprint (GP_LOG, "NaN "); + gprint (GP_LOG, " NaN "); else - gprint (GP_LOG, "%6.3f ", value); + gprint (GP_LOG, " %6.3f ", value); } void print_short (double value, short int ival) { if (ival == NAN_S_SHORT) - gprint (GP_LOG, "NaN "); + gprint (GP_LOG, " NaN "); else - gprint (GP_LOG, "%6.3f ", value); + gprint (GP_LOG, " %6.3f ", value); } @@ -370,6 +421,39 @@ case 3: /* filter names */ code = GetPhotcodebyNsec (seq); - gprint (GP_LOG, "%-6s ", code[0].name); - break; - } -} + gprint (GP_LOG, " %-6s ", code[0].name); + break; + + case 4: /* secfilt flags */ + if (seq == -1) { + gprint (GP_LOG, "%08x ", 0); + } else { + gprint (GP_LOG, "%08x ", secfilt[seq].flags); + } + break; + + case 5: /* average ap mags */ + if (seq == -1) { + print_double (NAN); + } else { + print_double (secfilt[seq].Map); + } + break; + + case 6: /* average ap mags */ + if (seq == -1) { + print_double (NAN); + } else { + print_short (0.001*secfilt[seq].M_20, secfilt[seq].M_20); + } + break; + + case 7: /* average ap mags */ + if (seq == -1) { + print_double (NAN); + } else { + print_short (0.001*secfilt[seq].M_20, secfilt[seq].M_80); + } + break; + + } +} Index: /trunk/Ohana/src/opihi/dvo/hosts.c =================================================================== --- /trunk/Ohana/src/opihi/dvo/hosts.c (revision 33662) +++ /trunk/Ohana/src/opihi/dvo/hosts.c (revision 33662) @@ -0,0 +1,103 @@ +# include "dvoshell.h" +# include +# define MAX_PATH_LENGTH 1024 + +// functions to manage the remote hosts +int hosts (int argc, char **argv) { + + if (argc < 2) { + gprint (GP_ERR, "USAGE: hosts (command) [options]\n"); + gprint (GP_ERR, " commands:\n"); + gprint (GP_ERR, " purge-temp : delete all tempfiles for this shell\n"); + return FALSE; + } + + if (!strncmp(argv[1], "purge-temp", MAX(strlen(argv[1]), 3))) { + glob_t pglob; + int PID = getpid(); + int N; + if ((N = get_argument (argc, argv, "-old-pid"))) { + remove_argument (N, &argc, argv); + PID = atoi(argv[N]); + remove_argument (N, &argc, argv); + } + + int VERBOSE = FALSE; + if ((N = get_argument (argc, argv, "-v"))) { + remove_argument (N, &argc, argv); + VERBOSE = TRUE; + } + if ((N = get_argument (argc, argv, "-verbose"))) { + remove_argument (N, &argc, argv); + VERBOSE = TRUE; + } + + int DRYRUN = TRUE; + if ((N = get_argument (argc, argv, "-commit"))) { + remove_argument (N, &argc, argv); + DRYRUN = FALSE; + } + + // XXX wrap this up in a function: + char *CATDIR = GetCATDIR(); + if (!CATDIR) { + gprint (GP_ERR, "CATDIR is not set\n"); + return FALSE; + } + SkyTable *sky = GetSkyTable(); + if (!sky) { + gprint (GP_ERR, "failed to load sky table for database\n"); + return FALSE; + } + HostTable *table = HostTableLoad (CATDIR, sky->hosts); + if (!table) { + gprint (GP_ERR, "ERROR: failure reading Host Table %s for database %s\n", sky->hosts, CATDIR); + return FALSE; + } + + int i; + for (i = 0; i < table->Nhosts; i++) { + pglob.gl_offs = 0; + char name[MAX_PATH_LENGTH]; + snprintf (name, MAX_PATH_LENGTH, "%s/dvo.results.%05d.*.fits", table->hosts[i].pathname, PID); + if (VERBOSE) gprint (GP_ERR, "checking %s\n", name); + glob (name, 0, NULL, &pglob); + int j; + for (j = 0; j < pglob.gl_pathc; j++) { + if (VERBOSE) gprint (GP_ERR, "unlink %s\n", pglob.gl_pathv[j]); + if (!DRYRUN) unlink (pglob.gl_pathv[j]); + } + globfree (&pglob); + } + return TRUE; + } + + if (!strncmp(argv[1], "get.results", MAX(strlen(argv[1]), 3))) { + int N; + + char *varname = NULL; + if ((N = get_argument (argc, argv, "-var"))) { + remove_argument (N, &argc, argv); + varname = strcreate(argv[N]); + remove_argument (N, &argc, argv); + } + + if (argc < 2) { + gprint (GP_ERR, "USAGE: hosts get.results [-var var]\n"); + return FALSE; + } + + if (varname) { + set_str_variable (varname, RESULT_FILE); + free (varname); + } else { + gprint (GP_LOG, "results: %s\n", RESULT_FILE); + } + return TRUE; + } + + + + gprint (GP_ERR, "error: unknown hosts command %s\n", argv[1]); + return FALSE; +} Index: /trunk/Ohana/src/opihi/dvo/init.c =================================================================== --- /trunk/Ohana/src/opihi/dvo/init.c (revision 33661) +++ /trunk/Ohana/src/opihi/dvo/init.c (revision 33662) @@ -23,8 +23,10 @@ int fitsed PROTO((int, char **)); int gcat PROTO((int, char **)); +int catlist PROTO((int, char **)); int getxtra PROTO((int, char **)); int gimages PROTO((int, char **)); int gstar PROTO((int, char **)); int gtypes PROTO((int, char **)); +int hosts PROTO((int, char **)); int images PROTO((int, char **)); int imbox PROTO((int, char **)); @@ -47,4 +49,5 @@ int paverage PROTO((int, char **)); int procks PROTO((int, char **)); +int remote PROTO((int, char **)); int showtile PROTO((int, char **)); int skycat PROTO((int, char **)); @@ -77,6 +80,8 @@ {1, "fitsed", fitsed, "fit stellar SEDs to objects"}, {1, "gcat", gcat, "get catalog at location"}, + {1, "catlist", catlist, "get list of catalogs for region / host"}, {1, "gimages", gimages, "get images at location"}, {1, "gstar", gstar, "get star statistics"}, + {1, "hosts", hosts, "remote host support functions"}, {1, "images", images, "plot image boxes"}, {1, "imbox", imbox, "plot expected image box"}, @@ -99,4 +104,5 @@ {1, "paverage", paverage, "plot average magnitude"}, {1, "procks", procks, "plot rocks"}, + {1, "remote", remote, "generic remote dvo client operation"}, {1, "showtile", showtile, "plot tile pattern"}, {1, "skycat", skycat, "show sky catalog boundaries"}, Index: /trunk/Ohana/src/opihi/dvo/mextract.c =================================================================== --- /trunk/Ohana/src/opihi/dvo/mextract.c (revision 33661) +++ /trunk/Ohana/src/opihi/dvo/mextract.c (revision 33662) @@ -25,4 +25,5 @@ selection = NULL; + // fprintf (stderr, "start..."); if ((N = get_argument (argc, argv, "-h"))) goto help; if ((N = get_argument (argc, argv, "--help"))) goto help; @@ -34,4 +35,20 @@ } + int PARALLEL = FALSE; + if ((N = get_argument (argc, argv, "-parallel"))) { + remove_argument (N, &argc, argv); + PARALLEL = TRUE; + } + + // this is used to NOT save the results in the results file + // use this option when mextract is used in a script which does its + // own job of packaging the results + int SKIP_RESULTS = FALSE; + if ((N = get_argument (argc, argv, "-skip-results"))) { + remove_argument (N, &argc, argv); + SKIP_RESULTS = TRUE; + } + + // init here so free in 'escape' block does not crash dvo_catalog_init (&catalog, TRUE); @@ -40,11 +57,16 @@ Nsecfilt = GetPhotcodeNsecfilt (); + // parse skyregion options. NOTE: this is stripped off in parallel operation and always + // defined for the client via the -skyregion option. The dvo_client parses this + // argument in the main program, before it is passed to the command (like mextract) + if ((selection = SetRegionSelection (&argc, argv)) == NULL) { + print_error(); + goto escape; + } + // init locally static variables (time refs) dbExtractMeasuresInit(); - // parse skyregion options - if ((selection = SetRegionSelection (&argc, argv)) == NULL) goto escape; - - // command-line is of the form: avextract field,field, field [where (field op value)...] + // command-line is of the form: mextract field,field, field [where (field op value)...] // parse the fields to be extracted and returned @@ -84,13 +106,31 @@ if ((skylist = SelectRegions (selection)) == NULL) goto escape; + // this does all the work of re-packaging the command, calling it on the remote machines, then loading in the results + if (PARALLEL && !HOST_ID) { + int status = HostTableParallelOps (argc, argv, RESULT_FILE, 0, VERBOSE); + + dbFreeFields (fields, Nfields); + dbFreeStack (stack, Nstack); + free (stack); + FreeSkyRegionSelection (selection); + dvo_catalog_free (&catalog); + + return status; + } + // load image data if needed (for fields listed below) loadImages = FALSE; mosaicMode = FALSE; for (i = 0; !loadImages && (i < Nfields); i++) { - if (fields[i].ID == MEAS_XCCD) loadImages = TRUE; // Are these still needed? Xccd and Yccd are in measurement - if (fields[i].ID == MEAS_YCCD) loadImages = TRUE; // and dbExtractMeasures doesn't look at the image for these - if (fields[i].ID == MEAS_XMOSAIC) loadImages = mosaicMode = TRUE; - if (fields[i].ID == MEAS_YMOSAIC) loadImages = mosaicMode = TRUE; + if (!MEASURE_HAS_XCCD) { + // I'm keeping this code because it gives a way of handling dvo dbs that don't have + // measure.xccd if we need it + if (fields[i].ID == MEAS_XCCD) loadImages = TRUE; + if (fields[i].ID == MEAS_YCCD) loadImages = TRUE; + } + if (fields[i].ID == MEAS_XMOSAIC) loadImages = mosaicMode = TRUE; + if (fields[i].ID == MEAS_YMOSAIC) loadImages = mosaicMode = TRUE; if (fields[i].ID == MEAS_EXTERN_ID) loadImages = mosaicMode = TRUE; + if (fields[i].ID == MEAS_FLAT) loadImages = mosaicMode = TRUE; } if (loadImages && !SetImageSelection (mosaicMode, selection)) goto escape; @@ -115,7 +155,15 @@ interrupt = FALSE; + // fprintf (stderr, "done setup..."); + for (i = 0; (i < skylist[0].Nregions) && !interrupt; i++) { + + // does this host ID match the desired location for the table? + if (!HostTableTestHost(skylist[0].regions[i], HOST_ID)) continue; + /* lock, load, unlock catalog */ - catalog.filename = skylist[0].filename[i]; + char hostfile[1024]; + snprintf (hostfile, 1024, "%s/%s.cpt", HOSTDIR, skylist[0].regions[i]->name); + catalog.filename = HOST_ID ? hostfile : skylist[0].filename[i]; catalog.catflags = LOAD_AVES | LOAD_MEAS | LOAD_SECF; catalog.Nsecfilt = Nsecfilt; @@ -126,5 +174,5 @@ if (!dvo_catalog_open (&catalog, NULL, FALSE, "r")) { gprint (GP_ERR, "ERROR: failure to open catalog file %s\n", catalog.filename); - exit (2); + return FALSE; } dvo_catalog_unlock (&catalog); @@ -132,12 +180,23 @@ /* XXX need to call dvo_catalog_chipcoords here passing the loaded images */ + // fprintf (stderr, "done read..."); + for (j = 0; (j < catalog.Naverage) && !interrupt; j++) { m = catalog.average[j].measureOffset; - dbExtractMeasuresInitAve (); // reset counters for saved fields + if (m > catalog.Nmeasure) { + gprint (GP_ERR, "ERROR: inconsistent average->measure offset. Unsorted database?\n"); + goto escape; + } + + dbExtractMeasuresInitAve (); // reset counters for saved fields (costs very little for (k = 0; (k < catalog.average[j].Nmeasure); k++, m++) { + if (catalog.measure[m].averef != j) { + gprint (GP_ERR, "ERROR: inconsistent measure->average link. Unsorted database?\n"); + goto escape; + } // extract the relevant values for this measurement - dbExtractMeasuresInitMeas (); // reset counters for saved fields + dbExtractMeasuresInitMeas (); // reset counters for saved fields (costs very little for (n = 0; n < Nfields; n++) { values[n] = dbExtractMeasures (&catalog.average[j], &catalog.secfilt[j*Nsecfilt], &catalog.measure[m], &fields[n]); @@ -173,6 +232,14 @@ interrupt = FALSE; + // fprintf (stderr, "done load..."); + for (n = 0; n < Nreturn; n++) { - ResetVector (vec[n], fields[n].type, MAX(1,Npts)); + ResetVector (vec[n], fields[n].type, Npts); + } + + // write vectors to a table (this is used by parallel dvo operations, but can be used elsewhere) + if (RESULT_FILE && !SKIP_RESULTS) { + int status = WriteVectorTableFITS (RESULT_FILE, "RESULT", vec, Nreturn, FALSE, NULL); + if (!status) goto escape; } @@ -185,4 +252,6 @@ SkyListFree (skylist); FreeSkyRegionSelection (selection); + + // fprintf (stderr, "done extr...\n"); return (TRUE); @@ -307,4 +376,7 @@ gprint (GP_ERR, " imageID : ID of source image (32 bit)\n"); gprint (GP_ERR, " externID : externID of source image (32 bit)\n"); + + gprint (GP_ERR, " Mcal_offset : difference wrt nominal zero point (clouds are positive)\n"); + gprint (GP_ERR, " flat : flat-field correction (measure.Mcal - image.Mcal)\n"); return (FALSE); } Index: /trunk/Ohana/src/opihi/dvo/region_list.c =================================================================== --- /trunk/Ohana/src/opihi/dvo/region_list.c (revision 33661) +++ /trunk/Ohana/src/opihi/dvo/region_list.c (revision 33662) @@ -113,5 +113,5 @@ /* check for Region selection from display */ if ((N = get_argument (*argc, argv, "-skyregion"))) { - if (N + 5 >= *argc) { + if (N + 4 >= *argc) { gprint (GP_ERR, "USAGE: -skyregion (RA) (RA) (DEC) (DEC)\n"); FreeSkyRegionSelection (selection); @@ -134,4 +134,47 @@ /* given possible options (by name, by list, by graph region), select SkyRegions */ +int SetSkyRegions (SkyRegionSelection *selection) { + + if (selection->name != NULL) { + gprint (GP_ERR, "name-based selection not yet implemented (in parallel mode)\n"); + return FALSE; + } + + if (selection->list != NULL) { + gprint (GP_ERR, "list-based selection not yet implemented (in parallel mode)\n"); + return FALSE; + } + + if (selection->useDisplay) { + double Rmin, Rmax, Dmin, Dmax, Radius; + Graphdata graphsky; + + if (!GetGraphdata (&graphsky, NULL, NULL)) { + gprint (GP_ERR, "region display not available\n"); + return FALSE; + } + Radius = MAX (fabs(graphsky.xmax), fabs(graphsky.ymax)); + Dmin = graphsky.coords.crval2 - Radius; + Dmax = graphsky.coords.crval2 + Radius; + + if ((Dmin <= -89) || (Dmax >= 89)) { + Rmin = 0; + Rmax = 360; + } else { + double Rmod = MAX (Radius / (cos(Dmin*RAD_DEG)), Radius / (cos(Dmax*RAD_DEG))); + Rmin = graphsky.coords.crval1 - Rmod; + Rmax = graphsky.coords.crval1 + Rmod; + } + + set_skyregion (Rmin, Rmax, Dmin, Dmax); + return TRUE; + } + if (selection->useSkyregion) { + return TRUE; + } + return FALSE; +} + +/* given possible options (by name, by list, by graph region), select SkyRegions */ SkyList *SelectRegions (SkyRegionSelection *selection) { @@ -140,16 +183,5 @@ /* determine region-file names */ if (selection->name != NULL) { - char filename[256]; - char *CATDIR = dvo_get_catdir(); - - ALLOCATE (skylist, SkyList, 1); - ALLOCATE (skylist[0].regions, SkyRegion *, 1); - ALLOCATE (skylist[0].regions[0], SkyRegion, 1); - ALLOCATE (skylist[0].filename, char *, 1); - skylist[0].ownElements = TRUE; // free these elements when freeing the list - - strcpy (skylist[0].regions[0][0].name, selection->name); - sprintf (filename, "%s/%s.cpt", CATDIR, selection->name); - skylist[0].filename[0] = strcreate (filename); + skylist = SkyListByName (sky, selection->name); return (skylist); } Index: /trunk/Ohana/src/opihi/dvo/remote.c =================================================================== --- /trunk/Ohana/src/opihi/dvo/remote.c (revision 33662) +++ /trunk/Ohana/src/opihi/dvo/remote.c (revision 33662) @@ -0,0 +1,34 @@ +# include "dvoshell.h" +# include +# define MAX_PATH_LENGTH 1024 + +// functions to manage the remote hosts +int remote (int argc, char **argv) { + + int N; + + if (argc < 2) { + gprint (GP_ERR, "USAGE: remote (command)\n"); + gprint (GP_ERR, " launch (command) on the parallel hosts\n"); + gprint (GP_ERR, " options:\n"); + gprint (GP_ERR, " -v : verbose mode:\n"); + return FALSE; + } + + // we can call any command remotely, but the collection of macros will + // not automatically be passed along. if we want to run a specific macro, + // need to point at the relevant input file and have that get loaded + + // if we specified a remote result file, the function above assumes that this is a FITS table + // with a set of vectors to load. + + int VERBOSE = FALSE; + if ((N = get_argument (argc, argv, "-v"))) { + remove_argument (N, &argc, argv); + VERBOSE = TRUE; + } + + // strip of the 'remote' and send the remaining arguments to the remote machine + int status = HostTableParallelOps (argc - 1, &argv[1], NULL, 0, VERBOSE); + return status; +} Index: /trunk/Ohana/src/opihi/dvo/skycoverage.c =================================================================== --- /trunk/Ohana/src/opihi/dvo/skycoverage.c (revision 33661) +++ /trunk/Ohana/src/opihi/dvo/skycoverage.c (revision 33662) @@ -1,9 +1,12 @@ # include "dvoshell.h" +// enum to define possible modes +enum {COVERAGE, DENSITY, MIN_UBERCAL, MIN_DMAG_SYS, MIN_MCAL, MAX_MCAL}; + int wordhash (char *word); int skycoverage (int argc, char **argv) { - int WITH_MOSAIC, SOLO_MOSAIC, ShowDensity; + int WITH_MOSAIC, SOLO_MOSAIC, mode; off_t i, Nimage; int N, status, TimeSelect, ByName, xs, ys; @@ -60,8 +63,24 @@ } - ShowDensity = FALSE; + mode = COVERAGE; if ((N = get_argument (argc, argv, "-density"))) { remove_argument (N, &argc, argv); - ShowDensity = TRUE; + mode = DENSITY; + } + if ((N = get_argument (argc, argv, "-min-ubercal"))) { + remove_argument (N, &argc, argv); + mode = MIN_UBERCAL; + } + if ((N = get_argument (argc, argv, "-min-dmag-sys"))) { + remove_argument (N, &argc, argv); + mode = MIN_DMAG_SYS; + } + if ((N = get_argument (argc, argv, "-min-mcal"))) { + remove_argument (N, &argc, argv); + mode = MIN_MCAL; + } + if ((N = get_argument (argc, argv, "-max-mcal"))) { + remove_argument (N, &argc, argv); + mode = MAX_MCAL; } @@ -194,4 +213,5 @@ V = (float *)buf[0].matrix.buffer; + // init the V array (depends on mode) for (ys = 0; ys < Ny; ys++) { for (xs = 0; xs < Nx; xs++) { @@ -199,8 +219,35 @@ status &= (r >= 0); status &= (r <= 360); + // are we in a part of the projection covering the sky or not? if (status) { - V[ys*Nx + xs] = ShowDensity ? 0 : 2; + switch (mode) { + case COVERAGE: + V[ys*Nx + xs] = 2; + break; + case DENSITY: + V[ys*Nx + xs] = 0; + break; + case MIN_UBERCAL: + case MIN_DMAG_SYS: + case MIN_MCAL: + V[ys*Nx + xs] = 1E6; + break; + case MAX_MCAL: + V[ys*Nx + xs] = -1E6; + break; + } } else { - V[ys*Nx + xs] = ShowDensity ? NAN : 0; + switch (mode) { + case COVERAGE: + V[ys*Nx + xs] = 0; + break; + case DENSITY: + case MIN_UBERCAL: + case MIN_DMAG_SYS: + case MIN_MCAL: + case MAX_MCAL: + V[ys*Nx + xs] = NAN; + break; + } } } @@ -257,8 +304,23 @@ xs = (int)Xs; ys = (int)Ys; - if (ShowDensity) { + switch (mode) { + case COVERAGE: + V[ys*Nx + xs] = 1; + break; + case DENSITY: V[ys*Nx + xs] += 1; - } else { - V[ys*Nx + xs] = 1; + break; + case MIN_UBERCAL: + V[ys*Nx + xs] = MIN(V[ys*Nx + xs], image[i].ubercalDist); + break; + case MIN_DMAG_SYS: + V[ys*Nx + xs] = MIN(V[ys*Nx + xs], image[i].dMagSys); + break; + case MIN_MCAL: + V[ys*Nx + xs] = MIN(V[ys*Nx + xs], image[i].Mcal); + break; + case MAX_MCAL: + V[ys*Nx + xs] = MAX(V[ys*Nx + xs], image[i].Mcal); + break; } } Index: /trunk/Ohana/src/opihi/dvo/skyregion.c =================================================================== --- /trunk/Ohana/src/opihi/dvo/skyregion.c (revision 33661) +++ /trunk/Ohana/src/opihi/dvo/skyregion.c (revision 33662) @@ -1,21 +1,31 @@ # include "dvoshell.h" - -#ifdef NO_MOVED_TO_DVO -static double RAs = 0.0; -static double RAe = 0.0; -static double DECs = 0.0; -static double DECe = 0.0; -#endif // define the sky region for which extractions are limited int skyregion (int argc, char **argv) { + int N; + + // dvo_client should have 2 standard arguments: -hostID and -hostdir + int SaveRegion = FALSE; + if ((N = get_argument (argc, argv, "-save"))) { + remove_argument (N, &argc, argv); + SaveRegion = TRUE; + } + if (argc == 1) { - double RAs, RAe, DECs, DECe; - get_skyregion(&RAs, &RAe, &DECs, &DECe); + double Rmin, Rmax, Dmin, Dmax; + get_skyregion(&Rmin, &Rmax, &Dmin, &Dmax); - gprint (GP_ERR, "current skyregion: %f - %f : %f - %f\n", RAs, RAe, DECs, DECe); - gprint (GP_ERR, "USAGE: skyregion (min RA) (max RA) (min DEC) (max DEC)\n"); - return (FALSE); + if (SaveRegion) { + set_variable ("Rmin", Rmin); + set_variable ("Rmax", Rmax); + set_variable ("Dmin", Dmin); + set_variable ("Dmax", Dmax); + return TRUE; + } else { + gprint (GP_ERR, "current skyregion: %f - %f : %f - %f\n", Rmin, Rmax, Dmin, Dmax); + gprint (GP_ERR, "USAGE: skyregion (min RA) (max RA) (min DEC) (max DEC)\n"); + return FALSE; + } } @@ -29,25 +39,4 @@ return (TRUE); } -#ifdef NOT_MOVED_TO_LIBDVO -int get_skyregion (double *Rs, double *Re, double *Ds, double *De) { - - *Rs = RAs; - *Re = RAe; - *Ds = DECs; - *De = DECe; - - return TRUE; -} - -int set_skyregion (double Rs, double Re, double Ds, double De) { - - RAs = Rs; - RAe = Re; - DECs = Ds; - DECe = De; - - return TRUE; -} -#endif /* find region which overlaps c at given depth (-1 : populated ) */ @@ -107,4 +96,6 @@ ALLOCATE (new[0].regions, SkyRegion *, 1); ALLOCATE (new[0].filename, char *, 1); + new[0].Nregions = 0; + new[0].ownElements = FALSE; // this list is only holding a view to the elements // output list Index: /trunk/Ohana/src/opihi/include/dvomath.h =================================================================== --- /trunk/Ohana/src/opihi/include/dvomath.h (revision 33661) +++ /trunk/Ohana/src/opihi/include/dvomath.h (revision 33662) @@ -111,6 +111,8 @@ /* vector handling */ +void InitVectors PROTO((void)); +void FreeVectors PROTO((void)); Vector *InitVector PROTO((void)); -void InitVectors PROTO((void)); +void FreeVectorArray PROTO((Vector **vec, int Nvec)); int CopyVector PROTO((Vector *out, Vector *in)); int ResetVector PROTO((Vector *vec, char type, int Nelements)); @@ -126,5 +128,13 @@ int IsVectorPtr PROTO((Vector *vec)); int ListVectors PROTO((void)); +int ListVectorsToList PROTO((char *name)); Vector *SelectVector PROTO((char *name, int mode, int verbose)); +int AssignVector PROTO((Vector *vec, char *name, int mode, int verbose)); +Vector **MergeVectors PROTO((Vector **vec, int *Nvec, Vector **invec, int Ninvec)); +Vector **MergeVectorsByIndex PROTO((Vector **vec, int *Nvec, Vector **invec, int Ninvec, int Nelements)); + +/* vector IO functions */ +int WriteVectorTableFITS PROTO((char *filename, char *extname, Vector **vec, int Nvec, int append, char *format)); +Vector **ReadVectorTableFITS PROTO((char *filename, char *extname, int *Nvec)); /* buffer handling */ @@ -140,4 +150,5 @@ int IsBufferPtr PROTO((Buffer *buf)); int PrintBuffers PROTO((int Long)); +int ListBuffersToList PROTO((char *name)); int CreateBuffer PROTO((Buffer *buf, int Nx, int Ny, int bitpix, float bzero, float bscale)); Buffer *SelectBuffer PROTO((char *name, int mode, int verbose)); Index: /trunk/Ohana/src/opihi/include/dvoshell.h =================================================================== --- /trunk/Ohana/src/opihi/include/dvoshell.h (revision 33661) +++ /trunk/Ohana/src/opihi/include/dvoshell.h (revision 33662) @@ -6,269 +6,4 @@ # ifndef DVOSHELL_H # define DVOSHELL_H - -#ifdef NOT_MOVED_TO_LIBDVO -typedef enum { - DB_STACK_NONE = 0, - DB_STACK_INT = 0x01, - DB_STACK_FIELD = 0x02, - DB_STACK_TEMP = 0x04, - DB_STACK_VALUE = 0x08, - DB_STACK_CLOSE_PAR = 0x10, - DB_STACK_OPEN_PAR, - DB_STACK_LOGIC, - DB_STACK_COMPARE, - DB_STACK_BITWISE, - DB_STACK_SUM, - DB_STACK_MULTIPLY, - DB_STACK_POWER, - DB_STACK_UNARY, -} dbStackTypes; - -/* magnitude types */ -enum {MAG_NONE, - MAG_INST, - MAG_CAT, - MAG_APER, - MAG_SYS, - MAG_REL, - MAG_CAL, - MAG_AVE, - MAG_REF, - MAG_ERR, - MAG_AVE_ERR, - MAG_PHOT_FLAGS, - MAG_CHISQ, - MAG_NCODE, - MAG_NPHOT, -}; - -/* measure fields */ -enum {MEAS_ZERO, - MEAS_GLON, - MEAS_GLAT, - MEAS_GLON_AVE, - MEAS_GLAT_AVE, - MEAS_ELON, - MEAS_ELAT, - MEAS_ELON_AVE, - MEAS_ELAT_AVE, - MEAS_RA, - MEAS_DEC, - MEAS_RA_AVE, - MEAS_DEC_AVE, - MEAS_RA_AVE_ERR, - MEAS_DEC_AVE_ERR, - MEAS_U_RA, - MEAS_U_DEC, - MEAS_U_RA_ERR, - MEAS_U_DEC_ERR, - MEAS_PAR, - MEAS_PAR_ERR, - MEAS_RA_OFFSET, - MEAS_DEC_OFFSET, - MEAS_RA_FIT_OFFSET, - MEAS_DEC_FIT_OFFSET, - MEAS_RA_OFFSET_ERR, - MEAS_DEC_OFFSET_ERR, - MEAS_CHISQ_POS, - MEAS_CHISQ_PM, - MEAS_CHISQ_PAR, - MEAS_TMEAN, - MEAS_TRANGE, - MEAS_NMEAS, - MEAS_NMISS, - MEAS_NPOS, - MEAS_OBJ_FLAGS, - MEAS_MAG, - MEAS_MINST, - MEAS_MCAT, - MEAS_MSYS, - MEAS_MREL, - MEAS_MCAL, - MEAS_EXPTIME, - MEAS_AIRMASS, - MEAS_ALT, - MEAS_AZ, - MEAS_PHOTCODE, - MEAS_PHOTCODE_EQUIV, - MEAS_TIME, - MEAS_FWHM, - MEAS_FWHM_MAJ, - MEAS_FWHM_MIN, - MEAS_THETA, - MEAS_POSANGLE, - MEAS_PLATESCALE, - MEAS_MXX, - MEAS_MXY, - MEAS_MYY, - MEAS_DOPHOT, - MEAS_DB_FLAGS, - MEAS_PHOT_FLAGS, - MEAS_XCCD, - MEAS_YCCD, - MEAS_XCCD_ERR, - MEAS_YCCD_ERR, - MEAS_POS_SYS_ERR, - MEAS_XMOSAIC, - MEAS_YMOSAIC, - MEAS_SKY, - MEAS_dSKY, - MEAS_DET_ID, - MEAS_OBJ_ID, - MEAS_CAT_ID, - MEAS_IMAGE_ID, - MEAS_PSF_QF, - MEAS_PSF_QF_PERFECT, - MEAS_PSF_CHISQ, - MEAS_PSF_NDOF, - MEAS_PSF_NPIX, - MEAS_CR_NSIGMA, - MEAS_EXT_NSIGMA, - MEAS_EXTERN_ID, -}; - -/* average fields */ -enum {AVE_ZERO, - AVE_RA, - AVE_DEC, - AVE_RA_ERR, - AVE_DEC_ERR, - AVE_GLON, - AVE_GLAT, - AVE_ELON, - AVE_ELAT, - AVE_U_RA, - AVE_U_DEC, - AVE_U_RA_ERR, - AVE_U_DEC_ERR, - AVE_PAR, - AVE_PAR_ERR, - AVE_CHISQ_POS, - AVE_CHISQ_PM, - AVE_CHISQ_PAR, - AVE_TMEAN, - AVE_TRANGE, - AVE_Xp, - AVE_NMEAS, - AVE_NMISS, - AVE_NPOS, - AVE_NPHOT, - AVE_NCODE, - AVE_MAG, - AVE_dMAG, - AVE_Xm, - AVE_OBJ_FLAGS, - AVE_TYPE, - AVE_TYPEFRAC, - AVE_OBJID, - AVE_CATID, - AVE_EXTID_HI, - AVE_EXTID_LO, -}; - -enum {IMAGE_ZERO, - IMAGE_RA, - IMAGE_DEC, - IMAGE_GLON, - IMAGE_GLAT, - IMAGE_ELON, - IMAGE_ELAT, - IMAGE_XM, - IMAGE_AIRMASS, - IMAGE_MCAL, - IMAGE_dMCAL, - IMAGE_PHOTCODE, - IMAGE_TIME, - IMAGE_FWHM, - IMAGE_FWHM_MEDIAN, - IMAGE_EXPTIME, - IMAGE_NSTAR, - IMAGE_NCAL, - IMAGE_SKY, - IMAGE_FLAGS, - IMAGE_CCDNUM, - IMAGE_NX_PIX, - IMAGE_NY_PIX, - IMAGE_THETA, - IMAGE_SKEW, - IMAGE_SCALE, - IMAGE_DSCALE, - IMAGE_APRESID, - IMAGE_DAPRESID, - IMAGE_SIDTIME, - IMAGE_LATITUDE, - IMAGE_DET_LIMIT, - IMAGE_SAT_LIMIT, - IMAGE_CERROR, - IMAGE_FWHM_MAJ, - IMAGE_FWHM_MIN, - IMAGE_FWHM_MAJ_MEDIAN, - IMAGE_FWHM_MIN_MEDIAN, - IMAGE_TRATE, - IMAGE_IMAGE_ID, - IMAGE_EXTERN_ID, - IMAGE_SOURCE_ID, - IMAGE_X_LL_CHIP, - IMAGE_X_LR_CHIP, - IMAGE_X_UL_CHIP, - IMAGE_X_UR_CHIP, - IMAGE_Y_LL_CHIP, - IMAGE_Y_LR_CHIP, - IMAGE_Y_UL_CHIP, - IMAGE_Y_UR_CHIP, - IMAGE_X_LL_FP, - IMAGE_X_LR_FP, - IMAGE_X_UL_FP, - IMAGE_X_UR_FP, - IMAGE_Y_LL_FP, - IMAGE_Y_LR_FP, - IMAGE_Y_UL_FP, - IMAGE_Y_UR_FP, - IMAGE_X_ERR_SYS, - IMAGE_Y_ERR_SYS, - IMAGE_MAG_ERR_SYS, - IMAGE_NFIT_PHOTOM, - IMAGE_NFIT_ASTROM, - IMAGE_NLINK_PHOTOM, - IMAGE_NLINK_ASTROM -}; - -enum {DVO_TABLE_AVERAGE, DVO_TABLE_MEASURE, DVO_TABLE_IMAGE}; -enum {DVO_DB_CMDLINE_ERROR, DVO_DB_CMDLINE_IS_END, DVO_DB_CMDLINE_IS_WHERE, DVO_DB_CMDLINE_IS_MATCH}; - -// options for selecting the ra,dec limits of the db selections -typedef struct { - char *name; - char *list; - int useDisplay; - int useSkyregion; -} SkyRegionSelection; - -// a single db field -typedef struct { - char *name; - int extract; - int table; - int ID; - int magMode; - char type; - PhotCode *photcode; -} dbField; - -// db boolean operations -typedef struct { - char *name; - char type; - int field; - opihi_flt FltValue; - opihi_int IntValue; -} dbStack; - -typedef struct { - opihi_flt Flt; - opihi_int Int; -} dbValue; - -#endif // notdef MOVED_TO_LIBDVO typedef struct { @@ -291,4 +26,9 @@ e_time t; } CMPstars; + +/** some globals used particularly by DVO_CLIENT **/ +int HOST_ID; +char *HOSTDIR; +char *RESULT_FILE; /*** dvo prototypes ***/ @@ -316,13 +56,9 @@ void FreeDVO PROTO((void)); int InitPhotcodes PROTO((void)); -#ifdef NOT_MOVED_TO_DVO -Image *LoadImages PROTO((off_t *Nimage)); -void FreeImages PROTO((Image *images)); -Image *MatchImage PROTO((unsigned int time, short int source, unsigned int imageID)); -#endif Coords *MatchMosaic PROTO((unsigned int time, short int source)); int Quality PROTO((Measure *measure, int IsDophot)); int SelectMags PROTO((int Nphot, int Tphot, int Ns, Average *average, Measure *measure, SecFilt *secfilt, int UL)); +int SetSkyRegions PROTO((SkyRegionSelection *selection)); SkyList *SelectRegions PROTO((SkyRegionSelection *selection)); SkyList *SkyListLoadFile PROTO((char *filename)); @@ -352,21 +88,5 @@ CMPstars *cmpReadText PROTO((FILE *f, off_t *nstars)); int RD_to_XYpic PROTO((double *x, double *y, double r, double d, Coords *coords, double Rmin, double Rmax, double Rmid, int *leftside)); -int wordhash (char *word); - -#ifdef NOT_MOVED_TO_LIBDVO -// dvo DB field functions -dbField *dbCmdlineFields PROTO((int argc, char **argv, int table, int *last, int *nfields)); -int dbCmdlineConditions PROTO((int argc, char **argv, int first, int *nextField)); -dbStack *dbRPN PROTO((int argc, char **argv, int *nstack)); -int dbCheckStack PROTO((dbStack *stack, int Nstack, int table, dbField **inFields, int *Nfields)); -int dbBooleanCond PROTO((dbStack *inStack, int NinStack, dbValue *fields)); -void dbInitStack PROTO((dbStack *stack)); -void dbFreeStack PROTO((dbStack *stack, int Nstack)); -void dbFreeEntry PROTO((dbStack *stack)); -void dbFreeTempEntry PROTO((dbStack *stack)); - -dbStack *dbBinary PROTO((dbStack *V1, dbStack *V2, char *op, dbValue *fields)); -dbStack *dbUnary PROTO((dbStack *V1, char *op, dbValue *fields)); -#endif +int wordhash PROTO((char *word)); int GetMagMode PROTO((char *string)); @@ -378,31 +98,6 @@ dbValue dbExtractMeasures PROTO((Average *average, SecFilt *secfilt, Measure *measure, dbField *field)); dbValue dbExtractImages PROTO((Image *image, off_t Nimage, off_t N, dbField *field)); -#ifdef NOT_MOVED_TO_LIBDVO -dbValue dbExtractAverages PROTO((Average *average, SecFilt *secfilt, Measure *measure, dbField *field)); -void dbInitField PROTO((dbField *field)); -void dbFreeFields PROTO((dbField *fields, int Nfields)); -int dbAstroRegionLimits PROTO((dbStack **stack, int *nstack, SkyRegionSelection *selection, int table)); -int get_skyregion (double *Rs, double *Re, double *Ds, double *De); -int set_skyregion (double Rs, double Re, double Ds, double De); -void FreeImageSelection (void); -char *strfloat PROTO((float value)); - - -void FreeSkyRegionSelection (SkyRegionSelection *selection); - - -int dbExtractMeasuresInitTransform (CoordTransformSystem target); -int dbExtractMeasuresInitAve (void); -int dbExtractMeasuresInitMeas (void); -int dbExtractMeasuresInit (void); - -int dbExtractAveragesInitTransform (CoordTransformSystem target); -int dbExtractAveragesInit (void); - -int dbExtractImagesInitTransform (CoordTransformSystem target); -int dbExtractImagesInit (void); -int dbExtractImagesReset (void); -#endif // NOT_MOVED_TO_LIBDVO - +int HostTableLaunchJobs PROTO((HostTable *table, char *basecmd, char *options)); +int HostTableParallelOps PROTO((int argc, char **argv, char *ResultFile, int Nelements, int VERBOSE)); # endif Index: /trunk/Ohana/src/opihi/include/pcontrol.h =================================================================== --- /trunk/Ohana/src/opihi/include/pcontrol.h (revision 33661) +++ /trunk/Ohana/src/opihi/include/pcontrol.h (revision 33662) @@ -252,5 +252,4 @@ int CheckDoneJob (Job *job, Host *host); int GetJobOutput (char *command, Host *host, JobOutput *output); -int rconnect (char *command, char *hostname, char *shell, int *stdio); int PclientCommand (Host *host, char *command, char *response, HostResp response_state); Index: /trunk/Ohana/src/opihi/include/shell.h =================================================================== --- /trunk/Ohana/src/opihi/include/shell.h (revision 33661) +++ /trunk/Ohana/src/opihi/include/shell.h (revision 33662) @@ -148,4 +148,5 @@ char *opihi_append PROTO((char *output, int *Noutput, char *start, char *stop)); void interpolate_slash PROTO((char *line)); +char *paste_args PROTO((int argc, char **argv)); /* macro functions (mapped to commands) */ Index: /trunk/Ohana/src/opihi/lib.shell/BufferOps.c =================================================================== --- /trunk/Ohana/src/opihi/lib.shell/BufferOps.c (revision 33661) +++ /trunk/Ohana/src/opihi/lib.shell/BufferOps.c (revision 33662) @@ -256,2 +256,16 @@ return (TRUE); } + +int ListBuffersToList (char *name) { + + int i; + char line[1024]; + + for (i = 0; i < Nbuffers; i++) { + sprintf (line, "%s:%d", name, i); + set_str_variable (line, buffers[i][0].name); + } + sprintf (line, "%s:n", name); + set_int_variable (line, Nbuffers); + return (Nbuffers); +} Index: /trunk/Ohana/src/opihi/lib.shell/ListOps.c =================================================================== --- /trunk/Ohana/src/opihi/lib.shell/ListOps.c (revision 33661) +++ /trunk/Ohana/src/opihi/lib.shell/ListOps.c (revision 33662) @@ -193,4 +193,6 @@ if (!strcmp (temp, "-add")) goto escape; if (!strcmp (temp, "-del")) goto escape; + if (!strcmp (temp, "-vectors")) goto escape; + if (!strcmp (temp, "-buffers")) goto escape; } Index: /trunk/Ohana/src/opihi/lib.shell/Makefile =================================================================== --- /trunk/Ohana/src/opihi/lib.shell/Makefile (revision 33661) +++ /trunk/Ohana/src/opihi/lib.shell/Makefile (revision 33662) @@ -24,4 +24,5 @@ $(SDIR)/VariableOps.$(ARCH).o \ $(SDIR)/VectorOps.$(ARCH).o \ +$(SDIR)/VectorIO.$(ARCH).o \ $(SDIR)/check_stack.$(ARCH).o \ $(SDIR)/command.$(ARCH).o \ Index: /trunk/Ohana/src/opihi/lib.shell/VectorIO.c =================================================================== --- /trunk/Ohana/src/opihi/lib.shell/VectorIO.c (revision 33662) +++ /trunk/Ohana/src/opihi/lib.shell/VectorIO.c (revision 33662) @@ -0,0 +1,233 @@ +# include "opihi.h" + +// write a set of vectors to a FITS file (vectors names become fits column names) +int WriteVectorTableFITS (char *filename, char *extname, Vector **vec, int Nvec, int append, char *format) { + + char *tformat = NULL; + Header header; + Matrix matrix; + Header theader; + FTable ftable; + + int j; + FILE *f = NULL; + + /* open file for outuput */ + if (append) { + f = fopen (filename, "a"); + } else { + f = fopen (filename, "w"); + } + if (f == (FILE *) NULL) { + gprint (GP_ERR, "can't open file for write\n"); + return (FALSE); + } + + if (!append) { + gfits_init_header (&header); + header.extend = TRUE; + gfits_create_header (&header); + gfits_create_matrix (&header, &matrix); + } + + gfits_create_table_header (&theader, "BINTABLE", extname); + + ALLOCATE (tformat, char, 2*Nvec); + if (format) { + // the bintable format string can only define the byte-width of each field. valid output columns are currently: + // B (char), I (16 bit short), J (32 bit int), E (32 bit float), D (64 bit double). + // the format string is just the sequence of types, eg: LIIJEED + // validate the format string + char *ptr = format; + for (j = 0; j < Nvec; j++) { + while (*ptr && OHANA_WHITESPACE (*ptr)) ptr++; + if (*ptr == 0) { + gprint (GP_ERR, "error in binary table format %s (insufficient format chars)\n", format); + goto escape; + } + if ((*ptr != 'B') && (*ptr != 'I') && (*ptr != 'J') && (*ptr != 'D') && (*ptr != 'E')) { + gprint (GP_ERR, "error in binary table format %s: invalid character %c\n", format, *ptr); + goto escape; + } + tformat[2*j + 0] = *ptr; + tformat[2*j + 1] = 0; // a bit sleazy : use a 2xN string to store N 1-byte strings + ptr ++; + } + while (*ptr && OHANA_WHITESPACE (*ptr)) ptr++; + if (*ptr) { + gprint (GP_ERR, "error in binary table format %s (extra characters in format)\n", format); + goto escape; + } + } else { + for (j = 0; j < Nvec; j++) { + // if the format is not defined, just use the native byte-widths + tformat[2*j + 0] = (vec[j][0].type == OPIHI_FLT) ? 'D' : 'J'; + tformat[2*j + 1] = 0; + } + } + + // define the columns of the table. XXX NOTE: we cannot have duplicate names in + // output table (because the data goes to the named column below). need to enforce + // this somehow + for (j = 0; j < Nvec; j++) { + gfits_define_bintable_column (&theader, &tformat[2*j], vec[j][0].name, NULL, NULL, 1.0, 0.0); + } + free (tformat); + + // generate the output array that carries the data + gfits_create_table (&theader, &ftable); + + // add the vectors to the output array + for (j = 0; j < Nvec; j++) { + if (vec[j][0].type == OPIHI_FLT) { + gfits_set_bintable_column_reformat (&theader, &ftable, vec[j][0].name, "double", vec[j][0].elements.Flt, vec[j][0].Nelements); + } else { + gfits_set_bintable_column_reformat (&theader, &ftable, vec[j][0].name, "int", vec[j][0].elements.Int, vec[j][0].Nelements); + } + } + + if (!append) { + gfits_fwrite_header (f, &header); + gfits_fwrite_matrix (f, &matrix); + gfits_free_header (&header); + gfits_free_matrix (&matrix); + } + gfits_fwrite_Theader (f, &theader); + gfits_fwrite_table (f, &ftable); + + gfits_free_header (&theader); + gfits_free_table (&ftable); + + fclose (f); + fflush (f); + return (TRUE); + + escape: + if (!append) { + gfits_free_header (&header); + gfits_free_matrix (&matrix); + } + gfits_free_header (&theader); + gfits_free_table (&ftable); + + if (tformat) free (tformat); + fclose (f); + fflush (f); + return (FALSE); +} + +// read the complete set of vectors from the given FITS file & extension +// XXX not quite right : I need to merge multiple vectors together: +// do not associate with an Opihi vector until later in mextract +Vector **ReadVectorTableFITS (char *filename, char *extname, int *nvec) { + + Header header; + FTable ftable; + + int i, j, k; + FILE *f = NULL; + + /* open file for input */ + f = fopen (filename, "r"); + if (f == (FILE *) NULL) { + gprint (GP_ERR, "can't open file for read\n"); + return NULL; + } + + ftable.header = &header; + + // read the full table data into a buffer + if (!gfits_fread_ftable (f, &ftable, extname)) { + fclose (f); + return NULL; + } + + // XXX handle binary and ascii tables (see read_vectors.c) + // find the columns in the table + int Nfields; + gfits_scan (&header, "TFIELDS", "%d", 1, &Nfields); + int Nrows = header.Naxis[1]; + + // how many output vectors do we need? depends on number of columns per field, but min is Nfields + int Nvec = 0; + int NVEC = Nfields; + Vector **vec = NULL; + ALLOCATE (vec, Vector *, NVEC); + + for (i = 0; i < Nfields; i++) { + int Nval; + char type[16], label[16], name[80]; + + // determine the column name, data type, and number of sub-fields + sprintf (label, "TTYPE%d", i + 1); + int status = gfits_scan (&header, label, "%s", 1, name); + assert (status); + + status = gfits_get_bintable_column_type_by_N (&header, i + 1, type, &Nval); + assert (status); + + if (Nvec + Nval >= NVEC) { + NVEC = Nvec + Nval + 16; + REALLOCATE (vec, Vector *, NVEC); + } + + int vecType = OPIHI_INT; + if (!strcmp (type, "double") || !strcmp (type, "float")) { + vecType = OPIHI_FLT; + } + + // generate the needed vectors + for (j = 0; j < Nval; j++) { + vec[Nvec + j] = InitVector(); + if (Nval == 1) { + strcpy (vec[Nvec + j]->name, name); + } else { + snprintf (vec[Nvec + j]->name, OPIHI_NAME_SIZE, "%s:%d", name, j); + } + ResetVector (vec[Nvec + j], vecType, Nrows); + } + + // read the actual table data into a column + void *data; + status = gfits_get_bintable_column (&header, &ftable, name, &data); + assert (status); + +# define ASSIGN_DATA(TYPE,OPTYPE) \ + /* assign the data to the actual vector */ \ + if (!strcmp (type, #TYPE)) { \ + TYPE *Ptr = data; \ + for (k = 0; k < Nrows; k++) { \ + for (j = 0; j < Nval; j++, Ptr++) { \ + vec[Nvec + j][0].elements.OPTYPE[k] = *Ptr; \ + } } } + + // assign the data to the actual vector + ASSIGN_DATA(char, Int); + ASSIGN_DATA(short, Int); + ASSIGN_DATA(int, Int); + ASSIGN_DATA(int64_t, Int); + ASSIGN_DATA(float, Flt); + ASSIGN_DATA(double, Flt); + + free (data); + Nvec += Nval; + } + + gfits_free_header (&header); + gfits_free_table (&ftable); + + fclose (f); + + *nvec = Nvec; + return vec; + +//escape: +// gfits_free_header (&header); +// gfits_free_matrix (&matrix); +// gfits_free_header (&theader); +// gfits_free_table (&ftable); +// +// fclose (f); +// fflush (f); +// return (FALSE); +} Index: /trunk/Ohana/src/opihi/lib.shell/VectorOps.c =================================================================== --- /trunk/Ohana/src/opihi/lib.shell/VectorOps.c (revision 33661) +++ /trunk/Ohana/src/opihi/lib.shell/VectorOps.c (revision 33662) @@ -15,15 +15,22 @@ // this function is NOT thread protected : it is only used in startup and/or shutdown +void FreeVectorArray (Vector **vec, int Nvec) { + + int i; + + if (!vec) return; + for (i = 0; i < Nvec; i++) { + if (!vec[i]) continue; + if (vec[i]->elements.Int) { + free (vec[i]->elements.Int); + } + free (vec[i]); + } + free (vec); +} + +// this function is NOT thread protected : it is only used in startup and/or shutdown void FreeVectors () { - - int i; - - for (i = 0; i < Nvectors; i++) { - if (vectors[i][0].elements.Int) { - free (vectors[i][0].elements.Int); - } - free (vectors[i]); - } - free (vectors); + FreeVectorArray (vectors, Nvectors); } @@ -93,4 +100,38 @@ } +// Assign the given Vector to the internal array of vectors by name +int AssignVector (Vector *vec, char *name, int mode, int verbose) { + + int i; + + if (name == NULL) goto error; + if (ISNUM(name[0])) goto error; + if (IsBuffer(name)) goto error; + + for (i = 0; (i < Nvectors) && (strcmp(vectors[i][0].name, name)); i++); + /* is a new vector */ + if (i == Nvectors) { + if (mode == OLDVECTOR) goto error; + pthread_mutex_lock (&mutex); + Nvectors ++; + REALLOCATE (vectors, Vector *, Nvectors); + vectors[i] = vec; + pthread_mutex_unlock (&mutex); + return TRUE; + } + /* is an old vector */ + if (mode == NEWVECTOR) goto error; + if (vectors[i]) { + if (vectors[i][0].elements.Ptr) free (vectors[i][0].elements.Ptr); + free (vectors[i]); + } + vectors[i] = vec; + return TRUE; + + error: + if (verbose) gprint (GP_ERR, "invalid vector %s\n", name); + return FALSE; +} + /* delete by pointer */ int DeleteVector (Vector *vec) { @@ -148,9 +189,9 @@ if (in[0].elements.Ptr) { if (in[0].type == OPIHI_FLT) { - ALLOCATE (out[0].elements.Flt, opihi_flt, out[0].Nelements); + ALLOCATE (out[0].elements.Flt, opihi_flt, MAX(1,out[0].Nelements)); memcpy (out[0].elements.Flt, in[0].elements.Flt, out[0].Nelements*sizeof(opihi_flt)); out[0].type = OPIHI_FLT; } else { - ALLOCATE (out[0].elements.Int, opihi_int, out[0].Nelements); + ALLOCATE (out[0].elements.Int, opihi_int, MAX(1,out[0].Nelements)); memcpy (out[0].elements.Int, in[0].elements.Int, out[0].Nelements*sizeof(opihi_int)); out[0].type = OPIHI_INT; @@ -164,8 +205,8 @@ out[0].Nelements = in[0].Nelements; if (type == OPIHI_FLT) { - ALLOCATE (out[0].elements.Flt, opihi_flt, out[0].Nelements); + ALLOCATE (out[0].elements.Flt, opihi_flt, MAX(1,out[0].Nelements)); out[0].type = OPIHI_FLT; } else { - ALLOCATE (out[0].elements.Int, opihi_int, out[0].Nelements); + ALLOCATE (out[0].elements.Int, opihi_int, MAX(1,out[0].Nelements)); out[0].type = OPIHI_INT; } @@ -173,8 +214,8 @@ } -// ResetVector (vecx, OPIHI_FLT, MAX (Npts, 1)); int ResetVector (Vector *vec, char type, int Nelements) { - vec[0].Nelements = Nelements; + // a vector can only have >= 0 elements + vec[0].Nelements = MAX(Nelements,0); if (type == OPIHI_FLT) { REALLOCATE (vec[0].elements.Flt, opihi_flt, MAX(1, Nelements)); @@ -190,10 +231,10 @@ int SetVector (Vector *vec, char type, int Nelements) { - vec[0].Nelements = Nelements; + vec[0].Nelements = MAX(Nelements,0); if (type == OPIHI_FLT) { - ALLOCATE (vec[0].elements.Flt, opihi_flt, Nelements); + ALLOCATE (vec[0].elements.Flt, opihi_flt, MAX(1,Nelements)); vec[0].type = OPIHI_FLT; } else { - ALLOCATE (vec[0].elements.Int, opihi_int, Nelements); + ALLOCATE (vec[0].elements.Int, opihi_int, MAX(1,Nelements)); vec[0].type = OPIHI_INT; } @@ -287,3 +328,108 @@ return (TRUE); } + +int ListVectorsToList (char *name) { + + int i; + char line[1024]; + + for (i = 0; i < Nvectors; i++) { + sprintf (line, "%s:%d", name, i); + set_str_variable (line, vectors[i][0].name); + } + sprintf (line, "%s:n", name); + set_int_variable (line, Nvectors); + return (Nvectors); +} + +// Take two arrays of vectors and merge equal named vectors. +// Output is a single array in (vec), with vectors lengths of len(vec) + len(invec) +// For ease, require that the order of the names match & number of vectors match +Vector **MergeVectors (Vector **vec, int *Nvec, Vector **invec, int Ninvec) { + + int i, j; + + if (vec == NULL) { + *Nvec = Ninvec; + return invec; + } + + myAssert (*Nvec == Ninvec, "programming error (1) %d vs %d", *Nvec, Ninvec); + + for (i = 0; i < Ninvec; i++) { + myAssert (!strcmp(vec[i]->name, invec[i]->name), "programming error (2) %s vs %s", vec[i]->name, invec[i]->name); + myAssert (vec[i]->type == invec[i]->type, "programming error (3), %d vs %d", vec[i]->type, invec[i]->type); + + int N = vec[i]->Nelements; + if (vec[i]->type == OPIHI_FLT) { + REALLOCATE (vec[i]->elements.Flt, opihi_flt, vec[i]->Nelements + invec[i]->Nelements); + for (j = 0; j < invec[i]->Nelements; j++) { + vec[i]->elements.Flt[N+j] = invec[i]->elements.Flt[j]; + } + } else { + REALLOCATE (vec[i]->elements.Int, opihi_int, vec[i]->Nelements + invec[i]->Nelements); + for (j = 0; j < invec[i]->Nelements; j++) { + vec[i]->elements.Int[N+j] = invec[i]->elements.Int[j]; + } + } + vec[i]->Nelements += invec[i]->Nelements; + } + return vec; +} + +// Take two arrays of vectors and merge equal named vectors, where the last vector +// specifies the element in the output vector. All input vectors must have a max sequence +// value of Nelements. Output is a single array in (vec), with vector lengths of +// Nelements for ease, require that the order of the names match & number of vectors match +Vector **MergeVectorsByIndex (Vector **vec, int *Nvec, Vector **invec, int Ninvec, int Nelements) { + + int i, j; + + // on first call, allocate a new vector, excluding the index + int newArray = FALSE; + if (vec == NULL) { + ALLOCATE (vec, Vector *, Ninvec - 1); + *Nvec = Ninvec - 1; + newArray = TRUE; + } + + myAssert (*Nvec == Ninvec - 1, "programming error (1)"); + + // find the index vector + int idx = Ninvec - 1; + myAssert (!strcmp(invec[idx]->name, "index"), "failed to find index vector"); + for (i = 0; i < Ninvec - 1; i++) { + // on first call, create the output vector + if (newArray) { + vec[i] = InitVector(); + strcpy (vec[i]->name, invec[i]->name); + ResetVector (vec[i], invec[i]->type, Nelements); + for (j = 0; j < Nelements; j++) { + if (vec[i][0].type == OPIHI_FLT) { + vec[i][0].elements.Flt[j] = NAN; + } else { + vec[i][0].elements.Int[j] = 0; // or NAN_INT? + } + } + } + + myAssert (!strcmp(vec[i]->name, invec[i]->name), "programming error (2)"); + myAssert (vec[i]->type == invec[i]->type, "programming error (2)"); + + // copy vector elements from input to output, matching location + if (vec[i]->type == OPIHI_FLT) { + for (j = 0; j < invec[i]->Nelements; j++) { + int seq = invec[idx]->elements.Int[j]; + vec[i]->elements.Flt[seq] = invec[i]->elements.Flt[j]; + } + } else { + for (j = 0; j < invec[i]->Nelements; j++) { + int seq = invec[idx]->elements.Int[j]; + vec[i]->elements.Int[seq] = invec[i]->elements.Int[j]; + } + } + } + return vec; +} + Index: /trunk/Ohana/src/opihi/lib.shell/parse.c =================================================================== --- /trunk/Ohana/src/opihi/lib.shell/parse.c (revision 33661) +++ /trunk/Ohana/src/opihi/lib.shell/parse.c (revision 33662) @@ -101,5 +101,9 @@ V1 ++; while (isspace (*V1)) V1++; - if (*V1 == 0) goto error; + if (*V1 == 0) { + // assign empty vector + set_str_variable (V0, ""); + goto escape; + } /* command replacement. execute the line, place answer in val. */ Index: /trunk/Ohana/src/opihi/lib.shell/parse.cleanup.c =================================================================== --- /trunk/Ohana/src/opihi/lib.shell/parse.cleanup.c (revision 33662) +++ /trunk/Ohana/src/opihi/lib.shell/parse.cleanup.c (revision 33662) @@ -0,0 +1,428 @@ +# include "opihi.h" + +enum {OP_SET_SCALAR, OP_INCR_SCALAR, OP_DECR_SCALAR, OP_ADD_SCALAR, OP_SUB_SCALAR, OP_SET_VECTOR, OP_INCR_VECTOR, OP_DECR_VECTOR, OP_ADD_VECTOR, OP_SUB_VECTOR}; + +char *parse (int *status, char *line) { + + double fval; + float *fptr; + char *newline, *N, *L, *val, *B, *V, *V0, *V1, *c1, *c2, *p, *px, *py; + int Nx, Ny, Nbytes, filestatus, size, NLINE, isBuffer, inRange; + FILE *f; + Vector *vec; + Buffer *buf; + + *status = TRUE; + + Ny = 0; + val = V0 = NULL; + + if (line == (char *) NULL) goto error; + + // what are we trying to do here? + /* case 1: $var = expression */ + if (line[0] == '$') { + /* find the target variable name and the rest of the line */ + interpolate_slash (line); + V0 = thisvar (line); + if (V0 == NULL) goto error; + V1 = aftervar (line); + if (V1 == NULL) goto error; + + if (!strcmp (V1, "++")) { mode = OP_INCR_SCALAR; } + if (!strcmp (V1, "--")) { mode = OP_DECR_SCALAR; } + if (!strncmp (V1, "+=", 2)) { mode = OP_ADD_SCALAR; } + if (!strncmp (V1, "-=", 2)) { mode = OP_SUB_SCALAR; } + if (!strcmp (V1, "=")) { mode = OP_SET_SCALAR; } + } + + /* case 2: vect[N] = value or buff[N][M] = value */ + V0 = thiscomm (line); + px = py = NULL; + if (strchr(V0, '[') != (char *) NULL) { + isBuffer = FALSE; + + N = strchr (V0, '['); + if (N == V0) goto error; + + // find end of 1st bracket + L = strchr (V0, ']'); + if (!L) goto error; + px = N + 1; + + // is there a second bracket? + N = L + 1; + if (*N == '[') { + isBuffer = TRUE; + + // find end of 2nd bracket + L = strchr (N, ']'); + if (!L) goto error; + py = N + 1; + } + + V1 = nextcomm (line); + + if (!strcmp (V1, "++")) { mode = OP_INCR_VECTOR; } + if (!strcmp (V1, "--")) { mode = OP_DECR_VECTOR; } + if (!strncmp (V1, "+=", 2)) { mode = OP_ADD_VECTOR; } + if (!strncmp (V1, "-=", 2)) { mode = OP_SUB_VECTOR; } + if (!strcmp (V1, "=")) { mode = OP_SET_VECTOR; } + } + free (V0); V0 = (char *) NULL; + + + + /* case 1: $var = expression */ + if (line[0] == '$') { + /* find the target variable name and the rest of the line */ + interpolate_slash (line); + V0 = thisvar (line); + if (V0 == NULL) goto error; + V1 = aftervar (line); + if (V1 == NULL) goto error; + + /* increment operator */ + if (!strcmp (V1, "++")) { + val = get_variable (V0); + if (val == NULL) { + fval = 1.0; + } else { + fval = atof (val) + 1.0; + } + set_variable (V0, fval); + goto escape; + } + + /* decrement operator */ + if (!strcmp (V1, "--")) { + val = get_variable (V0); + if (val == NULL) { + fval = -1.0; + } else { + fval = atof (val) - 1.0; + } + set_variable (V0, fval); + goto escape; + } + + if (!strncmp (V1, "+=", 2)) { + V1 ++; + if (*V1 == 0) goto error; + V1 ++; + if (*V1 == 0) goto error; + + val = get_variable (V0); + if (val == NULL) { + fval = 0.0; + } else { + fval = atof (val); + } + + /* dvomath returns a new string, or NULL, with the result of the expression */ + val = dvomath (1, &V1, &size, 0); + if (val == NULL) goto error; + fval += atof(val); + // save the result + set_variable (V0, fval); + goto escape; + } + + if (!strncmp (V1, "-=", 2)) { + V1 ++; + if (*V1 == 0) goto error; + V1 ++; + if (*V1 == 0) goto error; + + val = get_variable (V0); + if (val == NULL) { + fval = 0.0; + } else { + fval = atof (val); + } + + /* dvomath returns a new string, or NULL, with the result of the expression */ + val = dvomath (1, &V1, &size, 0); + if (val == NULL) goto error; + fval -= atof(val); + // save the result + set_variable (V0, fval); + goto escape; + } + + /* not an assignement, syntax error */ + if (*V1 != '=') goto error; + + /* find first non-WHITESPACE character after = */ + V1 ++; + while (isspace (*V1)) V1++; + if (*V1 == 0) { + // assign empty vector + set_str_variable (V0, ""); + goto escape; + } + + /* command replacement. execute the line, place answer in val. */ + if (*V1 == '`') { + + /* look for end of line, must be ` */ + B = V1 + strlen (V1) - 1; + if (*B != '`') goto error; + + /* val will hold the result */ + // XXX this is limiting!!! + ALLOCATE (val, char, 1024); + + /* B is the command to be executed */ + B = strncreate (V1 + 1, strlen(V1) - 2); + f = popen (B, "r"); + Nbytes = fread (val, 1, 1023, f); + val[Nbytes] = 0; + filestatus = pclose (f); + free (B); + + if (filestatus) gprint (GP_ERR, "warning: exit status of command %d\n", filestatus); + + /* convert all but last return to ' '. drop last return */ + for (B = val; *B != 0; B++) { + if (*B == '\n') { + if (B[1]) { + *B = ' '; + } else { + *B = 0; + } + } + } + // save the resulting value + set_str_variable (V0, val); + goto escape; /* frees temp variables */ + } + + /* simple variable assignment */ + /* dvomath returns a new string, or NULL, with the result of the expression */ + val = dvomath (1, &V1, &size, 0); + if (val == NULL) { + while (OHANA_WHITESPACE (*V1)) V1++; + val = strcreate (V1); + } + // save the result + set_str_variable (V0, val); + goto escape; /* frees temp variables */ + } + + /* case 2: vect[N] = value or buff[N][M] = value */ + V0 = thiscomm (line); + if (strchr(V0, '[') != (char *) NULL) { + /* get vector name (left in V0) */ + N = strchr (V0, '['); + if (N == V0) goto error; + + /* get bracket contents (left in V) */ + L = strchr (V0, ']'); + if (L == (char *) NULL) goto error; + + *N = 0; + V = strncreate (N+1, L - N - 1); + + /* expand value in brackets */ + val = dvomath (1, &V, &size, 0); + if (val == NULL) { + print_error (); + goto error; + } + Nx = atoi (val); + free (val); val = NULL; + free (V); + + isBuffer = FALSE; + N = L + 1; + if (*N == '[') { + isBuffer = TRUE; + L = strchr (N, ']'); + + V = strncreate (N+1, L - N - 1); + + /* expand value in brackets */ + val = dvomath (1, &V, &size, 0); + if (val == NULL) { + print_error (); + goto error; + } + Ny = atoi (val); + free (val); val = NULL; + free (V); + } + + /* find value for assignment */ + V1 = nextcomm (line); + if (V1 == (char *) NULL) goto error; + if (V1[0] != '=') goto error; + V1 ++; + + /*** assign vector element to value ***/ + val = dvomath (1, &V1, &size, 0); + if (val == (char *) NULL) { + print_error (); + goto error; + } + + /* find vector */ + if (isBuffer) { + if ((buf = SelectBuffer (V0, OLDBUFFER, TRUE)) == NULL) goto error; + free (V0); V0 = (char *) NULL; + + inRange = TRUE; + inRange &= (Nx < +1*buf[0].header.Naxis[0]); + inRange &= (Nx >= -1*buf[0].header.Naxis[0]); + inRange &= (Ny < +1*buf[0].header.Naxis[1]); + inRange &= (Ny >= -1*buf[0].header.Naxis[1]); + if (!inRange) { + gprint (GP_ERR, "no element %d,%d\n", Nx, Ny); + goto error; + } + if (Nx < 0) Nx += buf[0].header.Naxis[0]; + if (Ny < 0) Ny += buf[0].header.Naxis[1]; + + fptr = (float *) buf[0].matrix.buffer; + fptr[Nx + Ny*buf[0].header.Naxis[0]] = atof (val); + } else { + if ((vec = SelectVector (V0, OLDVECTOR, TRUE)) == NULL) goto error; + free (V0); V0 = (char *) NULL; + + inRange = TRUE; + inRange &= (Nx < +1*vec[0].Nelements); + inRange &= (Nx >= -1*vec[0].Nelements); + if (!inRange) { + gprint (GP_ERR, "no element %d\n", Nx); + goto error; + } + if (Nx < 0) Nx += vec[0].Nelements; + if (vec[0].type == OPIHI_FLT) { + vec[0].elements.Flt[Nx] = atof (val); + } else { + vec[0].elements.Int[Nx] = atol (val); + } + } + goto escape; + } + free (V0); V0 = (char *) NULL; + + /* case 3: {expression} */ + NLINE = MAX (128, strlen (line)); + ALLOCATE (newline, char, NLINE); + memset (newline, 0, NLINE); + for (L = line; *L != 0; ) { + + // copy elements from L (line) to newline up to first '{' + p = strchr (L, '{'); + newline = opihi_append (newline, &NLINE, L, p); + if (p == NULL) break; + L = p + 1; + + // check for \{ at this point, replace \ in newline with { + if ((p > line) && (*(p - 1) == 0x5c)) { + N = newline + strlen(newline) - 1; + *N = '{'; // replace \ with { + continue; + } + + /* check on the syntax of the line */ + L = p + 1; + c1 = strchr (L, '}'); + if (c1 == NULL) { + gprint (GP_ERR, "no close brackets!\n"); + goto error; + } + c2 = strchr (L, '{'); + if ((c2 != NULL) && (c2 < c1)) { + gprint (GP_ERR, "can't nest brackets!\n"); + goto error; + } + *c1 = 0; + + /* value in brackets must be a math expression of some sort. + isolated words are not valid here */ + val = dvomath (1, &L, &size, -1); + if (val == NULL) { + print_error (); + goto error; + } + + /* interpolate vector element into newline (being accumulated) */ + newline = opihi_append (newline, &NLINE, val, NULL); + + /* copy val to outline */ + L = c1 + 1; + free (val); val = (char *) NULL; + } + + free (line); line = (char *) NULL; + + /* \ protects next character */ + interpolate_slash (newline); + + REALLOCATE (newline, char, strlen (newline) + 1); + *status = TRUE; + return (newline); + +error: + gprint (GP_ERR, "syntax error\n"); + // assignment or increment operation: free line and return NULL + if (line != (char *) NULL) free (line); + if (val != (char *) NULL) free (val); + if (V0 != (char *) NULL) free (V0); + *status = FALSE; + return NULL; + +escape: + // assignment or increment operation: free line and return NULL + if (line != (char *) NULL) free (line); + if (val != (char *) NULL) free (val); + if (V0 != (char *) NULL) free (V0); + *status = TRUE; + return (NULL); +} + +/* this routine looks for math and/or logic expressions that + need to be evaluated and passes them on to "parenthesis", + which evaluates the expression */ + +/* There are two situations now which define an expression to + be evaluated: + 1) $var = expression -- everything after the = must be a math expression or a string + 2) {expression} -- everything within the {} must be a math expression +*/ + +/* case 1: $var = expression. test that + a) there is a $ as the first character + b) the variable defined by that $ is not null + c) there is an = sign following the variable + d) there is something following the = sign +*/ + +/* case 2: vect[N] = expression. check syntax: + a) last char must be ] + b) word must contain [ + c) between [ & ] must be an integer + d) vector must exist + e) there is an = sign following the first word + f) there is something following the = sign +*/ + + +/* case 3: we hunt for '{', and then the first '}' and pass everything + inbetween. no nested {{}}s are allowed! */ + +/* thisword ALLOCATES + thisvar ALLOCATES + thiscomm ALLOCATES + + nextword !ALLOCATE + nextcomm !ALLOCATE + + lastword !ALLOCATE + lastvar !ALLOCATE + + aftervar !ALLOCATE +*/ + Index: /trunk/Ohana/src/opihi/lib.shell/string.c =================================================================== --- /trunk/Ohana/src/opihi/lib.shell/string.c (revision 33661) +++ /trunk/Ohana/src/opihi/lib.shell/string.c (revision 33662) @@ -301,2 +301,23 @@ *out = *in; } + +// paste together argv[0] .. argv[N] into a single string +char *paste_args (int argc, char **argv) { + + int i; + + int length = 0; + for (i = 0; i < argc; i++) { + length += strlen(argv[i]) + 1; + } + + char *string = NULL; + ALLOCATE (string, char, length); + string[0] = 0; + for (i = 0; i < argc; i++) { + strcat (string, argv[i]); + if (i < argc - 1) strcat (string, " "); + } + return string; +} + Index: /trunk/Ohana/src/opihi/mana/findpeaks.c =================================================================== --- /trunk/Ohana/src/opihi/mana/findpeaks.c (revision 33661) +++ /trunk/Ohana/src/opihi/mana/findpeaks.c (revision 33662) @@ -76,7 +76,7 @@ } - ResetVector (vecx, OPIHI_FLT, MAX (Npeak, 1)); - ResetVector (vecy, OPIHI_FLT, MAX (Npeak, 1)); - ResetVector (vecz, OPIHI_FLT, MAX (Npeak, 1)); + ResetVector (vecx, OPIHI_FLT, Npeak); + ResetVector (vecy, OPIHI_FLT, Npeak); + ResetVector (vecz, OPIHI_FLT, Npeak); /* eliminate non-local peaks */ Index: /trunk/Ohana/src/opihi/mana/starcontour.c =================================================================== --- /trunk/Ohana/src/opihi/mana/starcontour.c (revision 33661) +++ /trunk/Ohana/src/opihi/mana/starcontour.c (revision 33662) @@ -22,6 +22,6 @@ N = 0; Npts = 100; - ResetVector (vecx, OPIHI_FLT, MAX (Npts, 1)); - ResetVector (vecy, OPIHI_FLT, MAX (Npts, 1)); + ResetVector (vecx, OPIHI_FLT, Npts); + ResetVector (vecy, OPIHI_FLT, Npts); Nx = buf[0].matrix.Naxis[0]; Index: /trunk/Ohana/src/opihi/pcontrol/Makefile =================================================================== --- /trunk/Ohana/src/opihi/pcontrol/Makefile (revision 33661) +++ /trunk/Ohana/src/opihi/pcontrol/Makefile (revision 33662) @@ -24,5 +24,4 @@ $(SRC)/init.$(ARCH).o \ $(SRC)/pcontrol.$(ARCH).o \ -$(SRC)/rconnect.$(ARCH).o \ $(SRC)/CheckBusyJob.$(ARCH).o \ $(SRC)/CheckDoneHost.$(ARCH).o \ Index: /trunk/Ohana/src/opihi/pcontrol/StartHost.c =================================================================== --- /trunk/Ohana/src/opihi/pcontrol/StartHost.c (revision 33661) +++ /trunk/Ohana/src/opihi/pcontrol/StartHost.c (revision 33662) @@ -14,10 +14,11 @@ if (VarConfig ("SHELL", "%s", shell) == NULL) strcpy (shell, "pclient"); - if (VerboseMode()) gprint (GP_ERR, "starting host within thread\n"); + if (VerboseMode()) gprint (GP_ERR, "starting remote connection to %s...", host[0].hostname); - pid = rconnect (command, host[0].hostname, shell, stdio); + int errorInfo; + pid = rconnect (command, host[0].hostname, shell, stdio, &errorInfo, TRUE); if (!pid) { /** failure to start: extend retry period **/ - if (VerboseMode()) gprint (GP_ERR, "failure to start %s\n", host[0].hostname); + if (VerboseMode()) gprint (GP_ERR, "failure to start %s (error %d)\n", host[0].hostname, errorInfo); gettimeofday (&now, (void *) NULL); if (ZTIME(host[0].next_start_try) || ZTIME(host[0].last_start_try)) {