## ===========================================================================
f.gam <- function()
{
  ## Purpose:
  ## -------------------------------------------------------------------------
  ## Arguments:
  ## -------------------------------------------------------------------------
  ## Author: Werner Stahel, Date:  6 Nov 2000, 09:00
}

# ============================================================================
f.spec1 <- function(a.sample)
{
  ## Purpose:  not working!!!
  ## -------------------------------------------------------------------------
  ## Arguments:
  ## -------------------------------------------------------------------------
  ## Author: Werner Stahel, Date: 27 Apr 2000, 09:27
  t.d <- f.getdata(a.sample, a.plot=F)
  t.d<<-f.baslcorr(t.d, a.function=f.basl.linls, a.plot=F) # f.baslcorr(t.d)
  r.pca<<-f.specpca(t.d, a.plot=F)
  g.summary(r.pca, a.data=t.d)
}
# ============================================================================
g.summary <- function(a.sq, a.data=NULL)
{
  ## Purpose:   plot summary of decomposition
  ## -------------------------------------------------------------------------
  ## Arguments:
  ## -------------------------------------------------------------------------
  ## Author: Werner Stahel, Date: 27 Apr 2000, 09:32
  on.exit(browser())
  par(mfrow=c(1,1))
  split.screen(c(2,1))
  screen(1)
  mtext(s.datatitle,3,1,cex=1.2)
  split.screen(c(1,2),1,erase=F)
  screen(3,F)
  g.plot2(a.y=a.sq$x[,1:2],a.ylab1="scores.1",a.ylab2="scores.2")
  screen(4,F)
  t.err <- f.resid(r.pca, a.data=a.data, a.plotextr=F)
  screen(1,F)
  g.dirspecs(a.sq$rotation[,1,drop=F], a.mf=NULL)
  close.screen(all=T)
  invisible(t.err)
}
# ============================================================================
f.savepca <- function()
{
  ## Purpose:   save the pca results
  ## -------------------------------------------------------------------------
  ## Arguments:
  ## -------------------------------------------------------------------------
  ## Author: Werner Stahel, Date: 25 Apr 2000, 15:54
  r.pca$x <- r.pca$x[,1:s.nc]
  r.pca$rotation <- r.pca$rotation[,1:s.nc]
  t.rw <- r.pca$rotwg
  if (!is.null(t.rw)) r.pca$rotwg <- t.rw[,1:s.nc]
  r.pca$errvar <- NULL
  assign(paste("r.pca",s.series,s.sample,sep="."),r.pca,pos=1)
}
# ============================================================================
g.finaldec <- function(a.sq=r.sq, a.nc=2, a.nspecs=a.nc, a.time=NULL)
{
  ## Purpose:   plot final decomposition
  ## -------------------------------------------------------------------------
  ## Arguments:
  ## -------------------------------------------------------------------------
  ## Author: Werner Stahel, Date: 14 Apr 2000, 10:00
  t.n <- nrow(a.sq$x)
  if (is.null(a.time)) a.time <- if(exists("s.time")) s.time else 1
  if (length(a.time)!=t.n) a.time <- 1:t.n
  close.screen(all=T)
  par(mfrow=c(1,1))
  par(cex=0.7,mex=0.7, mar=c(3,3,1,1))
  split.screen(c(a.nc+1,1))
  split.screen(c(1,2),1)
  screen(a.nc+2)
  g.scpairs(a.sq, a.nc=2)
  screen(a.nc+3)
  matplot(a.time,a.sq$x[,1:a.nc],type="l",xlab="time")
  screen(1,F)
  g.dirspecs(a.sq$rotation[,1:a.nspecs], a.mf=0)
  close.screen(all=T)
}
# ============================================================================
f.getbc <- function(a.sample=s.sample, a.series=s.series,
                    a.project=c.env$project, a.plot=s.plot)
{
  ## Purpose:
  ## -------------------------------------------------------------------------
  ## Arguments:
  ## -------------------------------------------------------------------------
  ## Author: Werner Stahel, Date: 11 Apr 2000, 10:51
  t.dbc <- get(paste(a.series,a.sample,"bc",sep="."))
  t.sd <<- attr(t.dbc,"sd")
  s.datatitle <<- paste(a.project,", ",a.series,".",a.sample,sep="")
  s.title <<- paste(s.datatitle,", baseline-corr.",sep="")
  if (a.plot) g.specsample(t.dbc, a.pldata=F, a.sd=t.sd)
  t.dbc
}
# ============================================================================
f.getdata <- function(a.sample=s.sample, a.series=s.series, a.sd=NULL,
                      a.project=c.env$project, a.plot=s.plot, a.plotc=a.plot,
                      a.mf=a.plot+a.plotc, ...)
{
  ## Purpose:   get stored data
  ## -------------------------------------------------------------------------
  ## Arguments:
  ## -------------------------------------------------------------------------
  ## Author: Werner Stahel, Date: 11 Apr 2000, 10:07
  on.exit(browser())
  t.s <- if (!(is.null(a.series)||a.series==""))
    paste(a.series,a.sample,sep=".") else a.sample
  t.d <- get(t.s)
  t.sd <<- if (is.null(a.sd)) f.ysd(apply(t.d,2,mean)) else
    rep(a.sd, length=ncol(t.d))
  s.datatitle <<- paste(a.project,", ",t.s,sep="")
  s.title <<- paste(s.datatitle,", raw spectra",sep="")
  if (a.plot|a.plotc)
    g.specsample(t.d, a.pldata=a.plot, a.plotc=a.plotc, a.sd=t.sd,
                 a.mf=a.mf, ...)
  if (!Browse) on.exit()
  t.d
}
# =======================================================================
f.basl.rs <- function(a.data)
{
  ## Purpose:   simple regression on index of a.data, used in f.baslcorr
  ## -------------------------------------------------------------------------
  ## Arguments: a.data    spectrum (after preliminary centering)
  ## -------------------------------------------------------------------------
  ## Author: Werner Stahel, Date: 20 Apr 2000, 09:47
  lqs(y~x, data=list(y=a.data, x=1:length(a.data)),method="lts")$coef
}
# =======================================================================
f.basl.quadls <- function(a.data, a.weights=1, a.fit=F)
{
  ## Purpose:   simple regression on index of a.data, used in f.baslcorr
  ## -------------------------------------------------------------------------
  ## Arguments: a.data    spectrum (after preliminary centering)
  ## -------------------------------------------------------------------------
  ## Author: Werner Stahel, Date: 20 Apr 2000, 09:47
  t.r <- lm(y~x+x^2, data=list(y=a.data, x=1:length(a.data)), weights=a.weights)
  if (a.fit) t.r$fit else t.r$coef
}
# =======================================================================
f.basl.linls <- function(a.data, a.weights=1, a.fit=F)
{
  ## Purpose:   simple regression on index of a.data, used in f.baslcorr
  ## -------------------------------------------------------------------------
  ## Arguments: a.data    spectrum (after preliminary centering)
  ## -------------------------------------------------------------------------
  ## Author: Werner Stahel, Date: 20 Apr 2000, 09:47
  t.r <- lsfit(1:length(a.data), a.data, a.weights)
  if (a.fit) a.data-t.r$resid else t.r$coef
}
# =======================================================================
#t.d <- t.db0 <- f.baslcorr(t.d, a.function=f.basl.rs, a.center=median) 
  ##   a.function   function to be used to determine baseline
  ##                median  a constant will be subtracted
  ##                f.basl.rs: a robustly estimated straight line
