SDDS ToolKit Programs and Libraries for C and Python
•All Classes Files Functions Variables Macros Pages
sddsoutlier.c
Go to the documentation of this file.
1/**
2 * @file sddsoutlier.c
3 * @brief Eliminates statistical outliers from SDDS data files.
4 *
5 * @details
6 * This program processes an input SDDS (Self Describing Data Set) file to identify and eliminate
7 * statistical outliers based on various criteria such as standard deviation limits, absolute
8 * limits, percentile limits, and more. Outliers can be removed, marked, or replaced with
9 * specified values. The processed data is then written to an output SDDS file.
10 *
11 * @section Usage
12 * ```
13 * sddsoutlier [<inputfile>] [<outputfile>]
14 * [-pipe=[input][,output]]
15 * [-verbose]
16 * [-noWarnings]
17 * -columns=<list-of-names>
18 * [-excludeColumns=<list-of-names>]
19 * -stDevLimit=<value>
20 * [-absLimit=<value>]
21 * [-absDeviationLimit=<value>[,neighbor=<number>]]
22 * [-maximumLimit=<value>]
23 * [-minimumLimit=<value>]
24 * [-chanceLimit=<minimumChance>]
25 * [-passes=<number>]
26 * [-percentileLimit=lower=<lowerPercent>,upper=<upperPercent>]
27 * [-unpopular=bins=<number>]
28 * [-invert]
29 * [-majorOrder]
30 * [-markOnly]
31 * [-replaceOnly={lastValue|nextValue|interpolatedValue|value=<number>}]
32 * ```
33 *
34 * @section Options
35 * | Required | Description |
36 * |---------------------------------------|---------------------------------------------------------------------------------------|
37 * | `-columns` | Specifies a list of column names to process. |
38 * | `-stDevLimit` | Sets the standard deviation limit for outlier detection. |
39 *
40 * | Optional | Description |
41 * |---------------------------------------|---------------------------------------------------------------------------------------|
42 * | `-pipe` | Use standard input/output as data streams. |
43 * | `-verbose` | Enables verbose output, displaying processing information. |
44 * | `-noWarnings` | Suppresses warning messages. |
45 * | `-excludeColumns` | Specifies a list of column names to exclude from processing. |
46 * | `-absLimit` | Sets an absolute value limit for outlier detection. |
47 * | `-absDeviationLimit` | Sets the absolute deviation limit from the mean with optional neighbors. |
48 * | `-maximumLimit` | Sets a maximum value for outlier detection. |
49 * | `-minimumLimit` | Sets a minimum value for outlier detection. |
50 * | `-chanceLimit` | Sets a minimum probability threshold for outlier detection (Gaussian statistics). |
51 * | `-passes` | Sets the number of passes for outlier detection. |
52 * | `-percentileLimit` | Sets percentile limits for outlier detection. |
53 * | `-unpopular` | Removes points not in the most populated bin of a histogram. |
54 * | `-invert` | Inverts the outlier selection criteria. |
55 * | `-majorOrder` | Specifies output file data ordering (row or column major). |
56 * | `-markOnly` | Marks identified outliers without removing them. |
57 * | `-replaceOnly` | Replaces outliers with specified values. |
58 *
59 * @subsection Incompatibilities
60 * - `-markOnly` is incompatible with:
61 * - `-replaceOnly`
62 * - For `-absDeviationLimit`:
63 * - Requires a positive value for the limit.
64 * - Optional `neighbor` parameter must be odd.
65 *
66 * @copyright
67 * - (c) 2002 The University of Chicago, as Operator of Argonne National Laboratory.
68 * - (c) 2002 The Regents of the University of California, as Operator of Los Alamos National Laboratory.
69 *
70 * @license
71 * This file is distributed under the terms of the Software License Agreement
72 * found in the file LICENSE included with this distribution.
73 *
74 * @authors
75 * M. Borland, C. Saunders, R. Soliday, H. Shang
76 */
77
78#include "mdb.h"
79#include "SDDS.h"
80#include "scan.h"
81#include "SDDSutils.h"
82#include <ctype.h>
83
84/* Enumeration for option types */
85enum option_type {
86 SET_COLUMNS,
87 SET_EXCLUDE,
88 SET_STDDEV_LIMIT,
89 SET_ABS_LIMIT,
90 SET_ABSDEV_LIMIT,
91 SET_VERBOSE,
92 SET_PIPE,
93 SET_NOWARNINGS,
94 SET_INVERT,
95 SET_MARKONLY,
96 SET_CHANCELIMIT,
97 SET_PASSES,
98 SET_REPLACE,
99 SET_MAXLIMIT,
100 SET_MINLIMIT,
101 SET_MAJOR_ORDER,
102 SET_PERCENTILE_LIMIT,
103 SET_UNPOPULAR,
104 N_OPTIONS
105};
106
107char *option[N_OPTIONS] = {
108 "columns",
109 "excludecolumns",
110 "stdevlimit",
111 "abslimit",
112 "absdeviationlimit",
113 "verbose",
114 "pipe",
115 "nowarnings",
116 "invert",
117 "markonly",
118 "chancelimit",
119 "passes",
120 "replaceonly",
121 "maximumlimit",
122 "minimumlimit",
123 "majororder",
124 "percentilelimit",
125 "unpopular",
126};
127
128#define USAGE "\n" \
129 "sddsoutlier [<inputfile>] [<outputfile>]\n" \
130 " [-pipe=[input][,output]]\n" \
131 " [-verbose] \n" \
132 " [-noWarnings] \n" \
133 " -columns=<list-of-names>\n" \
134 " [-excludeColumns=<list-of-names>]\n" \
135 " -stDevLimit=<value>\n" \
136 " [-absLimit=<value>] \n" \
137 " [-absDeviationLimit=<value>[,neighbor=<number>]]\n" \
138 " [-maximumLimit=<value>] \n" \
139 " [-minimumLimit=<value>]\n" \
140 " [-chanceLimit=<minimumChance>] \n" \
141 " [-passes=<number>]\n" \
142 " [-percentileLimit=lower=<lowerPercent>,upper=<upperPercent>]\n" \
143 " [-unpopular=bins=<number>]\n" \
144 " [-invert] \n" \
145 " [-majorOrder] \n" \
146 " [-markOnly]\n" \
147 " [-replaceOnly={lastValue|nextValue|interpolatedValue|value=<number>}]\n" \
148 "Options:\n" \
149 " -pipe=[input][,output]\n" \
150 " Use standard input and/or output as data streams.\n" \
151 " -verbose\n" \
152 " Enable verbose output, displaying processing information.\n" \
153 " -noWarnings\n" \
154 " Suppress warning messages.\n" \
155 " -columns=<list-of-names>\n" \
156 " Specify a comma-separated list of column names to process.\n" \
157 " -excludeColumns=<list-of-names>\n" \
158 " Specify a comma-separated list of column names to exclude from processing.\n" \
159 " -stDevLimit=<value>\n" \
160 " Point is an outlier if it is more than <value> standard deviations from the mean.\n" \
161 " -absLimit=<value>\n" \
162 " Point is an outlier if it has an absolute value greater than <value>.\n" \
163 " -absDeviationLimit=<value>[,neighbor=<number>]\n" \
164 " Point is an outlier if its absolute deviation from the mean exceeds <value>.\n" \
165 " If neighbor is provided, the mean is computed with the neighbors instead of the whole data.\n" \
166 " -minimumLimit=<value>\n" \
167 " Point is an outlier if it is less than <value>.\n" \
168 " -maximumLimit=<value>\n" \
169 " Point is an outlier if it is greater than <value>.\n" \
170 " -chanceLimit=<minimumChance>\n" \
171 " Point is an outlier if it has a probability less than <minimumChance> of occurring (Gaussian statistics).\n" \
172 " -percentileLimit=lower=<lowerPercent>,upper=<upperPercent>\n" \
173 " Point is an outlier if it is below the <lowerPercent> percentile or above the <upperPercent> percentile.\n" \
174 " -unpopular=bins=<number>\n" \
175 " Remove points that are not in the most populated bin based on a histogram with <number> bins.\n" \
176 " -invert\n" \
177 " Invert the outlier selection criteria.\n" \
178 " -majorOrder=row|column\n" \
179 " Specify output file data ordering as row or column major.\n" \
180 " -markOnly\n" \
181 " Mark identified outliers without removing them.\n" \
182 " -replaceOnly={lastValue|nextValue|interpolatedValue|value=<number>}\n" \
183 " Replace outliers with specified values or strategies.\n" \
184 " -passes=<number>\n" \
185 " Define the number of passes for outlier detection.\n\n" \
186 "Program by Michael Borland. (" __DATE__ " " __TIME__ ", SVN revision: " SVN_VERSION ")\n"
187
188#define OUTLIER_CONTROL_INVOKED 0x00001U
189#define OUTLIER_STDEV_GIVEN 0x00002U
190#define OUTLIER_FRACTION_GIVEN 0x00004U
191#define OUTLIER_STDEVLIMIT_GIVEN 0x00008U
192#define OUTLIER_UNPOPULAR_BINS 0x00010U
193#define OUTLIER_VERBOSE_GIVEN 0x00020U
194#define OUTLIER_ABSLIMIT_GIVEN 0x00040U
195#define OUTLIER_ABSDEVLIMIT_GIVEN 0x00080U
196#define OUTLIER_INVERT_GIVEN 0x00100U
197#define OUTLIER_MARKONLY 0x00200U
198#define OUTLIER_CHANCELIMIT_GIVEN 0x00400U
199#define OUTLIER_MAXLIMIT_GIVEN 0x00800U
200#define OUTLIER_MINLIMIT_GIVEN 0x01000U
201#define OUTLIER_REPLACELAST 0x02000U
202#define OUTLIER_REPLACENEXT 0x04000U
203#define OUTLIER_REPLACEINTERP 0x08000U
204#define OUTLIER_REPLACEVALUE 0x10000U
205#define OUTLIER_REPLACEFLAGS (OUTLIER_REPLACELAST | OUTLIER_REPLACENEXT | OUTLIER_REPLACEINTERP | OUTLIER_REPLACEVALUE)
206#define OUTLIER_PERCENTILE_LOWER 0x20000U
207#define OUTLIER_PERCENTILE_UPPER 0x40000U
208#define OUTLIER_PERCENTILE_FLAGS (OUTLIER_PERCENTILE_LOWER | OUTLIER_PERCENTILE_UPPER)
209
210typedef struct
211{
212 double stDevLimit, fractionLimit, absoluteLimit, absDevLimit;
213 double chanceLimit, replacementValue, maximumLimit, minimumLimit;
214 double percentilePoint[2];
215 long passes;
216 int32_t unpopularBins;
217 int32_t neighbors;
218 unsigned long flags;
219} OUTLIER_CONTROL;
220
221int64_t removeOutliers(SDDS_DATASET *SDDSin, int64_t rows, char **column, long columns, OUTLIER_CONTROL *outlierControl, int32_t *isOutlier);
222long meanStDevForFlaggedData(double *mean, double *stDev, double *data, int32_t *keep, int64_t rows);
223
224int main(int argc, char **argv) {
225 int iArg;
226 char **column, **excludeColumn;
227 long columns, excludeColumns;
228 char *input, *output;
229 SCANNED_ARG *scanned;
230 SDDS_DATASET SDDSin, SDDSout;
231 long readCode, dataLimitGiven, tmpfileUsed;
232 int64_t i, rows;
233 long noWarnings, isOutlierIndex;
234 int32_t *isOutlier;
235 OUTLIER_CONTROL outlierControl;
236 unsigned long pipeFlags, tmpFlags, majorOrderFlag, dummyFlags;
237 short columnMajorOrder = -1;
238
239 SDDS_RegisterProgramName(argv[0]);
240 argc = scanargs(&scanned, argc, argv);
241 if (argc < 3) {
242 bomb(NULL, USAGE);
243 }
244
245 output = input = NULL;
246 columns = excludeColumns = dataLimitGiven = 0;
247 column = excludeColumn = NULL;
248
249 outlierControl.flags = 0;
250 outlierControl.passes = 1;
251 outlierControl.neighbors = 0;
252 pipeFlags = tmpfileUsed = noWarnings = isOutlierIndex = 0;
253
254 for (iArg = 1; iArg < argc; iArg++) {
255 if (scanned[iArg].arg_type == OPTION) {
256 /* process options here */
257 switch (match_string(scanned[iArg].list[0], option, N_OPTIONS, 0)) {
258 case SET_MAJOR_ORDER:
259 majorOrderFlag = 0;
260 scanned[iArg].n_items--;
261 if (scanned[iArg].n_items > 0 &&
262 (!scanItemList(&majorOrderFlag, scanned[iArg].list + 1, &scanned[iArg].n_items, 0,
263 "row", -1, NULL, 0, SDDS_ROW_MAJOR_ORDER,
264 "column", -1, NULL, 0, SDDS_COLUMN_MAJOR_ORDER, NULL)))
265 SDDS_Bomb("invalid -majorOrder syntax/values");
266 if (majorOrderFlag & SDDS_COLUMN_MAJOR_ORDER)
267 columnMajorOrder = 1;
268 else if (majorOrderFlag & SDDS_ROW_MAJOR_ORDER)
269 columnMajorOrder = 0;
270 break;
271 case SET_COLUMNS:
272 if (columns)
273 SDDS_Bomb("only one -columns option may be given");
274 if (scanned[iArg].n_items < 2)
275 SDDS_Bomb("invalid -columns syntax");
276 column = tmalloc(sizeof(*column) * (columns = scanned[iArg].n_items - 1));
277 for (i = 0; i < columns; i++)
278 column[i] = scanned[iArg].list[i + 1];
279 break;
280 case SET_EXCLUDE:
281 if (excludeColumns)
282 SDDS_Bomb("only one -excludecolumns option may be given");
283 if (scanned[iArg].n_items < 2)
284 SDDS_Bomb("invalid -excludecolumns syntax");
285 excludeColumn = tmalloc(sizeof(*excludeColumn) * (excludeColumns = scanned[iArg].n_items - 1));
286 for (i = 0; i < excludeColumns; i++)
287 excludeColumn[i] = scanned[iArg].list[i + 1];
288 break;
289 case SET_STDDEV_LIMIT:
290 if (scanned[iArg].n_items != 2 || sscanf(scanned[iArg].list[1], "%lf", &outlierControl.stDevLimit) != 1 ||
291 outlierControl.stDevLimit <= 0)
292 SDDS_Bomb("invalid -stDevLimit syntax");
293 outlierControl.flags |= OUTLIER_CONTROL_INVOKED | OUTLIER_STDEV_GIVEN | OUTLIER_STDEVLIMIT_GIVEN;
294 break;
295 case SET_ABS_LIMIT:
296 if (scanned[iArg].n_items != 2 || sscanf(scanned[iArg].list[1], "%lf", &outlierControl.absoluteLimit) != 1 ||
297 outlierControl.absoluteLimit <= 0)
298 SDDS_Bomb("invalid -absLimit syntax");
299 outlierControl.flags |= OUTLIER_CONTROL_INVOKED | OUTLIER_ABSLIMIT_GIVEN;
300 break;
301 case SET_ABSDEV_LIMIT:
302 if (scanned[iArg].n_items < 2)
303 SDDS_Bomb("invalid -absDeviationLimit syntax");
304 if (scanned[iArg].n_items == 2) {
305 if (sscanf(scanned[iArg].list[1], "%lf", &outlierControl.absDevLimit) != 1 || outlierControl.absDevLimit <= 0)
306 SDDS_Bomb("invalid -absDeviationLimit syntax");
307 } else {
308 if (sscanf(scanned[iArg].list[1], "%lf", &outlierControl.absDevLimit) != 1 || outlierControl.absDevLimit <= 0)
309 SDDS_Bomb("invalid -absDeviationLimit syntax");
310 scanned[iArg].list += 2;
311 scanned[iArg].n_items -= 2;
312 if (scanned[iArg].n_items > 0 &&
313 (!scanItemList(&dummyFlags, scanned[iArg].list, &scanned[iArg].n_items, 0, "neighbors", SDDS_LONG, &(outlierControl.neighbors), 1, 0, NULL)))
314 SDDS_Bomb("invalid -absDeviationLimit syntax/value");
315 if (outlierControl.neighbors % 2 == 0)
316 outlierControl.neighbors += 1;
317 /* always make it an odd number */
318 scanned[iArg].list -= 2;
319 scanned[iArg].n_items += 2;
320 }
321 outlierControl.flags |= OUTLIER_CONTROL_INVOKED | OUTLIER_ABSDEVLIMIT_GIVEN;
322 break;
323 case SET_VERBOSE:
324 outlierControl.flags |= OUTLIER_VERBOSE_GIVEN;
325 break;
326 case SET_PIPE:
327 if (!processPipeOption(scanned[iArg].list + 1, scanned[iArg].n_items - 1, &pipeFlags))
328 SDDS_Bomb("invalid -pipe syntax");
329 break;
330 case SET_NOWARNINGS:
331 noWarnings = 1;
332 break;
333 case SET_INVERT:
334 outlierControl.flags |= OUTLIER_INVERT_GIVEN;
335 break;
336 case SET_MARKONLY:
337 outlierControl.flags |= OUTLIER_MARKONLY;
338 break;
339 case SET_CHANCELIMIT:
340 if (scanned[iArg].n_items != 2 ||
341 sscanf(scanned[iArg].list[1], "%lf", &outlierControl.chanceLimit) != 1 ||
342 outlierControl.chanceLimit <= 0)
343 SDDS_Bomb("invalid -chanceLimit syntax");
344 outlierControl.flags |= OUTLIER_CONTROL_INVOKED | OUTLIER_CHANCELIMIT_GIVEN;
345 break;
346 case SET_PASSES:
347 if (scanned[iArg].n_items != 2 ||
348 sscanf(scanned[iArg].list[1], "%ld", &outlierControl.passes) != 1 ||
349 outlierControl.passes < 1)
350 SDDS_Bomb("invalid -passes syntax");
351 break;
352 case SET_MAXLIMIT:
353 outlierControl.flags |= OUTLIER_MAXLIMIT_GIVEN | OUTLIER_CONTROL_INVOKED;
354 if (scanned[iArg].n_items != 2 ||
355 sscanf(scanned[iArg].list[1], "%lf", &outlierControl.maximumLimit) != 1)
356 SDDS_Bomb("invalid -maximumLimit syntax");
357 break;
358 case SET_MINLIMIT:
359 outlierControl.flags |= OUTLIER_MINLIMIT_GIVEN | OUTLIER_CONTROL_INVOKED;
360 if (scanned[iArg].n_items != 2 ||
361 sscanf(scanned[iArg].list[1], "%lf", &outlierControl.minimumLimit) != 1)
362 SDDS_Bomb("invalid -minimumLimit syntax");
363 break;
364 case SET_REPLACE:
365 if (scanned[iArg].n_items != 2)
366 SDDS_Bomb("invalid -replace syntax");
367 scanned[iArg].n_items -= 1;
368 if (!scanItemList(&tmpFlags, scanned[iArg].list + 1, &scanned[iArg].n_items, 0,
369 "lastvalue", -1, NULL, 0, OUTLIER_REPLACELAST,
370 "nextvalue", -1, NULL, 0, OUTLIER_REPLACENEXT,
371 "interpolatedvalue", -1, NULL, 0, OUTLIER_REPLACEINTERP,
372 "value", SDDS_DOUBLE, &outlierControl.replacementValue, 1, OUTLIER_REPLACEVALUE, NULL))
373 SDDS_Bomb("invalid -replace syntax/values");
374 outlierControl.flags |= tmpFlags | OUTLIER_CONTROL_INVOKED;
375 break;
376 case SET_PERCENTILE_LIMIT:
377 if (scanned[iArg].n_items < 3)
378 SDDS_Bomb("invalid -percentileLimit syntax");
379 scanned[iArg].n_items -= 1;
380 if (!scanItemList(&tmpFlags, scanned[iArg].list + 1, &scanned[iArg].n_items, 0,
381 "lower", SDDS_DOUBLE, &outlierControl.percentilePoint[0], 1, OUTLIER_PERCENTILE_LOWER,
382 "upper", SDDS_DOUBLE, &outlierControl.percentilePoint[1], 1, OUTLIER_PERCENTILE_UPPER, NULL) ||
383 !(tmpFlags & OUTLIER_PERCENTILE_LOWER) || !(tmpFlags & OUTLIER_PERCENTILE_UPPER) ||
384 outlierControl.percentilePoint[0] >= outlierControl.percentilePoint[1])
385 SDDS_Bomb("invalid -percentileLimit syntax");
386 outlierControl.flags |= tmpFlags | OUTLIER_CONTROL_INVOKED;
387 break;
388 case SET_UNPOPULAR:
389 if (scanned[iArg].n_items < 2)
390 SDDS_Bomb("invalid -unpopular syntax");
391 scanned[iArg].n_items -= 1;
392 if (!scanItemList(&tmpFlags, scanned[iArg].list + 1, &scanned[iArg].n_items, 0,
393 "bins", SDDS_LONG, &(outlierControl.unpopularBins), 1, OUTLIER_UNPOPULAR_BINS, NULL) ||
394 !(tmpFlags & OUTLIER_UNPOPULAR_BINS) || outlierControl.unpopularBins < 2)
395 SDDS_Bomb("invalid -unpopular syntax");
396 outlierControl.flags |= tmpFlags | OUTLIER_CONTROL_INVOKED;
397 break;
398 default:
399 fprintf(stderr, "Error: Unknown or ambiguous option: %s\n", scanned[iArg].list[0]);
400 exit(EXIT_FAILURE);
401 break;
402 }
403 } else {
404 if (!input)
405 input = scanned[iArg].list[0];
406 else if (!output)
407 output = scanned[iArg].list[0];
408 else
409 SDDS_Bomb("too many filenames seen");
410 }
411 }
412 if (outlierControl.flags & OUTLIER_REPLACEFLAGS && outlierControl.flags & OUTLIER_MARKONLY)
413 SDDS_Bomb("Cannot use -replaceOnly and -markOnly simultaneously.");
414
415 processFilenames("sddsoutlier", &input, &output, pipeFlags, noWarnings, &tmpfileUsed);
416
417 if (!(outlierControl.flags & OUTLIER_CONTROL_INVOKED)) {
418 outlierControl.flags |= OUTLIER_CONTROL_INVOKED | OUTLIER_STDEV_GIVEN | OUTLIER_STDEVLIMIT_GIVEN;
419 outlierControl.stDevLimit = 2;
420 }
421
422 if (!SDDS_InitializeInput(&SDDSin, input) ||
423 !SDDS_InitializeCopy(&SDDSout, &SDDSin, output, "w"))
424 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
425
426 if (columnMajorOrder != -1)
427 SDDSout.layout.data_mode.column_major = columnMajorOrder;
428 else
429 SDDSout.layout.data_mode.column_major = SDDSin.layout.data_mode.column_major;
430
431 if (outlierControl.flags & OUTLIER_MARKONLY &&
432 (isOutlierIndex = SDDS_GetColumnIndex(&SDDSout, "IsOutlier")) < 0) {
433 if (!SDDS_DefineColumn(&SDDSout, "IsOutlier", NULL, NULL, NULL, NULL, SDDS_SHORT, 0))
434 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
435 }
436
437 if (!SDDS_WriteLayout(&SDDSout))
438 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
439
440 if ((columns = expandColumnPairNames(&SDDSout, &column, NULL, columns, excludeColumn, excludeColumns, FIND_NUMERIC_TYPE, 0)) <= 0) {
441 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
442 SDDS_Bomb("No columns selected for outlier control.");
443 }
444
445 isOutlier = NULL;
446 while ((readCode = SDDS_ReadPage(&SDDSin)) > 0) {
447 if (!SDDS_CopyPage(&SDDSout, &SDDSin))
448 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
449
450 if ((rows = SDDS_CountRowsOfInterest(&SDDSout)) < 3) {
451 if (!SDDS_WritePage(&SDDSout))
452 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
453 continue;
454 }
455 if (outlierControl.flags & OUTLIER_MARKONLY) {
456 if (isOutlierIndex >= 0) {
457 if (!(isOutlier = SDDS_GetNumericColumn(&SDDSout, "IsOutlier", SDDS_LONG)))
458 SDDS_Bomb("Unable to retrieve 'IsOutlier' column from input file despite its existence.");
459 } else {
460 long i;
461 isOutlier = SDDS_Realloc(isOutlier, sizeof(*isOutlier) * rows);
462 if (!isOutlier)
463 SDDS_Bomb("Memory allocation failure.");
464 for (i = 0; i < rows; i++)
465 isOutlier[i] = 0;
466 }
467 }
468 if (outlierControl.flags & OUTLIER_VERBOSE_GIVEN)
469 fprintf(stderr, "%" PRId64 " rows in page %ld\n", rows, readCode);
470 if ((rows = removeOutliers(&SDDSout, rows, column, columns, &outlierControl, isOutlier)) == 0) {
471 if (!noWarnings)
472 fprintf(stderr, " No rows left after outlier control--skipping page.\n");
473 continue;
474 }
475 if (outlierControl.flags & OUTLIER_VERBOSE_GIVEN)
476 fprintf(stderr, "%" PRId64 " rows left after outlier control\n", rows);
477 if (rows != SDDS_CountRowsOfInterest(&SDDSout)) {
478 fprintf(stderr, "Problem with row selection:\n %" PRId64 " expected, %" PRId64 " counted\n", rows, SDDS_CountRowsOfInterest(&SDDSout));
479 exit(EXIT_FAILURE);
480 }
481 if ((isOutlier && !SDDS_SetColumnFromLongs(&SDDSout, SDDS_SET_BY_NAME, isOutlier, rows, "IsOutlier")) ||
482 !SDDS_WritePage(&SDDSout))
483 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
484 }
485 if (!SDDS_Terminate(&SDDSin) || !SDDS_Terminate(&SDDSout))
486 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
487 if (tmpfileUsed && !replaceFileAndBackUp(input, output))
488 exit(EXIT_FAILURE);
489 return EXIT_SUCCESS;
490}
491
492int64_t removeOutliers(SDDS_DATASET *dataset, int64_t rows, char **column, long columns, OUTLIER_CONTROL *outlierControl, int32_t *isOutlier) {
493 long icol, ipass;
494 int64_t irow, kept, killed, j, k, summed;
495 double *data, sum1, stDev, mean;
496 static int32_t *keep = NULL;
497 double lastGoodValue = 0;
498 int64_t irow0, irow1;
499
500 if (!SDDS_SetRowFlags(dataset, 1))
501 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
502
503 /* Allocate or reallocate the keep array */
504 keep = SDDS_Realloc(keep, sizeof(*keep) * rows);
505 if (!keep)
506 SDDS_Bomb("Memory allocation failure.");
507
508 if (!isOutlier) {
509 for (irow = 0; irow < rows; irow++)
510 keep[irow] = 1;
511 kept = rows;
512 } else {
513 for (irow = kept = 0; irow < rows; irow++)
514 if ((keep[irow] = !isOutlier[irow]))
515 kept++;
516 }
517
518 for (icol = 0; icol < columns; icol++) {
519 /* Loop over columns for which outlier control is to be done */
520 data = SDDS_GetColumnInDoubles(dataset, column[icol]);
521 if (!data)
522 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
523
524 for (ipass = 0; ipass < outlierControl->passes; ipass++) {
525 if (outlierControl->flags & OUTLIER_UNPOPULAR_BINS && rows > 1) {
526 double *hist, lo, hi, delta;
527 int64_t imin, imax, ih;
528 hist = tmalloc(sizeof(*hist) * outlierControl->unpopularBins);
529 find_min_max(&lo, &hi, data, rows);
530 make_histogram(hist, outlierControl->unpopularBins, lo, hi, data, rows, 1);
531 delta = (hi - lo) / outlierControl->unpopularBins; /* yes, not bins-1 */
532 index_min_max(&imin, &imax, hist, outlierControl->unpopularBins);
533 for (irow = killed = 0; irow < rows; irow++) {
534 ih = (data[irow] - lo) / delta;
535 if (ih != imax) {
536 killed++;
537 kept--;
538 keep[irow] = 0;
539 }
540 }
541 free(hist);
542 if (killed && (outlierControl->flags & OUTLIER_VERBOSE_GIVEN))
543 fprintf(stderr, "%" PRId64 " additional rows killed by column %s unpopular control\n", killed, column[icol]);
544 }
545
546 if (outlierControl->flags & OUTLIER_PERCENTILE_FLAGS) {
547 double percentileResult[2];
548 killed = 0;
549 if (compute_percentiles(percentileResult, outlierControl->percentilePoint, 2, data, rows)) {
550 for (irow = killed = 0; irow < rows; irow++) {
551 if ((data[irow] < percentileResult[0] || data[irow] > percentileResult[1]) && keep[irow]) {
552 killed++;
553 kept--;
554 keep[irow] = 0;
555 }
556 }
557 }
558 if (killed && (outlierControl->flags & OUTLIER_VERBOSE_GIVEN))
559 fprintf(stderr, "%" PRId64 " additional rows killed by column %s percentile outlier control\n", killed, column[icol]);
560 }
561
562 if (outlierControl->flags & OUTLIER_MINLIMIT_GIVEN) {
563 /* Limit values to exceed a minimum */
564 for (irow = killed = 0; irow < rows; irow++) {
565 if (keep[irow] && data[irow] < outlierControl->minimumLimit) {
566 kept--;
567 keep[irow] = 0;
568 killed++;
569 }
570 }
571 if (killed && (outlierControl->flags & OUTLIER_VERBOSE_GIVEN))
572 fprintf(stderr, "%" PRId64 " additional rows killed by column %s minimum value outlier control\n", killed, column[icol]);
573 }
574
575 if (outlierControl->flags & OUTLIER_MAXLIMIT_GIVEN) {
576 /* Limit values to exceed a maximum */
577 for (irow = killed = 0; irow < rows; irow++) {
578 if (keep[irow] && data[irow] > outlierControl->maximumLimit) {
579 kept--;
580 keep[irow] = 0;
581 killed++;
582 }
583 }
584 if (killed && (outlierControl->flags & OUTLIER_VERBOSE_GIVEN))
585 fprintf(stderr, "%" PRId64 " additional rows killed by column %s maximum value outlier control\n", killed, column[icol]);
586 }
587
588 if (outlierControl->flags & OUTLIER_ABSLIMIT_GIVEN) {
589 /* Limit absolute values */
590 for (irow = killed = 0; irow < rows; irow++) {
591 if (keep[irow] && fabs(data[irow]) > outlierControl->absoluteLimit) {
592 kept--;
593 keep[irow] = 0;
594 killed++;
595 }
596 }
597 if (killed && (outlierControl->flags & OUTLIER_VERBOSE_GIVEN))
598 fprintf(stderr, "%" PRId64 " additional rows killed by column %s absolute value outlier control\n", killed, column[icol]);
599 }
600
601 if (outlierControl->flags & OUTLIER_ABSDEVLIMIT_GIVEN) {
602 /* Limit absolute deviation from mean */
603 if (outlierControl->neighbors > 0) {
604 for (irow = killed = 0; irow < rows; irow++) {
605 if (!keep[irow])
606 continue;
607 mean = 0;
608 for (j = irow - outlierControl->neighbors / 2; j <= irow + outlierControl->neighbors / 2; j++) {
609 if (j < 0)
610 k = irow + outlierControl->neighbors / 2 - j;
611 else if (j > rows - 1)
612 k = irow - outlierControl->neighbors / 2 - (j - rows + 1);
613 else
614 k = j;
615 mean += fabs(data[k]);
616 }
617 mean = mean / outlierControl->neighbors;
618 if (keep[irow] && fabs(data[irow] - mean) > outlierControl->absDevLimit) {
619 keep[irow] = 0;
620 kept--;
621 killed++;
622 }
623 }
624 } else {
625 for (irow = sum1 = summed = 0; irow < rows; irow++) {
626 if (!keep[irow])
627 continue;
628 sum1 += data[irow];
629 summed++;
630 }
631 if (summed < 1)
632 continue;
633 mean = sum1 / summed;
634 for (irow = killed = 0; irow < rows; irow++)
635 if (keep[irow] && fabs(data[irow] - mean) > outlierControl->absDevLimit) {
636 keep[irow] = 0;
637 kept--;
638 killed++;
639 }
640 }
641 if (killed && (outlierControl->flags & OUTLIER_VERBOSE_GIVEN))
642 fprintf(stderr, "%" PRId64 " additional rows killed by column %s absolute deviation outlier control\n", killed, column[icol]);
643 }
644
645 if (outlierControl->flags & OUTLIER_STDEV_GIVEN && kept && meanStDevForFlaggedData(&mean, &stDev, data, keep, rows) && stDev) {
646 /* Limit deviation from mean in units of standard deviations */
647 for (irow = killed = 0; irow < rows; irow++)
648 if (keep[irow] && fabs(data[irow] - mean) > outlierControl->stDevLimit * stDev) {
649 keep[irow] = 0;
650 kept--;
651 killed++;
652 }
653 if (killed && (outlierControl->flags & OUTLIER_VERBOSE_GIVEN))
654 fprintf(stderr, "%" PRId64 " additional rows killed by column %s standard deviation outlier control\n", killed, column[icol]);
655 }
656
657 if (outlierControl->flags & OUTLIER_CHANCELIMIT_GIVEN) {
658 /* Limit improbability of a point occurring based on Gaussian statistics */
659 if (kept && meanStDevForFlaggedData(&mean, &stDev, data, keep, rows) && stDev) {
660 int64_t lastKept;
661 double gProb, probOfSeeing;
662 lastKept = kept;
663 for (irow = killed = 0; irow < rows; irow++) {
664 if (!keep[irow])
665 continue;
666 /* gProb = (probability of value >= x) */
667 /* (1-gProb)^n = probability of not seeing x in n trials */
668 gProb = normSigLevel((data[irow] - mean) / stDev, 2);
669 probOfSeeing = 1 - ipow(1 - gProb, lastKept);
670 if (probOfSeeing < outlierControl->chanceLimit) {
671 keep[irow] = 0;
672 kept--;
673 killed++;
674 }
675 }
676 if (killed && (outlierControl->flags & OUTLIER_VERBOSE_GIVEN))
677 fprintf(stderr, "%" PRId64 " additional rows killed by column %s chance limit outlier control\n", killed, column[icol]);
678 }
679 }
680 }
681
682 if (outlierControl->flags & OUTLIER_REPLACEFLAGS && (outlierControl->flags & OUTLIER_INVERT_GIVEN)) {
683 for (irow = 0; irow < rows; irow++)
684 keep[irow] = !keep[irow];
685 kept = rows - kept;
686 }
687
688 if (outlierControl->flags & OUTLIER_REPLACELAST) {
689 for (irow = 0; irow < rows; irow++) {
690 if (keep[irow]) {
691 lastGoodValue = data[irow];
692 break;
693 }
694 }
695 for (irow = 0; irow < rows; irow++) {
696 if (!keep[irow]) {
697 keep[irow] = 1;
698 data[irow] = lastGoodValue;
699 } else
700 lastGoodValue = data[irow];
701 }
702 kept = rows;
703 if (!SDDS_SetColumnFromDoubles(dataset, SDDS_SET_BY_NAME, data, rows, column[icol]))
704 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
705 } else if (outlierControl->flags & OUTLIER_REPLACENEXT) {
706 for (irow = rows - 1; irow >= 0; irow--) {
707 if (keep[irow]) {
708 lastGoodValue = data[irow];
709 break;
710 }
711 }
712 for (irow = rows - 1; irow >= 0; irow--) {
713 if (!keep[irow]) {
714 data[irow] = lastGoodValue;
715 keep[irow] = 1;
716 } else
717 lastGoodValue = data[irow];
718 }
719 kept = rows;
720 if (!SDDS_SetColumnFromDoubles(dataset, SDDS_SET_BY_NAME, data, rows, column[icol]))
721 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
722 } else if (outlierControl->flags & OUTLIER_REPLACEINTERP) {
723 irow0 = -1;
724 irow1 = -1;
725 for (irow = 0; irow < rows; irow++) {
726 if (!keep[irow]) {
727 if ((irow0 = irow - 1) >= 0) {
728 if ((irow1 = irow + 1) < rows) {
729 while (irow1 < rows && !keep[irow1])
730 irow1++;
731 }
732 if (irow1 < rows && keep[irow1]) {
733 /* irow is bracketed by irow0 and irow1, both of which have good data */
734 if (!keep[irow0])
735 SDDS_Bomb("Bracketing problem.");
736 for (; irow < irow1; irow++)
737 data[irow] = data[irow0] + (data[irow1] - data[irow0]) / (1.0 * irow1 - irow0) * (irow - irow0);
738 continue;
739 } else {
740 /* Ran off the end with bad points---just duplicate the last good point */
741 for (; irow < rows; irow++)
742 data[irow] = data[irow0];
743 continue;
744 }
745 } else {
746 /* No good point precedes this point---look for a good point following it */
747 for (irow1 = irow + 1; irow1 < rows; irow1++) {
748 if (keep[irow1])
749 break;
750 }
751 if (irow1 != rows) {
752 for (; irow < irow1; irow++)
753 data[irow] = data[irow1];
754 continue;
755 }
756 }
757 }
758 }
759 for (irow = 0; irow < rows; irow++)
760 keep[irow] = 1;
761 kept = rows;
762 if (!SDDS_SetColumnFromDoubles(dataset, SDDS_SET_BY_NAME, data, rows, column[icol]))
763 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
764 } else if (outlierControl->flags & OUTLIER_REPLACEVALUE) {
765 for (irow = 0; irow < rows; irow++) {
766 if (!keep[irow]) {
767 data[irow] = outlierControl->replacementValue;
768 keep[irow] = 1;
769 }
770 }
771 kept = rows;
772 if (!SDDS_SetColumnFromDoubles(dataset, SDDS_SET_BY_NAME, data, rows, column[icol]))
773 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
774 }
775 free(data);
776 }
777
778 if (outlierControl->flags & OUTLIER_INVERT_GIVEN) {
779 for (irow = 0; irow < rows; irow++)
780 keep[irow] = !keep[irow];
781 kept = rows - kept;
782 if (outlierControl->flags & OUTLIER_VERBOSE_GIVEN)
783 fprintf(stderr, "%" PRId64 " rows left after inversion\n", kept);
784 }
785
786 if (isOutlier && (outlierControl->flags & OUTLIER_MARKONLY)) {
787 for (irow = 0; irow < rows; irow++)
788 isOutlier[irow] = !keep[irow];
789 if (!SDDS_SetRowFlags(dataset, 1))
790 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
791 return rows;
792 }
793
794 if (!SDDS_AssertRowFlags(dataset, SDDS_FLAG_ARRAY, keep, rows))
795 SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
796 return kept;
797}
798
799long meanStDevForFlaggedData(double *mean, double *stDev, double *data, int32_t *keep, int64_t rows) {
800 int64_t irow, summed;
801 double sum1, sum2, value;
802
803 *mean = *stDev = 0;
804 for (irow = sum1 = summed = 0; irow < rows; irow++) {
805 if (!keep[irow])
806 continue;
807 sum1 += data[irow];
808 summed++;
809 }
810 if (summed < 2)
811 return 0;
812 *mean = sum1 / summed;
813 for (irow = sum2 = 0; irow < rows; irow++) {
814 if (keep[irow]) {
815 value = data[irow] - *mean;
816 sum2 += value * value;
817 }
818 }
819 *stDev = sqrt(sum2 / (summed - 1));
820 return 1;
821}
SDDS.h
SDDS (Self Describing Data Set) Data Types Definitions and Function Prototypes.
SDDS_InitializeCopy
int32_t SDDS_InitializeCopy(SDDS_DATASET *SDDS_target, SDDS_DATASET *SDDS_source, char *filename, char *filemode)
Definition SDDS_copy.c:40
SDDS_CopyPage
int32_t SDDS_CopyPage(SDDS_DATASET *SDDS_target, SDDS_DATASET *SDDS_source)
Definition SDDS_copy.c:578
SDDS_SetColumnFromDoubles
int32_t SDDS_SetColumnFromDoubles(SDDS_DATASET *SDDS_dataset, int32_t mode, double *data, int64_t rows,...)
Sets the values for a single data column using double-precision floating-point numbers.
Definition SDDS_dataprep.c:1642
SDDS_SetColumnFromLongs
int32_t SDDS_SetColumnFromLongs(SDDS_DATASET *SDDS_dataset, int32_t mode, int32_t *data, int64_t rows,...)
Sets the values for a single data column using long integer numbers.
Definition SDDS_dataprep.c:2045
SDDS_AssertRowFlags
int32_t SDDS_AssertRowFlags(SDDS_DATASET *SDDS_dataset, uint32_t mode,...)
Sets acceptance flags for rows based on specified criteria.
Definition SDDS_extract.c:148
SDDS_GetNumericColumn
void * SDDS_GetNumericColumn(SDDS_DATASET *SDDS_dataset, char *column_name, int32_t desiredType)
Retrieves the data of a specified numerical column as an array of a desired numerical type,...
Definition SDDS_extract.c:1617
SDDS_CountRowsOfInterest
int64_t SDDS_CountRowsOfInterest(SDDS_DATASET *SDDS_dataset)
Counts the number of rows marked as "of interest" in the current data table.
Definition SDDS_extract.c:364
SDDS_SetRowFlags
int32_t SDDS_SetRowFlags(SDDS_DATASET *SDDS_dataset, int32_t row_flag_value)
Sets the acceptance flags for all rows in the current data table of a data set.
Definition SDDS_extract.c:46
SDDS_GetColumnInDoubles
double * SDDS_GetColumnInDoubles(SDDS_DATASET *SDDS_dataset, char *column_name)
Retrieves the data of a specified numerical column as an array of doubles, considering only rows mark...
Definition SDDS_extract.c:872
SDDS_InitializeInput
int32_t SDDS_InitializeInput(SDDS_DATASET *SDDS_dataset, char *filename)
Definition SDDS_input.c:49
SDDS_Terminate
int32_t SDDS_Terminate(SDDS_DATASET *SDDS_dataset)
Definition SDDS_input.c:1448
SDDS_ReadPage
int32_t SDDS_ReadPage(SDDS_DATASET *SDDS_dataset)
Definition SDDS_input.c:1006
SDDS_WritePage
int32_t SDDS_WritePage(SDDS_DATASET *SDDS_dataset)
Writes the current data table to the output file.
Definition SDDS_output.c:1228
SDDS_DefineColumn
int32_t SDDS_DefineColumn(SDDS_DATASET *SDDS_dataset, const char *name, const char *symbol, const char *units, const char *description, const char *format_string, int32_t type, int32_t field_length)
Defines a data column within the SDDS dataset.
Definition SDDS_output.c:1715
SDDS_WriteLayout
int32_t SDDS_WriteLayout(SDDS_DATASET *SDDS_dataset)
Writes the SDDS layout header to the output file.
Definition SDDS_output.c:899
SDDS_GetColumnIndex
int32_t SDDS_GetColumnIndex(SDDS_DATASET *SDDS_dataset, char *name)
Retrieves the index of a named column in the SDDS dataset.
Definition SDDS_utils.c:1309
SDDS_PrintErrors
void SDDS_PrintErrors(FILE *fp, int32_t mode)
Prints recorded error messages to a specified file stream.
Definition SDDS_utils.c:432
SDDS_RegisterProgramName
void SDDS_RegisterProgramName(const char *name)
Registers the executable program name for use in error messages.
Definition SDDS_utils.c:288
SDDS_Bomb
void SDDS_Bomb(char *message)
Terminates the program after printing an error message and recorded errors.
Definition SDDS_utils.c:342
SDDS_Realloc
void * SDDS_Realloc(void *old_ptr, size_t new_size)
Reallocates memory to a new size.
Definition SDDS_utils.c:677
SDDS_LONG
#define SDDS_LONG
Identifier for the signed 32-bit integer data type.
Definition SDDStypes.h:61
SDDS_SHORT
#define SDDS_SHORT
Identifier for the signed short integer data type.
Definition SDDStypes.h:73
SDDS_DOUBLE
#define SDDS_DOUBLE
Identifier for the double data type.
Definition SDDStypes.h:37
SDDSutils.h
Utility functions for SDDS dataset manipulation and string array operations.
tmalloc
void * tmalloc(uint64_t size_of_block)
Allocates a memory block of the specified size with zero initialization.
Definition array.c:59
bomb
void bomb(char *error, char *usage)
Reports error messages to the terminal and aborts the program.
Definition bomb.c:26
index_min_max
int index_min_max(int64_t *imin, int64_t *imax, double *list, int64_t n)
Finds the indices of the minimum and maximum values in a list of doubles.
Definition findMinMax.c:116
find_min_max
int find_min_max(double *min, double *max, double *list, int64_t n)
Finds the minimum and maximum values in a list of doubles.
Definition findMinMax.c:33
ipow
double ipow(const double x, const int64_t p)
Compute x raised to the power p (x^p).
Definition ipow.c:33
make_histogram
long make_histogram(double *hist, long n_bins, double lo, double hi, double *data, int64_t n_pts, long new_start)
Compiles a histogram from data points.
Definition makeHistogram.c:33
match_string
long match_string(char *string, char **option, long n_options, long mode)
Matches a given string against an array of option strings based on specified modes.
Definition match_string.c:34
compute_percentiles
long compute_percentiles(double *position, double *percent, long positions, double *x, long n)
Computes multiple percentiles of an array of doubles.
Definition median.c:84
replaceFileAndBackUp
long replaceFileAndBackUp(char *file, char *replacement)
Replaces a file with a replacement file and creates a backup of the original.
Definition replacefile.c:75
scanargs
int scanargs(SCANNED_ARG **scanned, int argc, char **argv)
Definition scanargs.c:36
processPipeOption
long processPipeOption(char **item, long items, unsigned long *flags)
Definition scanargs.c:356
processFilenames
void processFilenames(char *programName, char **input, char **output, unsigned long pipeFlags, long noWarnings, long *tmpOutputUsed)
Definition scanargs.c:390
scanItemList
long scanItemList(unsigned long *flags, char **item, long *items, unsigned long mode,...)
Scans a list of items and assigns values based on provided keywords and types.
Definition scanitemlist.c:41
normSigLevel
double normSigLevel(double z0, long tails)
Computes the probability that a standard normal variable exceeds a given value.
Definition sigLevel.c:40
OUTLIER_CONTROL
Definition sddsoutlier.c:211
SDDS_DATASET
Definition SDDS.h:308
Generated by doxygen 1.12.0