[Rd] Minor bug with stats::isoreg

[Martin Maechler]

>>>>> Ivan Krylov 
>>>>>     on Thu, 28 Sep 2023 00:59:57 +0300 writes:

    > В Wed, 27 Sep 2023 13:49:58 -0700 Travers Ching
    > <traversc using gmail.com> пишет:

    >> Calling isoreg with an Inf value causes a segmentation
    >> fault, tested on R 4.3.1 and R 4.2. A reproducible
    >> example is: `isoreg(c(0,Inf))`

    > Indeed, the code in src/library/stats/src/isoreg.c
    > contains the following loop:

>     do {
> 	slope = R_PosInf;
> 	for (i = known + 1; i <= n; i++) {
> 	    tmp = (REAL(yc)[i] - REAL(yc)[known]) / (i - known);
>             // if `tmp` becomes +Inf or NaN...
>             // or both `tmp` and `slope` become -Inf...
> 	    if (tmp < slope) { // <-- then this is false
> 		slope = tmp;
> 		ip = i; // <-- so this assignment never happens
> 	    }
> 	}/* tmp := max{i= kn+1,.., n} slope(p[kn] -> p[i])  and
> 	  *  ip = argmax{...}... */
> 	INTEGER(iKnots)[n_ip++] = ip; // <-- heap overflow and crash // ...
>     } while ((known = ip) < n); // <-- this loop never terminates

> I'm not quite sure how to fix this. Checking for tmp <= slope would
> have been a one-character patch, but it changes the reference outputs
> and doesn't handle isnan(tmp), so it's probably not correct. The
> INTEGER(iKnots)[n_ip++] = ip; assignment should only be reached in case
> of knots, but since the `ip` index never progresses past the
> +/-infinity, the knot condition is triggered repeatedly.

> Least squares methods don't handle infinities well anyway, so maybe
> it's best to put the check in the R function instead:

> --- src/library/stats/R/isoreg.R	(revision 85226)
> +++ src/library/stats/R/isoreg.R	(working copy)
> @@ -22,8 +22,8 @@
>  {
>      xy <- xy.coords(x,y)
>      x <- xy$x
> -    if(anyNA(x) || anyNA(xy$y))
> -	stop("missing values not allowed")
> +    if(!all(is.finite(x)) || !all(is.finite(xy$y)))
> +	stop("missing and infinite values not allowed")
>      isOrd <- ((!is.null(xy$xlab) && xy$xlab == "Index")
>                || !is.unsorted(x, strictly = TRUE))
>      if(!isOrd) {

> -- 
> Best regards,
> Ivan


The above would not even be sufficient: 
It's the sum(y) really, because internally
  yc <- cumsum(c(0,y)) and actually  diff(yc)  is used
where you get to  Inf - Inf ==> NaN

> isoreg(c(5, 9, 1:2, 7e308, 5:8, 3, 8)))

 *** caught segfault ***
address 0x7e48000, cause 'memory not mapped'
/u/maechler/bin/R_arg: Zeile 160: 873336 Speicherzugriffsfehler  (Speicherabzug geschrieben) $exe $@

Also, the C code still does not work for long vectors,
so I want to change the C code anyway.


In any case:
  Thank you, Travers, Ben, and Ivan, for reporting and addressing
  the issue!


------

There is an interesting point here though:

For dealing with +/- Inf, we used to follow the following idea
in R quite keenly (and sometimes extremely): 

If 'Inf' leads a computation to "fail" (NB:  1/Inf  |-->  0 does *not* fail)
try to see what the mathematical *or* computational limit
 x --> Inf would be.
If that is easily defined, we use that.

So, often as a first step, look at what happens if you replace
Inf by 1e100  (and then also what happens if you are finite but
*close* to Inf, i.e. the 7e308 above).

Now here, at least in some cases, such a limit cases are clearly
detectable, e.g., when you let  y[2] --->  -Inf here

> n <- length(y0 <- c(5, 9, 1:2, 5:8, 3, 8))
> y2s <- c(10:0, -10, -20, -1000, -1e4, -1e10, -1e100, -1e200, -1e300)
> iSet <- vapply(y2s, function(y2) isoreg({y <- y0; y[2] <- y2; y})$yf, numeric(n))
> t(iSet) # *does* change as function of y2 *but* predictably
              [,1]         [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
 [1,]   4.5000e+00   4.5000e+00 4.50 4.50    5    6    6    6    6     8
 [2,]   4.2500e+00   4.2500e+00 4.25 4.25    5    6    6    6    6     8
 [3,]   4.0000e+00   4.0000e+00 4.00 4.00    5    6    6    6    6     8
 [4,]   3.7500e+00   3.7500e+00 3.75 3.75    5    6    6    6    6     8
 [5,]   3.5000e+00   3.5000e+00 3.50 3.50    5    6    6    6    6     8
 [6,]   3.2500e+00   3.2500e+00 3.25 3.25    5    6    6    6    6     8
 [7,]   3.0000e+00   3.0000e+00 3.00 3.00    5    6    6    6    6     8
 [8,]   2.7500e+00   2.7500e+00 2.75 2.75    5    6    6    6    6     8
 [9,]   2.5000e+00   2.5000e+00 2.50 2.50    5    6    6    6    6     8
[10,]   2.2500e+00   2.2500e+00 2.25 2.25    5    6    6    6    6     8
[11,]   2.0000e+00   2.0000e+00 2.00 2.00    5    6    6    6    6     8
[12,]  -2.5000e+00  -2.5000e+00 1.00 2.00    5    6    6    6    6     8
[13,]  -7.5000e+00  -7.5000e+00 1.00 2.00    5    6    6    6    6     8
[14,]  -4.9750e+02  -4.9750e+02 1.00 2.00    5    6    6    6    6     8
[15,]  -4.9975e+03  -4.9975e+03 1.00 2.00    5    6    6    6    6     8
[16,]  -5.0000e+09  -5.0000e+09 1.00 2.00    5    6    6    6    6     8
[17,]  -5.0000e+99  -5.0000e+99 0.00 0.00    0    0    0    0    0     0
[18,] -5.0000e+199 -5.0000e+199 0.00 0.00    0    0    0    0    0     0
[19,] -5.0000e+299 -5.0000e+299 0.00 0.00    0    0    0    0    0     0
>

so one could say that ideally,

      isoreg(c(5, -Inf, 1:2, 5:8, 3, 8))

should produce fitted values

             c(-Inf, -Inf, 0, 0, ..., 0)

and if someone has a +/-  elegant implementation
we could again allow  +/-Inf entries in  isoreg(), at least when
the Inf's have all the same sign.

Martin