f.baslcorr <- function(a.data=t.d, a.center=T, a.degree=1,
                       a.sd=t.sd, a.plot=s.plot, a.mf=2)
{
  ## Purpose:   baseline correction using quadratic function
  ## -------------------------------------------------------------------------
  ## Arguments:
  ##   a.data
  ## -------------------------------------------------------------------------
  ## Author: Werner Stahel, Date: 10 Apr 2000, 15:25
  on.exit(browser())
  c.env$step <<- "baseline corr."
  t.xx <- t.x <- 1:ncol(a.data)
  if (a.degree>1) t.xx <- cbind(t.x,t.x^2)
  t.wg <- 1/a.sd^2
  t.d <- scale(a.data,scale=F) 
  if (!a.center) {
  # subtract common baseline (for numerical reasons)
    t.mn <- apply(a.data,2,mean)
    t.d <- sweep(a.data,2,
          t.mn-lsfit(1:length(t.mn), t.mn, t.wg)$residuals)
  }  
  # individual baselines
##-   t.cf <- apply(t.d,1,a.function,a.weights=1/a.sd^2)
  if (a.degree==0) t.cf <- apply(t.wg*t(t.d),2,sum)/sum(t.wg)
    else  t.cf <- lsfit(t.xx, t(t.d), t.wg)$coef
  t.dbc <- if (is.matrix(t.cf)) {
    t.cf <- t(t.cf)
    t.c <- t.cf[,1]+outer(t.cf[,2],t.x)
    if (ncol(t.cf)>2) t.c <- t.c+outer(t.cf[,3],t.x^2)
    t.d-t.c } else  a.data-t.cf  # or sweep(a.data,1,t.cf)
  dimnames(t.dbc) <- dimnames(a.data)
#  dimnames(t.cf) <- list(dimnames(a.data)[[1]],names(t.r$coef))
  attr(t.dbc,"coef") <- t.cf
  t.sd <- a.sd[dimnames(t.dbc)[[2]]]
  attr(t.dbc,"sd") <- a.sd
  s.title <<- paste(s.datatitle, ", baseline deg.",a.degree,sep="")
  browser()
  if (a.plot) g.specsamplebc(a.data=t.dbc, a.sd=a.sd, a.cf=t.cf, a.mf=a.mf)
  if (!Browse) on.exit()
  t.dbc
}
# =======================================================================
g.plot2 <- function(a.x=NULL, a.y, a.xlab="", a.ylab1="", a.ylab2="",
                    a.mar=c(3.5,3.5,2,3.5), a.type="l", ...)
{
  ## Purpose:   plot 2 columns with different scale
  ## -------------------------------------------------------------------------
  ## Arguments:
  ## -------------------------------------------------------------------------
  ## Author: Werner Stahel, Date: 11 Apr 2000, 13:17
  on.exit(browser())
  par(mar=a.mar)
  if (is.null(a.x)) a.x <- 1:nrow(a.y)
  plot(a.x,a.y[,1],type=a.type,xlab=a.xlab, ylab=a.ylab1, ...)
  t.usr <- par("usr")
  t.rg <- range(a.y[,2])
  par(usr=c(t.usr[1:2],t.rg+diff(t.rg)*0.05*c(-1,1)))
  if (a.type=="l")  lines(a.x,a.y[,2],lty=2) else
    points(a.x,a.y[,2], pch=2)
  axis(4)
  mtext(a.ylab2,4,2.5)
  g.date(sure=F)
  if (!Browse) on.exit()
}
# =======================================================================
g.fan <- function(a.data=t.d)
{
  ## Purpose:   plot sample of spectra, shifted
  ## -------------------------------------------------------------------------
  ## Arguments:
  ## -------------------------------------------------------------------------
  ## Author: Werner Stahel, Date: 26 Apr 2000, 09:51
  t.wl <- as.numeric(dimnames(a.data)[[2]])
  g.manycurves(t.wl, t(scale(a.data,scale=F)),
               xlab="wavelength nm",ylab="centered abs.")
}
# =======================================================================
g.manycurves <- function(a.x=NULL, a.y=NULL, a.shift=NULL, 
                         a.mar=c(3.5,3.5,2,3.5), 
                         a.title=s.title, a.lwd=2, a.mf=1, a.cnm=NULL, ...)
{
  ## Purpose:   Plot many curves, shifted
  ## -------------------------------------------------------------------------
  ## Arguments:
  ##   a.y      matrix of data to be plotted like with  matplot
  ##   a.shift  space between "stories"
  ##   a.cnm    column names
  ## -------------------------------------------------------------------------
  ## Author: Werner Stahel, Date: 26 Apr 2000, 08:20
  on.exit(browser())
  if (is.null(a.y)) {
    if (is.null(a.x)) stop("g.manycurves: nothing to plot")
    a.y <- a.x
    a.x <- 1:nrow(a.y) }
  if (is.null(a.x)) a.x <- 1:nrow(a.y)
  t.nc <- ncol(a.y)
  if (is.null(a.cnm)) a.cnm <- dimnames(a.y)[[2]]
  if (is.null(a.cnm)) a.cnm <- as.character(1:t.nc)
  if (is.null(a.shift)) a.shift <- diff(range(a.y))/sqrt(t.nc)
  t.zeros <- seq(t.nc-1,0,-1)*a.shift
  t.y <- sweep(a.y,2,t.zeros,"+")
  t.ylim <- range(t.y)
  if (length(a.mf)>0) u.mf(a.mf)
  if (length(a.mar)) par(mar=rep(a.mar,length=4))
  matplot(a.x,t.y, axes=F, type="l", lwd=a.lwd, ylim=t.ylim, yaxs="i",
          main=a.title, ...)
  axis(1)
  axis(4,at=t.zeros, label=a.cnm, srt=0)
  box()
  t.xx <- par("usr")[1:2]
  matlines(t.xx,rbind(t.zeros,t.zeros))
  g.date(sure=F)
  if (!Browse) on.exit()
}
# =======================================================================
f.ysd <- function(a.y, a.prm=c.intvarpar)
{
  ## Purpose:    calculate standard dev. of absorbances
  ## -------------------------------------------------------------------------
  ## Arguments:
  ## -------------------------------------------------------------------------
  ## Author: Werner Stahel, Date: 10 Apr 2000, 11:41
  pmax(a.prm["sdmin"], a.prm["sdfac"]*sqrt(10^(2*a.y)+1) )
}
# =======================================================================
f.ysd.old <- function(a.y, a.prm=c.intvarpar)
{
  ## Purpose:    calculate standard dev. of absorbances
  ## -------------------------------------------------------------------------
  ## Arguments:
  ## -------------------------------------------------------------------------
  ## Author: Werner Stahel, Date: 10 Apr 2000, 11:41
  on.exit(browser())
  t.lis <- log10(10^(2*a.y)+1)
  t.lv <- a.prm["v0"]+(t.lis-a.prm["i0"])*
    ifelse(t.lis<a.prm["i0"],a.prm["s0"],a.prm["s1"])
  on.exit()
  10^(t.lv/2)
}
# =======================================================================
f.readspec <- function(a.filestring=s.filestring, a.series=s.series,
   a.sample=s.sample,a.project=c.env$project, a.wlnr=s.wlnr,
                       a.plot=s.plot, a.plotc=a.plot)
{
  ## Purpose:   read spectra
  ## -------------------------------------------------------------------------
  ## Arguments:
  ## -------------------------------------------------------------------------
  ## Author: Werner Stahel, Date:  4 Apr 2000, 13:34
  on.exit(browser())
  t.dd <- read.table(paste(a.filestring,a.sample,".txt",sep=""),header=T)
  t.specnames <<- names(t.dd)
  if (t.dd[1,1]>10) {
    t.wavel <- t.dd[a.wlnr,1]
    t.dd <- as.matrix(t.dd[,-1])
  } else {
    t.wavel <- as.numeric(row.names(t.dd))
    t.dd <- as.matrix(t.dd)
  }
  if (max(a.wlnr)>nrow(t.dd)+1) {
    a.wlnr <- a.wlnr[a.wlnr<=nrow(t.dd)+1]
    warning("not all wavelengths provided") }
  t.wldiff <- diff(a.wlnr[1:2])
  t.wlnr <- a.wlnr[-1]-round(t.wldiff/2)
  t.nwl <- length(t.wlnr)
  t.d <- apply(matrix(c(t.dd[a.wlnr[1]+seq(1,t.nwl*t.wldiff)-1,]),t.wldiff),2,mean)
  t.d <- matrix(t.d,ncol=t.nwl,byrow=T)
  s.wavel <<- t.wavel[t.wlnr]
  dimnames(t.d) <- list(paste("S",1:nrow(t.d),sep="."),format(round(s.wavel,2)))
  s.datatitle <<- paste(a.project,", ",a.series,".",a.sample,sep="")
  s.title <<- paste(s.datatitle,", raw spectra",sep="")
  t.sd <<- f.ysd(apply(t.d,2,mean))
  g.specsample(t.d, a.sd=NULL, a.pldata=a.plot, a.plotc=a.plotc)
  if (!Browse) on.exit()
t.d
}
# =======================================================================
g.specsample <- function(a.data=t.d, a.pldata=T, a.plotc=T, a.sd=t.sd, a.fac=0,
                         a.clr=1, a.title=s.title, a.mf=a.pldata+a.plotc, ...)
{
  ## Purpose:   plot a sample of spectra
  ## -------------------------------------------------------------------------
  ## Arguments:
  ## -------------------------------------------------------------------------
  ## Author: Werner Stahel, Date:  4 Apr 2000, 14:02
  on.exit(browser())
  t.wl <- as.numeric(dimnames(a.data)[[2]])
  par(lab=c(8,5,5))
  if (a.mf>0) u.mf(a.mf)
  if(a.pldata)  
    matplot(t.wl ,t(a.data), type="l", col=1, lab=c(10,5,5),
            xlab="wavelength nm", ylab="absorbance", main=a.title, ...)
  if (a.plotc) {
    if (length(a.sd)<=1) a.sd <- f.ysd(apply(a.data,2,mean))
    t.data <- scale(a.data,scale=F)
    t.tit <- ""
    if (a.fac==0) { t.data <- sweep(t.data,2,t.sd,"/")
                    t.tit <- " and scaled" }
    matplot(t.wl ,t(t.data), type="l", col=1, lab=c(10,5,5),
          lty=1+1:nrow(a.data), xlab="wavelength nm", ylab="absorbance",
          main=paste(a.title,"centered",t.tit), ...)
  if (a.fac>0) matlines(t.wl, a.fac*cbind(a.sd,-a.sd), col=a.clr, lty=1, lwd=2)
  }
  if (a.pldata|a.plotc)  g.date(sure=F)
  if (!Browse) on.exit()
  "done"
}
# =======================================================================
g.specsamplebc <- function(a.data=t.dbc, a.sd=t.sd, a.cf=t.cf, a.clr=1,
                           a.title=s.title, a.mf=2, ...)
{
  ## Purpose:   plot a sample of spectra
  ## -------------------------------------------------------------------------
  ## Arguments:
  ## -------------------------------------------------------------------------
  ## Author: Werner Stahel, Date:  4 Apr 2000, 14:02
  on.exit(browser())
  if (a.mf>0) u.mf(a.mf)
  t.wl <- as.numeric(dimnames(a.data)[[2]])
  par(lab=c(8,5,5))
  if (length(a.sd)<=1) a.sd <- attr(a.data,"sd")  
  matplot(t.wl ,t(sweep(a.data,2,a.sd,"/")), type="l", col=1, lab=c(10,5,5),
          lty=1+1:nrow(a.data), xlab="wavelength nm", ylab="absorbance",
          main=paste(a.title,"baseline corr. & scaled"), ...)
  if (length(a.cf)<=1) a.cf <- attr(a.data,"cf")  
  if (is.matrix(a.cf))
      g.plot2(a.y=a.cf, a.ylab1="baseline const",
              a.ylab2="baseline slope", a.xlab="index of spectrum",
              main="baseline coefficients")
      else  plot(a.cf, ylab="baseline const", xlab="index of spectrum",
                 main="baseline coefficients", type="l")
  g.date(sure=F)
  if (!Browse) on.exit()
  "done"
}
# =======================================================================
f.selwavel <- function(a.data=t.d, a.start=0, a.end=4000, a.plot=s.plot)
{
  ## Purpose:   select wavelengths
  ## -------------------------------------------------------------------------
  ## Arguments:
  ## -------------------------------------------------------------------------
  ## Author: Werner Stahel, Date:  4 Apr 2000, 14:07
  on.exit(browser())
  t.wl <- as.numeric(dimnames(a.data)[[2]])
  t.j <- t.wl>=a.start & t.wl<=a.end
  s.wavel <<- t.wl[t.j]
  t.sd <<- t.sd[t.j]
  t.d <- a.data[,t.j]
  if (a.plot) g.specsample(t.d, a.sd=t.sd, a.pldata=F)
  if (!Browse) on.exit()
  t.d
}

