samplegridicon: compute mean and std1.

src/Samplegridicon.c

@@ -2,13 +2,16 @@
 #include "cdo_int.h"
 #include "grid.h"
 
+
 typedef struct {
   int ncells;
   int *neighbor; // neighbor cell index
   int *parent;   // parent cell index
   int *child;    // child cell index
+  const char *filename;
 } cellindex_type;
 
+
 static
 void copy_data_to_index(int ncells, const double *restrict data, int *restrict cellindex)
 {
@@ -16,6 +19,15 @@ void copy_data_to_index(int ncells, const double *restrict data, int *restrict c
     cellindex[i] = (int) lround(data[i]);
 }
 
+static
+void free_cellindex(cellindex_type *cellindex)
+{
+  if ( cellindex->neighbor ) Free(cellindex->neighbor);
+  if ( cellindex->parent ) Free(cellindex->parent);
+  if ( cellindex->child ) Free(cellindex->child);
+  Free(cellindex);
+}
+
 static
 cellindex_type *read_cellindex(const char *filename)
 {
@@ -78,6 +90,10 @@ cellindex_type *read_cellindex(const char *filename)
         }
     }
 
+  // Fortran to C index
+  for ( int i = 0; i < 3*ncells; ++i ) cellindex->neighbor[i] -= 1;
+  for ( int i = 0; i < ncells; ++i ) cellindex->parent[i] -= 1;
+
   streamClose(streamID);
 
   Free(data);
@@ -113,114 +129,106 @@ int read_grid(const char *filename)
 }
 
 /**
-* Find the interval i-1 .. i in which an element x fits and return i, the 
-* bigger one of the interval borders or x itself if it is an interval border.
+* Return the first index of element x fits.
 *
-* If no interval can be found return the length of the array.
+* If no interval can be found return -1.
 
-* @param *array ascending or descending sorted list
-* @param nelem  length of the sorted list
-* @param x      the element to find a position for 
+* @param *array ascending sorted list
+* @param n      length of the sorted list
+* @param search the element to find a position for 
 */
 static
-int find_element(int x, int nelem, const int *restrict array)
+int find_index(int search, int n, const int *restrict array)
 {
-  int ii;
-  int mid = 0;
-  int first = 1;
-  int last = nelem;
-
-  if ( array[0] < array[nelem-1] ) // ascending order
+  int first = 0;
+  int last = n - 1;
+  int middle = (first+last)/2;
+ 
+  while ( first <= last )
     {
-      /* return the length of the array if x is out of bounds */
-      if ( x < array[0] || x > array[nelem-1] ) return nelem;
-
-      /* search for the interval in which x fits */
-      // implementation: binary search algorithm
-      for ( ii = 1; ii < nelem; ++ii )
-	{
-	  // binary search: divide search room in the middle
-	  // mid = first + ((last - first) >> 1);
-	  // faster!
-	  mid = (first + last) >> 1;
-      
-	  /* return the bigger interval border of the interval in which x fits */
-	  // if ( x >= array[mid-1] && x <= array[mid] ) break;
-	  // faster!
-	  if ( !(x < array[mid-1] || x > array[mid]) ) break;
-
-	  // binary search: ignore half of the search room
-	  if ( x > array[mid] )
-	    first = mid;
-	  else
-	    last = mid;
-	}
-    }
-  else
-    {
-      /* return the length of the array if x is out of bounds */
-      if ( x < array[nelem-1] || x > array[0] ) return nelem;
-
-      /* search for the interval in which x fits */
-      // implementation: binary search algorithm
-      for ( ii = 1; ii < nelem; ++ii )
-	{
-	  // binary search: divide search room in the middle
-	  // mid = first + ((last - first) >> 1);
-	  // faster!
-	  mid = (first + last) >> 1;
-      
-	  /* return the bigger interval border of the interval in which x fits */
-	  // if ( x >= array[mid] && x <= array[mid-1] ) break;
-	  // faster!
-	  if ( !(x < array[mid] || x > array[mid-1]) ) break;
-
-	  // binary search: ignore half of the search room
-	  if ( x < array[mid] )
-	    first = mid;
-	  else
-	    last = mid;
-	}
+      if ( array[middle] < search )
+        first = middle + 1;    
+      else if ( array[middle] == search )
+        {
+          for ( int i = middle; i >= 0; i-- )
+            {
+              if ( array[i] == search ) middle = i;
+              else break;
+            }
+          return middle;
+        }
+      else
+        last = middle - 1;
+ 
+      middle = (first + last)/2;
     }
 
-  if ( mid > 1 && IS_EQUAL(x,array[mid-1]) ) mid--;
-
-  return mid;
+  return -1;
 }
 
 static
