Merge branch 'qfix-v2' into main

src/cross-level/stats-cl.c

@@ -33,6 +33,8 @@ This file contains functions for statistics
 
 #include "stats-cl.h"
 
+#include <gsl/gsl_sort.h>         // sort data
+#include <gsl/gsl_statistics.h>   // statistical functions
 
 void quantile_swap(float *x, int from, int to);
 int comp(const void *a, const void *b);
@@ -657,66 +659,21 @@ float abst, tsq, p, z;
 
 /** Quantile
 +++ This function computes a quantile of an array. The quick select algo-
-+++ rithm is used for this purpose. Caution: the array will be screwed up.
++++ rithm is used for to sort the array. Caution: the array will be screwed up.
 +++ Copy the array before calling the quantile function.
 --- x:      array
 --- n:      length of array
 --- p:      probability [0,1]
 +++ Return: quantile
 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++**/
-float quantile(float *x, int n, float p){
-int left = 0, right = n - 1, r, w, k;
-float piv;
+float quantile(double *x, int n, float p){
+double q;
 
+  gsl_sort (x, 1, n);
 
-  // compute position from probability
-  k = (int) (n-1)/(1/p);
+  q = gsl_stats_quantile_from_sorted_data(x, 1, n, p);
 
-  // do until left and right converge at k
-  while (left < right){
-
-    r = left ;        // reader
-    w = right;        // writer
-    piv = x[(r+w)/2]; // pivot
-
-    // do until reader and writer are equal
-    while (r < w){
-
-      if (x[r] >= piv){
-        // if larger than pivot, put value to the end, decrease writer
-        quantile_swap(x, r, w);
-        w--;
-      } else {
-        // if smaller than pivot, skip, increase reader
-        r++;
-      }
-    }
-
-    // if reader was increased, decrese reader
-    if (x[r] > piv) r--;
-
-    // reader is on the end of the first k elements
-    if (k <= r){
-      right = r;
-    } else {
-      left = r + 1;
-    }
-
-  }
-
-  return x[k];
-}
-
-
-/** Function for swapping two array elements
-+++ This function is used in the quick select algorithm for quantiles.
-+++--------------------------------------------------------------------**/
-void quantile_swap(float *x, int from, int to){
-float tmp = x[from];
-
-  x[from] = x[to]; x[to] = tmp;
-
-  return;
+  return (float)q;
 }
 
 

src/cross-level/stats-cl.h

@@ -67,7 +67,7 @@ float tscore_Norm_z(float p);
 float tscore_Hills_inv_t(float p, int idf);
 float tscore_T_z(float t, int df);
 float tscore_T_p(float t, int df);
-float quantile(float *x, int n, float p);
+float quantile(double *x, int n, float p);
 int mode(int *x, int n);
 int n_uniq(int *x, int n);
 int **histogram(int *x, int n, int *n_uniq);

src/higher-level/cso-hl.c

@@ -214,7 +214,7 @@ short q25_, q75_;
 double mean, var;
 double skew, kurt;
 double skewscaled, kurtscaled;
-float *q_array = NULL;
+double *q_array = NULL; // need to be double for GSL quantile function
 bool alloc_q_array = false;
 
 
@@ -292,7 +292,7 @@ bool alloc_q_array = false;
   #pragma omp parallel private(o,t,w,minimum,maximum,q,q_array,mean,var,skew,kurt,n,k,skewscaled,kurtscaled,q25_,q75_,d_ce,ce,ce_left) shared(mask_,cs,nc,nw,nt,nodata,alloc_q_array,phl,t0,t1,nprod,ard) default(none)
   {
 
-    if (alloc_q_array) alloc((void**)&q_array, nt+1, sizeof(float));
+    if (alloc_q_array) alloc((void**)&q_array, nt+1, sizeof(double));
 
     #pragma omp for
     for (p=0; p<nc; p++){

src/higher-level/fold-hl.c

@@ -53,7 +53,7 @@ short minimum, maximum;
 double mean, var;
 double skew, kurt;
 double skewscaled, kurtscaled;
-float *q_array = NULL;
+double *q_array = NULL; // need to be double for GSL quantile function
 bool alloc_q_array = false;
 
 
@@ -67,7 +67,7 @@ bool alloc_q_array = false;
   {
 
     // initialize _STAts
-    if (alloc_q_array) alloc((void**)&q_array, ni, sizeof(float));
+    if (alloc_q_array) alloc((void**)&q_array, ni, sizeof(double));
 
     #pragma omp for
     for (p=0; p<nc; p++){

src/higher-level/stm-hl.c

@@ -48,7 +48,7 @@ short q25_, q75_;
 double mean, var;
 double skew, kurt;
 double skewscaled, kurtscaled;
-float *q_array = NULL;
+double *q_array = NULL; // need to be double for GSL quantile function
 bool alloc_q_array = false;
 
 
@@ -63,7 +63,7 @@ bool alloc_q_array = false;
   {
 
     // initialize stats
-    if (alloc_q_array) alloc((void**)&q_array, ni, sizeof(float));
+    if (alloc_q_array) alloc((void**)&q_array, ni, sizeof(double));
 
     #pragma omp for
     for (p=0; p<nc; p++){

src/lower-level/cloud-ll.c

@@ -431,8 +431,8 @@ float lowtemp = 0.0, hightemp = 0.0;
 float A, B;
 float tempprob;
 float *LANDPROB = NULL;
-float *CLEARLPROB = NULL;
-float *CLEARTEMP = NULL;
+double *CLEARLPROB = NULL; // need to be double for GSL quantile function
+double *CLEARTEMP = NULL;  // need to be double for GSL quantile function
 
 
   #ifdef FORCE_CLOCK
@@ -444,7 +444,7 @@ float *CLEARTEMP = NULL;
 
   if ((100.0*nland/(float)npix) >= 0.1){
 
-    alloc((void**)&CLEARTEMP, nland, sizeof(float));
+    alloc((void**)&CLEARTEMP, nland, sizeof(double));
 
     for (p=0, k=0; p<nc; p++){
       if (lnd_[p]) CLEARTEMP[k++] = temp_[p];
@@ -455,7 +455,7 @@ float *CLEARTEMP = NULL;
 
   } else {
 
-    alloc((void**)&CLEARTEMP, nclear, sizeof(float));
+    alloc((void**)&CLEARTEMP, nclear, sizeof(double));
 
     for (p=0, k=0; p<nc; p++){
       if (clr_[p]) CLEARTEMP[k++] = temp_[p];
@@ -492,13 +492,14 @@ float *CLEARTEMP = NULL;
 
 
   // dynamic cloud detection threshold over land
-  alloc((void**)&CLEARLPROB, nland, sizeof(float));
+  alloc((void**)&CLEARLPROB, nland, sizeof(double));
 
   for (p=0, k=0; p<nc; p++){
     if (lnd_[p]) CLEARLPROB[k++] = LANDPROB[p];
   }
 
   lclr_max = quantile(CLEARLPROB, nland, hi) + cldprob;
+
   free((void*)CLEARLPROB);
 
   #ifdef FORCE_CLOCK
@@ -531,8 +532,8 @@ float wtrtemp = 0.0;
 int nCLEARWTR = 0;
 bool  *CLEARWTR = NULL;
 float *WATERPROB = NULL;
-float *CLEARWPROB = NULL;
-float *CLEARTEMP = NULL;
+double *CLEARWPROB = NULL; // need to be double for GSL quantile function
+double *CLEARTEMP = NULL;  // need to be double for GSL quantile function
 
 
   #ifdef FORCE_CLOCK
@@ -551,7 +552,7 @@ float *CLEARTEMP = NULL;
   // test if there is clear water, if not skip and give water prob = -1
   if (nCLEARWTR > 0){
 
-    alloc((void**)&CLEARTEMP, nCLEARWTR, sizeof(float));
+    alloc((void**)&CLEARTEMP, nCLEARWTR, sizeof(double));
 
     for (p=0, k=0; p<nc; p++){
       if (CLEARWTR[p])CLEARTEMP[k++] = temp_[p];
@@ -575,7 +576,7 @@ float *CLEARTEMP = NULL;
     }
 
     // dynamic cloud detection threshold over water
-    alloc((void**)&CLEARWPROB, nCLEARWTR, sizeof(float));
+    alloc((void**)&CLEARWPROB, nCLEARWTR, sizeof(double));
 
     for (p=0, k=0; p<nc; p++){
       if (CLEARWTR[p]) CLEARWPROB[k++] = WATERPROB[p];
@@ -770,7 +771,7 @@ short *spr_  = NULL;
 short *toa_    = NULL;
 short *mask_   = NULL;
 short *marker_ = NULL;
-float *clear_  = NULL;
+double *clear_  = NULL; // need to be double for GSL quantile function
 ushort *dist_;
 char domains[2][NPOW_10] = { "NIR", "SWIR1" };
 
@@ -824,7 +825,7 @@ char domains[2][NPOW_10] = { "NIR", "SWIR1" };
 
       memmove(mask_, toa_, nc*sizeof(short));
 
-      alloc((void**)&clear_,  nland, sizeof(float));
+      alloc((void**)&clear_,  nland, sizeof(double));
       for (p=0, k=0; p<nc; p++){
         if (lnd_[p]) clear_[k++] = mask_[p];
       }
@@ -1408,7 +1409,7 @@ int i, j, p, nx, ny, nc, nf, ne, g, k;
 int skip, influence = 8;
 int nobj, id, size_max = 0, size;
 int *P = NULL, *best_P = NULL;
-float *qtemp = NULL;
+double *qtemp = NULL; // need to be double for GSL quantile function
 float res, radius, core, basetemp;
 float base_min, base_max, base_step, base;
 float best_match, match, shadow, total;
@@ -1546,7 +1547,7 @@ short  *temp_      = NULL;
   #pragma omp parallel private(k, g, x, y, p, size, radius, core, basetemp, base_min, base_max, base, height, shadow, total, match, best_match, qtemp, P, best_P) shared(nx, ny, res, nobj, influence, lowtemp, hightemp, dlapse, rlapse, wlapse, base_step, size_max, spr_, shdprob, cld_, slp_, atc, sun_, view_, QAI, CCL, shd_, array_x, array_y, array_temp, SIZE, temp_) default(none) 
   {
 
-    if (temp_ != NULL) alloc((void**)&qtemp, size_max, sizeof(float));
+    if (temp_ != NULL) alloc((void**)&qtemp, size_max, sizeof(double));
     alloc((void**)&P,       size_max, sizeof(int));
     alloc((void**)&best_P,  size_max, sizeof(int));