/* wcyav.c - weighted cyclic averages Public domain by various programmers. This program may be used by any person for any purpose. */ /* qsortd() and dorder() are translations of the FORTRAN subroutines QSORTD and DORDER by Robert Renka. Collected in the Naval Surface Warfare Center Library, 1992. For unlimited distribution. */ /************************************************************ Robert Renka Oak Ridge Natl. Lab. This subroutine uses an order N*LOG(N) quick sort to sort the real array x into increasing order. The algorithm is as follows. ind is initialized to the ordered sequence of indices 1,...,n, and all interchanges are applied to ind. x is divided into two portions by picking a central element t. The first and last elements are compared with t, and interchanges are applied as necessary so that the three values are in ascending order. Interchanges are then applied so that all elements greater than t are in the upper portion of the array and all elements less than t are in the lower portion. The upper and lower indices of one of the portions are saved in local arrays, and the process is repeated iteratively on the other portion. When a portion is completely sorted, the process begins again by retrieving the indices bounding another unsorted portion. Input parameters - n - Length of the array x. x - Vector of length n to be sorted. ind - Vector of length .GE. n. n and x are not altered by this routine. Output parameter - ind - sequence of indices 1,...,n permuted in the same fashion as x would be. Thus, the ordering on x is defined by y(i) = x(ind(i)). *********************************************************** NOTE -- iu and il must be dimensioned .GE. LOG(N) where LOG has base 2. ************************************************************/ void qsortd (double *x, int *ind, int *n) { /* Local variables */ int i, j, k, l, m; double r, t; int ij, il[31], it, iu[31], itt, indx; /* Local parameters - iu,il = Temporary storage for the upper and lower indices of portions of the array x m = Index for iu and il i,j = Lower and upper indices of a portion of x k,l = Indices in the range i,...,j ij = Randomly chosen index between i and j it,itt = Temporary storage for interchanges in ind indx = Temporary index for x r = Pseudo random number for generating ij t = Central element of x */ /* Parameter adjustments */ --ind; --x; /* Function Body */ if (*n <= 0) return; /* Initialize ind, m, i, j, and r */ for (i = 1; i <= *n; ++i) ind[i] = i; m = 1; i = 1; j = *n; r = .375; /* Top of loop */ L2: if (i >= j) goto L10; if (r > .5898437) goto L3; r += .0390625; goto L4; L3: r += -.21875; /* Initialize k */ L4: k = i; /* Select a central element of x and save it in t */ ij = i + (int) (r * (double) (j - i)); it = ind[ij]; t = x[it]; /* If the first element of the array is greater than t, interchange it with t */ indx = ind[i]; if (x[indx] <= t) goto L5; ind[ij] = indx; ind[i] = it; it = indx; t = x[it]; /* Initialize l */ L5: l = j; /* If the last element of the array is less than t, interchange it with t */ indx = ind[j]; if (x[indx] >= t) goto L7; ind[ij] = indx; ind[j] = it; it = indx; t = x[it]; /* If the first element of the array is greater than t, interchange it with t */ indx = ind[i]; if (x[indx] <= t) goto L7; ind[ij] = indx; ind[i] = it; it = indx; t = x[it]; goto L7; /* Interchange elements k and l */ L6: itt = ind[l]; ind[l] = ind[k]; ind[k] = itt; /* Find an element in the upper part of the array which is not larger than t */ L7: --l; indx = ind[l]; if (x[indx] > t) goto L7; /* Find an element in the lower part of the array which is not smaller than t */ L8: ++k; indx = ind[k]; if (x[indx] < t) goto L8; /* If k .LE. l, interchange elements k and l */ if (k <= l) goto L6; /* Save the upper and lower subscripts of the portion of the array yet to be sorted */ if (l - i <= j - k) goto L9; il[m-1] = i; iu[m-1] = l; i = k; ++m; goto L11; L9: il[m-1] = k; iu[m-1] = j; j = l; ++m; goto L11; /* Begin again on another unsorted portion of the array */ L10: --m; if (m == 0) return; i = il[m-1]; j = iu[m-1]; L11: if (j - i >= 11) goto L4; if (i == 1) goto L2; --i; /* Sort elements i+1,...,j. Note that 1 .LE. i .LT. j and j-i .LT. 11. */ L12: ++i; if (i == j) goto L10; indx = ind[i+1]; t = x[indx]; it = indx; indx = ind[i]; if (x[indx] <= t) goto L12; k = i; L13: ind[k+1] = ind[k]; --k; indx = ind[k]; if (t < x[indx]) goto L13; ind[k+1] = it; goto L12; } /* end qsortd */ /************************************************************ Robert Renka Oak Ridge Natl. Lab. This routine applies a set of permutations to a vector. Input parameters - n - Length of a and ip. ip - Vector containing the sequence of integers 1,...,n permuted in the same fashion that a is to be permuted. a - Vector to be permuted. n and ip are not altered by this routine. Output parameter - a - Reordered vector reflecting the permutations defined by ip. ************************************************************/ void dorder (double *a, int *ip, int *n) { /* Local variables */ int j, k, nn, ipj; double temp; /* Local parameters - nn = Local copy of n k = Index for ip and for the first element of a in a permutation j = Index for ip and a, j .GE. k ipj = ip(j) temp = Temporary storage for a(k) */ /* Parameter adjustments */ --ip; --a; /* Function Body */ nn = *n; if (nn < 2) return; k = 1; /* Loop on permutations */ L1: j = k; temp = a[k]; /* Apply permutation to a. ip(j) is marked (made negative) as being included in the permutation. */ L2: ipj = ip[j]; ip[j] = -ipj; if (ipj == k) goto L3; a[j] = a[ipj]; j = ipj; goto L2; L3: a[j] = temp; /* Search for an unmarked element of ip */ L4: ++k; if (k > nn) goto L5; if (ip[k] > 0) goto L1; goto L4; /* All permutations have been applied. Unmark ip. */ L5: for (k = 1; k <= nn; ++k) ip[k] = -ip[k]; return; } /* end dorder */ /* wcymen() and wcymed() are variations of an algorithm introduced by David Radcliffe in 2005. By Andy Allinger, 2016, placed in the public domain. */ #include #include /*!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! weighted cyclic mean for phase only __Name___Type_______________In/Out__Description_________________________ *x double array In data array *w double array In importance weights *ind int array Work permutation vector n int In length of arrays t double In period *ave double array Out average *ties pointer to int Out how many values in AVE[] *dev pointer to double Out standard deviation !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!*/ void wcymen(double *x, double *w, int *ind, int n, double t, double *ave, int *ties, double *dev) { /* Local variables */ int i; double r, v, vmin, t2, sw, swx, swx2, tol; /* begin */ if (1 == n) { *ties = 1; ave[0] = x[0]; *dev = 0.; return; } tol = sqrt(t * DBL_EPSILON); t2 = t * t; sw = 0.; swx = 0.; swx2 = 0.; /* move all data to the range 0...t, sort */ for (i = 0; i < n; ++i) { r = fmod(x[i], t); if (r < 0.) r += t; sw += w[i]; // accumulate sums swx += w[i] * r; swx2 += w[i] * r * r; x[i] = r; } (void) qsortd(x, ind, &n); (void) dorder(w, ind, &n); (void) dorder(x, ind, &n); /* find phase of seam of lowest weighted deviation */ vmin = DBL_MAX; *ties = 0; for (i = 0; i < n; ++i) { v = swx2 - swx * swx / sw; // sw @ weighted variance if (vmin - v > tol) { // new low deviation *ties = 0; vmin = v; } if (fabs(v - vmin) < tol) { //tie for new low deviation ave[*ties] = swx / sw; // weighted mean ++(*ties); } /* apply update formula to sum and sum of squares */ swx += w[i] * t; swx2 += w[i] * (2. * x[i] * t + t2); } for (i = 0; i < *ties; ++i) { r = ave[i]; if (r >= t) r -= t; // subtract if out of 0...T ave[i] = r; } *dev = sqrt(vmin / sw); // (sw @ variance / sw)^(1/2) return; } /* end wcymen */ /*!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! weighted cyclic median for phase only !NOTE THAT! unlike the unweighted cyclic median, there is not an update formula available. The mean deviation is computed N repetitions, in O(N^2) time ! __Name___Type_______________In/Out________Description___________________ *x double array In data array *w double array In importance weights *ind int array Work permuatation vector n int In length of X t double In period *ave double array Out average *ties pointer to int Out how many values in AVE[] *dev pointer to double Out mean absolute deviation !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!*/ void wcymed (double *x, double *w, int *ind, int n, double t, double *ave, int *ties, double *dev) { /* Local variables */ int i, j, k, j1; double r, s1, sw, tol, hsw, smin; /* begin */ if (1 == n) { *ties = 1; ave[0] = x[0]; return; } tol = t * DBL_EPSILON; /* move all data to the range 0...t, sort */ sw = 0.; r = 0.; // avoid compiler warning for (i = 0; i < n; ++i) { r = fmod(x[i], t); if (r < 0.) r += t; x[i] = r; sw += w[i]; } if (sw <= 0.) return; hsw = sw / 2.; (void) qsortd (x, ind, &n); (void) dorder (w, ind, &n); (void) dorder (x, ind, &n); /* initial value of the cyclical median */ sw = 0.; j1 = 0; // avoid compiler warning for (i = 0; i < n; ++i) { sw += w[i]; if (sw > hsw) { r = x[i]; j1 = i; break; // done seeking initial median } else if (sw == hsw) { r = 0.5 * (x[i] + x[i+1]); j1 = i; break; // done seeking initial median } } s1 = 0.; // sum of weighted absolute deviations for (i = 0; i < n; ++i) { s1 += w[i] * fabs(x[i] - r); } /* find phase of seam of lowest deviation */ smin = DBL_MAX; *ties = 0; for (i = 0; i < n; ++i) { if (smin - s1 > tol) { // new low deviation *ties = 0; smin = s1; } if (fabs(s1 - smin) < tol) { // tie for new low deviation ave[*ties] = r; ++(*ties); } /* recompute median and deviation */ x[i] += t; // advance seam sw -= w[i]; // don't include in total k = j1; // possible for median to stay same for (j = 0; j < n; ++j) { if (sw > hsw) { r = x[k]; j1 = k; break; } else if (sw == hsw) { if (k == n-1) { r = 0.5 * (x[n-1] + x[0]); } else { r = 0.5 * (x[k] + x[k+1]); } j1 = k; break; } ++k; if (k == n) k = 0; sw += w[k]; } // next j s1 = 0.; // sum of weighted absolute deviations for (j = 0; j < n; ++j) { s1 += w[j] * fabs(x[j] - r); } } // next i for (i = 0; i < *ties; ++i) { r = ave[i]; if (r >= t) r -= t; // result in 0 <= r < t ave[i] = r; } *dev = 0.5 * smin / hsw; return; } /* end wcymed */