-void samplegrid(cellindex_type *cellindex1, double *array1, cellindex_type *cellindex2, double *array2, int *samp2)
+void compute_child(cellindex_type *cellindex1, cellindex_type *cellindex2)
 {
   int ncells1 = cellindex1->ncells;
+  int *parent1 = cellindex1->parent;
   int ncells2 = cellindex2->ncells;
+  int *child2 = (int*) Malloc(4*ncells2*sizeof(int));
+  cellindex2->child = child2;
   for ( int i = 0; i< ncells2; ++i )
     {
-      int j = find_element(i+1, ncells1, cellindex1->parent);
-      if ( i%10000 == 0 )
-          printf("%d %d %d %d %d %d %d\n", i, j, cellindex1->parent[j-2], cellindex1->parent[j-1], cellindex1->parent[j],
-                 cellindex1->parent[j+1], cellindex1->parent[j+2]);
-      if ( j == i+1 )
+      for ( int k = 0; k < 4; ++k ) child2[i*4+k] = -1;
+      int j = find_index(i, ncells1, parent1);
+      if ( j < 0 ) continue;
+      for ( int k = 0; k < 4; ++k )
         {
+          if ( i != parent1[j+k] ) break;
+          child2[i*4+k] = j+k;
         }
+     // if ( i%10000 == 0 ) printf("%d %d %d %d %d %d\n", i, j, parent1[j], parent1[j+1], parent1[j+2], parent1[j+3]);
     }
 }
 
 static
-void compute_child(cellindex_type *cellindex1, cellindex_type *cellindex2)
+void compute_sum(int i, int *n, double *sum, double *sumq, int kci, cellindex_type **cellindex, double *array)
 {
-  int ncells2 = cellindex2->ncells;
-  cellindex2->child = (int*) Malloc(4*ncells2*sizeof(int));
+  // printf("compute: i, kci %d %d\n", i, kci);
+  int ncells2 = cellindex[kci]->ncells;
+  if ( i < 0 || i > ncells2 ) cdoAbort("Child grid cell index %d out of bounds %d!", i, ncells2);
+
+  for ( int k = 0; k < 4; ++k )
+    {
+      int index = cellindex[kci]->child[i*4+k];
+      if ( index == -1 ) break;
+      if ( kci == 1 )
+        {
+          *sum += array[index];
+          *sumq += array[index]*array[index];
+          *n += 1;
+        }
+      else compute_sum(index, n, sum, sumq, kci-1, cellindex, array);
+    }
 }
 
-/*
 static
-void map_index(int num_links, int dst_size, int *src_idx, int *dst_idx, double *src_array, double *dst_array, double missval)
+void samplegrid(double missval, int nci, cellindex_type **cellindex, double *array1, double *array2, double *array3)
 {
-  for ( int n = 0; n < dst_size; ++n ) dst_array[n] = missval;
-  for ( int n = 0; n < num_links; ++n ) dst_array[dst_idx[n]] = 0.;
-  for ( int n = 0; n < num_links; ++n ) dst_array[dst_idx[n]] += src_array[src_idx[n]];
+  int kci = nci-1;
+  int ncells2 = cellindex[kci]->ncells;
+#if defined(_OPENMP)
+#pragma omp parallel for default(none) shared(missval, ncells2, kci, cellindex, array1, array2, array3)
+#endif
+  for ( int i = 0; i < ncells2; ++i )
+    {
+      int n = 0;
+      double sum = 0, sumq = 0;
+      compute_sum(i, &n, &sum, &sumq, kci, cellindex, array1);
+      array2[i] = n ? sum/n : missval;  // mean
+      double var1 = (n*n > n) ? (sumq*n - sum*sum) / (n*n - n) : missval;
+      if ( var1 < 0 && var1 > -1.e-5 ) var1 = 0;
+      array3[i] = var_to_std(var1, missval); // std1
+    }
 }
-*/
+
 
 #include "pstream.h"
 
