Add 1st version of the experiment prbs-multiroom-experiment1

.gitignore

+prbs-multiroom-experiment1/cache/
+
+prbs-multiroom-experiment1/genfig/
+
+prbs-multiroom-experiment1/outchannels/
+
+prbs-multiroom-experiment1/experiment.html
\ No newline at end of file

data/outchannels.csv

+ResampleVarname,channelid
+Room_VVX02_Omsk_C,ltechsp::VVX02::Omsk_C
+Room_VVX02_PV01_K,ltechsp::VVX02::PV01_K
+Room_VVX02_SpTF01,ltechsp::VVX02::SpTF01
+Room_VVX03_Omsk_C,ltechsp::VVX03::Omsk_C
+Room_VVX03_PV01_K,ltechsp::VVX03::PV01_K
+Room_VVX03_SpTF01,ltechsp::VVX03::SpTF01
+Room_VEN16_SpTI01,ltechsp::VEN16::SpTI01
+Room_C0.11_RoomTemperatureSetpoint,northqsp::houseid_2021_roomid_19558::roomTemperatureSetpoint
+Room_C0.12_RoomTemperatureSetpoint,northqsp::houseid_2021_roomid_19557::roomTemperatureSetpoint
+Room_C0.15_RoomTemperatureSetpoint,northqsp::houseid_2021_roomid_19556::roomTemperatureSetpoint
+Room_C0.14_RoomTemperatureSetpoint,northqsp::houseid_2021_roomid_19563::roomTemperatureSetpoint
+Room_C0.13_RoomTemperatureSetpoint,northqsp::houseid_2021_roomid_19559::roomTemperatureSetpoint
+Room_C0.01C_RoomTemperatureSetpoint,northqsp::houseid_2021_roomid_19582::roomTemperatureSetpoint
+Room_C0.04_RoomTemperatureSetpoint,northqsp::houseid_2021_roomid_19574::roomTemperatureSetpoint
+Room_C0.03_RoomTemperatureSetpoint,northqsp::houseid_2021_roomid_19573::roomTemperatureSetpoint
+Room_C0.02_RoomTemperatureSetpoint,northqsp::houseid_2021_roomid_19572::roomTemperatureSetpoint
+Room_C0.01B_RoomTemperatureSetpoint,northqsp::houseid_2021_roomid_19571::roomTemperatureSetpoint
+Room_C0.01_RoomTemperatureSetpoint,northqsp::houseid_2021_roomid_19570::roomTemperatureSetpoint
+Room_C0.07A_RoomTemperatureSetpoint,northqsp::houseid_2021_roomid_19562::roomTemperatureSetpoint
+Room_C0.07_RoomTemperatureSetpoint,northqsp::houseid_2021_roomid_19561::roomTemperatureSetpoint
+Room_C0.06_RoomTemperatureSetpoint,northqsp::houseid_2021_roomid_19560::roomTemperatureSetpoint
+Room_C0.08_RoomTemperatureSetpoint,northqsp::houseid_2021_roomid_19554::roomTemperatureSetpoint
+Room_C1.08_RoomTemperatureSetpoint,northqsp::houseid_2034_roomid_19566::roomTemperatureSetpoint
+Room_C1.09_RoomTemperatureSetpoint,northqsp::houseid_2034_roomid_19256::roomTemperatureSetpoint
+Room_C1.07_RoomTemperatureSetpoint,northqsp::houseid_2034_roomid_19567::roomTemperatureSetpoint
+Room_C2.04_RoomTemperatureSetpoint,northqsp::houseid_2034_roomid_19580::roomTemperatureSetpoint
+Room_C2.03_RoomTemperatureSetpoint,northqsp::houseid_2034_roomid_19579::roomTemperatureSetpoint
+Room_C2.01_RoomTemperatureSetpoint,northqsp::houseid_2034_roomid_19578::roomTemperatureSetpoint
+Room_C1.04_RoomTemperatureSetpoint,northqsp::houseid_2034_roomid_19581::roomTemperatureSetpoint
+Room_C1.00_RoomTemperatureSetpoint,northqsp::houseid_2034_roomid_19575::roomTemperatureSetpoint
+Room_C1.03_RoomTemperatureSetpoint,northqsp::houseid_2034_roomid_19577::roomTemperatureSetpoint
+Room_C1.02_RoomTemperatureSetpoint,northqsp::houseid_2034_roomid_19576::roomTemperatureSetpoint

