sddsfindin2dgrid.c File Reference

Detailed Description

Searches a 2D grid to find locations based on specified SDDS input data.

This program processes 2D grid data from an SDDS input file to locate points where specific variables are closest to given target values. It supports both direct and file-based value inputs, optionally applies interpolation to refine the results, and writes output in SDDS format for further use or analysis.

Usage

sddsfindin2dgrid [<input>] [<output>] 
                 [-pipe=[input][,output]]
                 -gridVariables=<gridColumnName1>,<gridColumnName2>
                 -findLocationOf=<columnName1>,<columnName2>
                 {
                  [-valuesFile=<filename>] | 
                  [-atValues=<value1>,<value2>]
                 }
                 [-presorted]
                 [-interpolate]
                 [-mode={onePairPerPage|reuseFirstPage|all}]
                 [-inverse]
                 [-verbose]

Options

Required	Description
-gridVariables	Specifies the two columns representing the grid variables.
-findLocationOf	Specifies the two columns whose optimal grid values are to be found.
-valuesFile or -atValues	Specifies either a file with values or direct values for location finding.

Optional	Description
-pipe	Use SDDS pipe for input and/or output.
-presorted	Indicates that input data is pre-sorted, improving processing efficiency.
-interpolate	Enables linear interpolation for more precise location determination.
-mode	Sets the processing mode: onePairPerPage, reuseFirstPage, or all.
-inverse	Performs an inverse operation, finding data values for specified grid inputs.

Incompatibilities

• Only one of -valuesFile or -atValues may be specified.

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

Definition in file sddsfindin2dgrid.c.

#include "mdb.h"
#include "SDDS.h"
#include "scan.h"
#include <stdlib.h>
#include <stdio.h>
#include <inttypes.h>
#include <string.h>
#include <float.h>
#include <math.h>

Go to the source code of this file.

Functions
void	gridifyData (char *gridVariable[2], uint64_t gridPoints, short presorted)
int	findLocationInGrid (double at1, double at2, double *location, double *value, short interpolate)
uint64_t	findGridInterval (double value, short dimension)
int	main (int argc, char **argv)
double	distance (double *position, long *invalid)
int	compareGridLocations (const void *data1, const void *data2)

Function Documentation

compareGridLocations()

int compareGridLocations	(	const void *	data1,
		const void *	data2 )

Definition at line 565 of file sddsfindin2dgrid.c.

                                                               {
  uint64_t i1 = *((uint64_t *)data1);
  uint64_t i2 = *((uint64_t *)data2);
  double diff = -(gridValue[0][i1] - gridValue[0][i2]);
  if (diff == 0) {
    diff = -(gridValue[1][i1] - gridValue[1][i2]);
  }
 
  return (diff < 0 ? 1 : (diff > 0 ? -1 : 0));
}

distance()

double distance	(	double *	position,
		long *	invalid )

Definition at line 438 of file sddsfindin2dgrid.c.

                                                 {
  uint64_t ix, iy;
  double fx, fy;
  double v[2];
  long i;
 
  if (position[0] < 0 || position[0] >= ng[0] || position[1] < 0 || position[1] >= ng[1]) {
    *invalid = 1;
    printf("Invalid position: %.6le, %.6le for ng=%" PRIu64 ", %" PRIu64 "\n", position[0], position[1],
           ng[0], ng[1]);
    return DBL_MAX;
  }
  *invalid = 0;
 
  ix = position[0];
  iy = position[1];
  if (ix == (ng[0] - 1)) {
    ix--;
  }
  if (iy == (ng[1] - 1)) {
    iy--;
  }
 
  fx = position[0] - ix;
  fy = position[1] - iy;
 
  for (i = 0; i < 2; i++) {
    double v00, v10, v01, v11, v0, v1;
    v00 = valueAtLocation[i][(ix + 0) * ng[1] + (iy + 0)];
    v10 = valueAtLocation[i][(ix + 1) * ng[1] + (iy + 0)];
    v01 = valueAtLocation[i][(ix + 0) * ng[1] + (iy + 1)];
    v11 = valueAtLocation[i][(ix + 1) * ng[1] + (iy + 1)];
    v0 = v00 + (v10 - v00) * fx;
    v1 = v01 + (v11 - v01) * fx;
    v[i] = v0 + (v1 - v0) * fy;
    achieved[i] = v[i];
  }
 
  v[0] = target[0] - v[0];
  v[1] = target[1] - v[1];
  return v[0] * v[0] + v[1] * v[1];
}

