sddscorrelate.c File Reference

Detailed Description

Computes and evaluates correlations among columns of data in SDDS files.

This program reads an input SDDS file, calculates the correlation coefficients between specified or all numeric columns, and outputs the results to a new SDDS file. It supports options for rank-order correlation, standard deviation-based outlier removal, and column/row major data ordering.

Usage

sddscorrelate [<inputfile>] [<outputfile>]
              [-pipe=[input][,output]]
              [-columns=<list-of-names>]
              [-excludeColumns=<list-of-names>]
              [-withOnly=<name>]
              [-rankOrder]
              [-stDevOutlier[=limit=<factor>][,passes=<integer>]]
              [-majorOrder=row|column]
              [-threads=<number>]

Options

Option	Description
-pipe	Use standard input/output for input/output.
-columns	Specify columns to include in correlation analysis.
-excludeColumns	Specify columns to exclude from analysis.
-withOnly	Correlate only with the specified column.
-rankOrder	Use rank-order (Spearman) correlation instead of linear (Pearson).
-stDevOutlier	Remove outliers based on standard deviation.
-majorOrder	Set data ordering to row-major or column-major.
-threads	Set the number of threads for rank and correlation computation.

Incompatibilities

• -columns and -excludeColumns cannot be used together.
• -withOnly is mutually exclusive with -excludeColumns.

Features

• Linear (Pearson) or Rank-Order (Spearman) correlation calculation.
• Exclusion of specific columns or correlation with a single column.
• Outlier detection and removal using standard deviation thresholds.
• Customizable data ordering (row-major or column-major).

Output

The output SDDS file contains:

• Correlation coefficients and significance for column pairs.
• Number of points used in the correlation.
• Parameters summarizing the correlation analysis.

Copyright

• (c) 2002 The University of Chicago, as Operator of Argonne National Laboratory.
• (c) 2002 The Regents of the University of California, as Operator of Los Alamos National Laboratory.

License

This file is distributed under the terms of the Software License Agreement found in the file LICENSE included with this distribution.

Author

M. Borland, C. Saunders, R. Soliday

Definition in file sddscorrelate.c.

#include "mdb.h"
#include "SDDS.h"
#include "scan.h"
#include "SDDSutils.h"
#include <ctype.h>

Go to the source code of this file.

Functions
void	replaceWithRank (double *data, int64_t n)
double *	findRank (double *data, int64_t n)
void	markStDevOutliers (double *data, double limit, long passes, short *keep, int64_t n)
static long	correlationPairIndex (long i, long j, long columns)
int	main (int argc, char **argv)
int	compareData (const void *d1, const void *d2)

Function Documentation

compareData()

int compareData	(	const void *	d1,
		const void *	d2 )

Definition at line 450 of file sddscorrelate.c.

                                                {
  double diff = ((DATAnINDEX *)d1)->data - ((DATAnINDEX *)d2)->data;
  if (diff < 0)
    return -1;
  else if (diff > 0)
    return 1;
  else
    return 0;
}

correlationPairIndex()

static long correlationPairIndex ( long i, long j, long columns )

static

Definition at line 403 of file sddscorrelate.c.

                                                               {
  return i * (2 * columns - i - 1) / 2 + (j - i - 1);
}

findRank()

double * findRank	(	double *	data,
		int64_t	n )

Definition at line 460 of file sddscorrelate.c.

                                          {
  double *rank = malloc(sizeof(*rank) * n);
  if (!rank)
    return NULL;
  for (int64_t i = 0; i < n; i++)
    rank[i] = data[i];
  replaceWithRank(rank, n);
  return rank;
}

main()

int main	(	int	argc,
		char **	argv )