prbs-multiroom-experiment1/experiment.Rmd

+---
+title: "prbs-multiroom-experiment1"
+output:
+  html_document:
+    toc: true
+    toc_float: true
+---
+
+```{r setup, include=FALSE, cache=FALSE, purl=FALSE}
+knitr::opts_chunk$set(echo = TRUE, results = 'hide', message = FALSE, warning = FALSE, fig.height = 4, fig.width = 14, fig.path = 'genfig/', cache.path = 'cache/', cache = TRUE, cache.rebuild = FALSE, dpi=100)
+```
+
+```{r init, echo=FALSE, cache=FALSE}
+##----------------------------------------------------------------
+## Set the working directory
+setwd(".")
+## Source the files in the "functions" folder
+files <- dir("functions",full.names=TRUE)
+for(i in 1:length(files)) source(files[i])
+##----------------------------------------------------------------
+```
+
+# Introduction {.tabset}
+
+In this experiment the set points in each room are controlled following a Pseudo Random Binary Sequence (PRBS). This enables a controlled experiment where the control variables (room temperature setpoints) will be set independently (uncorrelated) to other variables (like time of day, external temperature, etc.). Further, by lagging the PRBS signals, each the change of set points in each room will be independent.
+
+
+## PRBS explanaition
+
+Here a small explanaition of how the PRBS works and a reference to a paper about it.
+
+```{r prbsacf}
+##----------------------------------------------------------------
+## Generate a PRBS signal
+## Use the function defined in the file "functions/prbs.R", which generates a PRBS signal
+## - n is the length of the register
+## - initReg just needs to be some initial value of 1,2,...
+##   it is the initial value of the registers and therefore only
+##   determines the start of the cycle
+## - lambda is the length of the smallest period in
+##   which the signal can change, given in samples
+x <- prbs(n=6)
+acf(x, lag.max=length(x))
+##  prbs
+plot(x,type="s")
+acf(x, lag.max=length(x))
+## See  to PRBS
+x <- prbs(n=6)
+## LagWithCycling to the right
+y <- lagWithCycling(x, lag=10)
+## cross cor.
+ccf(x,y, lag.max=length(x))
+```
+
+```{r prbshist}
+##----------------------------------------------------------------
+## Generate the PRBS signals for the PRBS1 experiment, where a single signal controls all heaters in the building
+## PRBS med n=6, lambda=4
+##   - Smallest period in one state for 5 minute sample period is then 4*5min=20min
+##   - Settling time of the system (T_s in the Godfrey 1980 paper) below the period (T_0 in the paper) in: lambda * (2^n-1) * 5 / 60 = 21 hours
+n <- 5
+lambda <- 4
+(x <- prbs(n, 37, lambda))
+length(x) / 24
+(xper <- c(0,cumsum(abs(diff(x)))))
+perlens <- unlist(lapply(split(xper,xper), length))
+table(perlens)
+hist(perlens, breaks=(0:n)*lambda+lambda*0.5)
+```
+
+```{r prbsplotts}
+plot(x, type="s")
+```
+
+
+## PRBS different levels
+
+Combine periods with different temperature levels.
+
+```{r prbsplotts2}
+## Create sequence
+testseq <- function(n, lambda, nhours, starthour, Tmin, Tmax, plotit=TRUE){
+    x <- prbs(n, 42, lambda)
+    t1 <- starthour + ((1:length(x))[abs(c(0,diff(x)))==1] / length(x)) * nhours
+    ## Round to 5 min steps
+    t1 <- round(t1*12) / 12
+    val1 <- rep(c(Tmin,Tmax), length(t1))[1:length(t1)]
+    ##
+    if(plotit){ plot(t1, val1, type="s") }
+    ##
+    return(data.frame(t=t1, Tset=val1))
+}
+## val <- testseq(n=5, lambda=1, nhours=96, Tmin=15, Tmax=24)
+## diff(val$t)
+
+## Four days between 10 and 25 degrees C
+val1 <- testseq(n=5, lambda=1, nhours=96, starthour=0, Tmin=10, Tmax=25)
+## Two days between 17 and 20 degrees C
+val2 <- testseq(n=5, lambda=1, nhours=48, starthour=96, Tmin=17, Tmax=20)
+## Two days between 20 and 22 degrees C
+val3 <- testseq(n=5, lambda=2, nhours=48, starthour=96+48, Tmin=20, Tmax=22)
+## One day between 21 and 22 degrees C
+val4 <- testseq(n=5, lambda=1, nhours=24, starthour=96+48+48, Tmin=21, Tmax=22)
+
+## Combined
+X <- rbind(val1, val2, val3, val4)
+plot(X$t/24, X$Tset, type="s")
+```
+
+
+## Lagged for each room
+
+Lag for each room to make them independent. Note that on each floor they are repeated, this should be such that each neighbouring room, also on top of each other, are lagged differently.
+
+```{r outputs}
+nms <- read.table("../data/outchannels.csv", header=TRUE, sep=",", as.is=TRUE)$ResampleVarname
+nms <- nms[grep("^Room_C", nms)]
+rooms <- sapply(strsplit(nms, "_"), function(x){x[2]})
+
+rooms <- rooms[order(as.numeric(gsub("[[:alpha:]]","",rooms)))]
+
+lagval <- numeric(length(rooms))
+lagval[rooms=="C0.00"] <- 0
+lagval[rooms=="C0.00B"] <- 0
+lagval[rooms=="C0.01"] <- 1
+lagval[rooms=="C0.02"] <- 1
+lagval[rooms=="C0.03"] <- 1
+lagval[rooms=="C0.04"] <- 1
+lagval[rooms=="C0.05"] <- 0#?
+lagval[rooms=="C0.06"] <- 2
+lagval[rooms=="C0.07"] <- 3
+lagval[rooms=="C0.07A"] <- 3
+lagval[rooms=="C0.08"] <- 4
+lagval[rooms=="C0.10"] <- 0#?
+lagval[rooms=="C0.11"] <- 5
+lagval[rooms=="C0.12"] <- 7
+lagval[rooms=="C0.13"] <- 6
+lagval[rooms=="C0.14"] <- 7
+lagval[rooms=="C0.15"] <- 8
+lagval[rooms=="C1.00"] <- 3
+lagval[rooms=="C1.02"] <- 4
+lagval[rooms=="C1.03"] <- 5
+lagval[rooms=="C1.04"] <- 6
+lagval[rooms=="C1.07"] <- 7
+lagval[rooms=="C1.08"] <- 0
+lagval[rooms=="C1.09"] <- 1
+lagval[rooms=="C1.12"] <- 2
+lagval[rooms=="C1.13"] <- 3
+lagval[rooms=="C2.01"] <- 6
+lagval[rooms=="C2.03"] <- 7
+lagval[rooms=="C2.04"] <- 8
+
+
+lagit <- function(val, nhours){
+    maxk <- max(lagval)
+    starthour <- floor(val$t[1] / 24) * 24
+    lapply(1:length(rooms), function(i){
+        k <- lagval[i]
+        val$t <- (val$t + k/maxk * nhours) %% nhours
+        val$t <- val$t + starthour
+        return(val[order(val$t), ])
+    })
+}
+
+L1 <- lagit(val1, 96)
+L2 <- lagit(val2, 48)
+L3 <- lagit(val3, 48)
+L4 <- lagit(val4, 24)
+
+
+L <- lapply(1:length(rooms), function(i){
+    rbind(L1[[i]], L2[[i]], L3[[i]], L4[[i]])
+})
+```
+
+```{r plotrooms}
+##setpar("ts", mfrow=c(5,5))
+lapply(1:length(L), function(i){
+    x <- L[[i]]
+    plot(x$t, x$Tset, type="s", ylab="")
+    title(main=rooms[i], line=-1)
+})
+##dev.copy2pdf(file="setpoints/setpoints.pdf")
+```
+
+
+## Write to files 
+
+```{r write-to-files}
+lapply(1:length(L), function(i){
+   x <- L[[i]]
+   ## Set starting date and time in hours
+   x$t <- asp("2019-04-13 00:00") + x$t * 3600
+   ## Write to separate .csv files
+   write.table(x, file=pst("outchannels/",rooms[i],".csv"), sep=",", row.names=FALSE)
+})
+```
+
+```{r write-the-R-code, purl=FALSE}
+##library(knitr)
+##purl("experiment.Rmd")
+```