# =======================================================================
f.errvar <- function(a.pc=r.pca, a.nc=s.nc, a.data=NULL, a.rerr=F)
{
  ## Purpose:   estimate error variances 
  ## -------------------------------------------------------------------------
  ## Arguments:
  ## -------------------------------------------------------------------------
  ## Author: Werner Stahel, Date:  3 Apr 2000, 11:28
  on.exit(browser())
  t.n <- nrow(a.pc$x)
  t.m <- nrow(a.pc$rotation)
  r.var <- matrix(NA,1+length(a.nc),t.m)
  t.dt <- if (is.null(a.data)) a.pc$x %*% t(a.pc$rotation) else a.data
  r.var <- rbind(apply(t.dt^2,2,sum))
  if (length(a.nc)>0) 
  for (t.nc in a.nc) {
    t.err <- t.dt - a.pc$x[,1:t.nc,drop=F] %*% t(a.pc$rotation[,1:t.nc,drop=F])
    r.var <- rbind(r.var, apply(t.err^2,2,sum)*(t.m/((t.m-t.nc))) )
  }
  else t.err <- t.dt
  r.var <- r.var/(t.n-1)
  dimnames(r.var) <- list(as.character(c(0,a.nc)),dimnames(a.pc$rotation)[[1]])
  on.exit()
  if(a.rerr) list(var=r.var, err=t.err) else r.var
}
# =======================================================================
g.errvar <- function(a.errvar=r.pca$errvar, a.nc=s.nc, a.prows=2,
                     a.title=s.title)
{
  ## Purpose:   plot residual error variation
  ## -------------------------------------------------------------------------
  ## Arguments:
  ## -------------------------------------------------------------------------
  ## Author: Werner Stahel, Date:  4 Apr 2000, 14:26
  if (!is.null(a.prows)&&a.prows>0) u.mf(a.prows)
  t.wl <- as.numeric(dimnames(a.errvar)[[2]])
  matplot(t.wl,sqrt(t(a.errvar[1:(a.nc+1),])),type="l",col=1, log="y",
          xlab="wavelength",ylab="residual st.dev.",
          main=paste(a.title," residual error standard dev."))
  abline(h=0)
  matplot(t.wl,t(a.errvar[1:a.nc,]-a.errvar[1+1:a.nc,]),type="l", 
          col=1, xlab="wavelength",ylab="change in res. st.dev.")
  abline(h=0)
  g.date(sure=F)
}