findGridInterval()

uint64_t findGridInterval	(	double	value,
		short	dimension )

Definition at line 543 of file sddsfindin2dgrid.c.

                                                         {
  uint64_t lower = 0, upper = ng[dimension] - 1, middle;
  double coordinate;
 
  coordinate = dimension == 0 ? gridValue[0][0] : gridValue[1][0];
  if (value <= coordinate)
    return 0;
  coordinate = dimension == 0 ? gridValue[0][upper * ng[1]] : gridValue[1][upper];
  if (value >= coordinate)
    return upper - 1;
 
  while (upper - lower > 1) {
    middle = lower + (upper - lower) / 2;
    coordinate = dimension == 0 ? gridValue[0][middle * ng[1]] : gridValue[1][middle];
    if (value < coordinate)
      upper = middle;
    else
      lower = middle;
  }
  return lower;
}

findLocationInGrid()

int findLocationInGrid	(	double	at1,
		double	at2,
		double *	location,
		double *	value,
		short	interpolate )

Definition at line 481 of file sddsfindin2dgrid.c.

                                                                                                   {
  uint64_t ix, iy, j;
  uint64_t ixBest = ng[0] / 2, iyBest = ng[1] / 2;
  double delta, bestDelta = DBL_MAX;
  long i;
 
  for (ix = 0; ix < ng[0]; ix++) {
    for (iy = 0; iy < ng[1]; iy++) {
      j = ix * ng[1] + iy;
      double delta0 = valueAtLocation[0][j] - at1;
      double delta1 = valueAtLocation[1][j] - at2;
      delta = delta0 * delta0 + delta1 * delta1;
      if (delta < bestDelta) {
        location[0] = gridValue[0][j];
        location[1] = gridValue[1][j];
        value[0] = valueAtLocation[0][j];
        value[1] = valueAtLocation[1][j];
        ixBest = ix;
        iyBest = iy;
        bestDelta = delta;
      }
    }
  }
 
  if (interpolate) {
    double result, start[2], step[2], lower[2], upper[2];
    start[0] = ixBest;
    start[1] = iyBest;
    step[0] = step[1] = 0.1;
    lower[0] = lower[1] = 0;
    upper[0] = ng[0] - 1;
    upper[1] = ng[1] - 1;
    target[0] = at1;
    target[1] = at2;
    if (simplexMin(&result, start, step, lower, upper, NULL, 2, 0, 1e-14, distance,
                   NULL, 1500, 3, 12, 3, 1, 0) >= 0) {
      double a00, a01, a10, a11, a0, a1, fx, fy;
      int64_t cellX = start[0] < 0 ? 0 : (int64_t)start[0];
      int64_t cellY = start[1] < 0 ? 0 : (int64_t)start[1];
      if (cellX >= (int64_t)(ng[0] - 1))
        cellX = ng[0] - 2;
      if (cellY >= (int64_t)(ng[1] - 1))
        cellY = ng[1] - 2;
      fx = start[0] - cellX;
      fy = start[1] - cellY;
 
      for (i = 0; i < 2; i++) {
        a00 = gridValue[i][(cellX + 0) * ng[1] + (cellY + 0)];
        a10 = gridValue[i][(cellX + 1) * ng[1] + (cellY + 0)];
        a01 = gridValue[i][(cellX + 0) * ng[1] + (cellY + 1)];
        a11 = gridValue[i][(cellX + 1) * ng[1] + (cellY + 1)];
        a0 = a00 + (a10 - a00) * fx;
        a1 = a01 + (a11 - a01) * fx;
        location[i] = a0 + (a1 - a0) * fy;
        value[i] = achieved[i];
      }
    }
  }
 
  return 1;
}