prbs-multiroom-experiment1/functions/asp.R

+## Define a generic method
+asp <- function(object, ...){
+    UseMethod("asp")
+}
+
+asp.default <- function(object){
+    return(object)
+}
+
+asp.character <- function(object, tz = "GMT", ...){
+    as.POSIXct(object, tz = tz, ...)
+}
+
+asp.POSIXct <- function(object){
+    object
+}
+
+asp.POSIXlt <- function(object){
+    as.POSIXct(object)
+}
+
+asp.numeric <- function(object){
+    ISOdate(1970, 1, 1, 0) + object
+}

prbs-multiroom-experiment1/functions/lagWithCycling.R

+## Shift the values in vector x, lag steps to the right
+lagWithCycling <- function(x, lag)
+{
+    if(lag == 0){ return(x) }
+    n <- length(x)
+    temp <- x[(lag+1):n]
+    x[(n-(lag-1)):n] <- x[1:lag]
+    x[1:(n-lag)] <- temp
+    x
+  }
+

prbs-multiroom-experiment1/functions/nams.R

+nams <- function(x) {
+    if(is.matrix(x)){
+        colnames(x)
+    } else {
+        names(x)
+    }
+}
+
+
+`nams<-` <- function(x, value) {
+    if(is.matrix(x)){
+        colnames(x) <- value
+    } else {
+        names(x) <- value
+    }
+  x
+}