# =======================================================================
f.specpca <- function(a.data=t.d, a.weights=1/t.sd, a.nc=s.nc0,
                      a.plot=s.plot, a.mf=2)
{
  ## Purpose:   PCA with weights
  ## -------------------------------------------------------------------------
  ## Arguments:
  ## -------------------------------------------------------------------------
  ## Author: Werner Stahel, Date:  4 Apr 2000, 14:16
  on.exit(browser())
  s.title <<- paste(s.datatitle, "PCA", sep=", ")
  t.dwg <- sweep(as.matrix(a.data), 2, a.weights, "*")
  r.pca <- f.prcomp2(t.dwg)
  # reverse directions if scores decrease
  t.dec <- r.pca$x[nrow(r.pca$x),]<r.pca$x[1,]
  if (any(t.dec)) {
    r.pca$x[,t.dec] <- -r.pca$x[,t.dec]
    r.pca$rotation[,t.dec] <- -r.pca$rotation[,t.dec] }    
  r.pca$rotwg <- sweep(r.pca$rotation,1,a.weights,"/")
  # error variances, weighted
  r.pca$errvar <- f.errvar(r.pca, 1:a.nc)
  if (a.plot)  g.specpca(r.pca, a.nc=a.nc, a.mf=a.mf)
  on.exit()
  r.pca
}
# =======================================================================
g.scpairs <- function(a.pca=r.pca, a.nc=s.nc, a.labels=NULL, a.lty=NULL,
                      a.lwd=0.5, a.cex=par("cex"), a.title=s.title)
{
  ## Purpose:   pairs of scores
  ## -------------------------------------------------------------------------
  ## Arguments:
  ## -------------------------------------------------------------------------
  ## Author: Werner Stahel, Date: 13 Apr 2000, 09:45
  on.exit(browser())
  a.nc <- min(a.nc,ncol(a.pca$x))
  t.x <- a.pca$x[,1:a.nc]
  if (is.logical(a.labels)) {
    llabels <- a.labels
    a.labels <- NULL } else llabels <- TRUE
  if (llabels) {
    if (is.null(a.labels)) a.labels <- dimnames(t.x)[[1]]
    if (is.null(a.labels)) a.labels <- as.character(1:nrow(t.x))
  }
  if (is.null(a.lty)) a.lty <- 1 + llabels
  if (a.nc>2) 
    pairs(t.x, panel = if (llabels) g.linesandtext else lines,
          main=a.title, a.cex=a.cex, lty=a.lty, a.labels=a.labels)
##-     g.pairs(t.x[,-1],t.x[,-ncol(t.x)], connect=TRUE, main=a.title,
##-             cex=a.cex, lty=a.lty)
  else {
    if (!is.numeric(split.screen())) u.mf(1)
    t.nm <- dimnames(t.x)[[2]]
    plot(t.x[,1],t.x[,2], type="b", pch="  ", cex=a.cex,
         lty=a.lty, lwd=a.lwd, xlab=t.nm[1], ylab=t.nm[2])
    text(t.x[,1],t.x[,2], a.labels, cex=a.cex) 
    title(a.title,outer=par("oma")[3]>0) }
  g.date(sure=F)
  if (!Browse) on.exit()
}
# =======================================================================
g.specpca <- function(a.pca=r.pca, a.nc=s.nc, a.clr=1, a.mf=2,
                      a.title=s.title)
{
  ## Purpose:   plots to evaluate pca
  ## -------------------------------------------------------------------------
  ## Arguments:
  ## -------------------------------------------------------------------------
  ## Author: Werner Stahel, Date: 13 Apr 2000, 09:39
  on.exit(browser())
  if (a.mf>0) u.mf(a.mf)
  par(mar=c(3.5,3.5,2,3.5))
  t.v <- c(a.pca$sdev[1:4]^2,rep(0,5))
  barplot(t.v/sum(a.pca$sdev^2),main=a.title, axes=F) 
  axis(2)
# ---
  t.sc <- a.pca$x[,1:a.nc]
  t.nr <- nrow(t.sc)
  par(usr=c(-t.nr,t.nr+1,1.05*min(t.sc),1.05*max(t.sc)))
  matlines(1:t.nr,t.sc, col=1)
  axis(4)
  box()
# ---
#  g.errvar(a.pca$errvar, a.nc=a.nc, a.prows=NULL)
  t.rot <- sweep(a.pca$rotation[,1:a.nc],2,a.pca$sdev[1:a.nc],"*")
  t.wl <- as.numeric(dimnames(t.rot)[[1]])
  matplot(t.wl, t.rot, type="l", col=a.clr, xlab="wavelength  nm",
          main="spectra of principal components")
  abline(h=0)
  g.legcorner(a.nc)
  g.date(sure=F)
  if (!Browse) on.exit()
  "g.specpca: done"
}
# =======================================================================
g.legcorner <- function(a.n, a.horfac=0.1, a.edge=0.03, a.text=1:a.n,
                        a.lty=1:a.n)
{
  ## Purpose:
  ## -------------------------------------------------------------------------
  ## Arguments:
  ## -------------------------------------------------------------------------
  ## Author: Werner Stahel, Date:  6 Apr 2000, 14:12
  t.us <- par("usr")
  legend(t.us[2]-a.horfac*diff(t.us[1:2]),t.us[4]-a.edge*diff(t.us[3:4]),
         as.character(a.text), lty=a.lty)
}
# =======================================================================
g.dirspecs <- function(a.dec=r.sq, a.mf=NA, a.title=s.title)
{
  ## Purpose:   plot rotation matrix as spectra
  ## -------------------------------------------------------------------------
  ## Arguments:
  ## -------------------------------------------------------------------------
  ## Author: Werner Stahel, Date:  5 Apr 2000, 08:13
  t.rot <- if(is.list(a.dec)) a.dec$rotwg else a.dec
  if (is.null(t.rot)) t.rot <- a.dec$rotation
  t.nc <- ncol(t.rot)
  # --- number of "rows" to plot
  if (is.numeric(a.mf)) {
    if (is.na(a.mf)) a.mf <- t.nc
    if(a.mf>0) u.mf(a.mf) }
  # ---
  t.wl <- as.numeric(dimnames(t.rot)[[1]])
  t.lab <- dimnames(t.rot)[[2]]
  for (l.k in 1:t.nc) {    
    if (is.numeric(split.screen())) screen(screen()+1)
    plot(t.wl, t.rot[,l.k], type="l",
         xlab="wavelength  nm", ylab=t.lab[l.k])
    abline(h=0,lty=4)
  }
  mtext(a.title,3,1,outer=T, cex=1.5)
  g.date(sure=F)
}