gridifyData()

void gridifyData	(	char *	gridVariable[2],
		uint64_t	gridPoints,
		short	presorted )

Definition at line 576 of file sddsfindin2dgrid.c.

                                                                              {
  long i;
  char s[256];
  uint64_t j, *index;
  double *buffer;
 
  for (i = 0; i < 2; i++) {
    double *copy = tmalloc(sizeof(double) * gridPoints);
    memcpy(copy, gridValue[i], gridPoints * sizeof(double));
    qsort((void *)copy, gridPoints, sizeof(double), double_cmpasc);
    ng[i] = 1;
    for (j = 1; j < gridPoints; j++) {
      if (copy[j - 1] != copy[j]) {
        ng[i] += 1;
      }
    }
    free(copy);
    if (ng[i] == gridPoints) {
      snprintf(s, sizeof(s), "Grid variable %s has only unique values.\n", gridVariable[i]);
      SDDS_Bomb(s);
    }
    if (ng[i] == 1) {
      snprintf(s, sizeof(s), "Grid variable %s has only one unique value.\n", gridVariable[i]);
      SDDS_Bomb(s);
    }
  }
  if (ng[0] * ng[1] != gridPoints) {
    snprintf(s, sizeof(s), "Input data does not form a grid (nx = %" PRIu64 ", ny = %" PRIu64 ", rows = %" PRIu64 ")\n",
             ng[0], ng[1], gridPoints);
    SDDS_Bomb(s);
  }
 
  if (!presorted) {
    /* Sort indices for points to place data in gridded order */
    index = tmalloc(sizeof(uint64_t) * gridPoints);
    for (j = 0; j < gridPoints; j++) {
      index[j] = j;
    }
 
    qsort((void *)index, gridPoints, sizeof(uint64_t), compareGridLocations);
 
    /* Copy data into sorted order in the global variables */
    for (i = 0; i < 2; i++) {
      buffer = tmalloc(sizeof(double) * gridPoints);
      for (j = 0; j < gridPoints; j++) {
        buffer[j] = gridValue[i][index[j]];
      }
      free(gridValue[i]);
      gridValue[i] = buffer;
 
      buffer = tmalloc(sizeof(double) * gridPoints);
      for (j = 0; j < gridPoints; j++) {
        buffer[j] = valueAtLocation[i][index[j]];
      }
      free(valueAtLocation[i]);
      valueAtLocation[i] = buffer;
    }
    free(index);
  }
}

main()

int main	(	int	argc,
		char **	argv )