prbs-multiroom-experiment1/functions/plot_ts.R

+## Define a generic method
+plot_ts <- function(object, ...){
+    UseMethod("plot_ts")
+}
+
+plot_ts.data.list <- function(object, patterns, tstart = NA, tend = NA, kseq = 0, colorpalette = NA, ...) {
+    DL <- object
+    ## Plot it for a period
+    oldpar <- setpar("ts")
+    on.exit(par(oldpar))
+    ##
+    if(is.na(tstart)) { tstart <- DL$t[1] - 1 }
+    if(is.na(tend)) { tend <- DL$t[length(DL$t)] }
+    DL <- subset(DL, period(tstart, DL$t, tend))
+    ##
+    ## Generate a data.frame with the series to be plotted
+    X <- lapply_cbind_df(patterns, function(pattern) {
+        ## Find the columns to plot
+        nms <- grep(pattern, names(DL), value = TRUE)
+        if(length(nms) == 0){
+            warning(pst("No names where found matching: ", pattern))
+            tmp <- as.data.frame(matrix(NA,nrow=length(DL$t),ncol=1))
+            names(tmp)[1] <- pattern
+            return(tmp)
+        }else{
+            ## Do the plotting
+            do.call("cbind", lapply(nms, function(nm){
+                if(is.null(dim(DL[[nm]]))) {
+                    ## It is a vector
+                    X <- data.frame(DL[[nm]])
+                    names(X) <- nm
+                    return(X)
+                } else {
+                    ## Its a matrix
+                    ## Find the columns with k and digits
+                    i <- grep("^k[[:digit:]]+$", nams(DL[[nm]]))
+                    if(length(i) > 0) {
+                        ## Try to return the kseq forecast columns
+                        ik <- which(pst("k",kseq) %in% nams(DL[[nm]]))
+                        X <- lag(DL[[nm]][ ,pst("k",kseq[ik])], lag = kseq[ik])
+                        if(is.null(dim(X))) {
+                            X <- as.data.frame(X)
+                            names(X) <- pst("k",kseq)
+                        }
+                    } else {
+                        ## Just return all
+                        X <- DL[[nm]]
+                    }
+                    nams(X) <- pst(nm, "_", nams(X))
+                    return(X)
+                }
+            }))
+        }
+    })
+    if(any(duplicated(nams(X)))){
+        X <- X[ ,unique(nams(X))]
+    }
+    X$t <- DL$t
+    ## Use the plot_ts function which takes the data.frame
+    nms <- unlist(getse(strsplit(nams(X), "_k"), 1))
+    plot_ts.data.frame(X, patterns, colorpalette = colorpalette, nms = nms, ...)
+}
+
+
+## Plot all with prefix
+plot_ts.data.frame <- function(object, patterns, xnm = "t", draw_grid = TRUE, xlabs = NA, ylabs = NA,
+    space_for_legend = 0.25, cex.legend = 1, xaxis_format = NA, colorpalette = NA,
+    ylim = NA, main = "", nms = NA, xat = NA, xnticks = 8, ...) {
+    ##
+    data <- object
+    patterns <- patterns[patterns != xnm]
+    ##
+    oldpar <- setpar("ts", mfrow = c(length(patterns), 1))
+    on.exit(par(oldpar))
+    ##
+    for (i in 1:length(patterns)) {
+        if (is.na(ylim[1])) {
+            ylim <- NA
+        } else {
+            ylim <- ylim[[i]]
+        }
+        if (class(colorpalette) == "list") {
+            colorpalette <- colorpalette[[i]]
+        }
+        ##
+        plot_ts_series(data, patterns[i], draw_grid, xnm, space_for_legend = space_for_legend,
+            cex.legend = cex.legend, colorpalette, ylim, main[i], nms = nms, ...)
+        if (!is.na(ylabs[1]))
+            title(ylab = ylabs[i], yaxt = "s")
+    }
+    if (any(nams(data) == xnm)) {
+        if (class(data[ ,xnm])[1] == "numeric") {
+            axis(1, data[ ,xnm], xaxt = "s")
+        } else {
+            ## makes too few ticks: axis.POSIXct(1, data[ ,xnm], format = xaxis_format, xaxt = "s")
+            if(is.na(xat[1])){ xat <- pretty(data[ ,xnm]) }
+            if(is.na(xaxis_format)){
+                if( all(as.numeric(xat,unit="secs") %% (24*3600) == 0) ){
+                    xaxis_format <- "%Y-%m-%d"
+                }else{
+                    xaxis_format <- "%Y-%m-%d %H:%M"
+                }
+            }
+            axis.POSIXct(1, data[ ,xnm], at = xat, format = xaxis_format, xaxt = "s")
+        }
+    } else {
+        axis(1, 1:nrow(data), xaxt = "s")
+    }
+}
+
+plot_ts.matrix <- plot_ts.data.frame
+
+## Plot all columns found with regex pattern
+plot_ts_series <- function(data, pattern, draw_grid, xnm, space_for_legend = 0.25, cex.legend = 1,
+                           colorpalette = NA, ylim = NA, main = "", nms = NA, ...) {
+    ## Use these names when finding columns to plot
+    if(is.na(nms[1])) {
+        nms <- nams(data)
+    }
+    iseq <- integer(0)
+    for (pf in strsplit(pattern, "\\|")[[1]]) {
+        iseq <- c(iseq, grep(pf, nms))
+    }
+    ## Check if xnm is in the data
+    if (any(nms == xnm)) {
+        ##iseq <- iseq[-which(nms[iseq] == xnm)]
+        x <- data[, xnm]
+    } else {
+        x <- 1:nrow(data)
+    }
+    ## Limits on y-axis
+    if (is.na(ylim[1])) {
+        ylim <- range(data[, iseq], na.rm = TRUE)
+    }
+    ##
+    if (any(is.infinite(ylim))){
+        ## If no values then plot empty
+        plot(x, x, type = "n")
+        legend("topright", paste0(nams(data)[iseq], ": NA"))
+    } else {
+        plot(x, x, type = "n", xlim = c(min(x), max(x) +
+                                                diff(range(x)) * space_for_legend), ylim = ylim, yaxt = "n", bty = "n",
+             xlab = "", ylab = "", ...)
+        yat <- pretty(scalerange(data[, iseq], 0.2))
+        axis(2, yat)
+        ## Grid
+        if(draw_grid){
+            abline(v = pretty(x), h = yat, col = "grey85", lty = 2)
+        }
+        box()
+        ##
+        if (is.na(colorpalette[1])) {
+            colorramp <- colorRampPalette(c("black","cyan","purple","blue","red","green"))
+            ##colorRampPalette(c("#00007F", "blue", "#007FFF", "cyan",
+            ## "#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000"))
+        } else {
+            colorramp <- colorRampPalette(colorpalette)
+        }
+        ##
+        colormap <- colorramp(length(iseq))
+        for (i in 1:length(iseq)) {
+            lines(x, data[, iseq[i]], col = colormap[i], ...)
+        }
+        rng <- do.call("rbind", lapply(1:length(iseq), function(i) {
+            paste(format(range(data[, iseq[i]], na.rm = TRUE), digits = 2), collapse = " to ")
+        }))
+        legend("topright", paste0(nams(data)[iseq], ": ", rng), lty = 1, col = colormap,
+               cex = cex.legend, bg="white", ...)
+    }
+    title(main = main, line = -1)
+    invisible(iseq)
+}
+
+
+## DL <- list(x=10:1, Y=data.frame(k1=c(1:10),k2=c(11:20)))
+## class(DL) <- "data.list"
+## plot_ts(DL)
+
+
+
+## Plot ts with ggplot2
+## ## Plot all with prefix plotmultigg <-
+## function(X,patterns,grid_delta_t=NA,xlabs=NA,ylabs=NA,...)  { ## Plot
+## require('ggplot2') require('reshape2') require('gridExtra')
+
+## ## Multiple with same y-axis ## Get the data to plot
+
+
+## length(patterns)
+
+## for(i in 1:length(patterns)) { ## Melt the data to plot, with 't' and the
+## columns tmp <- melt(X[ ,c('t',grep(patterns[i],nams(X)))], 't')
+
+## X <- Xplot() ## Only the selected boxes X <- X[X$id %in% selectedBoxIds(), ] ##
+## Lplot <- list() i <- 0 for(var in input$vars){ ## tmp <- X[X$variable==var, ]
+## if(nrow(tmp)>0){ ## By color i <- i + 1 Lplot[[i]] <- ggplot(tmp, aes(x=t,
+## y=value, color=id)) + labs(title=var) + geom_line() } } ml <-
+## marrangeGrob(Lplot, nrow=5, ncol=1, top='') print(ml)
+
+## setpar('ts',mfrow=c(length(patterns),1)) for(i in 1:length(patterns)) {
+## plotseries(X,patterns[i],grid_delta_t,...)  if(!is.na(ylabs[1])) title(ylab=ylabs[i],
+## yaxt='s') } axis.POSIXct(1, X$t,format=c('%m-%d %H:%M'), xaxt='s') }
+