@@ -229,7 +237,6 @@ void *Samplegridicon(void *argument)
 {
   int nrecs;
   int varID, levelID;
-  int nmiss;
 
   cdoInitialize(argument);
 
@@ -238,23 +245,31 @@ void *Samplegridicon(void *argument)
   int nsamplegrids = operatorArgc();
   if ( nsamplegrids < 2 ) cdoAbort("Parameter missing!");
 
-  cellindex_type *cellindex1 = read_cellindex(operatorArgv()[0]);
-  printf("ncells %d\n", cellindex1->ncells);
-  cellindex_type *cellindex2 = read_cellindex(operatorArgv()[1]);
-  printf("ncells %d\n", cellindex2->ncells);
+  cellindex_type *cellindex[nsamplegrids];
+  
+  for ( int i = 0; i < nsamplegrids; ++i )
+    {
+      cellindex[i] = read_cellindex(operatorArgv()[i]);
+      cellindex[i]->filename = operatorArgv()[i];
+      if ( cdoVerbose ) cdoPrint("Found %d grid cells in %s", cellindex[i]->ncells, cellindex[i]->filename);
+    }
 
-  compute_child(cellindex1, cellindex2);
+  for ( int i = 0; i < nsamplegrids-1; ++i )
+    compute_child(cellindex[i], cellindex[i+1]);
   
-  int gridID2 = read_grid(operatorArgv()[1]);
+  int gridID2 = read_grid(operatorArgv()[nsamplegrids-1]);
 
   int streamID1 = streamOpenRead(cdoStreamName(0));
 
   int vlistID1 = streamInqVlist(streamID1);
   int vlistID2 = vlistDuplicate(vlistID1);
+  int vlistID3 = vlistDuplicate(vlistID1);
 
   int taxisID1 = vlistInqTaxis(vlistID1);
   int taxisID2 = taxisDuplicate(taxisID1);
+  int taxisID3 = taxisDuplicate(taxisID1);
   vlistDefTaxis(vlistID2, taxisID2);
+  vlistDefTaxis(vlistID3, taxisID3);
 
   int ngrids = vlistNgrids(vlistID1);
   for ( int index = 0; index < ngrids; ++index )
@@ -265,73 +280,73 @@ void *Samplegridicon(void *argument)
         cdoAbort("Unsupported gridtype: %s with %d corners", gridNamePtr(gridtype), gridInqNvertex(gridID));
 
       vlistChangeGridIndex(vlistID2, index, gridID2);
+      vlistChangeGridIndex(vlistID3, index, gridID2);
     }
 
   int streamID2 = streamOpenWrite(cdoStreamName(1), cdoFiletype());
-
   streamDefVlist(streamID2, vlistID2);
 
+  int streamID3 = streamOpenWrite(cdoStreamName(2), cdoFiletype());
+  streamDefVlist(streamID3, vlistID3);
+
   int gridsize = vlistGridsizeMax(vlistID1);
   if ( vlistNumber(vlistID1) != CDI_REAL ) gridsize *= 2;
   double *array1 = (double *) Malloc(gridsize*sizeof(double));
 
   int gridsize2 = gridInqSize(gridID2);
+  if ( cdoVerbose ) cdoPrint("Target gridsize = %d", gridsize2);
   if ( vlistNumber(vlistID2) != CDI_REAL ) gridsize2 *= 2;
   double *array2 = (double *) Malloc(gridsize2*sizeof(double));