# =======================================================================
f.rotation <- function(a.pc=r.pca, a.dir=t.dir, a.plot=s.plot)
{
  ## Purpose:
  ## -------------------------------------------------------------------------
  ## Arguments:
  ## -------------------------------------------------------------------------
  ## Author: Werner Stahel, Date:  5 Apr 2000, 08:28
  on.exit(browser())
  s.title <<- paste(s.datatitle, "rotated directions", sep=", ")
  t.ip <- matrix( a.dir, ncol=2, byrow=T)
  t.nc <- nrow(t.ip)
  t.transf <- a.pc$x[t.ip[,2],1:t.nc]-a.pc$x[t.ip[,1],1:t.nc]
  dimnames(t.transf)[[1]] <- paste("dir",1:t.nc,sep=".")
  t.tinv <- solve(t.transf)
  t.sc <- a.pc$x[,1:t.nc]%*%t.tinv
  dimnames(t.sc)[[2]] <- paste("dir",1:t.nc,sep=".")
  t.rot <- a.pc$rotation[,1:t.nc]%*%t(t.transf)
  t.rotw <- if (is.null(a.pc$rotwg)) NULL else
    a.pc$rotwg[,1:t.nc]%*%t(t.transf)
  if (!Browse) on.exit()
  list(x=t.sc, rotation=t.rot, rotwg=t.rotw)
}
# =======================================================================
f.resid <- function(a.pc=r.pca, a.rmslim=c.rmslim, a.data=NULL, a.nc=s.nc,
                    a.relerr=F, 
                    a.plot=s.plot, a.plotextr=T, a.title=s.title)
{
  ## Purpose:   Analyze residuals
  ## -------------------------------------------------------------------------
  ## Arguments:
  ## -------------------------------------------------------------------------
  ## Author: Werner Stahel, Date:  5 Apr 2000, 09:11
  on.exit(browser())
  t.ev <- f.errvar(a.pc, a.nc=a.nc, a.rerr=T)
  t.errsd <- NULL
  t.err <- if (a.relerr) {
    t.errsd <- sqrt(t.ev$var[nrow(t.ev$var),])
    sweep(t.ev$err,2,pmax(t.errsd,median(t.errsd)),"/")
  } else t.ev$err
  t.m <- length(t.err)
  t.rms <- sqrt(apply(t.err^2,1,sum)/t.m)
  t.extr <- which(t.rms>a.rmslim)
  if (a.plot) {
#    u.mf(2)
    plot(t.rms, type="l", ylim=c(0,1.05*max(t.rms)),yaxs="i",
         xlab="sample",ylab="residual RMS", main="residual rms")
    abline(h=c(1,a.rmslim),lty=2)
    if (length(t.extr)>0 & a.plotextr) {
      t.wl <- as.numeric(dimnames(a.pc$rotation)[[1]])
      matplot(t.wl,t(t.ev$err[t.extr,,drop=F]),type="l",col=1,
              xlab="wavelength",ylab="residual")
      abline(h=0,col=2)
      if (length(t.errsd)>0) {
        lines(t.wl,t.errsd,lty=2,col=2)
        lines(t.wl,-t.errsd,lty=2,col=2)
      }
    }
    g.date(sure=F)
    mtext(a.title, 3,0.5,cex=1.5,outer=T)
  }
  if (!Browse) on.exit()
  invisible(list(var=t.ev$var, resid=t.ev$err, rms=t.rms))
}
c.rmslim <- 1.5
# =======================================================================
g.resid <-
function(a.pc=r.pca, a.rmslim=c.rmslim, a.data=NULL, a.title=s.title,
         a.plot=s.plot)
{
  ## Purpose:
  ## -------------------------------------------------------------------------
  ## Arguments:
  ## -------------------------------------------------------------------------
  ## Author: Werner Stahel, Date:  5 Apr 2000, 09:11
  on.exit(browser())
  t.err <- f.errvar(a.pc, a.nc=s.nc, a.data=a.data, a.rerr=T)
       t.errsd <- sqrt(t.err$var[2,])
  t.errst <- sweep(t.err$err,2,pmax(t.errsd,median(t.errsd)),"/")
  t.rms <- sqrt(apply(t.errst^2,1,sum)/length(t.errst))
  t.extr <- which(t.rms>a.rmslim)
  if (a.plot) {
    u.mf(2)
    plot(t.rms,xlab="sample",ylab="residual RMS",
         ylim=c(0,1.05*max(t.rms)),yaxs="i")
    abline(h=1,lty=2)
    t.wl <- as.numeric(dimnames(a.pc$rotation)[[1]])
    if (length(t.extr)>0) {
      matplot(t.wl,t(t.err$err[t.extr,,drop=F]),type="l",col=1,
              xlab="wavelength",ylab="residual")
      lines(t.wl,t.errsd,lty=2,col=2)
      lines(t.wl,-t.errsd,lty=2,col=2)
      abline(h=0,col=2)
    }
    g.date(sure=F)
    mtext(a.title, 3,0.5,cex=1.5,outer=T)
  }
  if (!Browse) on.exit()
  invisible(list(var=t.err$var, resid=t.err$err, st.resid=t.errst,
                 st.rms=t.rms))
}
# ======================================================
f.prcomp2 <- function(a.x, a.center=T, a.retx = TRUE)
{
  ## Purpose: PCA when number of variables >> number of observations
  ## -------------------------------------------------------------------------
  ## Arguments:
  ##   x   matrix (rows=units, columns=variables)
  ## -------------------------------------------------------------------------
  ## Author: Marcel Wolbers, Date: 17 Mar 2000, 11:45
  t.ct <- rep(0,ncol(a.x))
  # center data
  if (a.center) {
    t.ct <- apply(a.x,2,mean)
    a.x <- sweep(a.x,2,t.ct) }
  # svd
  t.d <- t(a.x)
  t.svd <- svd(t.d)
  # rearrange
  t.s <- list(d = t.svd$d, u = t.svd$v, v = t.svd$u)	# svd of original
  t.rank <- sum(t.s$d > 0)
  if(t.rank < nrow(t.d))
    t.s <- list(d=t.s$d[1:t.rank], u=t.s$u[,1:t.rank], v=t.s$v[,1:t.rank])
  t.cnm <- paste("comp",1:t.rank,sep=".")
  t.x <-  NULL
  if(a.retx) {
    t.x <- sweep(t.s$u,2,t.s$d,"*")
    dimnames(t.x) <- list(dimnames(a.x)[[1]],t.cnm) }
  t.s$d <- t.s$d/sqrt(max(1, t.rank - 1))
  names(t.s$d) <- t.cnm
  dimnames(t.s$v) <- list(dimnames(a.x)[[2]],t.cnm)
  list(sdev = t.s$d, rotation = t.s$v, x=t.x, center=t.ct)
}
# ======================================================
f.colmed <- function(data,window=5,fun=median)
{
  ## Purpose:  medians over  window  columns
  ## ----------------------------------------------------------------------
  ## Arguments:
  ## ----------------------------------------------------------------------
  ## Author: Werner Stahel, Date: 30 Apr 2003, 09:20
  on.exit(browser())
  lnc <- ncol(data)%/%window
  ldt <- matrix(t(as.matrix(data[,1:(lnc*window)])),5)
  lr <- matrix(apply(ldt,2,fun),nrow(data),byrow=TRUE)
  ldn <- dimnames(data)
  dimnames(lr) <- 
    list(ldn[[1]],ldn[[2]][seq(ceiling(window/2),by=window,length=ncol(lr))])
  if (!Browse) on.exit()
  lr
}
# ======================================================
g.linesandtext <- function(a.x, a.y, a.labels=as.character(1:length(a.x)),
  a.col=1, lty=2, a.lwd=0.5, a.cex=par("cex"))
{
  ## Purpose:   lines and text, to be used as panel function for pairs
  ## -------------------------------------------------------------------------
  ## Arguments:
  ## -------------------------------------------------------------------------
  ## Author: Werner Stahel, Date: 17 Mar 2000, 16:04
  lines(a.x, a.y, lty=lty, lwd=a.lwd, type="b", pch="  ", cex=a.cex)
  text(a.x, a.y, a.labels, cex=a.cex, col=a.col)
}