Definition at line 124 of file sddscorrelate.c.

                                {
  int iArg;
  char **column, **excludeColumn, *withOnly;
  long columns, excludeColumns;
  char *input, *output;
  SCANNED_ARG *scanned;
  SDDS_DATASET SDDSin, SDDSout;
  long i, j, row, readCode, rankOrder, iName1, iName2;
  long pairCount, pairIndex;
  int64_t rows;
  int32_t outlierStDevPasses;
  double **data, outlierStDevLimit;
  double *correlationValue, *significanceValue;
  long *correlationPoints;
  double **rank;
  short **accept;
  char s[SDDS_MAXLINE];
  unsigned long pipeFlags, dummyFlags, majorOrderFlag;
  short columnMajorOrder = -1;
  int threads = 1;
 
  SDDS_RegisterProgramName(argv[0]);
  argc = scanargs(&scanned, argc, argv);
  if (argc < 2)
    bomb(NULL, USAGE);
 
  output = input = withOnly = NULL;
  columns = excludeColumns = 0;
  column = excludeColumn = NULL;
  pipeFlags = 0;
  rankOrder = 0;
  outlierStDevPasses = 0;
  outlierStDevLimit = 1.0;
  rank = NULL;
  accept = NULL;
 
  for (iArg = 1; iArg < argc; iArg++) {
    if (scanned[iArg].arg_type == OPTION) {
      /* process options here */
      switch (match_string(scanned[iArg].list[0], option, N_OPTIONS, 0)) {
      case SET_MAJOR_ORDER:
        majorOrderFlag = 0;
        scanned[iArg].n_items--;
        if (scanned[iArg].n_items > 0 &&
            (!scanItemList(&majorOrderFlag, scanned[iArg].list + 1, &scanned[iArg].n_items, 0,
                           "row", -1, NULL, 0, SDDS_ROW_MAJOR_ORDER,
                           "column", -1, NULL, 0, SDDS_COLUMN_MAJOR_ORDER, NULL)))
          SDDS_Bomb("invalid -majorOrder syntax/values");
        if (majorOrderFlag & SDDS_COLUMN_MAJOR_ORDER)
          columnMajorOrder = 1;
        else if (majorOrderFlag & SDDS_ROW_MAJOR_ORDER)
          columnMajorOrder = 0;
        break;
      case SET_COLUMNS:
        if (columns)
          SDDS_Bomb("only one -columns option may be given");
        if (scanned[iArg].n_items < 2)
          SDDS_Bomb("invalid -columns syntax");
        column = tmalloc(sizeof(*column) * (columns = scanned[iArg].n_items - 1));
        for (i = 0; i < columns; i++)
          column[i] = scanned[iArg].list[i + 1];
        break;
      case SET_EXCLUDE:
        if (scanned[iArg].n_items < 2)
          SDDS_Bomb("invalid -excludeColumns syntax");
        moveToStringArray(&excludeColumn, &excludeColumns, scanned[iArg].list + 1, scanned[iArg].n_items - 1);
        break;
      case SET_WITHONLY:
        if (withOnly)
          SDDS_Bomb("only one -withOnly option may be given");
        if (scanned[iArg].n_items < 2)
          SDDS_Bomb("invalid -withOnly syntax");
        withOnly = scanned[iArg].list[1];
        break;
      case SET_PIPE:
        if (!processPipeOption(scanned[iArg].list + 1, scanned[iArg].n_items - 1, &pipeFlags))
          SDDS_Bomb("invalid -pipe syntax");
        break;
      case SET_RANKORDER:
        rankOrder = 1;
        break;
      case SET_STDEVOUTLIER:
        scanned[iArg].n_items--;
        outlierStDevPasses = 1;
        outlierStDevLimit = 1.0;
        if (!scanItemList(&dummyFlags, scanned[iArg].list + 1, &scanned[iArg].n_items, 0,
                          "limit", SDDS_DOUBLE, &outlierStDevLimit, 1, 0,
                          "passes", SDDS_LONG, &outlierStDevPasses, 1, 0, NULL) ||
            outlierStDevPasses <= 0 || outlierStDevLimit <= 0.0)
          SDDS_Bomb("invalid -stdevOutlier syntax/values");
        break;
      case SET_THREADS:
        if (scanned[iArg].n_items != 2 ||
            sscanf(scanned[iArg].list[1], "%d", &threads) != 1 || threads < 1)
          SDDS_Bomb("invalid -threads syntax");
        break;
      default:
        fprintf(stderr, "Error: unknown or ambiguous option: %s\n", scanned[iArg].list[0]);
        exit(EXIT_FAILURE);
        break;
      }
    } else {
      if (!input)
        input = scanned[iArg].list[0];
      else if (!output)
        output = scanned[iArg].list[0];
      else
        SDDS_Bomb("too many filenames seen");
    }
  }
 
  processFilenames("sddscorrelate", &input, &output, pipeFlags, 0, NULL);
 
  if (!SDDS_InitializeInput(&SDDSin, input))
    SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
 
  if (!columns)
    columns = appendToStringArray(&column, columns, "*");
  if (withOnly)
    columns = appendToStringArray(&column, columns, withOnly);
 
  if ((columns = expandColumnPairNames(&SDDSin, &column, NULL, columns, excludeColumn, excludeColumns, FIND_NUMERIC_TYPE, 0)) <= 0) {
    SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
    SDDS_Bomb("no columns selected for correlation analysis");
  }
 
  if (!SDDS_InitializeOutput(&SDDSout, SDDS_BINARY, 0, NULL, "sddscorrelate output", output) ||
      SDDS_DefineColumn(&SDDSout, "Correlate1Name", NULL, NULL, "Name of correlated quantity 1", NULL, SDDS_STRING, 0) < 0 ||
      SDDS_DefineColumn(&SDDSout, "Correlate2Name", NULL, NULL, "Name of correlated quantity 2", NULL, SDDS_STRING, 0) < 0 ||
      SDDS_DefineColumn(&SDDSout, "CorrelatePair", NULL, NULL, "Names of correlated quantities", NULL, SDDS_STRING, 0) < 0 ||
      SDDS_DefineColumn(&SDDSout, "CorrelationCoefficient", "r", NULL, "Linear correlation coefficient", NULL, SDDS_DOUBLE, 0) < 0 ||
      SDDS_DefineColumn(&SDDSout, "CorrelationSignificance", "P$br$n", NULL, "Linear correlation coefficient significance", NULL, SDDS_DOUBLE, 0) < 0 ||
      SDDS_DefineColumn(&SDDSout, "CorrelationPoints", NULL, NULL, "Number of points used for correlation", NULL, SDDS_LONG, 0) < 0 ||
      SDDS_DefineParameter(&SDDSout, "CorrelatedRows", NULL, NULL, "Number of data rows in correlation analysis", NULL, SDDS_LONG, NULL) < 0 ||
      SDDS_DefineParameter(&SDDSout, "sddscorrelateInputFile", NULL, NULL, "Data file processed by sddscorrelate", NULL, SDDS_STRING, input ? input : "stdin") < 0 ||
      SDDS_DefineParameter(&SDDSout, "sddscorrelateMode", NULL, NULL, NULL, NULL, SDDS_STRING, rankOrder ? "Rank-Order (Spearman)" : "Linear (Pearson)") < 0 ||
      SDDS_DefineParameter1(&SDDSout, "sddscorrelateStDevOutlierPasses", NULL, NULL, "Number of passes of standard-deviation outlier elimination applied", NULL, SDDS_LONG, &outlierStDevPasses) < 0 ||
      SDDS_DefineParameter1(&SDDSout, "sddscorrelateStDevOutlierLimit", NULL, NULL, "Standard-deviation outlier limit applied", NULL, SDDS_DOUBLE, &outlierStDevLimit) < 0) {
    SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
  }
 
  if (columnMajorOrder != -1)
    SDDSout.layout.data_mode.column_major = columnMajorOrder;
  else
    SDDSout.layout.data_mode.column_major = SDDSin.layout.data_mode.column_major;
 
  if (!SDDS_WriteLayout(&SDDSout))
    SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
 
  data = malloc(sizeof(*data) * columns);
  pairCount = columns * (columns - 1) / 2;
  correlationValue = pairCount ? malloc(sizeof(*correlationValue) * pairCount) : NULL;
  significanceValue = pairCount ? malloc(sizeof(*significanceValue) * pairCount) : NULL;
  correlationPoints = pairCount ? malloc(sizeof(*correlationPoints) * pairCount) : NULL;
  if (!data ||
      (pairCount && (!correlationValue || !significanceValue || !correlationPoints)) ||
      (rankOrder && !(rank = malloc(sizeof(*rank) * columns))) ||
      !(accept = malloc(sizeof(*accept) * columns))) {
    SDDS_Bomb("allocation failure");
  }
 
  while ((readCode = SDDS_ReadPage(&SDDSin)) > 0) {
    if ((rows = SDDS_CountRowsOfInterest(&SDDSin)) < 3)
      continue;
    if (!SDDS_StartPage(&SDDSout, columns * (columns - 1) / 2) ||
        !SDDS_SetParameters(&SDDSout, SDDS_SET_BY_NAME | SDDS_PASS_BY_VALUE, "CorrelatedRows", rows, NULL)) {
      SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
    }
    for (i = 0; i < columns; i++) {
      data[i] = SDDS_GetColumnInDoubles(&SDDSin, column[i]);
      if (!data[i])
        SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
      if (rankOrder)
        rank[i] = NULL;
      accept[i] = NULL;
      if (outlierStDevPasses) {
        accept[i] = malloc(sizeof(**accept) * rows);
        if (!accept[i])
          SDDS_Bomb("allocation failure");
      }
    }
#pragma omp parallel for if (threads > 1) num_threads(threads)
    for (i = 0; i < columns; i++) {
      if (rankOrder)
        rank[i] = findRank(data[i], rows);
      if (outlierStDevPasses)
        markStDevOutliers(data[i], outlierStDevLimit, outlierStDevPasses, accept[i], rows);
    }
#pragma omp parallel for private(j, iName1, iName2, pairIndex) if (threads > 1) num_threads(threads)
    for (i = 0; i < columns; i++) {
      for (j = i + 1; j < columns; j++) {
        long count;
        double correlation, significance;
        pairIndex = correlationPairIndex(i, j, columns);
        iName1 = i;
        iName2 = j;
        if (withOnly) {
          if (strcmp(withOnly, column[i]) == 0) {
            iName1 = j;
            iName2 = i;
          } else if (strcmp(withOnly, column[j]) == 0) {
            iName1 = i;
            iName2 = j;
          } else {
            correlationValue[pairIndex] = 0;
            significanceValue[pairIndex] = 0;
            correlationPoints[pairIndex] = 0;
            continue;
          }
        }
        correlation = linearCorrelationCoefficient(rankOrder ? rank[i] : data[i],
                                                   rankOrder ? rank[j] : data[j],
                                                   accept[i], accept[j], rows, &count);
        significance = linearCorrelationSignificance(correlation, count);
        correlationValue[pairIndex] = correlation;
        significanceValue[pairIndex] = significance;
        correlationPoints[pairIndex] = count;
      }
    }
    for (i = row = 0; i < columns; i++) {
      for (j = i + 1; j < columns; j++) {
        pairIndex = correlationPairIndex(i, j, columns);
        iName1 = i;
        iName2 = j;
        if (withOnly) {
          if (strcmp(withOnly, column[i]) == 0) {
            iName1 = j;
            iName2 = i;
          } else if (strcmp(withOnly, column[j]) == 0) {
            iName1 = i;
            iName2 = j;
          } else {
            continue;
          }
        }
        snprintf(s, sizeof(s), "%s.%s", column[iName1], column[iName2]);
        if (!SDDS_SetRowValues(&SDDSout, SDDS_SET_BY_INDEX | SDDS_PASS_BY_VALUE, row++,
                               0, column[iName1],
                               1, column[iName2],
                               2, s,
                               3, correlationValue[pairIndex],
                               4, significanceValue[pairIndex],
                               5, correlationPoints[pairIndex],
                               -1)) {
          SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
        }
      }
    }
    for (i = 0; i < columns; i++) {
      free(data[i]);
      if (rankOrder)
        free(rank[i]);
      if (accept[i])
        free(accept[i]);
    }
    if (!SDDS_WritePage(&SDDSout))
      SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
  }
 
  free(data);
  free(correlationValue);
  free(significanceValue);
  free(correlationPoints);
  if (rankOrder)
    free(rank);
  free(accept);
 
  if (!SDDS_Terminate(&SDDSin)) {
    SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
    exit(EXIT_FAILURE);
  }
  if (!SDDS_Terminate(&SDDSout)) {
    SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
    exit(EXIT_FAILURE);
  }
 
  return EXIT_SUCCESS;
}

