protorow <- list(testno=as.integer, info=as.integer, #
  drcond=as.double, scheme=as.integer, #
  xstatus=as.integer, zstatus=as.integer, #
  dx=as.double, dz=as.double, #
  prev_dx=as.double, prev_dz=as.double, #
  xbnd=as.double, zbnd=as.double, #
  oldcnt=as.integer, newcnt=as.integer, #
  dxrmax=as.double, dzrmax=as.double, #
  xerr=as.double, zerr=as.double, #
  berr=as.double, crmax=as.double, pgrowth=as.double, #
  seed0=as.integer,seed1=as.integer,seed2=as.integer,seed3=as.integer, #
  test_type=as.integer, #
  srcond=as.double, #
  dxbnd=as.double,dzbnd=as.double,dberr=as.double, #
  firstdx=as.double,firstdz=as.double, #
  prev_x=as.double, x=as.double, crmin=as.double, #
  srcnt=as.integer, time=as.double, xmin=as.double, #
  maxtdiff=as.double, maxtcdiff=as.double, maxtres=as.double, #
  residual=as.double, final.dx.x=as.double, final.dz.z=as.double, #
  xtmin=as.double, eqres=as.double, ymin=as.double, yerr=as.double, #
  y=as.double, xt=as.double, worst.x=as.double, worst.c=as.double, #
  x.stop.cnt=as.integer, z.stop.cnt=as.integer, #
  equed=as.integer, xmode=as.integer, xcond=as.double, #
  ytcond=as.double, srccond=as.double, srncond=as.double, #
  flags=as.integer, flags.equ.a=as.integer, flags.equ.b=as.integer, #
  flags.factor=as.integer, flags.soln=as.integer, #
  drccond=as.double, drncond=as.double,
  fact.time=as.double,
  init.zerr=as.double,
  init.xerr=as.double,
  equil.cnt=as.integer,
  total.time=as.double,
  itref.time=as.double,
  rcond.time=as.double,
  nrcond.time=as.double,
  crcond.time=as.double,
  berr.time=as.double,
  equil.time=as.double
  )
colnames <- names(protorow)
rowoff <- 1:length(colnames)
names(rowoff) <- colnames
savednames <- colnames[c(-rowoff["testno"])]
savedoff <- rowoff[savednames]
savedfn <- protorow[savednames]
alltests.names <- c("drcond", "crmax", "pgrowth", "seed0", "seed1", "seed2", "seed3", "dxbnd", "dzbnd", "dberr", "equed", "randx", "xcond")
join.names <- c("crmax", "seed0", "seed1", "seed2", "seed3", "randx", "xcond")
test.specific.names <-savednames[!(savednames %in% alltests.names)]
nuke.em.names <- c("prev_dx","prev_dz","dxbnd","dzbnd","dberr","firstdx","firstdz","prev_x","srcnt","time","maxtdiff","maxtcdiff")
nuke.em.off <- savedoff[nuke.em.names]
  
narow <- protorow[savednames]
for (i in 1:length(narow)) { narow[i] = protorow[[savedoff[i]]](NA); }

scheme.levels <- c(0,1,2,-1,-2)
scheme.labels <- c("single-res","double-res","double-x","working","wilk")
status.levels <- c(0,1,2,3)
status.labels <- c("unstable","working","conv","noprog")

s.workingeps <- 2^-24
d.workingeps <- 2^-53
s.rteps <- sqrt(s.workingeps)
d.rteps <- sqrt(d.workingeps)
workingeps <- function (prec="s") { return(switch(prec,s=s.workingeps,d=d.workingeps,-1)); }
rteps <- function(prec="s") { return(switch(prec,s=s.rteps,d=d.rteps,-1)); }
#tinybrk <- 2*sqrt(100)*2^-24
#strength.cuts <- c(0,tinybrk,rteps,Inf)
strength.lbls <- c("strong","weak","fail")
tiny.thresh <- function (n = 100, prec="s") { return(sqrt(n)*workingeps(prec)) }
strength.cuts <- function (n = 100, prec="s") { return(c(0, 2*tiny.thresh(n,prec), rteps(prec), Inf)) }
cthresh <- function (n = 100, prec="s") { return(max(10,sqrt(n)) * workingeps(prec)); }
condclass.cuts <- function (n = 100, prec="s") { return(c(0, cthresh(n,prec), Inf)); }
condclass.labels <- c("ill","well-enough")
condhard.cuts <- function (n = 100, prec="s") { return(c(0, cthresh(n,prec), rteps(prec), Inf)); }
condhard.labels <- c("hard","medium","easy")

#require(R.matlab)

categorize.estimates <- function (d, ...) {
  xcat <- cut(d$xbound / d$xerror, breaks=c(-Inf,0.01,0.1,10,100,Inf),labels=c("extreme under", "under", "ok", "over", "extreme over"))
  levels(xcat) <- c(levels(xcat),"hosed")
  xcat[d$xerr.cls=="fail" & d$xbnd.cls=="fail"] <- factor("hosed")
  d$xcat <- xcat

  zcat <- cut(d$zbound / d$zerror, breaks=c(-Inf,0.01,0.1,10,100,Inf),labels=c("extreme under", "under", "ok", "over", "extreme over"))
  levels(zcat) <- c(levels(zcat),"hosed")
  zcat[d$zerr.cls=="fail" & d$zbnd.cls=="fail"] <- factor("hosed")
  d$zcat <- zcat
  return(d)
}