prbs-multiroom-experiment1/functions/prbs.R

+prbs <- function(n, initReg=666, lambda=1)
+  {
+    ## Function for generating PRBS sequences.
+    ## - n is the length of the register
+    ## - initReg just needs to be some initial value of 1,2,...
+    ##   it is the initial value of the registers and therefore only
+    ##   determines the start of the cycle
+    ## - lambda is the length of the smallest period in
+    ##   which the signal can change, given in samples
+    ##
+    ## Check the input
+    if(n < 2 | n > 11){ stop("n must be between 1 and 11") }
+    ## Do init
+    reg <- intToBits(as.integer(initReg))
+    print(reg)
+    x <- vector()
+    N <- 2^n - 1
+    ## Do the shift according to the value of n
+    for(i in 1:N)
+      {
+        ## Make the xor operation according to Godfrey80
+        if(n <= 4 | n == 6){ reg[n+1] <- xor(reg[1],reg[2]) }
+        else if(n == 5 | n == 11){ reg[n+1] <- xor(reg[1],reg[3]) }
+        else if(n == 7 | n == 10){ reg[n+1] <- xor(reg[1],reg[4]) }
+        else if(n == 8){ reg[n+1] <- as.raw(sum(as.integer(c(reg[1],reg[5],reg[6],reg[7])))%%2) }
+        else if(n == 9){ reg[n+1] <- xor(reg[1],reg[5]) }
+        ## Keep the value of the first position in the register
+        x[i] <- as.integer(reg[1])
+        ## Shift the register
+        reg[1:n] <- reg[2:(n+1)]
+      }
+    ## Return x with each element repeated lambda times
+    rep(x, rep(lambda,N))
+  }

