1#### misc-goodies.R 2#### ~~~~~~~~~~~~~~ SfS - R - goodies that are NOT in 3#### "/u/sfs/R/SfS/R/u.goodies.R" 4#### "/u/sfs/R/SfS/R/p.goodies.R" 5 6###--- Original: From 'S' in /u/sfs/S/misc-goodies.S 7###--- ======== But start doing *less* here ! 8 9###================================================================== 10### Functions <<<<<<<< Please use a few subsections like "Plotting"... 11###================================================================== 12 13### ___Note___ we have some of these headers __MESS__ 14### But we leave it because of RCS {rather dismantle everything into 4-6 pieces 15 16##-#### Vector, Matrix (or higher Array) stuff ######## 17##-### -------------------------------------- ######## 18 19last <- function(x, length.out = 1, na.rm = FALSE) 20{ 21 ## Purpose: last element(s) of a vector 22 ## Author: Werner Stahel, Date: Tue Jan 21 17:29:42 1992 23 ## ---------------------------------------------------------------- 24 ## Arguments: 25 ## x: vector 26 ## length.out: if positive, return the length.out last elements of x, 27 ## if negative, the last length.out elements are dropped 28 ## ---------------------------------------------------------------- 29 if (na.rm) 30 x <- x[!is.na(x)] 31 n <- length(x) 32 x[sign(length.out)*(n-abs(length.out)+1):n] 33} 34 35empty.dimnames <- function(a) 36{ 37 ## 'Remove' all dimension names from an array for compact printing. 38 n <- length(da <- dim(a)) 39 if(n == 0) return(a) 40 dimnames(a) <- lapply(1:n, function(i) rep.int("", da[i])) 41 a 42} 43 44 45##-#### Plot / Devices related stuff ######## 46##-### ----------------------------- ######## 47##-### >>>>> "p.goodies.S" or "ps.goodies.S" ######## 48 49errbar <- function(x, y, yplus, yminus, cap = 0.015, 50 ylim = range(y, yplus, yminus), 51 xlab = deparse(substitute(x)), 52 ylab = deparse(substitute(y)), ... ) 53{ 54 ## Purpose: Makes a plot with error bars 55 ## Authors: Charles Geyer, Statistics, U. Chicago, geyer@galton.uchicago.edu 56 ## Martin Maechler, Date: 11 Apr 91 and Mar 27 1992, 12:32 57 ## ---------------------------------------------------------------- 58 ## Arguments: --- see help(..) page ---> ?errbar 59 ## ----------------------------------------======= 60 61 plot( x, y, ylim=ylim, xlab=xlab, ylab=ylab, ... ) 62 xcoord <- par()$usr[1:2] 63 segments( x, yminus, x, yplus ) 64 smidge <- cap * ( xcoord[2] - xcoord[1] ) / 2 65 segments( x - smidge, yminus, x + smidge, yminus ) 66 segments( x - smidge, yplus, x + smidge, yplus ) 67} 68## C.Monatsname , etc.. sind jetzt bei der zugehoerigen Funktion 69## u.Datumvonheute in /u/sfs/S/u.goodies.S 70 71cum.Vert.funkt <- function(x, Quartile = TRUE, titel = TRUE, Datum = TRUE, 72 rang.axis = n <= 20, xlab = "", main = "", ...) 73{ 74 ## Ziel: Kumulative Verteilung von x aufzeichnen, auf Wunsch auch Median 75 ## und Quartile 76 op <- par(xaxs = "r", yaxs = "r", las = 1)# the default anyway 77 on.exit(par(op)) 78 r <- plotStep(x, xlab = xlab, main = main, ...) 79 #### FIXME : stepfun() / ecdf() instead 80 n <- length(x) 81 if(rang.axis) 82 axis(4, at = (0:n)/n, labels = 0:n, pos = par("usr")[1])#, las = 1) 83 if(titel) mtext("Kumulative Verteilungsfunktion", 3, line = 0.5) 84 if(Quartile) for(i in 1:3) abline(h = i/4, lty = 2) 85 if(Datum) p.datum() 86 invisible(r) 87} 88 89 90## This was "plot.step()" but that's in conflict with S3 methods 91plotStep <- function(ti, y, 92 cad.lag = TRUE, 93 verticals = !cad.lag, 94 left.points = cad.lag, 95 right.points = FALSE, 96 end.points = FALSE, 97 98 add = FALSE, 99 100 pch = par('pch'), 101 xlab = deparse(substitute(ti)), 102 ylab = deparse(substitute(y)), 103 main = NULL, 104 ...) 105 106#####- FIXME ----------- use stepfun(), plot.stepfun() etc !!! ---------------- 107 108{ 109 ## Purpose: plot step-function f(x)= sum{ y[i] * 1_[ t[i-1], t[i] ] (x) } 110 ## ------------------------------------------------------------------------- 111 ## Arguments: for missing 'y', do empirical CDF; ==> ON-LINE Help "?plot.step" 112 ## ------------------------------------------------------------------------- 113 ## Author: Martin Maechler, 1990, U.Washington, Seattle; improved -- Dec.1993 114 ## 115 ## EXAMPLE: ##-- Plot empirical cdf Fn(x) for a small n: 116 ## xx <- runif(20); plot.step(xx); plot.step( xx, cad.lag = F ) 117 ## plot.step( runif(20), add=T, cad.lag=F) 118 xlab 119 ylab 120 if(missing(y)) { 121 if(is.vector(ti) && is.numeric(ti)) { # -- Do empirical CDF -- 122 nt <- length(ti) 123 ti <- sort(ti) 124 dt <- (ti[nt] - ti[1])/20 125 ti <- c(ti[1] - dt, ti, ti[nt] + dt) 126 n <- nt + 1 127 y <- (0:nt)/nt 128 } else { 129 xy <- xy.coords(ti,NULL,NULL,NULL) 130 ti <- c(xy$x[1], xy$x) 131 y <- xy$y 132 n <- length(y) 133 } 134 } else { 135 n <- length(y) 136 if(length(ti) != (n + 1)) stop("length(ti) MUST == length(y) + 1") 137 } 138 if(length(ti) != (n + 1) || length(y) != n) 139 stop("NEVER CALLED! --length(ti) MUST == length(y) + 1") 140 if(missing(main)) main <- deparse(sys.call()) 141 142 n1 <- n+1 143 ##-- horizontal segments: 144 if (add) segments(ti[-n1], y, ti[-1], y, ...) 145 else { 146 plot(ti, c(y[1],y), type = 'n', xlab = xlab, ylab = ylab, main = main, ...) 147 segments(ti[-n1], y, ti[-1], y) 148 } 149 if(left.points) points(ti[-n1],y, pch = pch) 150 if(right.points) points(ti[-1], y, pch = pch) 151 ##-- col=0 <==> "erase" : 152 if(! end.points) points(ti[c(1,n1)], y[c(1,n)], pch = pch, col = 0) 153 if(verticals) { 154 if (add) segments(ti[2:n], y[-n], ti[2:n], y[-1], ...) 155 else segments(ti[2:n], y[-n], ti[2:n], y[-1]) 156 } 157 invisible(list(t = ti, y = y)) 158} 159 160histBxp <- 161 function(x, nclass, breaks, probability = FALSE, include.lowest = TRUE, 162 xlab = deparse(substitute(x)), ..., width = 0.2, 163 boxcol = 3, medcol = 2, medlwd = 5, whisklty = 2, staplelty = 1) 164{ 165 ## Purpose: Plot a histogram and a boxplot 166 ## ------------------------------------------------------------------------- 167 ## Arguments: ---> see help(hist.bxp) ! 168 ## ------------------------------------------------------------------------- 169 ## Authors: Christian Keller, Date: 10 Nov 95, (Martin Maechler, Jan 96) 170 ## calls p.hboxp(.) ! 171 172 ## determine the height of the plot 173 h <- 174 if(missing(breaks)){ 175 if(missing(nclass)) 176 hist(x, include.lowest = include.lowest, plot = FALSE) 177 else 178 hist(x, nclass = nclass, 179 include.lowest = include.lowest, plot = FALSE) 180 } 181 else 182 hist(x, breaks = breaks, 183 include.lowest = include.lowest, plot = FALSE) 184 ymax <- max(h$counts) 185 ymin <- - ymax * width # range: (-w,1)*ymax instead of (0,1)*ymax 186 187 ##------- drawing the histogram ------------- 188 hist(x, breaks = h$breaks, probability = probability, 189 include.lowest = include.lowest, plot = TRUE, xlab = xlab, 190 ylim = c(ymin, ymax), axes = FALSE, ...) 191 axis(1) 192 axis(2, at = pretty(c(0,ymax), n = 5), srt = 90) ## ph, 8.5.00: n instead of nint 193 abline(h = 0) # 194 ##-------- drawing the boxplot -------------- 195 196 ##-- scale a range 197 scale.r <- function(x1,x2, fact = 1.1) 198 (x1+x2)/2 + c(-fact,fact) * (x2-x1)/2 199 200 ##-- since 4% extra space above x-axis (just below ymin): 201 ##- cat("par$usr[3:4]:", par("usr")[3:4], 202 ##- " ymin -.04 *(ymax-ymin)",ymin -.04 *(ymax-ymin),"\n") 203 ##-- NOTE: Always have (seemingly): par("usr")[3] == ymin -.04 *(ymax-ymin) 204 205##-O- ORIGINAL VERSION (Keller & Keller) : 206##-O- p.hboxp(x, ymin, -.04 *(ymax-ymin), 207##-O- boxcol=boxcol, medcol=medcol, 208##-O- medlwd=medlwd, whisklty=whisklty, staplelty=staplelty) 209 210 ##---- This is much better for width <=.1 or so... 211 ##-- but you should leave some white space -> scale down 212 ##-- The scaling factor is really a KLUDGE but works for a wide range! 213 p.hboxp(x, scale.r(par("usr")[3], 0, ## ph, 8.5.00: changed f=.9 to f=.8 214 f = .8 - max(0, .15 - width)*(1+(par("mfg")[3] >= 3))), 215 boxcol = boxcol, medcol = medcol, 216 medlwd = medlwd, whisklty = whisklty, staplelty = staplelty) 217} 218 219 220##-#### Print & Strings ######## 221##-### =============== ######## 222 223ccat <- ## character 'concat' 224 function(...) paste0(..., collapse = "") 225vcat <- ## (numeric) vector 'concat' 226 function(vec, sep = " ") paste(vec, collapse = sep) 227 228paste.vec <- function(name, digits = options()$digits) 229{ 230 ## Purpose: Utility for "showing vectors" 231 ## ------------------------------------------------------------------------- 232 ## Example: x <- 1:4; paste.vec(x) ##-> "x = 1 2 3 4" 233 paste(paste(deparse(substitute(name))), "=", 234 paste(format(name, digits = digits), collapse = " ")) 235} 236signi <- function(x, digits = 6) round(x, digits - trunc(log10(abs(x)))) 237 238##' NB: Since ~ R 3.3.0 (May 2016), use base R's "new" strrep(x, times) instead 239 240repChar <- function(char, no) paste(rep.int(char, no), collapse = "") 241## correct, but slower than the next one: 242bl.string <- function(no) repChar(" ", no) 243## faster: 244bl.string <- function(no) sprintf("%*s", no, "") 245 246### symnum : standard R function !! 247 248wrapFormula <- function(f, data, wrapString = "s(*)") 249{ 250 ## Purpose: Mainly: Construct a useful gam() formula from "Y ~ ." 251 ## ---------------------------------------------------------------------- 252 ## Arguments: f : the initial formula; typically something like "Y ~ ." 253 ## data: data.frame to which the formula applies 254 ## ---------------------------------------------------------------------- 255 ## Author: Martin Maechler, Date: 22 May 2007, 18:03 256 257 form <- formula(terms(f, data = data)) 258 if(length(form) != 3) 259 stop("invalid formula; need something like 'Y ~ .'") 260 wrapS <- strsplit(wrapString, "\\*")[[1]] 261 stopifnot(length(wrapS) == 2) 262 cc <- gsub("([^+ ]+)", paste0(wrapS[1], "\\1", wrapS[2]), 263 format(form[[3]])) 264 form[[3]] <- parse(text = cc, srcfile = NULL)[[1]] 265 form 266} 267 268##' Capture Output and print first and last parts, eliding middle parts. 269##' Particularly useful for teaching purposes, and e.g., in Sweave 270##' 271##' @title Capture output and Write / Print First and Last Parts 272##' @param EXPR the (literal) expression the output is to be captured 273##' @param first integer: how many lines should be printed at beginning 274##' @param last integer: how many lines should be printed at the end. 275##' @param middle numeric (or NA logical): 276##' @param i.middle index start of middle part 277##' @param dotdots string to be used for elided lines 278##' @param n.dots number of \code{dotdots} ....{FIXME} 279##' @return return value of \code{\link{capture.output}(EXPR)} 280##' @author Martin Maechler 281## -> ../man/capture-n-write.Rd 282capture.and.write <- function(EXPR, first, last = 2, 283 middle = NA, i.middle, 284 dotdots = " ....... ", n.dots = 2) { 285 co <- capture.output(EXPR) 286 writeLines(head(co, first)) 287 catDots <- function(M) cat(rep.int(paste0(dotdots,"\n"), M), sep="") 288 catDots(n.dots) 289 if(is.numeric(middle)) { 290 stopifnot(length(middle) == 1, middle >= 0, middle == round(middle)) 291 i0 <- first+2 292 if(missing(i.middle)) { 293 i.middle <- max(i0, length(co) %/% 2 - middle %/% 2) 294 } else { ## !missing(i.middle) 295 if(i.middle < i0) 296 stop("'i.middle' is too small, should be at least ", i0) 297 } 298 writeLines(co[i.middle-1 + seq_len(middle)]) 299 catDots(n.dots) 300 } 301 writeLines(tail(co, last)) 302 invisible(co) 303} 304 305 306 307##-#### "Calculus" Mathematical stuff ######## 308##-### ----------------------------- ######## 309 310polyn.eval <- function(coef, x) 311{ 312 ## Purpose: compute coef[1] + coef[2]*x + ... + coef[p+1]* x^p 313 ## if coef is vector, x can be any array; result : of same dim. as x 314 ## if coef is matrix, x must be vector; dim(result) = len(x) * nrow(coef) 315 ## coef = matrix: evaluate SEVERAL polynomials (of same degree) 316 ## ---- contains coefficient-vectors as ROWS ==> coef[,i] <-> x^{i-1} 317 ## Author: Martin Maechler <maechler@stat.math.ethz.ch> 318 if(is.null(dim(coef))) { 319 lc <- length(coef) 320 if (lc == 0) 0 else { 321 r <- coef[lc] 322 if (lc > 1) 323 for (i in (lc-1):1) r <- coef[i] + r*x 324 r 325 } 326 } else { #-- coef is MATRIX -- 327 dc <- dim(coef) 328 lc <- dc[2]; dc <- dc[1] 329 n <- length(x) 330 if (lc == 0) matrix(0, n, dc) else { 331 r <- matrix(coef[,lc], n, dc, byrow = TRUE) 332 if (lc > 1) 333 for (i in (lc-1):1) r <- r*x + matrix(coef[,i], n, dc, byrow = TRUE) 334 r 335 } 336 } 337} 338 339## negative x .. may make sense in some cases,.... but not yet : 340##digitsBase <- function(x, base = 2, ndigits = 1 + floor(log(max(abs(x)),base))) 341digitsBase <- function(x, base = 2, ndigits = 1 + floor(1e-9 + log(max(x,1), base))) 342{ 343 ## Purpose: Give the vector A of the base-_base_ representation of _n_: 344 ## ------- n = sum_{k=0}^M A_{M-k} base ^ k , where M = length(a) - 1 345 ## Value: MATRIX M where M[,i] corresponds to x[i] 346 ## Author: Martin Maechler, Date: Dec 4, 1991 347 ## ---------------------------------------------------------------- 348 ## ----> help(digitsBase) ! 349 ## ------------------------------ 350 if(any(x < 0)) 351 stop("'x' must be non-negative integers") 352 if(any(x != trunc(x))) 353 stop("'x' must be integer-valued") 354 r <- matrix(0, nrow = ndigits, ncol = length(x)) 355 if(ndigits >= 1) for (i in ndigits:1) { 356 r[i,] <- x %% base 357 if (i > 1) x <- x %/% base 358 } 359 class(r) <- "basedInt" 360 attr(r, "base") <- base 361 r 362} 363 364bi2int <- function(xlist, base) 365 vapply(xlist, function(u) polyn.eval(rev(u), base), numeric(1)) 366 367as.intBase <- function(x, base = 2) { 368 xl <- if(is.character(x)) lapply(strsplit(x,""), as.integer) 369 else if(is.numeric(x) && is.matrix(x)) tapply(x, col(x), c) 370 else if(!is.list(x)) 371 stop("'x' must be character, list or a digitsBase() like matrix") 372 bi2int(xl, base) 373} 374 375as.integer.basedInt <- function(x, ...) 376 as.intBase(x, base = attr(x, "base")) 377 378print.basedInt <- function (x, ...) { 379 cat(sprintf("Class 'basedInt'(base = %d) [1:%d]\n", 380 attr(x,"base"), ncol(x))) 381 cx <- x 382 attr(cx,"base") <- NULL 383 print(unclass(cx), ...) 384 invisible(x) 385} 386 387sHalton <- function(n.max, n.min = 1, base = 2, leap = 1) 388{ 389 ## Purpose: Halton sequence H(k,b) for k=n.min:n.max -- for Quasi Monte Carlo 390 ## ---------------------------------------------------------------------- 391 ## Author: Martin Maechler, Date: 29 Jul 2004, 21:34 392 393 stopifnot((leap <- as.integer(leap)) >= 1) 394 ## now do this via digitsBase(), later go directly 395 nd <- as.integer(1 + log(n.max, base)) 396 dB <- digitsBase(if(leap == 1) n.min:n.max else seq(n.min, n.max, by=leap), 397 base = base, ndigits = nd) 398 colSums(dB/base^(nd:1)) 399} 400 401QUnif <- function(n, min = 0, max = 1, n.min = 1, p, leap = 1, silent = FALSE) 402{ 403 ## Purpose: p-dimensional ''Quasi Random'' sample in [min,max]^p 404 ## ---------------------------------------------------------------------- 405 ## Author: Martin Maechler, Date: 29 Jul 2004, 21:43 406 ## Example: plot(QUnif(1000, 2), cex=.6, pch=20, xaxs='i', yaxs='i') 407 stopifnot(1 <= (n <- as.integer(n)), length(n) == 1, 408 1 <= (p <- as.integer(p)), length(p) == 1, 409 length(min) == p || length(min) == 1, 410 length(max) == p || length(max) == 1, 411 1 <= (n.min <- as.integer(n.min)), 412 1 <= (leap <- as.integer(leap)), 413 (n.max <- n.min + (n - 1:1)*leap) < .Machine$integer.max) 414 pr. <- c(2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83, 415 89,97,101,103,107,109,113,127,131,137,139,149,151,157,163,167,173, 416 179,181,191, 193,197,199,211,223,227,229,233,239,241,251,257,263, 417 269,271,277,281,283,293,307,311,313,317,331,337,347,349,353,359, 418 367,373,379,383,389,397,401,409,419,421,431,433,439,443,449,457) 419 if(length(pr.) < p) { 420 if(!silent) 421 message("enlarging internal prime table for \"large\" p=",p) 422 Lp <- log(p) 423 pr. <- primes(p*(Lp + log(Lp))) ## using p_n/n < log n + log log n 424 } 425 pr <- pr.[1:p] 426 if(leap > 1 && any(leap == pr) && length(pr.) >= p+1) # take a non-leap pr 427 pr <- c(pr[leap != pr], pr.[p+1]) 428 max <- rep.int(max, p) 429 min <- rep.int(min, p) 430 dU <- max - min 431 r <- matrix(0., n, p) 432 for(j in 1:p) 433 r[,j] <- min[j] + dU[j] * 434 sHalton(n.max, n.min, base = pr[j], leap = leap) 435 r 436} 437 438 439 440chars8bit <- function(i = 1:255) 441{ 442 ## Purpose: Compute a character vector from its "ASCII" codes. 443 ## We seem to have to use this complicated way thru text and parse. 444 445 ## Author: Martin Maechler, Original date: Wed Dec 4, 1991 446 ## this is an improved version of make.ASCII() from ~/S/Good-string.S ! 447 ## ---------------------------------------------------------------- 448 i <- as.integer(i) 449 if(any(i < 0 | i > 255)) stop("'i' must be in 0:255") 450 if(any(i == 0)) 451 warning("\\000 (= 'nul') is no longer allowed in R strings") 452 i8 <- apply(digitsBase(i, base = 8), 2, paste, collapse="") 453 c8 <- paste0('"\\', i8, '"') 454 eval(parse(text = paste0("c(",paste(c8, collapse=","),")"))) 455} 456 457strcodes <- function(x, table = chars8bit(1:255)) 458{ 459 ## Purpose: R (code) implementation of old S's ichar() 460 ## ---------------------------------------------------------------------- 461 ## Arguments: x: character vector 462 ## ---------------------------------------------------------------------- 463 ## Author: Martin Maechler, Date: 23 Oct 2003, 12:42 464 465 lapply(strsplit(x, ""), match, table = table) 466} 467 468## S-PLUS has AsciiToInt() officially, and ichar() in library(examples): 469AsciiToInt <- ichar <- function(strings) unname(unlist(strcodes(strings))) 470 471helppdf <- function(topic, viewer = getOption("pdfviewer"), quiet = !interactive(), 472 ...) { 473 if(!tryCatch(is.character(topic) && length(topic) == 1L, 474 error = function(e) FALSE)) 475 topic <- deparse1(substitute(topic)) 476 hh <- help(topic, help_type = "pdf", ...) 477 pdfile <- with(attributes(hh), paste(topic, type, sep=".")) 478 ## almost all rendering & pdf creation happens here: 479 print(hh, msg=!quiet)# utils:::print.help_files_with_topic() |--> .show_help_on_topic_offline() 480 if(length(viewer)) 481 system(paste(viewer, pdfile), wait = FALSE) 482 ans <- file.path(getwd(), pdfile) 483 if(quiet) invisible(ans) else ans 484} 485 486 487 488 489##-#### "Miscellaneous" (not any other category) ######## 490##-### ============= ------------------------- ######## 491 492uniqueL <- function(x, isuniq = !duplicated(x), need.sort = is.unsorted(x)) 493{ 494 ## return list(ix, uniq) 495 ## such that all(x == uniq[ix]) and (of course) uniq == x[isuniq] 496 if(need.sort) { 497 xs <- sort(x, index.return = TRUE) 498 ixS <- xs $ ix 499 isuniq <- isuniq[ixS] 500 x <- xs$x 501 } 502 ix <- as.integer(cumsum(isuniq)) 503 if(need.sort) 504 ix <- ix[sort.list(ixS)] 505 list(ix = ix, xU = x[isuniq]) 506} 507 508 509##' Constructor of a "list" (really an environment) of functions (and more) 510##' which all *share* the same environment in which they exist 511##' --> ../man/funEnv.Rd 512##' ~~~~~~~~~~~~~~~~ 513funEnv <- function(..., envir = NULL, parent = parent.frame(), 514 hash = (...length() > 100), size = max(29L, ...length())) { 515 e <- list2env(list(...), envir=envir, parent=parent, hash=hash, size=size) 516 for(n in names(e)) ## iff function or formula, set environment to 'e': 517 if(is.function(e[[n]]) || (is.call(e[[n]]) && 518 inherits(e[[n]], "formula"))) 519 environment(e[[n]]) <- e 520 e 521} 522 523 524is.whole <- function(x, tolerance = sqrt(.Machine$double.eps)) 525{ 526 ## Tests if a numeric scalar (or vector, matrix or array) is a whole 527 ## number; returns an boolean object of the same dimension as x, each entry 528 ## indicating whether the corresponding entry in x is whole. 529 is.whole.scalar <- 530 if (is.integer(x)) { 531 function(x) TRUE 532 } else if (is.numeric(x)) { 533 function(x) isTRUE(all.equal(x, round(x), tolerance = tolerance)) 534 } else if (is.complex(x)) { 535 function(x) 536 isTRUE(all.equal(Re(x), round(Re(x)), tolerance = tolerance)) && 537 isTRUE(all.equal(Im(x), round(Im(x)), tolerance = tolerance)) 538 } else stop("Input must be of type integer, numeric or complex.") 539 540 if (is.null(dim(x))) 541 vapply(x, is.whole.scalar, NA) 542 else 543 apply(x, seq_along(dim(x)), is.whole.scalar) 544} 545 546##' 547##' @title Generate Random Date/Time Sequences 548##' @param n number of entries to generate 549##' @param min, max character strings or \R objects inheriting from \code{"POSIXt"}. 550##' @return vector 551##' @author Martin Maechler 552## 553## __ NOT YET EXPORTED 554## FIXME: consider 'mean = Sys.time(), delta.tim = "1 month"' 555## ----- ==> min = mean - as.difftime(delta.tim), 556## max = mean - as.difftime(delta.tim) 557## now <- Sys.time(); del <- as.difftime(100, units="weeks") 558## rDatetime(100, now-del, now+del) 559rDatetime <- function(n, min = "1900-01-01", max = "2100-12-31") { 560 if(is.character(min) || inherits(min, "POSIXt")) 561 min <- as.POSIXct(min) 562 else stop("'min' must be string (coercable to \"POSIXct\") or \"POSIXt\" object") 563 if(is.character(max) || inherits(max, "POSIXt")) 564 max <- as.POSIXct(max) 565 else stop("'max' must be string (coercable to \"POSIXct\") or \"POSIXt\" object") 566 stopifnot(length(min) == 1, length(max) == 1) 567 structure(runif(n, as.numeric(min), as.numeric(max)), 568 class = c("POSIXct", "POSIXt"), tzone = "") 569} 570 571### 572### autoreg(), mean.cor() etc ... not yet 573### 574### if we take them, use different file !! 575 576 577 578####========== This is from /u/maechler/S/Good.S ============= 579####========== --------------------------------- ============= 580 581##-#### Plot / Devices related stuff ######## 582##-### ----------------------------- ######## 583 584mpl <- function(mat, ...) { 585 matplot(1:nrow(mat), mat, xaxt = 'n',...) 586 if(0 == length(dn <- dimnames(mat)[[1]])) 587 axis(1) else 588 axis(1, at = 1:nrow(mat), labels = dn) 589} 590 591roundfixS <- function(x, method = c("offset-round", "round+fix", "1greedy")) 592{ 593 ## Purpose: y := r2i(x) with integer y *and* sum(y) == sum(x) 594 ## Author: Martin Maechler, 28 Nov 2007 595 n <- length(x) 596 x0 <- floor(x) 597 e <- x - x0 ## == (x %% 1) in [0, 1) 598 S. <- sum(e) 599 stopifnot(all.equal(S., (S <- round(S.)))) 600 method <- match.arg(method) 601 602 ## The problem is equivalent to transforming 603 ## e[] \in [0,1) into f[] \in {0,1}, with sum(e) == sum(f) 604 ## Goal: transform e[] into f[] gradually, by "shifting" mass 605 ## such that the sum() remains constant 606 607 switch(method, 608 "offset-round" = { 609 ## This is going to be equivalent to 610 ## r := round(x + f) with the correct f \in [-1/2, 1/2], or 611 ## r == floor(x + f + 1/2) = floor(x + g), g \in [0, 1] 612 ## 613 ## Need sum(floor(e + g)) = S; 614 ## since sum(floor(e)) == 0, sum(floor(e+1)) == n, 615 ## we just need to floor(.) the S smallest, and ceiling(.) the others 616 if(S > 0) { 617 r <- numeric(n) # all 0; set to 1 those corresponding to large e: 618 r[sort.list(e, decreasing=TRUE)[1:S]] <- 1 619 x0 + r 620 } else x 621 }, ## end{offset-round} 622 623 "round+fix" = { 624 r <- round(e) 625 if((del <- S - sum(r)) != 0) { # need to add +/- 1 to 'del' entries 626 s <- sign(del) ## +1 or -1: add +1 only to r < x entries, 627 aD <- abs(del) ## and -1 only to r > x entries, 628 ## those with the "worst" rounding are made a bit worse 629 if(del > 0) { 630 iCand <- e > r 631 dx <- (e - r)[iCand] # > 0 632 } else { ## del < 0 633 iCand <- e < r 634 dx <- (e - x)[iCand] # > 0 635 } 636 ii <- sort.list(dx, decreasing = TRUE)[1:aD] 637 r[iCand][ii] <- r[iCand][ii] + s 638 } 639 640 return(x0 + r) 641 642 }, ## end{round+fix} 643 644 "1greedy" = { 645 ii <- e != 0 646 while(any(ii)) { 647 ci <- cumsum(ii) # used to revert u[ii] subsetting 648 m <- length(e. <- e[ii]) 649 ie <- sort.list(e.) # both ends are relevant 650 left <- e.[ie[1]] < 1 - e.[ie[m]] 651 iThis <- if(left) 1 else m 652 iother <- if(left) m else 1 653 J <- which.max(ci == ie[iThis]) ## which(.)[1] but faster 654 I <- which.max(ci == ie[iother]) 655 r <- x[J] 656 x[J] <- k <- if(left) floor(r) else ceiling(r) 657 mass <- r - k # if(left) > 0 else < 0 658 if(m <= 2) { # short cut and evade rounding error 659 if(m == 1) { # should happen **rarely** 660 if(!(min(abs(mass), abs(1-mass)) < 1e-10)) 661 warning('m==1 in "1greedy" w/ mass not close to {0,1}') 662 } else { ## m==2 663 x[I] <- round(x[I] + mass) 664 } 665 break ## ii <- FALSE 666 } 667 else { ## m >= 3 668 e[J] <- if(left) 0 else 1 669 ii[J] <- FALSE 670 ## and move it's mass to the other end: 671 e.new <- e[I] + mass 672 if(e.new > 1) 673 stop("e[I] would be > 1 -- internal error") 674 else if(e.new < 0) 675 stop("e[I] would be < 0 -- internal error") 676 x[I] <- x[I] + mass 677 e[I] <- e.new 678 } ## m >= 3 679 } ## end{while} 680 x 681 682 }) # end{switch} 683}## roundfixS 684 685 686seqXtend <- function(x, length., method = c("simple","aim","interpolate"), 687 from = NULL, to = NULL) 688{ 689 ## Purpose: produce a seq(.) covering the range of 'x' and INCLUDING x 690 ## Author: Martin Maechler, Date: 28 Nov 2007 =======> ../man/seqXtend.Rd 691 x <- unique(sort(x)) 692 n <- length(x) 693 method <- match.arg(method) 694 if(length. > n) { 695 if((from_is1 <- is.null(from))) from <- x[1] 696 if((from_isL <- is.null(to))) to <- x[n] 697 if(method == "interpolate") { 698 if(!from_is1) { 699 if(from > x[1]) 700 stop("'from' > min(x) not allowed for method ", method) 701 x <- c(from, x) 702 } 703 if(!from_isL) { 704 if(to < x[n]) 705 stop("'to' < max(x) not allowed for method ", method) 706 x <- c(x, to) 707 } 708 n <- length(x) 709 dx <- x[-1] - x[-n] ## == diff(x) 710 w <- as.numeric(x[n] - x[1]) ## == sum(dx); 711 ## as.n..(.) -> works with "Date" etc 712 nn <- length. - n ## need 'nn' new points in 'n - 1' intervals 713 ## how many in each? 714 ## Want them approximately equidistant, ie. of width ~= w / (nn + 1) 715 ## but do this smartly such that dx[i] / (k1[i] + 1) {= stepsize in interval i} 716 ## is approximately constant 717 k1 <- (nn + n-1) * dx / w - 1 ## ==> sum(k1) == nn 718 ## now "round" the k1[] such that sum(.) remains == nn 719 k <- roundfixS(k1) ## keep the right border, drop the left 720 seqI <- function(i) seq(x[i], x[i+1], length.out=k[i]+2)[-1] 721 l.seq <- lapply(1:(n-1), seqI) 722 ## do.call(c, *), e.g. for new (R-devel 4.1.x) c.Date() [KH]: 723 c(x[1], if(is.object(x)) do.call(c, l.seq) else unlist(l.seq)) 724 } else { 725 nn <- switch(method, "simple" = length., 726 "aim" = length. - n + from_is1 + from_isL) 727 ## a more sophisticated 'method' would have to use iteration, *or* 728 ## interpolate between the 'x' values instead 729 ## which might be considered to be too far from seq() 730 unique(sort(c(x, seq(from, to, length.out = nn)))) 731 } 732 } else x 733}## {seqXtnd} 734 735plotDS <- 736function(x, yd, ys, xlab = "", ylab = "", ylim = rrange(c(yd, ys)), 737 xpd = TRUE, do.seg = TRUE, seg.p = .95, 738 segP = list(lty = 2, lwd = 1, col = 2), 739 linP = list(lty = 1, lwd = 2.5, col = 3), ...) 740{ 741 ## Purpose: Plot Data & Smooth 742 ## ------------------------------------------------------------------------- 743 ## Arguments: do.seg: logical, plot "residual segments" iff T (= default). 744 ## ------------------------------------------------------------------------- 745 ## Author: Martin Maechler, 1990-1994 746 ## 2007: allow ys to be a (xs,ys)-xycoords structure, where {x[] \in xs[]} 747 if((hasMoreSmooth <- !is.numeric(ys))) { 748 ysl <- xy.coords(ys) 749 ixs <- match(x, ysl$x) 750 if(any(is.na(ixs))) 751 stop("'x' inside the 'ys' structure must contain all the observational 'x'") 752 ys <- ysl$y[ixs] 753 } 754 if(is.unsorted(x)) { 755 i <- sort.list(x) 756 x <- x[i] 757 yd <- yd[i] 758 ys <- ys[i] 759 } 760 addDefaults <- function(listArg) { 761 ## trick such that user can call 'segP = list(col = "pink")' : 762 nam <- deparse(substitute(listArg)) 763 P <- as.list(formals(sys.function(sys.parent()))[[nam]])[-1] # w/o "list" 764 for(n in names(listArg)) P[[n]] <- listArg[[n]] 765 P 766 } 767 768 plot(x, yd, xlab = xlab, ylab = ylab, ylim = ylim, ...) #pch = pch, 769 if(!missing(linP)) 770 linP <- addDefaults(linP) 771 if(hasMoreSmooth) 772 lines(ysl, xpd = xpd, lty = linP$lty, lwd = linP$lwd, col = linP$col) 773 else lines(x, ys, xpd = xpd, lty = linP$lty, lwd = linP$lwd, col = linP$col) 774 if(do.seg) { 775 if(!missing(segP)) 776 segP <- addDefaults(segP) 777 segments(x, seg.p*ys + (1-seg.p)*yd, x, yd, 778 xpd = xpd, lty = segP$lty, lwd = segP$lwd, col = segP$col) 779 } 780 invisible() 781} 782 783 784 785##-#### Matrix (or higher Array) stuff ######## 786##-### ------------------------------ ######## 787 788colcenter <- function(mat) sweep(mat,2, apply(mat,2,mean)) 789 790col01scale <- function(mat, scale.func = function(x) diff(range(x)), 791 location.func = mean) 792{ 793 ##-- See also 'scale' (std. S func) -- 794 mat <- sweep(mat,2, apply(mat,2, location.func)) 795 sweep( mat, 2, apply(mat,2, scale.func), "/") 796} 797 798## diag.ex <- function(n) --- now renamed : 799diagX <- function(n) 800{ 801 ## Purpose: Returns "the other diagonal" matrix 802 ## Author: Martin Maechler, Date: Tue Jan 14 1992; Nov.2002 803 ## ---------------------------------------------------------------- 804 ## Arguments: n: integer dimension of matrix 805 ## ---------------------------------------------------------------- 806 m <- numeric(n * n) 807 m[1L+ (n-1L)* seq_len(n)] <- 1 808 dim(m) <- c(n,n) 809 m 810} 811 812xy.grid <- function(x,y) 813{ 814 ## Purpose: Produce the grid used by persp, contour, .. as N x 2 matrix 815 nx <- length(x) 816 ny <- length(y) 817 cbind(rep.int(x,rep.int(ny,nx)), rep.int(y,nx)) 818} 819 820rot2 <- function(xy, phi) 821{ 822 ## Purpose: rotate xy-points by angle 'phi' (in radians) 823 ## ------------------------------------------------------------------------- 824 ## Arguments: xy : n x 2 matrix; phi: angle (in [0, 2pi]) 825 ## ------------------------------------------------------------------------- 826 ## Author: Martin Maechler, Date: 26 Oct 94, 22:16 827 co <- cos(phi); s <- sin(phi) 828 xy %*% t( matrix(c(co,s, -s, co), 2,2) ) 829} 830 831tapplySimpl <- function(X, INDEX, FUN, ...) 832{ 833 ## Purpose: Nicer result for tapply(..) when Function returns 834 ## vector AND there is >= 2 "INDEX", i.e., categories. 835 ## ------------------------------------------------------------------------- 836 ## Arguments: as for tapply, 837 ## FUN: Must return [named, if possible] FIXED length vector 838 ## [num/char] EVEN for NULL and NA ! 839 ## ------------------------------------------------------------------------- 840 ## Author: Martin Maechler, Date: 14 Jun 93, 17:34 841 rl <- tapply(X, INDEX, FUN, ..., simplify = TRUE) 842 if (is.list(rl)) { #-- when >=2 indices AND length(FUN(x)) > 1 --- 843 if(any(Nas <- unlist(lapply(rl, is.null)))) 844 rl[Nas] <- list(FUN(NULL)) 845 array(unlist(rl), 846 dim = c(length(rl[[1]]), dim(rl)), 847 dimnames = c(list(names(rl[[1]])), dimnames(rl)) ) 848 } else rl 849} 850 851 852##-#### "Calculus" Mathematical stuff ######## 853##-### ----------------------------- ######## 854 855u.log <- function(x, c = 1) 856{ 857 ## Purpose: log(.) only for high x- values ... identity for low ones 858 ## This f(x) is continuously differentiable (once). 859 ## f(x) = x for |x| <= c 860 ## f(x) = sign(x)*c*(1 + log(|x|/c)) for |x| >= c 861 ## ------------------------------------------------------------------------- 862 ## Arguments: x: numeric vector; c: scalar > 0 863 ## ------------------------------------------------------------------------- 864 ## Author: Martin Maechler, Date: 24 Jan 95, 17:28 865 if(!is.numeric(c)|| c < 0) stop("'c' must be positive number") 866 r <- x 867 to.log <- abs(x) > c ; x <- x[to.log] 868 r[to.log] <- sign(x) * c * (1 + log(abs(x/c))) 869 r 870} 871 872xy.unique.x <- function(x, y, w, fun.mean = mean, ...) 873{ 874 ## Purpose: given 'smoother data' (x_i, y_i) [and maybe weight w_i] 875 ## with multiple x_i, use unique x's, replacing y's by their mean 876 ## ------------------------------------------------------------------------- 877 ## Author: Martin Maechler, Date: 8 Mar 93, 16:36 878 ##--*--*--*--*--*--*--*--*--*-- x,y,w treatment --*--*--*--*--*--*--*--*-- 879 if(missing(x)) x <- time(y) else 880 if(missing(y)) { 881 if(is.list(x)) { 882 if(any(is.na(match(c("x", "y"), names(x))))) 883 stop("cannot find x and y in list") 884 y <- x$y; x <- x$x; if(!is.null(x$w)) w <- x$w 885 } else if(is.complex(x)) { 886 y <- Im(x); x <- Re(x) 887 } else if(is.matrix(x) && ncol(x) == 2) { 888 y <- x[, 2]; x <- x[, 1] 889 } else if(is.matrix(x) && ncol(x) == 3) { 890 y <- x[, 2]; w <- x[, 3]; x <- x[, 1] 891 } else { 892 y <- x; x <- time(x) 893 } 894 } 895 n <- length(x) 896 if(n != length(y)) stop("lengths of x and y must match") 897 if(missing(w)) w <- rep.int(1,n) 898 else if(n != length(w)) stop("lengths of x and w must match") 899 ##--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*-- 900 gr <- match(x, ux <- unique(x, ...)) 901 cbind(x = ux, 902 y = tapply(y, gr, FUN = fun.mean), 903 w = tapply(w, gr, FUN = sum)) 904} 905 906 907 908##-#### Non-calculus ("Discrete") Mathematical stuff ######## 909##-### -------------------------------------------- ######## 910 911lseq <- function(from, to, length) 912{ 913 ## Purpose: seq(.) : equidistant on log scale 914 ## ---------------------------------------------------------------------- 915 ## Author: Martin Maechler, Date: 3 Feb 2005, 08:34 916 stopifnot(from > 0) 917 exp(seq(log(from), log(to), length.out = length)) 918} 919 920inv.seq <- function(i) { 921 ## Purpose: 'Inverse seq': Return a short expression for the 'index' 'i' 922 ## -------------------------------------------------------------------- 923 ## Arguments: i: vector of (usually increasing) integers. 924 ## -------------------------------------------------------------------- 925 ## Author: Martin Maechler, Date: 3 Oct 95, 18:08 926 ## -------------------------------------------------------------------- 927 ## EXAMPLES: cat(rr <- inv.seq(c(3:12, 20:24, 27, 30:33)),"\n"); eval(rr) 928 ## r2 <- inv.seq(c(20:13, 3:12, -1:-4, 27, 30:31)); eval(r2); r2 929 li <- length(i <- as.integer(i)) 930 if(li == 0) return(expression(NULL)) 931 else if(li == 1) return(as.expression(i)) 932 ##-- now have: length(i) >= 2 933 di1 <- abs(diff(i)) == 1 #-- those are just simple sequences n1:n2 ! 934 i <- i + 0 # coercion to "double", so result has no 'L' appended integers. 935 s1 <- i[!c(FALSE,di1)] # beginnings 936 s2 <- i[!c(di1,FALSE)] # endings 937 mkseq <- function(i, j) if (i==j) i else call(':', i, j) 938 as.call(c(list(as.name('c')), 939 mapply(s1, s2, FUN=mkseq, SIMPLIFY=FALSE, USE.NAMES=FALSE))) 940} 941 942iterate.lin.recursion <- function(x, coeff, delta = 0, nr.it) 943{ 944 r <- c(x, numeric(nr.it)) 945 n <- length(x) 946 ic <- length(coeff):1 947 for(i in 1:nr.it) 948 r[n + i] <- delta + c(coeff %*% r[n + i - ic]) 949 r 950} 951 952quadrant <- function(x,y=NULL) { 953 xy <- xy.coords(x,y); x <- xy$x; y <- xy$y 954 Sgn <- function(u) ifelse(u >= 0, 1, -1) 955 y <- Sgn(y); 2 - y + (y != Sgn(x)) 956} 957 958n.code <- function(n, ndig = 1, dec.codes = c("","d","c","k")) 959{ 960 ##-- convert "round integers" to short char.strings 961 ##-- useful to build-up variable names in simulations 962 ##-- e.g., 963 nd <- length(dec.codes) 964 e10 <- pmin(floor(log10(n) + 1e-12), nd - 1) 965 if (any(e10 < 0)) { 966 e10 <- pmax(0, e10) ; warning("some 'n' too small") 967 } 968 ## IDEA: Things like 969 ## ---- n.code(c(2000,1e4,5e4,6e5,7e6,8e7), 970 ## dec. = c("","d","c","k","-","-","M")) 971 ## could work; (not quite yet, see ex. above) 972##- if(any(id <- is.na(dec.codes) | dec.codes == "-")) { 973##- ## then use previous code for these (things like "20k", "300k") 974##- ## sequentially from the left: 975##- for(k in which(id)) { 976##- dec.codes[k] <- dec.codes[k - 1] 977##- ii <- 1+e10 == k 978##- e10[ii] <- e10[ii] - 1 979##- } 980##- } 981 paste0(round(n/ 10^(e10 + 1 - ndig)), dec.codes[1 + e10]) 982} 983 984code2n <- function(ncod, ndig = 1, dec.codes = c("","d","c","k")) 985{ 986 ## The inverse function to n.code 987 le <- nchar(ncod) 988 cod <- substring(ncod, le, le) 989 as.integer(substring(ncod, 1, le-1)) * 10^(match(cod, dec.codes)-1) 990} 991 992nr.sign.chg <- function(y) 993{ 994 ## Purpose: Compute number of sign changes in sequence 995 ## Be careful with y[i] that were 0 !! 996 y <- sign(c(y)) 997 y <- y[y != 0] 998 sum(y[-1] != y[-length(y)]) 999} 1000 1001unif <- function(n, round.dig = 1 + trunc(log10(n))) 1002{ 1003 ## Purpose: Give regular points on [-c,c] with mean 0 and variance ~= 1 1004 if(n %% 2 == 0) { 1005 if(n > 0) round((2 * 1:n - (n + 1)) * sqrt(3/(n * (n + 1))), round.dig) 1006 } else { 1007 m <- n %/% 2 #--> m+1 = (n+1)/2 1008 ( - m:m) * round(sqrt(6/((m + 1) * n)), round.dig) 1009 } 1010} 1011 1012prt.DEBUG <- function(..., LEVEL = 1) { 1013 stop("prt.DEBUG() is defunct: use a 'verbose' argument or options(verbose=.) instead") 1014 ## if (exists("DEBUG", where = 1) && DEBUG >= LEVEL )# 1015 ## ## 1016 ## cat(paste0("in '", sys.call(sys.nframe()-1)[1], "':"), ..., "\n") 1017} 1018 1019 1020 1021##' @title Read an Emacs Org Table by read.table() 1022## --> ../man/read.org.table.Rd 1023read.org.table <- function(file, header = TRUE, skip = 0, 1024 encoding = "native", fileEncoding = "", 1025 text, ...) { 1026 ## file - text handling --- cut'n'paste from read.table()'s header 1027 if (missing(file) && !missing(text)) { 1028 file <- textConnection(text, encoding = "UTF-8") 1029 on.exit(close(file)) 1030 } 1031 if(is.character(file)) { 1032 file <- if(nzchar(fileEncoding)) 1033 file(file, "rt", encoding = fileEncoding) else file(file, "rt") 1034 on.exit(close(file)) 1035 } 1036 if(!inherits(file, "connection")) 1037 stop("'file' must be a character string or connection") 1038 if(!isOpen(file, "rt")) { 1039 open(file, "rt") 1040 on.exit(close(file)) 1041 } 1042 if(skip > 0L) readLines(file, skip) 1043 ll <- readLines(file, encoding=encoding) 1044 close(file); on.exit() 1045 ## drop |--------+---------+--------+--| : 1046 if(length(i <- grep("---+\\+--", ll[1:3]))) ll <- ll[-i] 1047 ## drop beginning and ending "|" : 1048 ll <- sub("^ *\\|", "", 1049 sub("\\| *$", "", ll)) 1050 if(header) { ## assume header in first 2 lines 1051 ii <- if(nchar(ll[1]) < 2) 2 else 1 1052 ## header line 1053 hl <- ll[ii] 1054 ## drop header line(s) 1055 ll <- ll[-seq_len(ii)] 1056 ## split the header lines into column names 1057 col.names <- sub("^ +", "", sub(" +$", "", strsplit(hl, " *\\| *") [[1L]])) 1058 } 1059 ## drop empty lines at end only 1060 while(grepl("^ *$", tail(ll, 1L))) ll <- ll[-length(ll)] 1061 f.ll <- textConnection(ll)# , encoding = "UTF-8" is *NOT* good 1062 on.exit(close(f.ll)) 1063 read.table(f.ll, header=FALSE, sep = "|", 1064 col.names = col.names, ...) # , encoding = "UTF-8" *not* good 1065} 1066