classify.conds <- function(d, n=100, prec="s") {
  wellill.cuts.v <- condclass.cuts(n=n,prec=prec)
  difficulty.cuts.v <- condhard.cuts(n=n,prec=prec)

  d$xcond.cls <- cut(d$srncond, wellill.cuts.v, labels=condclass.labels)
  d$zcond.cls <- cut(d$srccond, wellill.cuts.v, labels=condclass.labels)
  d$xhard.cls <- cut(d$srncond, difficulty.cuts.v, labels=condhard.labels)
  d$zhard.cls <- cut(d$srccond, difficulty.cuts.v, labels=condhard.labels)
  d$hard.cls <- cut(pmin(d$srncond, d$srccond), difficulty.cuts.v, labels=condhard.labels)

  return (d);
}

just.classify.conds <- function(d, n=100, prec="s", fudge=1.0) {
  wellill.cuts.v <- fudge*condclass.cuts(n,prec)
  difficulty.cuts.v <- fudge*condhard.cuts(n,prec)

  xcond.cls <- cut(d$srncond, wellill.cuts.v, labels=condclass.labels)
  zcond.cls <- cut(d$srccond, wellill.cuts.v, labels=condclass.labels)
  xhard.cls <- cut(d$srncond, difficulty.cuts.v, labels=condhard.labels)
  zhard.cls <- cut(d$srccond, difficulty.cuts.v, labels=condhard.labels)
  hard.cls <- cut(pmin(d$srncond, d$srccond), difficulty.cuts.v, labels=condhard.labels)

  return (list(xcond.cls=xcond.cls,
               zcond.cls=zcond.cls,
               xhard.cls=xhard.cls,
               hard.cls=hard.cls));
}

classify.errs <- function(d, n=100, prec="s") {
  strength.cuts.v <- strength.cuts(n,prec)

  d$zerr.cls <- cut(d$zerr,strength.cuts.v,labels=strength.lbls)
  d$zbnd.cls <- cut(d$zbnd,strength.cuts.v,labels=strength.lbls)
  d$xerr.cls <- cut(d$xerr,strength.cuts.v,labels=strength.lbls)
  d$xbnd.cls <- cut(d$xbnd,strength.cuts.v,labels=strength.lbls)

  return(d)
}

sanitize.file <- function (d, n=100, prec="s") {

  d$scheme <- factor(d$scheme, levels=scheme.levels, labels=scheme.labels)
  d$xstatus <- ordered(d$xstatus, levels=status.levels, labels=status.labels)
  d$zstatus <- ordered(d$zstatus, levels=status.levels, labels=status.labels)

  d$zerror <- pmin(pmax(d$zerr, workingeps(prec)),1.0)
  d$zbound <- pmin(pmax(d$zbnd, tiny.thresh(n,prec)),1.0)
  d$xerror <- pmin(pmax(d$xerr, workingeps(prec)),1.0)
  d$xbound <- pmin(pmax(d$xbnd, tiny.thresh(n,prec)),1.0)

  d$prec <- factor(prec)

  d <- categorize.estimates(d, n, prec)
  d <- classify.conds(d, n, prec)
  d <- classify.errs(d, n, prec)

  return(d) 
}

## readmatlab <- function (con, ntests.max = 100000, n=100) {
##   nfields <- length(protorow)
##   openstate <- isOpen(con)
##   if (!openstate) { open(con,"rb") }
##   tmp <- readMat(con)
##   if (!openstate) { close(con) }
##   ntests <- dim(tmp$data)[2]
##   if (ntests > ntests.max) {
##     tmp$data = tmp$data[,1:ntests.max]
##   }
##   d <- data.frame(t(tmp$data[savedoff,]))
##   names(d) <- savednames

##   for (i in 1:length(savednames)) {
##     d[,i] <- savedfn[[i]](d[,i])
##   }

##   return (sanitize.file(d,n=n))
## }

readraw <- function (con, ntests.max = 100000, n=100, prec="s") {
  nfields <- length(protorow)
  if (is.character(con)) { con <- file(con, "rb"); }
  openstate <- isOpen(con)
  if (!openstate) { open(con,"rb") }
  if (!isOpen(con)) { return(data.frame(narow)[FALSE,]) }
  v <- readBin(con,double(),n=nfields*ntests.max,size=4,endian="little")
  if (!openstate) { close(con) }

  ntests <- as.integer(length(v) / nfields)
  v <- matrix(data=v, nrow=ntests, ncol=nfields, byrow=TRUE)
  v <- v[,savedoff]
  d <- data.frame(v)
  names(d) <- savednames

  for (i in 1:length(savednames)) {
    d[,i] <- savedfn[[i]](d[,i])
  }

  d$n <- n
  d <- sanitize.file(d, n=n, prec=prec)
  
  return(d)
}

