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