-  int *samp2 = (int *) Malloc(gridsize2*sizeof(int));
-  for ( int i = 0; i < gridsize2; ++i ) samp2[i] = 0;
-  for ( int i = 0; i < gridsize2; ++i ) array2[i] = 0;
+  double *array3 = (double *) Malloc(gridsize2*sizeof(double));
 
   int tsID = 0;
   while ( (nrecs = streamInqTimestep(streamID1, tsID)) )
     {
       taxisCopyTimestep(taxisID2, taxisID1);
-
+      taxisCopyTimestep(taxisID3, taxisID1);
       streamDefTimestep(streamID2, tsID);
+      streamDefTimestep(streamID3, tsID);
 
       for ( int recID = 0; recID < nrecs; recID++ )
         {
+          int nmiss;
           streamInqRecord(streamID1, &varID, &levelID);
           streamReadRecord(streamID1, array1, &nmiss);
 
-          streamDefRecord(streamID2, varID, levelID);
+          double missval = vlistInqVarMissval(vlistID1, varID);
 
-          int gridID = vlistInqVarGrid(vlistID1, varID);
+          samplegrid(missval, nsamplegrids, cellindex, array1, array2, array3);
 
-          //for ( int i = 0; i < gridsize2; ++i ) array2[i] = array1[i];
-          samplegrid(cellindex1, array1, cellindex2, array2, samp2);
-
-          if ( nmiss )
-            {
-              nmiss = 0;
-              double missval = vlistInqVarMissval(vlistID2, varID);
-              for ( int i = 0; i < gridsize2; i++ )
-                if ( DBL_IS_EQUAL(array2[i], missval) ) nmiss++;
-            }
+          nmiss = 0;
+          for ( int i = 0; i < gridsize2; ++i )
+            if ( DBL_IS_EQUAL(array2[i], missval) ) nmiss++;
 
+          streamDefRecord(streamID2, varID, levelID);
           streamWriteRecord(streamID2, array2, nmiss);
+
+          nmiss = 0;
+          for ( int i = 0; i < gridsize2; ++i )
+            if ( DBL_IS_EQUAL(array3[i], missval) ) nmiss++;
+
+          streamDefRecord(streamID3, varID, levelID);
+          streamWriteRecord(streamID3, array3, nmiss);
         }
 
       tsID++;
     }
 
+  streamClose(streamID3);
   streamClose(streamID2);
   streamClose(streamID1);
 
   vlistDestroy(vlistID2);
   gridDestroy(gridID2);
 
-  if ( samp2 ) Free(samp2);
+  if ( array3 ) Free(array3);
   if ( array2 ) Free(array2);
   if ( array1 ) Free(array1);
 
-  if ( cellindex1->neighbor ) Free(cellindex1->neighbor);
-  if ( cellindex1->parent ) Free(cellindex1->parent);
-  if ( cellindex1->child ) Free(cellindex1->child);
-
-  if ( cellindex2->neighbor ) Free(cellindex2->neighbor);
-  if ( cellindex2->parent ) Free(cellindex2->parent);
-  if ( cellindex2->child ) Free(cellindex2->child);
+  for ( int i = 0; i < nsamplegrids; ++i ) free_cellindex(cellindex[i]);
 
   cdoFinish();
 

src/modules.c

@@ -699,7 +699,7 @@ static modules_t Modules[] =
   { Rotuv,          RotuvbHelp,        RotuvOperators,         1,   CDI_REAL,  1,  1 },
   { Runpctl,        RunpctlHelp,       RunpctlOperators,       1,   CDI_REAL,  1,  1 },
   { Runstat,        RunstatHelp,       RunstatOperators,       1,   CDI_REAL,  1,  1 },
-  { Samplegridicon, NULL,              SamplegridiconOperators,1,   CDI_REAL,  1,  1 },
+  { Samplegridicon, NULL,              SamplegridiconOperators,1,   CDI_REAL,  1,  2 },
   { Seascount,      NULL,              SeascountOperators,     1,   CDI_BOTH,  1,  1 },
   { Seaspctl,       SeaspctlHelp,      SeaspctlOperators,      1,   CDI_REAL,  3,  1 },
   { Seasstat,       SeasstatHelp,      SeasstatOperators,      1,   CDI_REAL,  1,  1 },