readraw8 <- function (con, ntests.max = 100000, n=100) {
  nfields <- length(protorow)
  openstate <- isOpen(con)
  if (!openstate) { open(con,"rb") }
  if (!isOpen(con)) { return(data.frame(narow)[FALSE,]) }
  v <- readBin(con,double(),n=nfields*ntests.max,size=8,endian="little")
  if (!openstate) { close(con) }

  ntests <- as.integer(length(v) / nfields)
  v <- matrix(data=v, nrow=ntests, ncol=nfields, byrow=TRUE)
  v <- v[,savedoff]
  d <- data.frame(v)
  names(d) <- savednames

#   for (i in 1:length(savednames)) {
#     d[,i] <- savedfn[[i]](d[,i])
#   }

  d = sanitize.file(d, n=n)
  
  return(d)
}

library(lattice)
#library(sm)

getcounts <- function(ds, var="x", ...) {
  err <- ds[,paste(var,"err.cls",sep="")]
  bnd <- ds[,paste(var,"bnd.cls",sep="")]
  M <- table(bnd,err,dnn=c("bound","error"))
  return(M)
}

errbrk.width <- function(prec="s") { return(switch(prec,s=0.25,d=0.5,NA)); }
errbrk.steps <- function(prec="s") {
  ew <- errbrk.width(prec);
  return (seq(log2(workingeps(prec))-4*ew, 4*ew, ew));
}
errbrks <- function(prec="s") { return (2^errbrk.steps(prec)); }
errbrk.corner <- function(prec="s") { es <- errbrk.steps(prec); return (es[-length(es)]); }
errbrks.corner <- function(prec="s") { eb <- errbrks(prec); return (eb[-length(eb)]); }
err.colorscale = c(  1,gray( seq(0.5,0,by=-0.01) )  )
err.ncolors = length(err.colorscale)
err.regions=hsv(seq(.95,1,length=100),seq(1,.3,length=100),seq(.3,.9,length=100))

panel.tbl <- function(x, y, z, breaks, logcnt=T, prec="s", ...) {
  if (!is.list(breaks)) {
    xcut <- cut(x, breaks=breaks)
    ycut <- cut(y, breaks=breaks)
  } else {
    xcut <- cut(x, breaks=breaks$x)
    ycut <- cut(y, breaks=breaks$y)
  }
  tbl <- table(xcut,ycut)
  if (logcnt) { 
    errbnd.tbl <- log10(errbnd.tbl)
    errbnd.tbl[errbnd.tbl<0] <- NA
  }
  ec <- errbrks.corner(prec)
  grid <- expand.grid(error=ec, bound=ec)
  grid$count <- as.vector(errbnd.tbl)
  levelplot(count ~ error * bound, data=grid, col.regions=col.regions, scales=list(log=T),
            panel=panel.errors, ... )
  
}

panel.errors <- function(x,y,z,zcol,...) {
#panel.levelplot(...,draw=FALSE)
  for (crd in c(-6, -4, -2, 0)) {
    panel.abline(h=crd,col=gray(0.8),lty=3,lwd=0.5)
    panel.abline(v=crd,col=gray(0.8),lty=3,lwd=0.5)
  }
  panel.abline(0,1,col=col)
  for (crd in c(2*2^-24*sqrt(100), 2^-12)) {
    panel.abline(h=log10(crd),col=gray(0.5))
    panel.abline(v=log10(crd),col=gray(0.5))
  }
  panel.levelplot(x=x,y=y,z=z,zcol=zcol,...)
#  lns <- contourLines(x,y,matrix(z,nrow=,nlevels=5)
#  llines(lns)
}

ploterrors <- function(ds, var="x", prec="s", full=T, col.regions=hsv(seq(0.9,1,length=100),seq(0.9,0.2,length=100),seq(0.2,0.9,length=100)), cuts=10, col="black", ...) {
  eb <- errbrks(prec)
  if (full) {
    err.cut <- cut(ds[,paste(var,"error",sep="")], breaks=eb)
    bnd.cut <- cut(ds[,paste(var,"bound",sep="")], breaks=eb)
  }
  else {
    err.cut <- cut(ds[,paste(var,"err",sep="")], breaks=eb)
    bnd.cut <- cut(ds[,paste(var,"bnd",sep="")], breaks=eb)
  }
  errbnd.tbl <- table(err.cut,bnd.cut)
  errbnd.tbl <- log10(errbnd.tbl)
  errbnd.tbl[errbnd.tbl<0] <- NA
  ec <- errbrks.corner(prec)
  grid <- expand.grid(error=ec, bound=ec)
  grid$count <- as.vector(errbnd.tbl)
  levelplot(count ~ error * bound, data=grid, col.regions=col.regions, scales=list(log=T),
            panel=panel.errors, prec=prec, ... )
}


# comp on x, y is xtrue ytrue ycomp xcomp and c