markStDevOutliers()

void markStDevOutliers	(	double *	data,
		double	limit,
		long	passes,
		short *	keep,
		int64_t	n )

Definition at line 407 of file sddscorrelate.c.

                                                                                        {
  double sum1, sum2, variance, mean, absLimit;
  long pass;
  int64_t i, summed, kept;
 
  for (i = 0; i < n; i++)
    keep[i] = 1;
  kept = n;
  for (pass = 0; pass < passes && kept; pass++) {
    sum1 = 0.0;
    summed = 0;
    for (i = 0; i < n; i++) {
      if (keep[i]) {
        sum1 += data[i];
        summed += 1;
      }
    }
    if (summed < 2)
      return;
    mean = sum1 / summed;
    sum2 = 0.0;
    for (i = 0; i < n; i++) {
      if (keep[i])
        sum2 += sqr(data[i] - mean);
    }
    variance = sum2 / summed;
    if (variance > 0.0) {
      absLimit = limit * sqrt(variance);
      for (i = 0; i < n; i++) {
        if (keep[i] && fabs(data[i] - mean) > absLimit) {
          keep[i] = 0;
          kept--;
        }
      }
    }
  }
}

replaceWithRank()

void replaceWithRank	(	double *	data,
		int64_t	n )

Definition at line 470 of file sddscorrelate.c.

                                              {
  DATAnINDEX *indexedData = NULL;
  int64_t i, j, iStart, iEnd;
 
  indexedData = SDDS_Malloc(sizeof(*indexedData) * n);
  for (i = 0; i < n; i++) {
    indexedData[i].data = data[i];
    indexedData[i].originalIndex = i;
  }
  qsort(indexedData, n, sizeof(*indexedData), compareData);
  for (i = 0; i < n; i++)
    data[indexedData[i].originalIndex] = (double)i;
  for (i = 0; i < n - 1; i++) {
    if (data[i] == data[i + 1]) {
      iStart = i;
      for (j = i + 2; j < n; j++) {
        if (data[j] != data[i])
          break;
      }
      iEnd = j - 1;
      double averageRank = (iStart + iEnd) / 2.0;
      for (j = iStart; j <= iEnd; j++)
        data[indexedData[j].originalIndex] = averageRank;
      i = iEnd;
    }
  }
  free(indexedData);
}