[Bug] r.stats.quantile produces segfault

[tpilz]

Describe the bug
r.stats.quantile terminates with a segfault during quantile calculation.

To Reproduce
Using North Carolina sample dataset: r.stats.quantile base=zipcodes cover=elevation quantiles=3 -p

Screenshots
Valgrind gives:

valgrind -v --leak-check=yes r.stats.quantile base=zipcodes cover=elevation quantiles=3 -p
[...]
==1793== Invalid read of size 4
==1793==    at 0x402924: get_quantile (main.c:95)
==1793==    by 0x402924: compute_quantiles (main.c:427)
==1793==    by 0x402924: main (main.c:695)
==1793==  Address 0xc is not stack'd, malloc'd or (recently) free'd
[...]

System description (please complete the following information):

• Operating System: Linux openSUSE 15.3
• GRASS GIS version:

g.version -rge
version=8.0.1
date=2022
revision=bf433781b
build_date=2022-03-11
build_platform=x86_64-pc-linux-gnu
build_off_t_size=8
libgis_revision=bf433781b
libgis_date=2022-02-24T20:37:11+00:00
proj=8.2.1
gdal=3.4.1
geos=3.10.2
sqlite=3.36.

[tpilz]

... maybe caused by c700a8a? Just noticed that with GRASS 8.0.0 it works fine.

[nilason]

quants is array of size num_quants:

quants = G_calloc(num_quants, sizeof(DCELL));

grass/raster/r.stats.quantile/main.c

Lines 90 to 95 in fdff46c

	if (n >= num_quants) {
	/* stop condition for initialize_bins() */
	return (double)bc->total + bc->total;
	}
	rnk = quants[n] * (bc->total - 1);

