funny axis ranges; GPretty(.) vs. pretty(.) and all that...

[Martin Maechler]

[This is something like a bug report;
 maybe somewhat longish & technical ..]

As an introduction, just try the following code
(it should work both in R and S).
I think it screws up the postscript() driver both for S and R, but this is
not the issue here.


is.R <- function() { ## returns  'TRUE'  iff we are using 'R'
  exists("version") && !is.null(vl <- version$language) && vl == "R"
}

p.axis.range <- function(k.set = c(-100,-9:-6,-3,3,5,20),
                         y.k.add = 1)
{
  ## Purpose: Plot axis scaling & labeling with funny range(.)
  ## -------------------------------------------------------------------------
  ## Arguments: k.set: (0, 10^k  will be the ranges
  ## -------------------------------------------------------------------------
  ## Author: Martin Maechler, Date: 12 Jan 98, 18:15

  lk <- length(k.set); slk <- ceiling(sqrt(lk))
  op <- par(mfrow=c(slk,slk), mar = .2+c(3,2,1,1), oma=rep(1,4))
  on.exit(par(op))
  ##        oma does not really work
  if(is.R()) { colx <- rainbow(lk); coly <- terrain.colors(lk)
             } else colx <- lk + (coly <- 1:lk)
  for(i in 1:lk) {
    k <- k.set[i]
    if(is.R()) {
      plot.new()
      plot.window(xlim= c(0,10^ k),
                  ylim= c(0,10^(k+y.k.add)),"")
      box(lty='dashed')
    } else { ## S4 and S-plus
      plot.default(0,0,xlim= c(0,10^ k), ylim= c(0,10^(k+1)),
                   type = 'n', xlab='', ylab='')
    }
    cat("\n--------------------------------\n\n",
        "xlim= c(0, 10^", k,")    ==> par('usr')[1:2]=", par("usr")[1:2],"\n")
    par(col=colx[i]); if(is.R()) axis(1, col.axis= colx[i]) else axis(1)
    mtext(paste("xlim= (0, 10^",k,")",sep=""),	       line=-4)
    par(col=coly[i]); if(is.R()) axis(2, col.axis= coly[i]) else axis(2)
    mtext(paste("ylim= 10^",y.k.add," * xlim",sep=""), line=-6)
    par(col=if(is.R())'black' else 1)
    mtext(paste("par(\"usr\")= 10 ^", k," *"), line= -9, cex= .4)
    mtext(paste(format(par("usr")*10^-k),collapse=",  "), line= -11, cex= .3)
  }
  invisible(NULL)
}

p.axis.range()

##-----------------------------------------------------------------------------

##>> Now, what's going on? [in R, don't know for S; this all hidden in 
##>> 			calls such as 
##>> 			.Internal(.Cur.picture(), "S_cur_pic", T, 1) ]

## The main problem is  GPretty() in  $RHOME/src/graphics/graphics.c !

##R: plot.default()  -->   plot.new();  plot.window(xlim, ylim, log)

##C: 1:  do_plot_window()   in  $RHOME/src/main/plot.c

##		checks its arguments and then
##		GCheckState(); GScale(xmin, xmax, 1); GScale(ymin, ymax, 2);

##   2:  GScale(.)	    in $RHOME/src/graphics/graphics.c

##		checks args; sets up 'xmin', 'xmax' depending on 'log',..axt'
##		then  GPretty(&xmin, &xmax, &n);
##		      ~~~~~~~
##		and puts its result  in  GP->[x/y]axp  -- i.e. par("xaxp") or y.

##   3:  GPretty(.)   [also in ...graphics.c  is SIMILAR to but different from
##		       		$RHOME/src/appl/pretty.c
##		       which is called from
##		       		$RHOME/src/library/base/R/pretty

##	uses  FLT_EPSILON  (/usr/include/float.h)  in 3 places, where the above
##	pretty(.)  uses FLT_EPSILON only once ...
##
## Why use it at all (the way it is done)?
##
## S-plus's  pretty(.) is quite different;   also the  p.axis.range(.) below:
## ---> Comparison of  pretty(.) and Gpretty(.)  : see below.

So it's the  FLT_EPSILON which is about 1e-7,
and leads to the funny  -1e-7, +1e-7  axis ranges.

%%>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>
%%>=>  Questions/Remarks						     %%
%%>=>  									     %%
%%>=>  1) Why (and how if..) should GPretty(.) be different from pretty?     %%
%%>=>  									     %%
%%>=>  2) Couldn't we at least replace FLT_EPSILON  by DOUBLE_EPSILON ?	     %%
%%>=>  									     %%
%%>=>  3) axis(.) still has another bug: it draws the line and writes labels %%
%%>=>  	   OUTSIDE the "usr" coordinates.				     %%
%%>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>



###------ Try the following code to find out about the differences
###------ of the 2 versions of  pretty(.) : ------------------------

Cpretty <- Cpretty1 <- function (x, n = 5)
  unlist(.C("pretty", l = min(x), u = max(x), n = as.integer(n)))
Cpretty2 <- function (x, n = 5) {
  ##-- return what  .C("pretty... doe
  n <- length(p <- pretty(x,n))
  c(l=p[1], u=p[n], n=n-1)
}

Cpretty(pi)
all(Cpretty1(pi) == Cpretty2(pi))
all(Cpretty1(c(0,pi)) == Cpretty2(c(0,pi)))
for(i in 1:100) {
  x_c(runif(1),pi)
  for(k in -14:14) {
    ok <- all(Cpretty1(x*10^k) == Cpretty2(x*10^k))
    if(!ok) cat("NOT ok:  x=",deparse(x)," k=",k,"\n")
  }
}

Gpretty <- function (x, n = 5)
{
  plot.new(); op <- par(); on.exit(par(op))
  par(lab=c(n, par("lab")[-1]))
  plot.window(xlim=c(min(x),max(x)), ylim=0:1, "")
  par("xaxp")
}

##--- Gpretty() and Cpretty() are DIFFERENT much more than the same !!! -------

vcat <- function(vec, sep = " ", dig=2)
  paste(formatC(vec, dig=dig), collapse = sep)

show.ok <- TRUE
show.nok <- TRUE
show.nok <- !show.ok
for(i in 1:100) {
  x <<- sort(c(runif(1),4*rexp(1)))
  n <- max(1, rpois(1,lam=3))
  cat("\n")
  for(k in -14:14) { cat(".")
    ok <- all((cp <- Cpretty(x*10^k,n=n)) == (gp <- Gpretty(x*10^k,n=n)))
    if(!ok && show.nok) cat("NOT the same:  n=",n,"x=",vcat(x)," k=",k,
                "; Cp=",vcat(cp), " Gp=",vcat(gp),"\n")
    else if(ok && show.ok) cat("the same:  n=",n,"x=",deparse(x)," k=",k,
                              ";  Cp=Gp=",vcat(cp),"\n")
  }
}
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._