u.mf <- function(row=1, col=1, mar=c(3.5,3.5,2,1), oma=c(1,1,2,1),
                 mgp=c(2,0.7,0))
{
  ## Purpose:   multiple frames, with defaults
  ## -------------------------------------------------------------------------
  ## Arguments:
  ## -------------------------------------------------------------------------
  ## Author: Werner Stahel, Date: 17 Mar 2000, 14:47
  par(mfrow=c(row,col), mar=mar, oma=oma, mgp=mgp)
}


## ===========================================================================
g.date <-
function(sure=T, outer = NULL, project=c.env$project, step=c.env$step,
         stamp=c.env$stamp, ...)
{
## Purpose:   plot date and project information
## -------------------------------------------------------------------------
## Arguments:
##   sure     if F, the function only plots its thing if  c.goptions$date>0
##   outer    if T, the date is written in the outer margin
##   project  project title
##   step     title of step of data analysis
##   ...      arguments to  mtext , e.g., line=3
## -------------------------------------------------------------------------
## Author: Werner Stahel, Date: 13 Aug 96, 09:00
#  on.exit(browser())
  if (is.null(outer)) outer <- par('oma')[4]>0 
  t.txt <- u.date()
    if (!is.null(project)) t.txt <- paste(t.txt,project,sep=" | ")
    if (!is.null(step)) t.txt <- paste(t.txt,step,sep=" | ")
  if( sure | stamp==2 | ( stamp==1 &
     # last figure on page
     { t.tfg <- par("mfg") ; all(t.tfg[1:2]==t.tfg[3:4]) } )  ) 
       mtext(t.txt, 4, cex = 0.6, adj = 0, outer = outer, ...) 
  invisible(t.txt)
}

if (!exists("c.env")) c.env <- list(project="",step="",stamp=1,out=1)
s.prt <- T
s.plot <- T
##- Project <- list(short="Spectra")
##- Step <- list(short="")
##- c.goptions <- list(date=1)
c.thisisS <- F
Browse <- F

c.intvarpar <- c(sdfac=5e-06, sdmin=1.5e-04,i0=2.5,v0=-7.7,s0=2/30,s1=0.881)

## ===========================================================================
warnings <-
function (...) 
{
    if (!(n <- length(last.warning))) 
        return()
    tb <- table(names(last.warning))
    cat(paste("Warning: ",names(tb)," (",tb,")\n",sep=""))
    invisible(tb)
  }