}

  } else {

    repmeans<-matrix(data=as.numeric(NA), ncol=NCOL(x), nrow=ncol(wts))

    for(i in 1:ncol(wts)){

      repmeans[i,]<-t(colSums(wts[,i]*x*pw)/sum(pw*wts[,i]))

    }

  }

  repmeans<-drop(repmeans)

  v <- svrVar(repmeans, scale, rscales,mse=design$mse, coef=rval)

}

# this is easy to fix, but for now...

return(v)

  attr(rval,"var") <-v

  attr(rval, "statistic")<-"mean"

  if (return.replicates){

    attr(repmeans,"scale")<-design$scale

    attr(repmeans,"rscales")<-design$rscales

    attr(repmeans,"mse")<-design$mse

    rval<-list(mean=rval, replicates=repmeans)

  }

  if (is.character(deff) || deff){

      nobs<-length(design$pweights)

      npop<-sum(design$pweights)

      vsrs<-unclass(svyvar(x,design,na.rm=na.rm, return.replicates=FALSE,estimate.only=TRUE))/length(design$pweights)

      if (deff!="replace")

        vsrs<-vsrs*(npop-nobs)/npop

      attr(rval,"deff") <- v/vsrs

  }

  class(rval)<-"svrepstat"

  rval

}

source("harness.R")

for (s in c("alabama", "california", "acs3yr" )) { 

  svydata <- readRDS(paste0(s,".rds"))

  names(svydata) <- tolower(names(svydata))

  print("loaded")

  svydsgn <- svrepdesign(

      weight = ~pwgtp ,

      repweights = 'pwgtp[1-9]' ,

      scale = 4 / 80 ,

      rscales = rep( 1 , 80 ) ,

      mse = TRUE ,

      data = svydata) 

  print("initialized")

  for (r in 1:5) {

    timing <- system.time({

      agep <- svymean(~agep, svydsgn, se=TRUE)

      print(agep)

      print(relp)

    })[[3]]

    log.result("survey", sys, conf, s, r, timing)

    clearResultRecycler()

  }

}