I believe the fix for this is to change line 90 to:

    if (n > num_quants) {

[nilason]

I believe the fix for this is to change line 90 to:

    if (n > num_quants) {

Sorry, that can't be right.

[nilason]

crash log:

0   r.stats.quantile              	       0x100082d4c compute_quantiles + 288 (main.c:431)
1   r.stats.quantile              	       0x100082d44 compute_quantiles + 280 (main.c:427)
2   r.stats.quantile              	       0x100082104 main + 1916 (main.c:695)
3   dyld                          	       0x1004190f4 start + 520

[nilason]

Some debug info:

* thread #1, queue = 'com.apple.main-thread', stop reason = EXC_BAD_ACCESS (code=1, address=0x0)
    frame #0: 0x0000000100002d4c r.stats.quantile`compute_quantiles at main.c:431:13
   428 			double k, v;
   429 			int i0, i1;
   430 	
-> 431 			while (b->origin + b->count < next)
   432 			    b++;
   433 	
   434 			k = next - b->origin;

│ (int) cat = 92                                                                                                                  ││                               │
│ ◆─(basecat *) bc = 0x0000000103815a80                                                                                           ││                               │
│ ├─◆─(size_t *) slots = 0x0000600000004410                                                                                       ││                               │
│ ├─(size_t) total = 29535                                                                                                        ││                               │
│ ├─(size_t) num_values = 0                                                                                                       ││                               │
│ ├─(DCELL) min = 89.626907348632813                                                                                              ││                               │
│ ├─(DCELL) max = 134.30082702636719                                                                                              ││                               │
│ ├─(DCELL) slot_size = +Inf                                                                                                      ││                               │
│ ├─(int) num_slots = 0                                                                                                           ││                               │
│ ├─◆─(unsigned char *) slot_bins = 0x0000000000000000                                                                            ││                               │
│ ├─(int) num_bins_alloc = 0                                                                                                      ││                               │
│ ├─(int) num_bins_used = 0                                                                                                       ││                               │
│ ├─◆─(bin *) bins = 0x0000000000000000                                                                                           ││                               │
│ ├─◆─(DCELL *) values = 0x0000000000000000                                                                                       ││                               │
│ └─◆─(DCELL *) quants = 0x000060000000c000                                                                                       ││                               │
│ (int) quant = 0                                                                                                                 ││                               │
│ ◆─(bin *) b = 0x0000000000000000                                                                                                ││                               │
│ (double) next = 9844.6666666666661                                                                                              ││                               │
│ (double) k = 3.0493302575189637E-314                                                                                            ││                               │
│ (double) v = 4.9406564584124654E-322                                                                                            ││                               │
│ (int) i0 = 0                                                                                                                    ││                               │
│ (int) i1 = 100                                                                                                                  ││                               │
│ (DCELL) min = 55.578792572021484                                                                                                ││                               │
│ (int) cols = 179                                                                                                                ││                               │
│ (int) rows = 165                                                                                                                ││                               │
│ (int) num_cats = 100                                                                                                            ││                               │
│ (int) num_slots = 1000                                                                                                          ││                               │
│ (CELL) cmax = 27610                                                                                                             ││                               │
│ (CELL) cmin = 27511                                                                                                             ││                               │
│ ◆─(DCELL *) quants = 0x0000600000004190                                                                                         ││                               │
│ ◆─(basecat *) basecats = 0x0000000103813800                                                                                     ││                               │
│ (int) num_quants = 2                                                                                                            ││                               │
│ (DCELL) max = 156.32986450195313                                                                                                ││                               │

[ecodiv]

I can confirm this bug for GRASS GIS 8.1.dev (9c4ded9) on Ubuntu 20.04
In GRASS 7.8.6 (2021) on Ubuntu 20.04 the function works

[ecodiv]

@wenzeslaus is this a difficult one to fix? Waiting till 8.2.0 means this function won't be available for a while. Don't want to push, just to know whether I and my students are better off reverting to an older version of grass for the time being.

[nilason]

I tracked the immediate cause of the segfault. In get_slot_counts() the bc->num_slots is set to 0, hence:

grass/raster/r.stats.quantile/main.c

Line 170 in b95f3ee

bc->slots = G_calloc(bc->num_slots, sizeof(unsigned int));

allocates zero objects for bc->slots.

and bc->bins never gets allocated:

grass/raster/r.stats.quantile/main.c

Lines 218 to 223 in b95f3ee

	if (bc->num_slots == 0)
	continue;
	bc->num_bins_alloc = num_quants * 2;
	bc->bins = G_calloc(bc->num_bins_alloc, sizeof(struct bin));
	bc->slot_bins = G_calloc(bc->num_slots, sizeof(unsigned char));

This is faulted in compute_quantiles() (line 431) :

grass/raster/r.stats.quantile/main.c

Lines 424 to 433 in b95f3ee

	struct bin *b = &bc->bins[0];
	for (quant = 0; quant < num_quants; quant++) {
	double next = get_quantile(bc, quant);
	double k, v;
	int i0, i1;
	while (b->origin + b->count < next)
	b++;

I don't know for sure what would be the correct fix for this, maybe @metzm have?

[wenzeslaus]

@ecodiv I marked this as 8.2.0 because I though c700a8a (#2108 fix quantile algorithm) is only in the main branch and from the comments it seemed it is the likely cause and 8.0.0 is fine. However, now I see it was backported to the 8.0 branch in 44949fa.

As for the timelines, the planned release dates for 8.2.0 and 8.0.2 are actually similar (see https://github.com/OSGeo/grass/milestones).

Notably, r.stats.quantile doesn't have any tests, so it is not checked by the automated checks in CI. (same for and r.quantile but that one was tested manually according to #2108; r.neighbors has tests)

[ecodiv]

Thanks for the feedback @wenzeslaus. I wish I knew how to write such tests. On a different note, not sure it is useful, but tested it on different versions in Windows.

Version 8.01 on Windows, installed using OSGeo4W runs without error (also when running with debug level 5), but also without showing the results.

C:\Windows>r.stats.quantile base=zipcodes cover=elevation quantiles=3 -p
Computing histograms
 100%
 100%
Computing bins
Binning data
 100%
Sorting bins
 100%
Computing quantiles

Version 8.1.dev on Windows (code revision d21dd0a), installed using standalone installer runs (somewhat unexpectedly given the above) fine.

C:\Windows>r.stats.quantile base=zipcodes cover=elevation quantiles=3 -p
Computing histograms
 100%
 100%
Computing bins
Binning data
 100%
Sorting bins
 100%
Computing quantiles
D3/5: G_option_to_separator(): key = separator -> sep = ':'
Printing quantiles
27511:0:33.333333:134.703466
27511:1:66.666667:143.974879
27513:0:33.333333:140.805303
27513:1:66.666667:146.319600
27518:0:33.333333:115.158768
27518:1:66.666667:129.913274
27529:0:33.333333:96.420019
27529:1:66.666667:103.635132
27539:0:33.333333:123.106649
27539:1:66.666667:133.501918
27601:0:33.333333:89.867775
27601:1:66.666667:98.400887
27603:0:33.333333:91.626602
27603:1:66.666667:103.526230
27604:0:33.333333:75.851888
27604:1:66.666667:87.067993
27605:0:33.333333:102.307660
27605:1:66.666667:112.607625
27606:0:33.333333:108.138123
27606:1:66.666667:124.474116
27607:0:33.333333:124.764328
27607:1:66.666667:135.923869
27608:0:33.333333:91.529701
27608:1:66.666667:104.507629
27610:0:33.333333:81.774063
27610:1:66.666667:92.784810

[metzm]

I think I found the change in #2108 that causes this segfault, it could be fixed with a single line. A corresponding PR should come today or tomorrow.