Definition at line 137 of file sddsfindin2dgrid.c.

                                {
  long iArg;
  SDDS_DATASET SDDSin, SDDSout, SDDSvalues;
  SCANNED_ARG *scanned;
  unsigned long pipeFlags = 0;
  uint64_t gridPoints = 0, atValues = 0, iv = 0, irow = 0;
  char *input = NULL, *output = NULL, *fileForValues = NULL;
  char *findLocationOf[2] = {NULL, NULL}, *gridVariable[2] = {NULL, NULL};
  double *atValue[2] = {NULL, NULL}, value[2] = {0.0, 0.0};
  short interpolate = 0, restarted = 0, needPage = 0, presorted = 0, inverse = 0, verbose = 0;
  unsigned long mode = MODE_ALL;
 
  SDDS_RegisterProgramName(argv[0]);
 
  argc = scanargs(&scanned, argc, argv);
  if (argc == 1) {
    fprintf(stderr, "%s", USAGE);
    exit(EXIT_FAILURE);
  }
 
  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 CLO_FINDLOCATIONOF:
        if (scanned[iArg].n_items != 3 ||
            !strlen(findLocationOf[0] = scanned[iArg].list[1]) ||
            !strlen(findLocationOf[1] = scanned[iArg].list[2])) {
          SDDS_Bomb("Invalid -findLocationOf syntax.\n");
        }
        if (strcmp(findLocationOf[0], findLocationOf[1]) == 0) {
          SDDS_Bomb("Invalid -findLocationOf values: two variables are the same.\n");
        }
        break;
      case CLO_GRIDVARIABLES:
        if (scanned[iArg].n_items != 3 ||
            !strlen(gridVariable[0] = scanned[iArg].list[1]) ||
            !strlen(gridVariable[1] = scanned[iArg].list[2])) {
          SDDS_Bomb("Invalid -gridVariables syntax.\n");
        }
        if (strcmp(gridVariable[0], gridVariable[1]) == 0) {
          SDDS_Bomb("Invalid -gridVariables values: two variables are the same.\n");
        }
        break;
      case CLO_VALUESFILE:
        if (scanned[iArg].n_items != 2 ||
            !strlen(fileForValues = scanned[iArg].list[1])) {
          SDDS_Bomb("Invalid -valuesFile syntax.\n");
        }
        if (atValues > 0) {
          SDDS_Bomb("Cannot use -valuesFile and -atValues together.\n");
        }
        break;
      case CLO_ATVALUES:
        atValue[0] = SDDS_Realloc(atValue[0], sizeof(double) * (atValues + 1));
        atValue[1] = SDDS_Realloc(atValue[1], sizeof(double) * (atValues + 1));
        if (scanned[iArg].n_items != 3 ||
            sscanf(scanned[iArg].list[1], "%le", &atValue[0][atValues]) != 1 ||
            sscanf(scanned[iArg].list[2], "%le", &atValue[1][atValues]) != 1) {
          SDDS_Bomb("Invalid -atValues syntax.\n");
        }
        if (fileForValues) {
          SDDS_Bomb("Cannot use -valuesFile and -atValues together.\n");
        }
        atValues++;
        break;
      case CLO_PIPE:
        if (!processPipeOption(scanned[iArg].list + 1, scanned[iArg].n_items - 1, &pipeFlags)) {
          SDDS_Bomb("Invalid -pipe syntax.\n");
        }
        break;
      case CLO_INTERPOLATE:
        interpolate = 1;
        break;
      case CLO_PRESORTED:
        presorted = 1;
        break;
      case CLO_MODE:
        mode = 0;
        if ((scanned[iArg].n_items -= 1) != 1 ||
            !scanItemList(&mode, scanned[iArg].list + 1, &scanned[iArg].n_items, 0,
                          "onepairperpage", -1, NULL, 0, MODE_ONEPAIRPERPAGE,
                          "reusefirstpage", -1, NULL, 0, MODE_REUSEFIRSTPAGE,
                          "all", -1, NULL, 0, MODE_ALL,
                          NULL) ||
            bitsSet(mode) != 1) {
          SDDS_Bomb("Invalid -mode syntax.\n");
        }
        break;
      case CLO_INVERSE:
        inverse = 1;
        break;
      case CLO_VERBOSE:
        verbose = 1;
        break;
      default:
        fprintf(stderr, "Invalid option: %s\n", scanned[iArg].list[0]);
        fprintf(stderr, "%s", USAGE);
        exit(EXIT_FAILURE);
      }
    } else {
      if (!input) {
        input = scanned[iArg].list[0];
      } else if (!output) {
        output = scanned[iArg].list[0];
      } else {
        SDDS_Bomb("Too many filenames provided.\n");
      }
    }
  }
 
  if (!findLocationOf[0] || !findLocationOf[1]) {
    SDDS_Bomb("Must provide -findLocationOf option.\n");
  }
  if (!gridVariable[0] || !gridVariable[1]) {
    SDDS_Bomb("Must provide -gridVariables option.\n");
  }
  if (!atValues && !fileForValues) {
    SDDS_Bomb("Must provide either -atValues or -valuesFile option.\n");
  }
 
  processFilenames("sddsfindin2dgrid", &input, &output, pipeFlags, 0, NULL);
 
  if (fileForValues) {
    if (!SDDS_InitializeInput(&SDDSvalues, fileForValues)) {
      SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
    }
    if (SDDS_ReadPage(&SDDSvalues) <= 0) {
      SDDS_Bomb("Unable to read values file.\n");
    }
    if ((atValues = SDDS_RowCount(&SDDSvalues)) > 0) {
      if (verbose)
    fprintf(stderr, "%" PRIu64 " values in values file\n", atValues);
      if (inverse) {
        if (!(atValue[0] = SDDS_GetColumnInDoubles(&SDDSvalues, gridVariable[0]))) {
          SDDS_Bomb("Unable to retrieve values of first grid variable in values file.\n");
        }
        if (!(atValue[1] = SDDS_GetColumnInDoubles(&SDDSvalues, gridVariable[1]))) {
          SDDS_Bomb("Unable to retrieve values of second grid variable in values file.\n");
        }
      } else {
        if (!(atValue[0] = SDDS_GetColumnInDoubles(&SDDSvalues, findLocationOf[0]))) {
          SDDS_Bomb("Unable to retrieve values of first findLocationOf variable in values file.\n");
        }
        if (!(atValue[1] = SDDS_GetColumnInDoubles(&SDDSvalues, findLocationOf[1]))) {
          SDDS_Bomb("Unable to retrieve values of second findLocationOf variable in values file.\n");
        }
      }
    }
    if (SDDS_ReadPage(&SDDSvalues) > 0) {
      SDDS_Bomb("Values file contains multiple pages, which is not supported.\n");
    }
    SDDS_Terminate(&SDDSvalues);
    SDDS_ClearErrors();
  }
 
  if (!SDDS_InitializeInput(&SDDSin, input)) {
    SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
  }
 
  if (!SDDS_InitializeOutput(&SDDSout, SDDS_BINARY, 0, NULL, NULL, output) ||
      !SDDS_TransferAllParameterDefinitions(&SDDSout, &SDDSin, 0) ||
      !SDDS_TransferColumnDefinition(&SDDSout, &SDDSin, gridVariable[0], NULL) ||
      !SDDS_TransferColumnDefinition(&SDDSout, &SDDSin, gridVariable[1], NULL) ||
      !SDDS_TransferColumnDefinition(&SDDSout, &SDDSin, findLocationOf[0], NULL) ||
      !SDDS_TransferColumnDefinition(&SDDSout, &SDDSin, findLocationOf[1], NULL) ||
      !SDDS_WriteLayout(&SDDSout) ||
      !SDDS_StartPage(&SDDSout, atValues * 1000)) {
    SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
  }
  SDDS_ClearErrors();
 
  irow = 0;
  iv = 0;
  restarted = 0;
  while (1) {
    double location[2];
    needPage = (irow == 0) ? 1 : 0;
    if (mode == MODE_ONEPAIRPERPAGE) {
      if (iv == atValues) {
        break;
      }
      needPage = 1;
    } else if (mode == MODE_REUSEFIRSTPAGE) {
      if (iv == atValues) {
        break;
      }
      if (iv == 0) {
        needPage = 1;
      }
    } else if (mode == MODE_ALL) {
      if (iv == atValues) {
        needPage = 1;
        iv = 0;
        restarted = 1; /* Signals that SDDS_ReadPage <= 0 is acceptable */
      }
    }
 
    if (needPage) {
      if (SDDS_ReadPage(&SDDSin) <= 0) {
        if (!restarted) {
          SDDS_Bomb("Too few pages in input file for number of location requests.\n");
        }
        break;
      }
      if (gridValue[0]) {
        free(gridValue[0]);
        free(gridValue[1]);
        free(valueAtLocation[0]);
        free(valueAtLocation[1]);
        gridValue[0] = gridValue[1] = valueAtLocation[0] = valueAtLocation[1] = NULL;
      }
      if ((gridPoints = SDDS_RowCount(&SDDSin)) <= 0) {
        SDDS_Bomb("First page of input file is empty.\n");
      }
      if (!(gridValue[0] = SDDS_GetColumnInDoubles(&SDDSin, gridVariable[0])) ||
          !(gridValue[1] = SDDS_GetColumnInDoubles(&SDDSin, gridVariable[1]))) {
        SDDS_Bomb("Grid variables are missing from input file.\n");
      }
      if (!(valueAtLocation[0] = SDDS_GetColumnInDoubles(&SDDSin, findLocationOf[0])) ||
          !(valueAtLocation[1] = SDDS_GetColumnInDoubles(&SDDSin, findLocationOf[1]))) {
        SDDS_Bomb("Location variables are missing from input file.\n");
      }
 
      gridifyData(gridVariable, gridPoints, presorted);
    }
 
    if (inverse) {
      /* Perform ordinary 2D interpolation on the grid */
      double x0, x1, y0, y1, v1, v2, fx, fy;
      double x, y;
      uint64_t ix, iy, ig;
 
      x = atValue[0][iv];
      y = atValue[1][iv];
      ix = findGridInterval(x, 0);
      iy = findGridInterval(y, 1);
      x0 = gridValue[0][ix * ng[1]];
      x1 = gridValue[0][(ix + 1) * ng[1]];
      y0 = gridValue[1][iy];
      y1 = gridValue[1][iy + 1];
      fx = (x - x0) / (x1 - x0);
      fy = (y - y0) / (y1 - y0);
 
      ig = ix * ng[1] + iy;
      v1 = valueAtLocation[0][ig] * (1 - fx) + valueAtLocation[0][ig + ng[1]] * fx;
      ig = ix * ng[1] + iy + 1;
      v2 = valueAtLocation[0][ig] * (1 - fx) + valueAtLocation[0][ig + ng[1]] * fx;
      location[0] = v1 * (1 - fy) + v2 * fy;
      value[0] = x;
 
      ig = ix * ng[1] + iy;
      v1 = valueAtLocation[1][ig] * (1 - fx) + valueAtLocation[1][ig + ng[1]] * fx;
      ig = ix * ng[1] + iy + 1;
      v2 = valueAtLocation[1][ig] * (1 - fx) + valueAtLocation[1][ig + ng[1]] * fx;
      location[1] = v1 * (1 - fy) + v2 * fy;
      value[1] = y;
      if (!SDDS_SetRowValues(&SDDSout, SDDS_SET_BY_INDEX | SDDS_PASS_BY_VALUE, irow++,
                             0, value[0], 1, value[1], 2, location[0], 3, location[1],
                             -1)) {
        SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
      }
    } else {
      if (verbose)
    fprintf(stderr, "Finding location for values (%le, %le)\n", atValue[0][iv], atValue[1][iv]);
      if (!findLocationInGrid(atValue[0][iv], atValue[1][iv], &location[0], &value[0], interpolate)) {
        fprintf(stderr, "Couldn't find location for %s=%.6le, %s=%.6le\n",
                findLocationOf[0], atValue[0][iv],
                findLocationOf[1], atValue[1][iv]);
        exit(EXIT_FAILURE);
      }
      if (verbose)
    fprintf(stderr, "Location is (%le, %le) with values (%le, %le)\n", location[0], location[1], value[0], value[1]);
      if (!SDDS_SetRowValues(&SDDSout, SDDS_SET_BY_INDEX | SDDS_PASS_BY_VALUE, irow++,
                             0, location[0], 1, location[1], 2, value[0], 3, value[1],
                             -1)) {
        SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
      }
    }
    iv++;
  }
 
  if (!SDDS_WritePage(&SDDSout) || !SDDS_Terminate(&SDDSout) ||
      !SDDS_Terminate(&SDDSin)) {
    SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
  }
  if (gridValue[0]) {
    free(gridValue[0]);
    free(gridValue[1]);
    free(valueAtLocation[0]);
    free(valueAtLocation[1]);
    gridValue[0] = gridValue[1] = valueAtLocation[0] = valueAtLocation[1] = NULL;
  }
  free(atValue[0]);
  free(atValue[1]);
  free_scanargs(&scanned, argc);
  return EXIT_SUCCESS;
}