prbs-multiroom-experiment1/functions/pst.R

+## Paste helper functions
+pst <- function(...) {
+    paste0(...)
+}

prbs-multiroom-experiment1/functions/setpar.R

+setpar <- function(tmpl = NA, ...) {
+    ## Get par list
+    p <- par(no.readonly = TRUE)
+    ## Templates
+    if (!is.na(tmpl)) {
+        if (tmpl == "ts") {
+            par(mfrow = c(3, 1), oma = c(3, 0, 2, 2), mar = c(0, 4, 0.5, 0), xaxt = "n", 
+                mgp = c(2.2, 0.4, 0), tcl = -0.4)
+        }
+        if (tmpl == "pdf") {
+            par(mar = c(4, 4, 1, 1), mgp = c(2.2, 0.7, 0), tcl = -0.4)
+        }
+    }
+    ## Replace all the parameters given in prm Get only the ... parameters
+    i <- which(!nams(match.call()) %in% nams(match.call(expand.dots = FALSE)))
+    if (length(i) > 0) {
+        par(...)
+        ## prm <- as.list(match.call()[i]) p <- list() for(i in 1:length(prm)) { p$new <-
+        ## eval(prm[[i]]) nams(p)[i] <- nams(prm)[i] } par(p)
+    }
+    ## Set par and return the original par options(warn = (-1)) options(warn = 1)
+    invisible(p)
+}