diff --git a/DESCRIPTION b/DESCRIPTION
new file mode 100644
index 0000000..5f9b53f
--- /dev/null
+++ b/DESCRIPTION
@@ -0,0 +1,15 @@
+Package: plrs
+Version: 0.99.0
+Type: Package
+Title: Piecewise Linear Regression Splines (PLRS) for the association
+        between DNA copy number and gene expression
+Author: Gwenael G.R. Leday
+Maintainer: Gwenael G.R. Leday to <gleday@few.vu.nl>
+Depends: R (>= 2.10), Biobase
+Imports: BiocGenerics, CGHbase, graphics, grDevices, ic.infer, marray,
+         methods, quadprog, Rcsdp, stats, stats4, utils
+Suggests: mvtnorm, methods
+Description: The present package implements a flexible framework for
+    modeling the relationship between DNA copy number and gene
+    expression data using Piecewise Linear Regression Splines (PLRS).
+License: GPL (>=2.0)
diff --git a/NAMESPACE b/NAMESPACE
new file mode 100644
index 0000000..f2b7dcc
--- /dev/null
+++ b/NAMESPACE
@@ -0,0 +1,43 @@
+import(methods)
+
+importMethodsFrom(marray, lines, points)
+
+importFrom(BiocGenerics, residuals)
+importFrom(Biobase, exprs)
+importFrom(CGHbase, segmented, probloss, probnorm, probgain, probamp)
+importFrom(graphics, abline, box, polygon)
+importFrom(grDevices, dev.off, jpeg, pdf)
+importFrom(ic.infer, ic.weights, pbetabar)
+importFrom(quadprog, solve.QP)
+importFrom(Rcsdp, csdp)
+importFrom(stats, lm, p.adjust, pf, quantile, rnorm, runif, vcov)
+importFrom(stats4, coef, summary, plot)
+importFrom(utils, combn, head)
+
+export("plrs.sim","plrs.cb","plrs.series")
+exportClasses(
+     "plrs",
+     "plrs.select",
+     "plrs.series" 
+)
+exportMethods(
+    "coef",
+    "fitted",
+    "effects",
+    "knots",
+    "model.matrix",
+    "plot",
+    "predict",
+    "print",
+    "residuals",
+    "show",
+    "summary"
+)
+export(
+	"criteria",
+	"modify.conf",
+	"plrs",
+	"plrs.select",
+	"plrs.test"
+)
+
diff --git a/R/AllClasses.r b/R/AllClasses.r
new file mode 100644
index 0000000..25cbfa6
--- /dev/null
+++ b/R/AllClasses.r
@@ -0,0 +1,33 @@
+################
+## All Classes
+################
+
+setClass("plrs",
+         representation(coefficients = "numeric",
+                        fitted.values = "numeric",
+                        residuals = "numeric",
+                        X = "matrix",
+                        data = "list",
+                        mdata = "list",
+                        QP = "list",
+                        test = "list",
+                        cb = "list",
+                        selected = "logical",
+                        type = "character",
+                        call.arg = "list")
+)
+setClass("plrs.select",
+         representation(table = "matrix",
+                        model = "plrs",
+                        crit = "character")
+)
+setClass("plrs.series",
+         representation(coefficients = "matrix",
+                        effects = "list",
+                        test = "matrix",
+                        general = "matrix",
+                        modelsType = "list",
+                        call.arg = "list")
+)
+
+
diff --git a/R/AllGenerics.r b/R/AllGenerics.r
new file mode 100644
index 0000000..c601d33
--- /dev/null
+++ b/R/AllGenerics.r
@@ -0,0 +1,6 @@
+################
+## All Generics
+################
+
+setGeneric("fitted",    function(object,...) standardGeneric("fitted"))
+setGeneric("predict",   function(object,...) standardGeneric("predict"))
diff --git a/R/Bmat.r b/R/Bmat.r
new file mode 100644
index 0000000..8fc10ff
--- /dev/null
+++ b/R/Bmat.r
@@ -0,0 +1,47 @@
+# Internal function
+# Build design matrix with truncated power basis functions
+
+# Author: Gwenael G.R. Leday
+
+.Bmat <- function(cghseg, knots, continuous=TRUE, general.intercept=TRUE){
+
+	nKnots <- length(knots)
+	nval <- nKnots + 1
+	xx <- cghseg
+	mat <- matrix(cghseg)
+	if(general.intercept) mat <- cbind(rep(1, length(xx)) ,mat)
+
+	# Truncated p-th power function
+	tpower <- function(x, t, p) (x - t) ^ p * (x > t)
+
+	if(continuous){
+		if(nKnots>=1) mat <- cbind(mat, tpower(xx, knots[1], 1))
+		if(nKnots>=2) mat <- cbind(mat, tpower(xx, knots[2], 1))
+		if(nKnots>=3) mat <- cbind(mat, tpower(xx, knots[3], 1))
+		nval2 <- 0
+	}
+	else{
+		if(nKnots>=1){
+			mat <- cbind(mat, tpower(xx, knots[1], 0))
+			mat <- cbind(mat, tpower(xx, knots[1], 1))
+		}
+		if(nKnots>=2){
+			mat <- cbind(mat, tpower(xx, knots[2], 0))
+			mat <- cbind(mat, tpower(xx, knots[2], 1))
+		}
+		if(nKnots>=3){
+			mat <- cbind(mat, tpower(xx, knots[3], 0))
+			mat <- cbind(mat, tpower(xx, knots[3], 1))
+		}
+		nval2 <- nval - 1
+	}
+
+	# Labels
+	theta.1 <- paste(paste(paste("theta",0:(nval-1),sep=""),".",sep=""),1,sep="")
+	if(!general.intercept & nval2==0) theta.0 <- NULL
+	else{ theta.0 <- paste(paste(paste("theta",0:nval2,sep=""),".",sep=""),0,sep="")}
+	colnames(mat) <- sort(c(theta.0, theta.1))
+
+	return(mat)
+}
+
diff --git a/R/convertToTime.r b/R/convertToTime.r
new file mode 100644
index 0000000..b776416
--- /dev/null
+++ b/R/convertToTime.r
@@ -0,0 +1,13 @@
+# Internal function
+# Convert seconds into a "HH:MM:SS" format
+
+# Author: Gwenael G.R. Leday
+
+.convertToTime <- function(x){
+	h <- as.character(x%/%3600)
+	m <- as.character((x%%3600)%/%60)
+	s <- as.character(round((x%%3600)%%60))
+	if(nchar(m)==1) m <- paste(0,m,sep="")
+	if(nchar(s)==1) s <- paste(0,s,sep="")
+	return(paste(h,m,s,sep=":"))
+}
\ No newline at end of file
diff --git a/R/criteria.r b/R/criteria.r
new file mode 100644
index 0000000..d26507e
--- /dev/null
+++ b/R/criteria.r
@@ -0,0 +1,77 @@
+# Extract AIC, AICC, BIC and OSAIC from a 'plrs' object
+
+# Author: Gwenael G.R. Leday
+
+criteria <- function(obj, crit = "all"){
+
+	# Error handling
+	if(class(obj)!="plrs") stop("An object of class \"plrs\" is required")
+	if(length(crit)==0) stop("No criterion")
+	if(length(crit)==1) if(crit=="all") crit <- c("aic","aicc","bic","osaic")
+	if(!is.character(crit)) stop("'crit' has to be a character")
+	if(!"all"%in%crit & !"aic"%in%crit & !"aicc"%in%crit &
+         !"bic"%in%crit & !"osaic"%in%crit) stop("Wrong criterion")
+	bool <- TRUE
+	if("osaic"%in%crit){
+		if(!obj@call.arg$constr) bool <- FALSE
+		else{
+			if(length(coef(obj))!=1){
+				if(obj@QP$meq == nrow(t(obj@QP$Amat))) bool <- FALSE
+			}
+			else{
+				bool <- FALSE
+			}
+		}
+	}
+
+	rss <- sum(residuals(obj)^2)
+	n <- length(obj@data$expr)
+
+	# Compute criteria for the segmented model
+	if(rss != 0){
+		if(obj@call.arg$constr & length(coef(obj))>1) p <- length(coef(obj)) - obj@QP$meq
+		else p <- length(coef(obj))
+
+		if("aic" %in% crit || "aicc" %in% crit || !bool) aic <- n*log(rss/n) + 2*p
+		if("aicc" %in% crit) aicc <- aic + ((2*p)*(p+1))/(n-p-1)
+		if("bic" %in% crit) bic <- n*log(rss/n) + log(n)*p
+		if("osaic" %in% crit){ # As defined in Hughes et al (2003)
+			if(bool){
+				library(mvtnorm)
+				# Calculate weights
+				obj.unconstr <- .plrs.fit(X=obj@X, matconstr=NULL, mytype=obj@type, modelFull=obj, callArg=obj@call.arg)
+				cov.mat <- obj.unconstr@QP$vcov/obj.unconstr@QP$S2
+				uiw <- t(obj@QP$Amat)
+
+				# Aknowledgement: Taken from package ic.infer
+				if(obj@QP$meq==0){
+					#wt.bar <- ic.weights(uiw %*% cov.mat %*% t(uiw))
+					wt.bar <- ic.weights(crossprod(t(crossprod(t(uiw),cov.mat)),t(uiw)))
+				}
+          	 	 	else{
+					wt.bar <- ic.weights(solve(solve(crossprod(t(crossprod(t(uiw),cov.mat)),t(uiw)))[-(1:obj@QP$meq), -(1:obj@QP$meq)]))
+				}
+				wt.bar <- wt.bar[sort(names(wt.bar))]
+				r <- (length(wt.bar)-1)
+				osaic.npar <- sum(wt.bar*(p-r+0:r))
+				osaic <- n*log(rss/n) + 2*osaic.npar
+			}
+			else{
+				osaic <- aic
+			}
+		}
+	}
+	else{
+		if("aic" %in% crit) aic <- -Inf
+		if("aicc" %in% crit) aicc <- -Inf
+		if("bic" %in% crit) bic <- -Inf
+		if("osaic" %in% crit) osaic <- -Inf
+	}
+
+	out <- NULL
+	if("aic" %in% crit) out <- c(out, list("aic" = aic))
+	if("aicc" %in% crit) out <- c(out, list("aicc" = aicc))
+	if("bic" %in% crit) out <- c(out, list("bic" = bic))
+	if("osaic" %in% crit) out <- c(out, list("osaic" = osaic))
+	return(out)
+}
diff --git a/R/estimateKnots.r b/R/estimateKnots.r
new file mode 100644
index 0000000..14873c4
--- /dev/null
+++ b/R/estimateKnots.r
@@ -0,0 +1,46 @@
+# Estimate knots
+
+# Author: Gwenael G.R. Leday
+
+.estimateKnots <- function(cghseg, cghcall, probloss=NULL, probnorm=NULL, probgain=NULL, probamp=NULL){
+
+	val <- sort(unique(cghcall))
+	nval <- length(val)
+	knots <- NULL
+	
+	# Midpoints of intermediate intervals
+	if(is.null(probloss) || is.null(probnorm) || is.null(probgain)){
+		if(nval>1){
+			for(k in 1:(nval-1)){
+				knots <- c(knots,(max(cghseg[cghcall==val[k]]) +  min(cghseg[cghcall==val[k+1]]))/2)
+			}
+		}
+	}
+	# Estimate using call probabilities
+	else{
+		p <- function(alpha, cghseg, leftState, rightState){
+			ind1 <- as.numeric(cghseg <= alpha)
+			ind2 <- 1-ind1
+			as.numeric(ind1%*%leftState + ind2%*%rightState)
+		}
+		for(j in 1:(length(val)-1)){
+			# Reduce search to adjacent states range
+			lb <- min(cghseg[cghcall==val[j]])
+			ub <- max(cghseg[cghcall==val[j+1]])
+
+			# Search maximum
+			xx <- seq(from=lb, to=ub, length.out=1000)
+			if(val[j]==-1) tomax <- sapply(xx, p, cghseg=cghseg, leftState=probloss, rightState=probnorm)
+			if(val[j]==0) tomax <- sapply(xx, p, cghseg=cghseg, leftState=probnorm, rightState=probgain)
+			if(val[j]==1) tomax <- sapply(xx, p, cghseg=cghseg, leftState=probgain, rightState=probamp)
+
+			# Take midpoints of interval
+			optvals <- xx[which(tomax==max(tomax))]
+			leftpts <- max(cghseg[cghseg<=min(optvals)])
+			rightpts <- min(cghseg[cghseg>min(optvals)])
+
+			knots <- c(knots, mean(c(leftpts, rightpts)))
+		}
+	}
+	return(knots)
+}
\ No newline at end of file
diff --git a/R/make.mat.constr.r b/R/make.mat.constr.r
new file mode 100644
index 0000000..e63dad2
--- /dev/null
+++ b/R/make.mat.constr.r
@@ -0,0 +1,190 @@
+# Internal function
+# Build matrix of constraints
+
+# Author: Gwenael G.R. Leday
+
+.make.mat.constr <- function(constr.slopes,
+                             constr.intercepts,
+                             continuous,
+                             val){
+
+	nval <- length(val)
+	
+	if(nval==2){
+		if(all(val==c(-1,1))|all(val==c(-1,2))|all(val==c(1,2))) constr.slopes <- 1 
+		if(constr.slopes==2){
+			if(all(val==c(-1,0))){
+				if(!continuous){
+					if(constr.intercepts){
+						mat <- matrix(c(0,1,0,1,
+	                                              0,0,0,-1,
+	                                              0,0,1,0),3,4, byrow=TRUE)
+					}
+					else{
+						mat <- matrix(c(0,1,0,1,
+	                                              0,0,0,-1),2,4, byrow=TRUE)
+					}
+				}
+				else{
+						mat <- matrix(c(0,1,1,
+	                                              0,0,-1),2,3, byrow=TRUE)
+				}
+			}
+			if(all(val==c(0,1))|all(val==c(0,2))){
+				if(!continuous){
+					if(constr.intercepts){
+						mat <- matrix(c(0,1,0,0,
+	                                              0,0,0,1,
+	                                              0,0,1,0),3,4, byrow=TRUE)
+					}
+					else{
+						mat <- matrix(c(0,1,0,0,
+	                                              0,0,0,1),2,4, byrow=TRUE)
+					}
+				}
+				else{
+						mat <- matrix(c(0,1,0,
+	                                              0,0,1),2,3, byrow=TRUE)
+				}
+			}
+		}
+		if(constr.slopes==1){
+			if(!continuous){
+				if(constr.intercepts){
+					mat <- matrix(c(0,1,0,0,
+                                              0,1,0,1,
+                                              0,0,1,0),3,4, byrow=TRUE)
+				}
+				else{
+					mat <- matrix(c(0,1,0,0,
+                                              0,1,0,1),2,4, byrow=TRUE)
+				}
+			}
+			else{
+				mat <- matrix(c(0,1,0,
+                                        0,1,1),2,3, byrow=TRUE)
+			}
+		}
+	}
+	if(nval==3){
+		if(all(val==c(-1,1,2))) constr.slopes <- 1 
+		if(constr.slopes==2){
+			if(all(val==c(-1,0,1))|all(val==c(-1,0,2))){
+				if(!continuous){
+					if(constr.intercepts){
+						mat <- matrix(c(0,1,0,1,0,0,
+	                                              0,0,0,0,0,1,
+	                                              0,0,0,-1,0,0,
+	                                              0,0,1,0,0,0,
+	                                              0,0,0,0,1,0),5,6, byrow=TRUE)
+					}
+					else{
+						mat <- matrix(c(0,1,0,1,0,0,
+	                                              0,0,0,0,0,1,
+	                                              0,0,0,-1,0,0),3,6, byrow=TRUE)
+					}
+				}
+				else{
+						mat <- matrix(c(0,1,1,0,
+	                                              0,0,0,1,
+	                                              0,0,-1,0),3,4, byrow=TRUE)
+				}
+			}
+			if(all(val==c(0,1,2))){
+				if(!continuous){
+					if(constr.intercepts){
+						mat <- matrix(c(0,1,0,0,0,0,
+	                                              0,0,0,1,0,0,
+	                                              0,0,0,1,0,1,
+	                                              0,0,1,0,0,0,
+	                                              0,0,0,0,1,0),5,6, byrow=TRUE)
+					}
+					else{
+						mat <- matrix(c(0,1,0,0,0,0,
+	                                              0,0,0,1,0,0,
+	                                              0,0,0,1,0,1),3,6, byrow=TRUE)
+					}
+				}
+				else{
+						mat <- matrix(c(0,1,0,0,
+	                                              0,0,1,0,
+	                                              0,0,1,1),3,4, byrow=TRUE)
+				}
+			}
+		}
+		if(constr.slopes==1){
+			if(!continuous){
+				if(constr.intercepts){
+					mat <- matrix(c(0,1,0,0,0,0,
+                                              0,1,0,1,0,0,
+                                              0,1,0,1,0,1,
+                                              0,0,1,0,0,0,
+                                              0,0,0,0,1,0),5,6, byrow=TRUE)
+				}
+				else{
+					mat <- matrix(c(0,1,0,0,0,0,
+                                              0,1,0,1,0,0,
+                                              0,1,0,1,0,1),3,6, byrow=TRUE)
+				}
+			}
+			else{
+					mat <- matrix(c(0,1,0,0,
+                                              0,1,1,0,
+                                              0,1,1,1),3,4, byrow=TRUE)
+			}
+		}
+	}
+	if(nval==4){
+		if(constr.slopes==2){
+			if(!continuous){
+				if(constr.intercepts){
+					mat <- matrix(c(0,1,0,1,0,0,0,0,
+                                              0,0,0,0,0,1,0,0,
+                                              0,0,0,0,0,1,0,1,
+                                              0,0,0,-1,0,0,0,0,
+                                              0,0,1,0,0,0,0,0,
+                                              0,0,0,0,1,0,0,0,
+                                              0,0,0,0,0,0,1,0),7,8, byrow=TRUE)
+				}
+				else{
+					mat <- matrix(c(0,1,0,1,0,0,0,0,
+                                              0,0,0,0,0,1,0,0,
+                                              0,0,0,0,0,1,0,1,
+                                              0,0,0,-1,0,0,0,0),4,8, byrow=TRUE)
+				}
+			}
+			else{
+					mat <- matrix(c(0,1,1,0,0,
+                                              0,0,0,1,0,
+                                              0,0,0,1,1,
+                                              0,0,-1,0,0),4,5, byrow=TRUE)
+			}
+		}
+		if(constr.slopes==1){
+			if(!continuous){
+				if(constr.intercepts){
+					mat <- matrix(c(0,1,0,0,0,0,0,0,
+                                              0,1,0,1,0,0,0,0,
+                                              0,1,0,1,0,1,0,0,
+                                              0,1,0,1,0,1,0,1,
+                                              0,0,1,0,0,0,0,0,
+                                              0,0,0,0,1,0,0,0,
+                                              0,0,0,0,0,0,1,0),7,8, byrow=TRUE)
+				}
+				else{
+					mat <- matrix(c(0,1,0,0,0,0,0,0,
+                                              0,1,0,1,0,0,0,0,
+                                              0,1,0,1,0,1,0,0,
+                                              0,1,0,1,0,1,0,1),4,8, byrow=TRUE)
+				}
+			}
+			else{
+					mat <- matrix(c(0,1,0,0,0,
+                                              0,1,1,0,0,
+                                              0,1,1,1,0,
+                                              0,1,1,1,1),4,5, byrow=TRUE)
+			}
+		}
+	}
+	return(mat)
+}
\ No newline at end of file
diff --git a/R/modify.conf.r b/R/modify.conf.r
new file mode 100644
index 0000000..47ada01
--- /dev/null
+++ b/R/modify.conf.r
@@ -0,0 +1,76 @@
+# Rearrange calls to avoid low number of measurement for each state
+
+# Author: Gwenael G.R. Leday
+
+################## Details
+## discard=FALSE:
+# 1. Reference label is normals (coded 0)
+# 2. If not enough observations:
+# - 'losses' merge to 'normals'
+# - 'gains' merge to 'normals'
+# - 'amplifications' merge to 'gains' or 'normals' (if no 'gains'
+#   or not enough observations as well)
+#
+## discard=TRUE:
+# If number of observations is lower than min.obs,
+# these observations are discarded (replaced by NAs)
+
+modify.conf <- function(cghcall, min.obs=3, discard=TRUE){
+
+	val <- sort(unique(cghcall))
+	nval <- length(val)
+	distr <- as.matrix(table(cghcall))
+	val.ref <- which(val==0)
+	out <- cghcall
+
+	if(nval!=1){
+		if(discard){
+			ind <- distr<min.obs
+			if(any(ind)){
+				xx <- val[which(ind)]
+				for(i in 1:length(xx)) out[out==xx[i]] <- NA
+			}
+		}
+		else{
+			if(-1%in%val){
+				if(distr[val==-1]<min.obs){
+					out[out==-1] <- 0
+					distr <- as.matrix(table(out))
+					val <- sort(unique(out))
+				}
+			}
+			if(2%in%val){
+				if(distr[val==2]<min.obs){
+					out[out==2] <- 1
+					distr <- as.matrix(table(out))
+					val <- sort(unique(out))
+				}
+			}
+			if(1%in%val){
+				if(distr[val==1]<min.obs){
+					if(2%in%out) out[out==2] <- 1
+					else{
+						out[out==1] <- 0
+					}	
+					distr <- as.matrix(table(out))
+					val <- sort(unique(out))
+				}
+			}
+			if(0%in%val){
+				if(distr[val==0]<min.obs){
+					if(1%in%out) out[out==0] <- 1
+					else{
+						if(-1%in%out) out[out==0] <- -1
+					}
+					distr <- as.matrix(table(out))
+					val <- sort(unique(out))
+				}
+			}
+			if(length(val)==1){
+				out <- rep(0,length(cghcall))
+				names(out) <- names(cghcall)
+			}
+		}
+	}
+	out
+}
diff --git a/R/plot.plrs.r b/R/plot.plrs.r
new file mode 100644
index 0000000..43f8ae8
--- /dev/null
+++ b/R/plot.plrs.r
@@ -0,0 +1,143 @@
+# Plot data points, model fit and confidence bands
+
+# Author: Gwenael G.R. Leday
+
+plot.plrs <- function(x,
+                      col.line = "black",
+                      col.pts = c("red", "blue", "green2", "green4"),
+                      col.cb = "yellow",
+                      xlim = c(floor(min(x@data$cghseg)),ceiling(max(x@data$cghseg))),
+                      ylim = c(floor(min(x@data$expr)),ceiling(max(x@data$expr))),
+                      pch = 16, lwd=4, cex = 1.2,
+                      xlab="", ylab="", main = "",                      
+                      add = FALSE,
+                      lty = 1,
+                      lin = FALSE,...){
+
+	# Error handling
+	if(class(x)!="plrs") stop("An object of class \"plrs\" is required")
+
+	# Get data and set colors per state
+	seg <- x@data$cghseg
+	rna <- x@data$expr
+	called <- x@data$cghcall
+	val  <- x@data$conf
+	nval <- length(val)
+	pl <- FALSE
+	if(x@type=="Linear" || x@type=="Intercept") pl <- TRUE
+	if(!is.vector(col.pts)) stop("col.pts is not correct")
+	if(length(col.pts)!=1 & length(col.pts)!=4){
+		warning("col.pts length wrong: set to default")
+		col.pts <- c("red", "blue", "green2", "green4")
+	}
+	colr <- rep(col.pts[1],length(called))
+	if(length(col.pts)>1){
+		for(i in val) colr[called==i] <- col.pts[c(-1,0,1,2)==i]
+	}
+	if(length(pch)>1){
+		if(length(pch)!=4) stop("Either one or four values as input")
+		pchtp <- rep(pch[1],length(called))
+		for(i in val) pchtp[called==i] <- pch[c(-1,0,1,2)==i]
+	}
+	else{
+		pchtp <- rep(pch,length(called))
+	}
+	if(length(x@cb)!=0){
+		xCB <- x@cb$x
+		ylbCB <- x@cb$inf
+		yubCB <- x@cb$sup
+	}
+
+	# Main plot
+	if(!add){
+		plot(seg, rna,
+                 main = main,
+                 xlab = xlab, ylab = ylab,
+                 ylim = ylim, xlim = xlim,
+                 pch = pchtp, cex = cex, lwd=lwd,
+                 col = colr)
+	}
+
+	# Add model
+	if(!pl){
+		xx <- NULL
+		if(!is.null(knots(x))) xx <- c(knots(x)+1e-5,knots(x)-1e-5)
+		xx <- c(xx, seq(from=xlim[1]-0.5, to=xlim[2]+0.5, length=50-length(knots(x))))
+		if(!is.null(knots(x))) xx[xx%in%knots(x)] <- xx[xx%in%knots(x)] + 1e-4
+		xx <- sort(xx)
+		yy <- predict(x, newcghseg=xx)
+		nknots <- length(knots(x))
+		xLine <- xx
+		yLine <- yy
+		xPts <- seg
+		yPts <- rna
+		ct <- 0
+		bool <- TRUE
+		while(bool){
+			ct <- ct + 1
+			if(ct<=nknots){
+				# cb
+				if(length(x@cb)!=0){
+					indCB <- (xCB<knots(x)[ct]) & (!is.na(xCB))
+					polygon(c(xCB[indCB],xCB[indCB][length(xCB[indCB]):1]), c(ylbCB[indCB], yubCB[indCB][length(yubCB[indCB]):1]), col = col.cb, border = NA)
+					xCB[indCB] <- NA
+					ylbCB[indCB] <- NA
+					yubCB[indCB] <- NA
+				}
+
+				# Data points
+				indPts <- (xPts<knots(x)[ct]) & (!is.na(xPts))
+				points(xPts[indPts], yPts[indPts], pch = pchtp[indPts], cex = cex, lwd=lwd, col = colr[indPts])
+				xPts[indPts] <- NA
+				yPts[indPts] <- NA
+
+				# Line
+				indLine <- (xLine<knots(x)[ct]) & (!is.na(xLine))
+				lines(xLine[indLine], yLine[indLine], col=col.line, lwd=lwd, lty=lty)
+				xLine[indLine] <- NA
+				yLine[indLine] <- NA
+			}
+			else{
+				# cb
+				if(length(x@cb)!=0){
+					indCB <- !is.na(xCB)
+					polygon(c(xCB[indCB],xCB[indCB][length(xCB[indCB]):1]), c(ylbCB[indCB], yubCB[indCB][length(yubCB[indCB]):1]), col = col.cb, border = NA)
+					xCB[indCB] <- NA
+					ylbCB[indCB] <- NA
+					yubCB[indCB] <- NA
+				}
+
+				# Data points
+				indPts <- !is.na(xPts)
+				points(xPts[indPts], yPts[indPts], pch = pchtp[indPts], cex = cex, lwd=lwd, col = colr[indPts])
+				xPts[indPts] <- NA
+				yPts[indPts] <- NA
+
+				# Line
+				indLine <- !is.na(xLine)
+				lines(xLine[indLine], yLine[indLine], col=col.line, lwd=lwd, lty=lty)
+				xLine[indLine] <- NA
+				yLine[indLine] <- NA
+			}
+			if(all(is.na(xLine))){
+				bool <- FALSE
+				box()
+			}
+		}
+		if(lin){
+			# Plot simple line
+			obj <- plrs(expr=rna, cghseg=seg, cghcall=NULL,
+			            constr = x@call.arg$constr)
+			plot(obj, add=TRUE, col.line=col.line, lty=2, lwd=lwd)
+		}
+	}
+	else{
+		if(length(x@cb)!=0){
+			polygon(c(xCB,xCB[length(xCB):1]), c(ylbCB, yubCB[length(yubCB):1]), col = col.cb, border = NA)
+			points(seg, rna, pch = pchtp, cex = cex, lwd=lwd, col = colr)
+			box()
+		}
+		if(length(coef(x))==2) abline(a=coef(x)[1] , b=coef(x)[2], lwd=lwd, col=col.line, lty=lty)
+		if(length(coef(x))==1) abline(a=coef(x)[1] , b=0, lwd=lwd, col=col.line, lty=lty)
+	}
+}
diff --git a/R/plrs-methods.r b/R/plrs-methods.r
new file mode 100644
index 0000000..c10b3b4
--- /dev/null
+++ b/R/plrs-methods.r
@@ -0,0 +1,128 @@
+######################
+## Methods for "plrs"
+######################
+
+setMethod(
+	f = "coef",
+	signature = "plrs",
+	definition = function(object,...) object@coefficients
+)
+setMethod(
+	f = "fitted",
+	signature = "plrs",
+	definition = function(object,...) object@fitted.values
+)
+setMethod(
+	f = "residuals",
+	signature = "plrs",
+	definition = function(object,...) object@residuals
+)
+setMethod(
+	f = "model.matrix",
+	signature = "plrs",
+	definition = function(object,...) object@X
+)
+setMethod(
+	f = "plot",
+	signature = "plrs",
+	definition = function(x,y,...) plot.plrs(x,...)
+)
+setMethod(
+	f = "knots",
+	signature = "plrs",
+	definition = function(Fn,...) Fn@mdata$mknots
+)
+setMethod(
+	f = "print",
+	signature = "plrs",
+	definition = function(x,...){
+		# Display
+		cat("\n")
+		cat("Object of class \"plrs\"\n\n")
+		if(x@selected){
+			cat("Selected spline coefficients:\n")
+		}
+		else{
+			cat("Spline coefficients:\n")
+		}
+		print(round(coef(x),5))
+		cat("\n")
+
+		if(x@call.arg$constr){
+			cat("Model is constrained:\n")
+			cat("constr.slopes =", x@call.arg$constr.slopes,"\n")
+			if(!is.null(x@call.arg$constr.slopes.to.zero)) cat("constr.slopes.to.zero =", x@call.arg$constr.slopes.to.zero,"\n")
+			cat("constr.intercepts =", x@call.arg$constr.intercepts,"\n\n")
+		}
+		if(length(x@test)!=0){
+			cat("Testing:\n")
+			cat("stat =", x@test$stat, "\n")
+			cat("quantile = ", x@test$qbetabar,"  (alpha = ", x@test$alpha, ")","\n", sep="")
+			cat("p-value = ", format.pval(x@test$pvalue),"\n\n")
+		}
+	}
+)
+setMethod(
+	f = "show",
+	signature = "plrs",
+	definition = function(object) print(object)
+)
+setMethod(
+	f = "summary",
+	signature = "plrs",
+	definition = function(object,...){
+		print(object)
+		cat("Configuration:\n")
+		print(effects(object))
+	}
+)
+setMethod(
+	f = "predict",
+	signature = "plrs",
+	definition = function(object, ...) predict.plrs(object, ...)
+)
+
+setMethod(
+	f = "effects",
+	signature = "plrs",
+	definition = function(object,...){
+		mat <- table(object@mdata$mcghcall)
+		confMat <- matrix(NA,3,4)
+		colnames(confMat) <- c("Loss","Normal","Gain","Ampl.")
+		rownames(confMat) <- c("n.obs","I.effect", "S.effect")
+		labs <- as.numeric(dimnames(mat)[[1]])
+		vals <- as.numeric(mat)
+		confMat[1,c(-1,0,1,2)%in%labs] <- vals[labs%in%c(-1,0,1,2)]
+		if(length(labs)!=1){
+			if(-1%in%labs){
+				confMat[2,1] <- round(coef(object)["theta1.0"],3)
+				confMat[3,1] <- round(coef(object)["theta0.1"],3)
+			}
+			if(1%in%labs){
+				if(-1%in%labs){
+					confMat[2,3] <- round(coef(object)["theta2.0"],3)
+					confMat[3,3] <- round(coef(object)["theta2.1"],3)
+				}
+				else{	
+					confMat[2,3] <- round(coef(object)["theta1.0"],3)
+					confMat[3,3] <- round(coef(object)["theta1.1"],3)
+				}
+			}
+			if(2%in%labs){
+				if(-1%in%labs){
+					confMat[2,4] <- round(coef(object)["theta3.0"],3)
+					confMat[3,4] <- round(coef(object)["theta3.1"],3)
+				}
+				else{
+					confMat[2,4] <- round(coef(object)["theta2.0"],3)
+					confMat[3,4] <- round(coef(object)["theta2.1"],3)
+				}
+			}
+		}
+		else{
+			confMat[2,2] <- 0
+			confMat[3,2] <- round(coef(object)["theta0.1"],3)
+		}
+		return(t(confMat))
+	}
+)
diff --git a/R/plrs.cb.r b/R/plrs.cb.r
new file mode 100644
index 0000000..b6e87da
--- /dev/null
+++ b/R/plrs.cb.r
@@ -0,0 +1,37 @@
+# Likelihood-based confidence bands for (constrained) plrs
+
+# Author: Gwenael G.R. Leday
+
+plrs.cb <- function(object, alpha=0.05, newcgh=NULL){
+
+	if(!class(object)%in%c("plrs","plrs.select")) stop("An object of class \"plrs\" or \"plrs.select\" is required")
+	if(length(object@test)==0 & !object@selected) stop("Likelihood ratio test using plrs.test() is required before computing bounds")
+	if(is.null(newcgh)){
+		if(!is.null(knots(object))) newcgh <- sort(c(knots(object)+1e-5,knots(object)-1e-5))
+		xx <- seq(from=floor(min(object@data$cghseg))-0.5, to=ceiling(max(object@data$cghseg))+0.5, length=100-length(knots(object)))
+		newcgh <- sort(c(newcgh, xx))
+	}
+	if(!object@call.arg$constr || !object@call.arg$constr.intercepts) stop("Model has to be fully constrained")
+	if(object@selected){
+		obj <- plrs(expr=object@mdata$mexpr, cghseg=object@mdata$mcghseg, cghcall=object@mdata$mcghcall,
+				probloss = object@mdata$mprobloss,
+				probnorm = object@mdata$mprobnorm,
+				probgain = object@mdata$mprobgain,
+				probamp = object@mdata$mprobamp,
+				knots = NULL,
+				continuous = object@call.arg$continuous,
+				constr = object@call.arg$constr,
+				constr.slopes = object@call.arg$constr.slopes,
+				constr.intercepts = object@call.arg$constr.intercepts)
+		obj <- plrs.test(obj, alpha=alpha)
+	}else{
+		obj <- object
+	}
+
+	allCI <- sapply(newcgh, .respCI, object=obj, alpha=alpha)
+
+	object@cb <- list(inf=allCI[1,], sup=allCI[2,], x=newcgh)
+
+	return(object)
+}
+
diff --git a/R/plrs.fit.r b/R/plrs.fit.r
new file mode 100644
index 0000000..6a75cd3
--- /dev/null
+++ b/R/plrs.fit.r
@@ -0,0 +1,132 @@
+# Internal function
+# Fit a (constrained) piecewise linear regression spline
+
+# Author: Gwenael G.R. Leday
+
+.plrs.fit <- function(X, matconstr, mytype, modelFull, callArg){
+
+	nval <- length(modelFull@mdata$mconf)
+	constr <- ifelse(is.null(matconstr), FALSE, TRUE)
+	expr <- modelFull@mdata$mexpr
+
+	### Fit (constrained) plrs
+	if(nval!=1){
+		# Unconstrained model
+		if(!constr){
+			lm.obj <- lm(expr~.-1, data=as.data.frame(cbind(expr,X)))
+			sol <- coef(lm.obj)
+			cov.mat <- vcov(lm.obj)
+			S2 <- summary(lm.obj)$sigma^2
+			df.error <- lm.obj$df.res
+		}
+		# Constrained model
+		else{
+			# Input solve.QP
+			Dmat <- crossprod(X,X)
+			dvec <- crossprod(X,expr)
+
+			# Input
+			Amat <- t(matconstr)
+			bvec <- rep(0, nrow(matconstr))
+
+			## Quadratic programming
+			sol <- solve.QP(Dmat, as.vector(dvec), Amat, bvec, meq = 0)$sol
+			names(sol) <- colnames(X)
+		}
+	}
+	### Fit simple linear model
+	else{
+		# Unconstrained
+		if(!constr){
+			lm.obj <- lm(expr~.-1, data=as.data.frame(cbind(expr,X)))
+			sol <- coef(lm.obj)
+			cov.mat <- vcov(lm.obj)
+			S2 <- summary(lm.obj)$sigma^2
+			df.error <- lm.obj$df.res
+		}
+		# Constrain slope to be positive or null
+		else{
+			# Input
+			Dmat <- crossprod(X,X)
+			dvec <- crossprod(X,expr)
+			Amat <- t(t(c(0,1)))
+			bvec <- 0
+
+			## Quadratic programming
+			sol <- solve.QP(Dmat, as.vector(dvec), Amat, bvec, meq = 0)$sol
+
+			names(sol) <- colnames(X)
+		}
+		chg.pt <- numeric()
+	}
+
+	### Output plrs object
+	confModel <- function(modeFull, mylabs){
+
+		theconf <- modelFull@mdata$mconf
+		theknots <- modeFull@mdata$mknots
+		idx1 <- substr(mylabs, 6,6)
+		idx2 <- substr(mylabs, 8,8)
+		myidx <- as.numeric(unique(idx1))
+
+		if(sum(idx2=="0")==1){
+			cont <- TRUE
+		}
+		else{
+			cont <- FALSE
+		}
+		if(length(mylabs)==1){
+			myconf <- 0
+			myknots <- numeric()
+		}
+		else{
+			myconf <- theconf[myidx+1]
+			deleted.knots <- !(1:length(theknots))%in%myidx
+			if(all(deleted.knots)){
+				myknots <- numeric()
+			}
+			else{
+				myknots <- theknots[!deleted.knots]
+			}			
+		}
+		return(list(conf=myconf, knots=myknots, conti=cont))
+	}
+	tp <- confModel(modelFull, colnames(X))
+	chg.pt <- tp$knots
+	continuous <- tp$conti
+	mdata <- modelFull@mdata
+	mdata$mconf <- tp$conf
+	mdata$mknots <- tp$knots
+
+	fit.val <- as.vector(crossprod(t(X),sol)) #as.vector(X %*% sol)
+	names(fit.val) <- names(expr)
+	if(constr){
+		qp <- list("Dmat" = Dmat, "dvec" = dvec, "Amat" = Amat,
+					"bvec" = bvec, "meq" = 0)
+	}
+	else{
+		qp <- list("vcov" = cov.mat, "S2" = S2, "df.error" = df.error)
+	}
+
+	call.arg <- list("continuous" = callArg$continuous,
+			     "constr" = callArg$constr,
+			     "constr.intercepts" = modelFull@call.arg$constr.intercepts,
+			     "constr.slopes" = modelFull@call.arg$constr.slopes,
+			     "config" = paste(as.character(mdata$mconf),collapse=","))
+
+	out <- new("plrs",
+			coefficients = sol,
+			fitted.values = fit.val,
+			residuals = expr - fit.val,
+			X = X,
+			data = modelFull@data,
+			mdata = mdata,
+			QP = qp,
+			test = list(),
+			cb = list(),
+			selected = TRUE,
+			type = mytype,
+			call.arg = call.arg)
+
+	return(out)
+}
diff --git a/R/plrs.r b/R/plrs.r
new file mode 100644
index 0000000..d324b2a
--- /dev/null
+++ b/R/plrs.r
@@ -0,0 +1,324 @@
+# Fit a (constrained) piecewise linear regression spline
+
+# Author: Gwenael G.R. Leday
+
+plrs <- function(expr, cghseg, cghcall=NULL,
+                 probloss = NULL,
+                 probnorm = NULL,
+                 probgain = NULL,
+                 probamp = NULL,
+                 knots = NULL,
+                 continuous = FALSE,
+                 constr = TRUE,
+                 constr.slopes = 2,
+                 constr.intercepts = TRUE,
+                 min.obs=3,
+                 discard.obs=TRUE){
+
+	#### Error handling
+	## check input data
+	if(!is.vector(expr)) stop("expr must be a vector")
+	if(!is.vector(cghseg)) stop("cghseg must be a vector")
+	if(any(is.na(expr))) stop("Missing values are not allowed in expr")
+	if(any(is.na(cghseg))) stop("Missing values are not allowed in cghseg")
+	if(is.null(cghcall)) cghcall <- rep(0,length(cghseg))
+	else{
+		if(!is.vector(cghcall)){
+			stop("cghcall must be a vector")
+		}
+		else{
+			if(any(is.na(cghcall))) stop("Missing values are not allowed in cghcall")
+			if(2%in%cghcall & !1%in%cghcall){
+				cghcall[cghcall==2] <- 1
+				warning("No gains: observations labelled amplifications are changed to gains")
+			}
+		}
+	}
+	if(length(expr)!=length(cghseg))	stop("expr and cghseg are incompatible")
+	if(length(expr)!=length(cghcall)) 	stop("expr and cghcall are incompatible")
+	if(length(cghseg)!=length(cghcall)) stop("cghseg and cghcall are incompatible")
+	if(!is.null(probloss) & !is.vector(probloss)) 	stop("probloss must be a vector")
+	if(!is.null(probnorm) & !is.vector(probnorm)) 	stop("probnorm must be a vector")
+	if(!is.null(probgain) & !is.vector(probgain)) 	stop("probgain must be a vector")
+	if(!is.null(probamp) & !is.vector(probamp)) 	stop("probamp must be a vector")
+
+	val <- sort(unique(cghcall))
+	conf <- val
+	nval <- length(val)
+	if(nval>4) stop("Configurations with more than 4 states are not allowed")
+	out <- NULL
+
+	## Check configuration
+	# min.obs
+	if(!is.numeric(min.obs)) stop("'min.obs' must be a numeric")
+	min.obs <- round(min.obs)
+	if(min.obs<1) stop("'min.obs' must be a positive integer")
+	# discard.obs
+	if(!is.logical(discard.obs)) stop("'discard.obs' must be a logical")
+	# Nb of obs
+	mcghcall <- modify.conf(cghcall, min.obs=min.obs, discard=discard.obs)
+
+	mexpr <- expr
+	mcghseg <- cghseg
+	mprobloss <- probloss
+	mprobnorm <- probnorm
+	mprobgain <- probgain
+	mprobamp <- probamp
+	if(any(is.na(mcghcall))){
+		mexpr <- mexpr[!is.na(mcghcall)]
+		mcghseg <- mcghseg[!is.na(mcghcall)]
+		if(!is.null(mprobloss)) mprobloss <- probloss[!is.na(mcghcall)]
+		if(!is.null(mprobnorm)) mprobnorm <- mprobnorm[!is.na(mcghcall)]
+		if(!is.null(mprobgain)) mprobgain <- mprobgain[!is.na(mcghcall)]
+		if(!is.null(mprobamp)) mprobamp <- mprobamp[!is.na(mcghcall)]
+		mcghcall <- mcghcall[!is.na(mcghcall)]
+	}
+	val <- sort(unique(mcghcall))
+	nval <- length(val)
+
+	if(nval!=1){
+		# knots
+		if(is.null(knots)){
+			chg.pt <- .estimateKnots(cghseg = mcghseg, cghcall = mcghcall,
+     	                                     probloss = mprobloss, probnorm = mprobnorm,
+                                           probgain = mprobgain, probamp = mprobamp)
+			if(sum(order(chg.pt)==1:length(chg.pt))!=length(chg.pt)){
+				warning("Cannot fit a plrs model: estimated knots are non-ordered (this may be caused by strong outliers).")
+			}
+			# Recalling
+			cl <- rep(val[1],length(mcghcall))
+			if(length(chg.pt)>=1) cl[mcghseg>=chg.pt[1]] <- val[2]
+			if(length(chg.pt)>=2) cl[mcghseg>=chg.pt[2]] <- val[3]
+			if(length(chg.pt)>=3) cl[mcghseg>=chg.pt[3]] <- val[4]
+			names(cl) <- names(mcghcall)
+			# Reconfigure
+			tp <- modify.conf(cl, min.obs=min.obs, discard=discard.obs)
+			# If not identical, this means that min obs is not fulfilled for a state after recalling, hence: 
+			if(!identical(cl,tp)){
+				mcghcall <- tp
+				if(any(is.na(mcghcall))){
+					mexpr <- mexpr[!is.na(mcghcall)]
+					mcghseg <- mcghseg[!is.na(mcghcall)]
+					if(!is.null(mprobloss)) mprobloss <- probloss[!is.na(mcghcall)]
+					if(!is.null(mprobnorm)) mprobnorm <- mprobnorm[!is.na(mcghcall)]
+					if(!is.null(mprobgain)) mprobgain <- mprobgain[!is.na(mcghcall)]
+					if(!is.null(mprobamp)) mprobamp <- mprobamp[!is.na(mcghcall)]
+					mcghcall <- mcghcall[!is.na(mcghcall)]
+				}
+				val <- sort(unique(mcghcall))
+				nval <- length(val)
+				chg.pt <- .estimateKnots(cghseg = mcghseg, cghcall = mcghcall,
+     	                                     probloss = mprobloss, probnorm = mprobnorm,
+                                           probgain = mprobgain, probamp = mprobamp)
+				if(sum(order(chg.pt)==1:length(chg.pt))!=length(chg.pt)){
+					warning("Cannot fit a plrs model: estimated knots are non-ordered (this may be caused by strong outliers).")
+				}
+				# Recalling
+				cl <- rep(val[1],length(mcghcall))
+				if(length(chg.pt)>=1) cl[mcghseg>=chg.pt[1]] <- val[2]
+				if(length(chg.pt)>=2) cl[mcghseg>=chg.pt[2]] <- val[3]
+				if(length(chg.pt)>=3) cl[mcghseg>=chg.pt[3]] <- val[4]
+				names(cl) <- names(mcghcall)
+			}
+			mcghcall <- cl
+			rm(cl)
+		}
+		else{
+			if(!is.numeric(knots)) stop("'knots' is not correct")
+			if(length(knots)!=(nval-1))  stop("wrong number of knots")
+			chg.pt <- knots
+		}
+
+		# continuous
+		if(is.null(continuous) || !is.logical(continuous) ||
+              (is.vector(continuous) & length(continuous)>1)) stop("'continuous' is not correct")
+
+		# constr
+		if(!(is.logical(constr) & length(constr)==1)) stop("'constr' must be a logical")
+		else{
+			# constr.slopes
+			if(!is.null(constr.slopes)){
+				if(!(is.numeric(constr.slopes) & length(constr.slopes)==1))
+					stop("'constr.slopes' must be a numeric")
+				if(sum(constr.slopes%in%c(1,2))!=length(constr.slopes))
+					stop("wrong type of constraints for slopes")
+			}
+			else{
+				constr.slopes <- 2
+			}
+
+			# constr.intercepts
+			if(!is.null(constr.intercepts)){
+				if(!(is.logical(constr.intercepts) & length(constr.intercepts)==1))
+					stop("'constr.intercepts' must be a logical")
+			}
+			else{
+				constr.intercepts <- FALSE
+			}
+			if(constr & !constr.intercepts & is.null(constr.slopes)) constr <- FALSE
+		}
+	}
+
+
+	### Fit (constrained) plrs
+	if(nval!=1){
+		# Unconstrained model
+		if(!constr){
+			# Design matrix
+			X <- .Bmat(mcghseg, knots = chg.pt, 
+				    continuous = continuous,
+				    general.intercept = TRUE)
+
+			# fit
+			#sol <- lm.fit(x=X, y=em)$coef
+			lm.obj <- lm(mexpr~.-1, data=as.data.frame(cbind(mexpr,X)))
+			sol <- coef(lm.obj)
+			cov.mat <- vcov(lm.obj)
+			S2 <- summary(lm.obj)$sigma^2
+			df.error <- lm.obj$df.res
+		}
+		# Constrained model
+		else{
+			# Design matrix
+			X <- .Bmat(mcghseg, knots=chg.pt, 
+				    continuous=continuous,
+				    general.intercept=TRUE)
+
+			# Input solve.QP
+			Dmat <- crossprod(X,X)
+			dvec <- crossprod(X,mexpr)
+
+			# Get matrix of constraints
+			get.mat <- .make.mat.constr(constr.intercepts=constr.intercepts, constr.slopes=constr.slopes, continuous=continuous, val=val)
+
+			# Input
+			Amat <- t(get.mat)
+			bvec <- rep(0, nrow(get.mat))
+
+			## Quadratic programming
+			sol <- solve.QP(Dmat, as.vector(dvec), Amat, bvec, meq = 0)$sol
+			names(sol) <- colnames(X)
+		}
+	}
+	### Fit simple linear model
+	else{
+		# Design matrix
+		X <- cbind(1,mcghseg)
+		colnames(X) <- c("theta0.0","theta0.1")
+
+		# Unconstrained
+		if(!constr){
+			#sol <- lm.fit(x=X, y=mexpr)$coef
+			lm.obj <- lm(mexpr~.-1, data=as.data.frame(cbind(mexpr,X)))
+			sol <- coef(lm.obj)
+			cov.mat <- vcov(lm.obj)
+			S2 <- summary(lm.obj)$sigma^2
+			df.error <- lm.obj$df.res
+		}
+		# Constrain slope to be positive or null
+		else{
+			# Input
+			Dmat <- crossprod(X,X)
+			dvec <- crossprod(X,mexpr)
+			Amat <- t(t(c(0,1)))
+			bvec <- 0
+			meq <- 0
+			
+			## Quadratic programming
+			sol <- solve.QP(Dmat, as.vector(dvec), Amat, bvec, meq = meq)$sol
+
+			names(sol) <- colnames(X)
+		}
+		chg.pt <- numeric()
+	}
+
+
+	### Output plrs object
+	fit.val <- as.vector(crossprod(t(X),sol))
+	names(fit.val) <- names(mexpr)
+	dat <- list("expr" = expr, "cghseg" = cghseg, "cghcall" = cghcall,
+                  "probloss" = probloss, "probnorm" = probnorm, "probgain" = "probgain", "probamp" = probamp,
+                  "conf" = conf, "knots" = chg.pt)
+
+	mdat <- list("mexpr" = mexpr, "mcghseg" = mcghseg, "mcghcall" = mcghcall,
+                   "mprobloss" = mprobloss, "mprobnorm" = mprobnorm, "mprobgain" = mprobgain, "mprobamp" = mprobamp,
+                   "mconf" = val, "mknots" = chg.pt)
+
+	if(constr){
+		qp <- list("Dmat" = Dmat, "dvec" = dvec, "Amat" = Amat,
+					"bvec" = bvec, "meq" = 0)
+	}
+	else{
+		qp <- list("vcov" = cov.mat, "S2" = S2, "df.error" = df.error)
+	}
+
+	## Type of model: Functions needed
+	# Identify type of model based on estimated parameters.
+	is.intercept <- function(x){
+		if(("theta0.0" %in% names(x)) & length(x)==1){
+			return(TRUE)
+		}
+		else{
+			return(FALSE)
+		}
+	}
+	is.linear <- function(x){
+		if(all(c("theta0.0","theta0.1") %in% names(x)) & length(x)==2){
+			return(TRUE)
+		}
+		else{
+			return(FALSE)
+		}
+	}
+	is.piecewise <- function(x){
+		if(!is.intercept(x) & !is.linear(x)){
+			return(TRUE)
+		}
+		else{
+			return(FALSE)
+		}
+	}
+	is.piecewise.level <- function(x){
+		IDintercpts <- as.numeric(substr(names(x),8,9))
+		if(is.piecewise(x) & sum(IDintercpts)==0){
+			return(TRUE)
+		}
+		else{
+			return(FALSE)
+		}
+	}
+	# Def: is piecewise linear if it is piecewise and not level
+	is.piecewise.linear <- function(x){
+		if(is.piecewise(x) & !is.piecewise.level(x)){
+			return(TRUE)
+		}
+		else{
+			return(FALSE)
+		}
+	}
+	if(is.intercept(sol)) model.type <- "Intercept"
+	if(is.linear(sol)) model.type <- "Linear"
+	if(is.piecewise.level(sol)) model.type <- "Piecewise level"
+	if(is.piecewise.linear(sol)) model.type <- "Piecewise linear"
+
+	call.arg <- list("continuous" = continuous,
+			     "constr" = constr,
+			     "constr.intercepts" = constr.intercepts,
+			     "constr.slopes" = constr.slopes,
+			     "config" = paste(as.character(val),collapse=","))
+	out <- new("plrs",
+			coefficients = sol,
+			fitted.values = fit.val,
+			residuals = mexpr - fit.val,
+			X = X,
+			data = dat,
+			mdata = mdat,
+			QP = qp,
+			test = list(),
+			cb = list(),
+			selected = FALSE,
+			type = model.type,
+			call.arg = call.arg)
+
+	return(out)
+}
diff --git a/R/plrs.select-methods.r b/R/plrs.select-methods.r
new file mode 100644
index 0000000..8e1efd9
--- /dev/null
+++ b/R/plrs.select-methods.r
@@ -0,0 +1,51 @@
+############################
+## Methods for "plrs.select"
+############################
+
+setMethod(
+	f = "plot",
+	signature = "plrs.select",
+	definition = function(x,y,...) plot.plrs(x@model,...)
+)
+setMethod(
+	f = "print",
+	signature = "plrs.select",
+	definition = function(x,...){
+		# Display
+		cat("\n")
+		cat("Object of class \"plrs.select\"\n\n")
+		cat(toupper(x@crit), "model selection procedure", "\n\n")
+
+		cat("Coefficients of selected spline:\n")
+		print(round(coef(x@model),5))
+		cat("\n")
+		if(x@model@call.arg$constr){
+			cat("Model is constrained:\n")
+			cat("constr.slopes =", x@model@call.arg$constr.slopes,"\n")
+			if(!is.null(x@model@call.arg$constr.slopes.to.zero)) cat("constr.slopes.to.zero =", x@model@call.arg$constr.slopes.to.zero,"\n")
+			cat("constr.intercepts =", x@model@call.arg$constr.intercepts,"\n\n")
+		}
+	}
+)
+setMethod(
+	f = "show",
+	signature = "plrs.select",
+	definition = function(object) print(object)
+)
+setMethod(
+	f = "summary",
+	signature = "plrs.select",
+	definition = function(object,...){
+		# Display
+		print(object)
+		cat("Model selection table:\n")
+		if(nrow(object@table)>5){
+			print(head(object@table))
+			cat("...\n\n")
+		}
+		else{
+			print(object@table)
+			cat("\n")
+		}
+	}
+)
\ No newline at end of file
diff --git a/R/plrs.select.r b/R/plrs.select.r
new file mode 100644
index 0000000..58552ba
--- /dev/null
+++ b/R/plrs.select.r
@@ -0,0 +1,49 @@
+# Select the best model based on information criteria
+# AIC, AICC, BIC or OSAIC
+
+# Author: Gwenael G.R. Leday
+
+plrs.select <- function(object, crit = ifelse(object@call.arg$constr,"osaic","aic")){
+
+	criter <- crit
+	if(is.null(criter)) stop("No criterion")
+	if(length(criter)>1) stop("Multiple criteria are not allowed")
+	if(!is.character(criter)) stop("Criterion has to be a character")
+	if(!criter%in%c("aic","aicc","bic","osaic")) stop("Wrong criterion")
+
+	# Generate all design and constraint matrices
+	allsubmodels <- .plrs.submodels(object)
+
+	# Fit all submodels
+	allModels <- rep(list(NULL),length(allsubmodels$allMats))
+	for(k in 1:length(allsubmodels$allMats)){
+		allModels[[k]] <- .plrs.fit(X = allsubmodels$allMats[[k]],
+                                        matconstr = allsubmodels$allConstr[[k]],
+                                        mytype = allsubmodels$allTypes[k],
+                                        modelFull = object, callArg=object@call.arg)
+	}
+	allModels <- c(allModels, object)
+	allTypes <- c(allsubmodels$allTypes, "Piecewise linear")
+
+	# Compute criterion
+	f <- function(x) criteria(x, criter)[criter]
+	crits <- unlist(lapply(allModels, f))
+	crits <- matrix(crits, length(crits), 1)
+	rownames(crits) <- paste("model",1:length(allModels), sep="")
+	colnames(crits) <- criter
+
+	# Output
+	best <- which.min(crits)
+	bmodel <- allModels[[best]]
+	bmodel@selected <- TRUE
+	bmodel@type <- allTypes[best]
+
+	out <- new("plrs.select",
+                 "table" = crits,
+                 "model" = bmodel,
+                 "crit" = criter)
+	return(out)
+}
+
+
+
diff --git a/R/plrs.series-methods.r b/R/plrs.series-methods.r
new file mode 100644
index 0000000..a0545df
--- /dev/null
+++ b/R/plrs.series-methods.r
@@ -0,0 +1,62 @@
+################################
+## Methods for "plrs.series"
+################################
+
+setMethod(
+	f = "print",
+	signature = "plrs.series",
+	definition = function(x,...){
+		cat("\n")
+		cat("Object of class \"plrs.series\"\n\n")
+		cat(nrow(x@coefficients),"genes\n")
+		if(is.null(x@call.arg$control.select)){
+			cat("\n")
+			cat("Type of fitted model:\n")
+			cat("continuous =", x@call.arg$control.model$continuous,"\n")
+			if(x@call.arg$control.model$constr){
+				cat("constr.slopes =", x@call.arg$control.model$constr.slopes,"\n")
+				cat("constr.intercepts =", x@call.arg$control.model$constr.intercepts,"\n\n")
+			}
+			else{
+				cat("constr =", x@call.arg$control.model$constr,"\n\n")
+			}
+		}else{
+			cat("Models selected with", toupper(x@call.arg$control.select$crit), "\n\n")
+		}
+	}
+)
+setMethod(
+	f = "show",
+	signature = "plrs.series",
+	definition = function(object){
+		print(object)
+	}
+)
+setMethod(
+	f = "summary",
+	signature = "plrs.series",
+	definition = function(object,...){
+		# Display
+		print(object)
+		cat("Configuration:\n")
+		print(object@general)
+		cat("\n")
+		cat("Nb of models regarding their types:\n")
+		cat("Intercept        ",object@modelsType$summary[1],"\n")
+		cat("Simple linear    ",object@modelsType$summary[2],"\n")
+		if(nrow(object@general)>1){
+			cat("Piecewise level  ",object@modelsType$summary[3],"\n")
+			cat("Piecewise linear ",object@modelsType$summary[4],"\n\n")
+		}
+		else{
+			cat("\n")
+		}
+		if(length(object@test)!=0){
+			cat("Testing:\n")
+			cat("Number of rejected null hypothesis:", sum(object@test[,3]<.1, na.rm=T), "genes\n")
+			cat("(at 0.1 significance level based on \n Benjamini-Hochberg corrected p-values)\n\n")
+		}
+	}
+)
+
+
diff --git a/R/plrs.series.r b/R/plrs.series.r
new file mode 100644
index 0000000..b91c40a
--- /dev/null
+++ b/R/plrs.series.r
@@ -0,0 +1,228 @@
+# Fit plrs models for a series of arrays
+
+# Author: Gwenael G.R. Leday
+
+plrs.series <- function(expr, cghseg, cghcall=NULL,
+				probloss = NULL,
+				probnorm = NULL,
+				probgain = NULL,
+				probamp = NULL,
+				control.model = list(continuous=FALSE, constr=TRUE, constr.slopes=2, constr.intercepts=TRUE, min.obs=3, discard.obs=TRUE),
+				control.select = list(crit = ifelse(control.model$constr,"osaic","aic")),
+				control.test = list(testing=TRUE, cb=FALSE, alpha=0.05),
+				control.output = list(save.models = FALSE, save.plots = FALSE, plot.lin = FALSE, type = "jpeg")){
+
+	#### Check input
+	if(class(expr)=="ExpressionSet"){
+		expr <- exprs(expr)
+	}
+	if(class(cghseg)=="cghSeg"){
+		cghseg <- segmented(cghseg)
+		cghcall <- NULL
+	}
+	if(class(cghseg)=="cghCall"){
+		cghcall <- CGHbase:::calls(cghseg)
+		if(sum(dim(cghcall)==c(0,0))!=0){
+			cghcall <- NULL
+			probloss <- probnorm <- probgain <- probamp <- NULL
+			cghseg <- segmented(cghseg)
+		}
+		else{
+			probloss <- probloss(cghseg)
+			probnorm <- probnorm(cghseg)
+			probgain <- probgain(cghseg)
+			probamp <- probamp(cghseg)
+			cghseg <- segmented(cghseg)
+		}
+	}
+	if(!(is.matrix(expr) & is.matrix(cghseg))){
+		stop("input data have to be in matrix format")
+	}
+	else{
+		if(any(is.na(expr))) stop("Missing values are not allowed in expr")
+		if(!sum(dim(expr)==dim(cghseg))==2)
+			stop("dimensions of input data are not the same")
+	}
+	if(is.null(cghcall)) cghcall <- matrix(0, nrow(expr), ncol(expr), dimnames=dimnames(expr))
+	if(!is.matrix(cghcall)){
+		stop("input data have to be in matrix format")
+	}
+	else{
+		if(!sum(dim(expr)==dim(cghcall))==2)
+			stop("dimensions of input data are not the same")
+	}
+	# control.model
+	if(!is.list(control.model)) stop("control.model is not a list")
+	if(!"continuous"%in%names(control.model)) control.model$continuous <- FALSE
+	if(!"constr"%in%names(control.model)) control.model$constr <- TRUE
+	if(!"constr.slopes"%in%names(control.model)) control.model$constr.slopes <- 2
+	if(!"constr.intercepts"%in%names(control.model)) control.model$constr.intercepts <- TRUE
+	if(!"min.obs"%in%names(control.model)) control.model$min.obs <- 3
+	if(!"discard.obs"%in%names(control.model)) control.model$discard.obs <- TRUE
+	# control.select
+	if(!is.list(control.select) & !is.null(control.select)) stop("control.select is not a list")
+	if(length(control.select)>0){
+		if(!"crit"%in%names(control.select)) control.select$crit <- ifelse(control.model$constr,"osaic","aic")
+	}
+	# control.test
+	if(!is.list(control.test)) stop("control.test is not a list")
+	if(!"testing"%in%names(control.test)) control.test$testing <- TRUE
+	if(!"cb"%in%names(control.test)) control.test$cb <- FALSE
+	if(!"alpha"%in%names(control.test)) control.test$alpha <- FALSE
+	if(!control.test$testing & control.test$cb) control.test$testing <- TRUE
+	# control.output
+	if(!is.list(control.output)) stop("control.output is not a list")
+	if(!"save.models"%in%names(control.output)) control.output$save.models <- FALSE
+	if(!"save.plots"%in%names(control.output)) control.output$save.plots <- FALSE
+	if(!"plot.lin"%in%names(control.output)) control.output$plot.lin <- FALSE
+	if(!"type"%in%names(control.output)) control.output$type <- "jpeg"
+
+
+	# Matrix of coefficients
+	ngenes <- nrow(expr)
+	est <- matrix(NA, ngenes, 8)
+	nam1 <- colnames(est) <- c("theta0.0", "theta0.1", "theta1.0", "theta1.1", "theta2.0", "theta2.1", "theta3.0", "theta3.1")
+	rownames(est) <- paste("gene",1:ngenes,sep="")
+
+	# Matrix of model types
+	modelTypes <- matrix(F, ngenes, 4)
+	colnames(modelTypes) <- c("Intercept", "Linear", "Piecewise level", "Piecewise linear")
+	rownames(modelTypes) <- rownames(est)
+
+	# Matrix of effects
+	Ieff <- Seff <- matrix(F, ngenes, 4)
+	nam3 <- colnames(Ieff) <- colnames(Seff) <- c("Loss", "Norm.", "Gain", "Ampl.")
+	rownames(modelTypes) <- rownames(est)
+
+	# Testing
+	if(control.test$testing){
+		res.test <- matrix(NA, ngenes, 3, dimnames=list(rownames(est), c("stat","raw.pval","BH.adj.pval")))
+	}
+
+	# Summary matrix on aberrations
+	config <- matrix(NA, ngenes, 4)
+	colnames(config) <- as.character(-1:2)
+	rownames(config) <- rownames(est)
+
+	# Create directories for output
+	if(control.output$save.plots) dir.create("plrsSeriesPlots")
+	if(control.output$save.models) dir.create("plrsSeriesObjects")
+	if(is.null(rownames(expr))){
+		lab <- 1:nrow(expr)
+	}
+	else{
+		lab <- rownames(expr)
+	}
+
+	# Relabel rownames duplicates
+	ind <- duplicated(lab)
+	cp <- 1
+	while(sum(ind)!=0){
+		lab[ind] <- paste(lab[ind], "_", cp, sep="")
+		ind <- duplicated(lab)
+		if(sum(ind)!=0) lab[ind] <- substr(lab[ind], 1, nchar(lab[ind])-2)
+		cp <- cp + 1
+	}
+
+	if(ngenes>100){
+		cat("In progress... \n\n")
+		toprint <- round(quantile(1:ngenes, probs=seq(0, 1, 0.10))[-1])
+		t <- proc.time()
+		ct <- 1
+	}
+
+	for(x in 1:ngenes){
+		if(ngenes>100){
+			if(x==toprint[ct]){
+				xx <- (proc.time() - t)[3]
+				cat(names(toprint[ct]), " done (", toprint[ct]," genes),\t  ",
+					"time elapsed = ", .convertToTime(xx),"\n", sep="")
+				ct <- ct + 1
+			}
+		}
+#cat("i = ", x, "\n")
+		error <- try(
+		model <- plrs(expr=expr[x,], cghseg=cghseg[x,], cghcall=cghcall[x,],
+                          probloss = probloss[x,],
+                          probnorm = probnorm[x,],
+                          probgain = probgain[x,],
+                          probamp = probamp[x,],
+                          continuous = control.model$continuous,
+                          constr = control.model$constr,
+                          constr.slopes = control.model$constr.slopes,
+                          constr.intercepts = control.model$constr.intercepts,
+                          min.obs = control.model$min.obs,
+                          discard.obs = control.model$discard.obs)
+		,T)
+		if(!inherits(error, "try-error")){
+			vec <- as.vector(table(model@data$cghcall))
+			config[x,colnames(config)%in%as.character(model@data$conf)] <- vec
+			if(!is.null(control.select)) model <- plrs.select(model, crit = control.select$crit)@model
+			est[x,nam1%in%names(coef(model))] <- coef(model)
+			modelTypes[x,colnames(modelTypes)%in%model@type] <- TRUE
+			tabb <- effects(model)
+			indeff <- !is.na(tabb[,2])
+			Ieff[x,indeff] <- tabb[indeff,2]
+			Seff[x,indeff] <- tabb[indeff,3]
+
+			if(control.test$testing || control.test$cb){
+				model <- plrs.test(model)
+				if(length(model@test)!=0){
+					res.test[x,1] <- model@test$stat
+					res.test[x,2] <- model@test$pvalue
+					if(control.test$cb){
+						model <- plrs.cb(model, alpha=0.05)
+					}
+				}
+			}
+			if(control.output$save.models){
+				save(model, file = paste("plrsSeriesObjects/gene_",lab[x],".RData",sep=""))
+			}
+			if(control.output$plot.lin || control.output$save.plots){
+				if(control.output$save.plots){
+					if(control.output$type=="pdf") pdf(paste("plrsSeriesPlots/gene",lab[x],".pdf",sep=""))
+					if(control.output$type=="jpeg") jpeg(paste("plrsSeriesPlots/gene",lab[x],".jpeg",sep=""))
+				}
+				plot(model, lin=control.output$plot.lin)
+				if(control.output$save.plots) dev.off()
+			}
+		}
+	}
+	est <- round(est,8)
+
+	### Distribution of nb observations
+	general <- apply(config, 2, function(x)sum(!is.na(x)))
+	tp <- matrix(NA, 4, 4)
+	for(d in 1:ncol(config)){
+		if(general[d]!=0) tp[,d] <- summary(config[,d])[c(1,3,4,6)]
+	}
+	general <- rbind(general, tp)
+	colnames(general) <- colnames(config)
+	rownames(general) <- c("Nb genes","Min. obs","Median. obs","Mean. obs","Max. obs")
+	if(sum(general[1,]!=0)==1){
+		general <- t(as.matrix(general[1,]))
+		rownames(general) <- "Nb genes"
+	}
+
+	### Type of models
+	type.model <- as.matrix(colSums(modelTypes))
+	rownames(type.model) <- c("Intercept", "Simple linear", "Piecewise level", "Piecewise linear")
+	colnames(type.model) <- "NbModels"
+
+	### Testing
+	if(control.test$testing) res.test[,3] <- p.adjust(res.test[,2], method = "BH")
+	else{
+		res.test <- matrix(nrow=0,ncol=0)
+	}
+
+	cat("\n")
+	out <- new("plrs.series",
+			coefficients = est,
+			effects = list(I=Ieff, S=Seff),
+			test = res.test,
+			general = general,
+			modelsType = list(summary=type.model, all=modelTypes),
+			call.arg = list(control.model=control.model, control.select=control.select))
+	return(out)
+}
+
diff --git a/R/plrs.sim.r b/R/plrs.sim.r
new file mode 100644
index 0000000..479b650
--- /dev/null
+++ b/R/plrs.sim.r
@@ -0,0 +1,100 @@
+# Generate PLRS relationships
+
+plrs.sim <- function(n=80, states=4, sigma=1, x=NULL){
+
+	if(states==1){
+		int <- 1
+		slp <- 1
+		knot <- NULL
+		if(is.null(x)){
+			x <- runif(round(n), min = 0, max = 1)
+			x <- x[sample(1:length(x))]
+		}
+		y <- int[1] + slp[1]*x + rnorm(n, mean = 0, sd = sigma)
+		cal <- rep(0,length(x))
+
+		xtrue <- seq(0, 1, length=100)
+		ytrue <- int[1] + slp[1]*xtrue
+	}
+	if(states==2){
+		int <- c(1,1.5)
+		slp <- c(0,1)
+		knot <- 0.5
+		if(is.null(x)){
+			x1 <- runif(round(n/2), min = 0, max = 0.5)
+			x2 <- runif(round(n/2), min = 0.5, max = 1)
+			x <- c(x1,x2)
+		}
+		y <- x
+		y[x<=0.5] <- slp[1]*x[x<=0.5] + int[1]
+		y[x>0.5] <- slp[2]*x[x>0.5] + int[2] - slp[2]*0.5
+		y <- y + rnorm(n, mean = 0, sd = sigma)
+		cal <- rep(0,n)
+		cal[x>=0.5] <- 1
+		
+		xtrue <- seq(0, 1, length=100)
+		ytrue <- y
+		ytrue[xtrue<=0.5] <- slp[1]*xtrue[xtrue<=0.5] + int[1]
+		ytrue[xtrue>0.5] <- slp[2]*xtrue[xtrue>0.5] + int[2] - slp[2]*0.5
+	}
+	if(states==3){
+		int <- c(1,1.5,1.5)
+		slp <- c(1,0,1)
+		knot <- c(0.5,1)
+		if(is.null(x)){
+			x1 <- runif(round(n/3), min = 0, max = 0.5)
+			x2 <- runif(round(n/3), min = 0.5, max = 1)
+			x3 <- runif(round(n/3), min = 1, max = 1.5)
+			x <- c(x1,x2,x3)
+		}
+
+		y <- x
+		y[x<=0.5] <- slp[1]*x[x<=0.5] + int[1]
+		y[x>0.5&x<=1] <- slp[2]*x[x>0.5&x<=1] + int[2] - slp[2]*0.5
+		y[x>1] <- slp[3]*x[x>1] + int[3] - slp[3]*0.5
+		y <- y + rnorm(n, mean = 0, sd = sigma)
+		cal <- rep(-1,n)
+		cal[x>0.5&x<=1] <- 0
+		cal[x>1] <- 1
+
+		xtrue <- seq(0, 1.5, length=100)
+		ytrue <- y
+		ytrue[xtrue<=0.5] <- slp[1]*xtrue[xtrue<=0.5] + int[1]
+		ytrue[xtrue>0.5&xtrue<=1] <- slp[2]*xtrue[xtrue>0.5&xtrue<=1] + int[2] - slp[2]*0.5
+		ytrue[xtrue>1] <- slp[3]*xtrue[xtrue>1] + int[3] - slp[3]*0.5
+	}
+	if(states==4){
+		int <- c(1,1.5,1.5,1.75)
+		slp <- c(1,0,1,2)
+		knot <- c(0.5,1,1.5)
+		if(is.null(x)){
+			x1 <- runif(round(n/4), min = 0, max = 0.5)
+			x2 <- runif(round(n/4), min = 0.5, max = 1)
+			x3 <- runif(round(n/4), min = 1, max = 1.5)
+			x4 <- runif(round(n/4), min = 1.5, max = 2)
+			x <- c(x1,x2,x3,x4)
+		}
+
+		y <- x
+		y[x<=0.5] <- slp[1]*x[x<=0.5] + int[1]
+		y[x>0.5&x<=1] <- slp[2]*x[x>0.5&x<=1] + int[2]- slp[2]*0.5
+		y[x>1&x<=1.5] <- slp[3]*x[x>1&x<=1.5] + int[3] - slp[3]*0.5
+		y[x>1.5] <- slp[4]*x[x>1.5] + int[4] - slp[4]*1
+		y <- y + rnorm(n, mean = 0, sd = sigma)
+		cal <- rep(-1,n)
+		cal[x>0.5&x<=1] <- 0
+		cal[x>1&x<=1.5] <- 1
+		cal[x>1.5] <- 2
+
+		xtrue <- seq(0, 2, length=100) #c(seq(from=0.3, to=0.7, length=10),seq(from=0.8, to=1.2, length=10),seq(from=1.3, to=1.7, length=10))
+
+		ytrue <- rep(0, length(xtrue))
+		ytrue[xtrue<=0.5] <- slp[1]*xtrue[xtrue<=0.5] + int[1]
+		ytrue[xtrue>0.5&xtrue<=1] <- slp[2]*xtrue[xtrue>0.5&xtrue<=1] + int[2] - slp[2]*0.5
+		ytrue[xtrue>1&xtrue<=1.5] <- slp[3]*xtrue[xtrue>1&xtrue<=1.5] + int[3] - slp[3]*0.5
+		ytrue[xtrue>1.5] <- slp[4]*xtrue[xtrue>1.5] + int[4] - slp[4]*1
+	}
+
+	out <- list(seg=x, expr=y, cal=cal, knots=knot, xtrue=xtrue, ytrue=ytrue)
+	return(out)
+}
\ No newline at end of file
diff --git a/R/plrs.submodels.r b/R/plrs.submodels.r
new file mode 100644
index 0000000..fdeb227
--- /dev/null
+++ b/R/plrs.submodels.r
@@ -0,0 +1,940 @@
+# Internal function
+# Generate all submodels (all design and constraint matrices)
+
+# Author: Gwenael G.R. Leday
+
+.plrs.submodels <- function(object){
+
+	constr.slopes <- object@call.arg$constr.slopes
+	constr.intercepts <- object@call.arg$constr.intercepts
+	continuous <- object@call.arg$continuous
+	val <- object@mdata$mconf
+	nval <- length(val)
+
+	if(!object@call.arg$constr){ # If no constraints make standard search
+		constr.slopes <- 1
+		allConstr <- NULL
+	}
+	if(nval==1){
+		allLabs <- list("theta0.0")
+		allConstr <- list(NULL)
+		allTypes <- "Intercept"
+	}
+	if(nval==2){
+		if(all(val==c(-1,1))|all(val==c(-1,2))|all(val==c(1,2))) constr.slopes <- 1
+		if(constr.slopes==2){
+			if(all(val==c(-1,0))){
+				if(!continuous){
+					allLabs <- list("theta0.0", c("theta0.0", "theta0.1"), c("theta0.0", "theta1.0"),
+                                            c("theta0.0", "theta0.1", "theta1.0"), c("theta0.0", "theta0.1", "theta1.1"))
+					allTypes <- c("Intercept", "Linear", "Piecewise level", "Piecewise linear", "Piecewise linear")
+					if(constr.intercepts){
+						# loss-normal: 5 submodels
+						allConstr <- list(NULL,
+									matrix(c(0,1),1,2),
+									matrix(c(0,1),1,2),
+									matrix(c(0,1,0,0,0,1),2,3, byrow=TRUE),
+									matrix(c(0,1,1,0,0,-1),2,3, byrow=TRUE))
+					}
+					else{
+						# loss-normal: 5 submodels
+						allConstr <- list(NULL,
+									matrix(c(0,1),1,2),
+									NULL,
+									matrix(c(0,1,0),1,3, byrow=TRUE),
+									matrix(c(0,1,1,0,0,-1),2,3, byrow=TRUE))
+					}
+				}
+				else{
+					# loss-normal: 3 submodels
+					allLabs <- list("theta0.0", c("theta0.0", "theta0.1"), c("theta0.0", "theta0.1", "theta1.1"))
+					allConstr <- list(NULL,
+								matrix(c(0,1),1,2),
+								matrix(c(0,1,1,0,0,-1),2,3, byrow=TRUE))
+					allTypes <- c("Intercept", "Linear", "Piecewise linear")
+				}
+			}
+			if(all(val==c(0,1))|all(val==c(0,2))){
+				if(!continuous){
+					allLabs <- list("theta0.0", c("theta0.0", "theta0.1"), c("theta0.0", "theta1.0"),
+                                            c("theta0.0", "theta1.1"), c("theta0.0", "theta0.1", "theta1.0"),
+                                            c("theta0.0", "theta0.1", "theta1.1"), c("theta0.0", "theta1.0", "theta1.1"))
+					allTypes <- c("Intercept", "Linear", "Piecewise level", "Piecewise linear",
+                                            "Piecewise linear", "Piecewise linear", "Piecewise linear")
+					if(constr.intercepts){
+						# normal-gain: 7 submodels
+						allConstr <- list(NULL,
+									matrix(c(0,1),1,2),
+									matrix(c(0,1),1,2),
+									matrix(c(0,1),1,2),
+									matrix(c(0,1,0,0,0,1),2,3, byrow=TRUE),
+									matrix(c(0,1,0,0,0,1),2,3, byrow=TRUE),
+									matrix(c(0,1,0,0,0,1),2,3, byrow=TRUE))
+					}
+					else{
+						# normal-gain: 7 submodels
+						allConstr <- list(NULL,
+									matrix(c(0,1),1,2),
+									NULL,
+									matrix(c(0,1),1,2),
+									matrix(c(0,1,0),1,3, byrow=TRUE),
+									matrix(c(0,1,0,0,0,1),2,3, byrow=TRUE),
+									matrix(c(0,0,1),1,3, byrow=TRUE))
+					}
+				}
+				else{
+					# normal-gain: 4 submodels
+					allLabs <- list("theta0.0", c("theta0.0", "theta0.1"), c("theta0.0", "theta1.1"), c("theta0.0", "theta0.1", "theta1.1"))
+					allConstr <- list(NULL,
+								matrix(c(0,1),1,2),
+								matrix(c(0,1),1,2),
+								matrix(c(0,1,0,0,0,1),2,3, byrow=TRUE))
+					allTypes <- c("Intercept", "Linear", "Piecewise linear", "Piecewise linear")
+				}
+			}
+		}
+		if(constr.slopes==1){
+			if(!continuous){
+				allLabs <- list("theta0.0", c("theta0.0", "theta0.1"), c("theta0.0", "theta1.0"),
+                                      c("theta0.0", "theta1.1"), c("theta0.0", "theta0.1", "theta1.0"),
+                                      c("theta0.0", "theta0.1", "theta1.1"), c("theta0.0", "theta1.0", "theta1.1"))
+				allTypes <- c("Intercept", "Linear", "Piecewise level", "Piecewise linear",
+                                      "Piecewise linear", "Piecewise linear", "Piecewise linear")
+				if(object@call.arg$constr){
+					if(constr.intercepts){
+						# gain-amp: 7 submodels
+						allConstr <- list(NULL,
+									matrix(c(0,1),1,2),
+									matrix(c(0,1),1,2),
+									matrix(c(0,1),1,2),
+									matrix(c(0,1,0,0,0,1),2,3, byrow=TRUE),
+									matrix(c(0,1,0,0,1,1),2,3, byrow=TRUE),
+									matrix(c(0,1,0,0,0,1),2,3, byrow=TRUE))
+					}
+					else{
+						# gain-amp: 7 submodels
+						allConstr <- list(NULL,
+									matrix(c(0,1),1,2),
+									NULL,
+									matrix(c(0,1),1,2),
+									matrix(c(0,1,0),1,3, byrow=TRUE),
+									matrix(c(0,1,0,0,1,1),2,3, byrow=TRUE),
+									matrix(c(0,0,1),1,3, byrow=TRUE))
+					}
+				}
+			}
+			else{
+				# gain-amp: 4 submodels
+				allLabs <- list("theta0.0", c("theta0.0", "theta0.1"), c("theta0.0", "theta1.1"), c("theta0.0", "theta0.1", "theta1.1"))
+				if(object@call.arg$constr) allConstr <- list(NULL,
+							matrix(c(0,1),1,2),
+							matrix(c(0,1),1,2),
+							matrix(c(0,1,0,0,1,1),2,3, byrow=TRUE))
+				allTypes <- c("Intercept", "Linear", "Piecewise linear", "Piecewise linear", "Piecewise linear")
+			}
+		}
+	}
+	if(nval==3){
+		if(all(val==c(-1,1,2))) constr.slopes <- 1 
+		if(constr.slopes==2){
+			if(all(val==c(-1,0,1))|all(val==c(-1,0,2))){
+				if(!continuous){
+					# loss-normal-gain: 23 models
+					allLabs <- list("theta0.0",
+						c("theta0.0", "theta0.1"),
+						c("theta0.0", "theta1.0"),
+						c("theta0.0", "theta2.0"),
+						c("theta0.0", "theta2.1"),
+						c("theta0.0", "theta0.1", "theta1.0"),
+						c("theta0.0", "theta0.1", "theta1.1"),
+						c("theta0.0", "theta0.1", "theta2.0"),
+						c("theta0.0", "theta0.1", "theta2.1"),
+						c("theta0.0", "theta1.0", "theta2.0"),
+						c("theta0.0", "theta1.0", "theta2.1"),
+						c("theta0.0", "theta2.0", "theta2.1"),
+						c("theta0.0", "theta0.1", "theta1.0", "theta1.1"),
+						c("theta0.0", "theta0.1", "theta1.0", "theta2.0"),
+						c("theta0.0", "theta0.1", "theta1.0", "theta2.1"),
+						c("theta0.0", "theta0.1", "theta1.1", "theta2.0"),
+						c("theta0.0", "theta0.1", "theta1.1", "theta2.1"),
+						c("theta0.0", "theta0.1", "theta2.0", "theta2.1"),
+						c("theta0.0", "theta1.0", "theta2.0", "theta2.1"),
+						c("theta0.0", "theta0.1", "theta1.0", "theta1.1", "theta2.0"),
+						c("theta0.0", "theta0.1", "theta1.0", "theta1.1", "theta2.1"),
+						c("theta0.0", "theta0.1", "theta1.0", "theta2.0", "theta2.1"),
+						c("theta0.0", "theta0.1", "theta1.1", "theta2.0", "theta2.1"))
+					allTypes <- c("Intercept", "Linear", "Piecewise level", "Piecewise level", "Piecewise linear",
+                                               "Piecewise linear", "Piecewise linear", "Piecewise linear", "Piecewise linear",
+                                               "Piecewise level", "Piecewise linear", "Piecewise linear", "Piecewise linear",
+                                               "Piecewise linear", "Piecewise linear", "Piecewise linear", "Piecewise linear",
+                                               "Piecewise linear", "Piecewise linear", "Piecewise linear", "Piecewise linear",
+                                               "Piecewise linear", "Piecewise linear")
+					if(constr.intercepts){
+						allConstr <- list(NULL,
+							matrix(c(0,1),1,2),
+							matrix(c(0,1),1,2),
+							matrix(c(0,1),1,2),
+							matrix(c(0,1),1,2),
+							matrix(c(0,1,0,0,0,1),2,3, byrow=TRUE),
+							matrix(c(0,1,1,0,0,-1),2,3, byrow=TRUE),
+							matrix(c(0,1,0,0,0,1),2,3, byrow=TRUE),
+							matrix(c(0,1,0,0,0,1),2,3, byrow=TRUE),
+							matrix(c(0,1,0,0,0,1),2,3, byrow=TRUE),
+							matrix(c(0,1,0,0,0,1),2,3, byrow=TRUE),
+							matrix(c(0,1,0,0,0,1),2,3, byrow=TRUE),
+							matrix(c(0,1,0,1,0,0,0,-1,0,0,1,0),3,4, byrow=TRUE),
+							matrix(c(0,1,0,0,0,0,1,0,0,0,0,1),3,4, byrow=TRUE),
+							matrix(c(0,1,0,0,0,0,1,0,0,0,0,1),3,4, byrow=TRUE),
+							matrix(c(0,1,1,0,0,0,-1,0,0,0,0,1),3,4, byrow=TRUE),
+							matrix(c(0,1,1,0,0,0,-1,0,0,0,0,1),3,4, byrow=TRUE),
+							matrix(c(0,1,0,0,0,0,1,0,0,0,0,1),3,4, byrow=TRUE),
+							matrix(c(0,1,0,0,0,0,1,0,0,0,0,1),3,4, byrow=TRUE),
+							matrix(c(0,1,0,1,0, 0,0,0,-1,0, 0,0,1,0,0, 0,0,0,0,1),4,5, byrow=TRUE),
+							matrix(c(0,1,0,1,0, 0,0,0,-1,0, 0,0,1,0,0, 0,0,0,0,1),4,5, byrow=TRUE),
+							matrix(c(0,1,0,0,0, 0,0,1,0,0, 0,0,0,1,0, 0,0,0,0,1),4,5, byrow=TRUE),		
+							matrix(c(0,1,1,0,0, 0,0,-1,0,0, 0,0,0,1,0, 0,0,0,0,1),4,5, byrow=TRUE))
+					}
+					else{
+						allConstr <- list(NULL,
+							matrix(c(0,1),1,2),
+							NULL,
+							NULL,
+							matrix(c(0,1),1,2),
+							matrix(c(0,1,0),1,3, byrow=TRUE),
+							matrix(c(0,1,1,0,0,-1),2,3, byrow=TRUE),
+							matrix(c(0,1,0),1,3, byrow=TRUE),
+							matrix(c(0,1,0,0,0,1),2,3, byrow=TRUE),
+							NULL,
+							matrix(c(0,0,1),1,3, byrow=TRUE),
+							matrix(c(0,0,1),1,3, byrow=TRUE),
+							matrix(c(0,1,0,1,0,0,0,-1),2,4, byrow=TRUE),
+							matrix(c(0,1,0,0),1,4, byrow=TRUE),
+							matrix(c(0,1,0,0,0,0,0,1),2,4, byrow=TRUE),
+							matrix(c(0,1,1,0,0,0,-1,0),2,4, byrow=TRUE),
+							matrix(c(0,1,1,0,0,0,-1,0,0,0,0,1),3,4, byrow=TRUE),
+							matrix(c(0,1,0,0,0,0,0,1),2,4, byrow=TRUE),
+							matrix(c(0,0,0,1),2,4, byrow=TRUE),
+							matrix(c(0,1,0,1,0, 0,0,0,-1,0),2,5, byrow=TRUE),
+							matrix(c(0,1,0,1,0, 0,0,0,-1,0, 0,0,1,0,0),3,5, byrow=TRUE),
+							matrix(c(0,1,0,0,0, 0,0,0,0,1),2,5, byrow=TRUE),		
+							matrix(c(0,1,1,0,0, 0,0,-1,0,0, 0,0,0,0,1),3,5, byrow=TRUE))
+					}
+				}
+				else{
+					# loss-normal-gain: 6 models
+					allLabs <- list("theta0.0",
+                                            c("theta0.0", "theta0.1"),
+                                            c("theta0.0", "theta2.1"),
+                                            c("theta0.0", "theta0.1", "theta1.1"),
+                                            c("theta0.0", "theta0.1", "theta2.1"),
+                                            c("theta0.0", "theta0.1", "theta1.1", "theta2.1"))
+					allTypes <- c("Intercept", "Linear", "Piecewise linear", "Piecewise linear", "Piecewise linear", "Piecewise linear")
+					allConstr <- list(NULL,
+                                                matrix(c(0,1),1,2),
+                                                matrix(c(0,1),1,2),
+                                                matrix(c(0,1,1,0,0,-1),2,3, byrow=TRUE),
+                                                matrix(c(0,1,0,0,0,1),2,3, byrow=TRUE),
+                                                matrix(c(0,1,1,0,0,0,-1,0,0,0,0,1),3,4, byrow=TRUE))
+				}
+			}
+			if(all(val==c(0,1,2))){
+				if(!continuous){
+					# normal-gain-amp: 31 models
+					allLabs <- list("theta0.0",
+                                              c("theta0.0", "theta0.1"),
+                                              c("theta0.0", "theta1.0"),
+                                              c("theta0.0", "theta1.1"),
+                                              c("theta0.0", "theta2.0"),
+                                              c("theta0.0", "theta2.1"),
+                                              c("theta0.0", "theta0.1", "theta1.0"),
+                                              c("theta0.0", "theta0.1", "theta1.1"),
+                                              c("theta0.0", "theta0.1", "theta2.0"),
+                                              c("theta0.0", "theta0.1", "theta2.1"),
+                                              c("theta0.0", "theta1.0", "theta1.1"),
+                                              c("theta0.0", "theta1.0", "theta2.0"),
+                                              c("theta0.0", "theta1.0", "theta2.1"),
+                                              c("theta0.0", "theta1.1", "theta2.0"),
+                                              c("theta0.0", "theta1.1", "theta2.1"),
+                                              c("theta0.0", "theta2.0", "theta2.1"),
+                                              c("theta0.0", "theta0.1", "theta1.0", "theta1.1"),
+                                              c("theta0.0", "theta0.1", "theta1.0", "theta2.0"),
+                                              c("theta0.0", "theta0.1", "theta1.0", "theta2.1"),
+                                              c("theta0.0", "theta0.1", "theta1.1", "theta2.0"),
+                                              c("theta0.0", "theta0.1", "theta1.1", "theta2.1"),
+                                              c("theta0.0", "theta0.1", "theta2.0", "theta2.1"),
+                                              c("theta0.0", "theta1.0", "theta1.1", "theta2.0"),
+                                              c("theta0.0", "theta1.0", "theta1.1", "theta2.1"),
+                                              c("theta0.0", "theta1.0", "theta2.0", "theta2.1"),
+                                              c("theta0.0", "theta1.1", "theta2.0", "theta2.1"),
+                                              c("theta0.0", "theta0.1", "theta1.0", "theta1.1", "theta2.0"),
+                                              c("theta0.0", "theta0.1", "theta1.0", "theta1.1", "theta2.1"),
+                                              c("theta0.0", "theta0.1", "theta1.0", "theta2.0", "theta2.1"),
+                                              c("theta0.0", "theta0.1", "theta1.1", "theta2.0", "theta2.1"),
+                                              c("theta0.0", "theta1.0", "theta1.1", "theta2.0", "theta2.1"))
+					allTypes <- c("Intercept", "Linear", "Piecewise level", "Piecewise linear", "Piecewise level",
+                                            "Piecewise linear", "Piecewise linear", "Piecewise linear", "Piecewise linear",
+                                            "Piecewise linear", "Piecewise linear", "Piecewise level", "Piecewise linear",
+                                            "Piecewise linear", "Piecewise linear", "Piecewise linear", "Piecewise linear",
+                                            "Piecewise linear", "Piecewise linear", "Piecewise linear", "Piecewise linear",
+                                            "Piecewise linear", "Piecewise linear", "Piecewise linear", "Piecewise linear",
+                                            "Piecewise linear", "Piecewise linear", "Piecewise linear", "Piecewise linear",
+                                            "Piecewise linear", "Piecewise linear")
+					if(constr.intercepts){
+						allConstr <- list(NULL,
+			                                    matrix(c(0,1),1,2),
+			                                    matrix(c(0,1),1,2),
+			                                    matrix(c(0,1),1,2),
+			                                    matrix(c(0,1),1,2),
+			                                    matrix(c(0,1),1,2),
+			                                    matrix(c(0,1,0,0,0,1),2,3, byrow=TRUE),
+			                                    matrix(c(0,1,0,0,0,1),2,3, byrow=TRUE),
+			                                    matrix(c(0,1,0,0,0,1),2,3, byrow=TRUE),
+			                                    matrix(c(0,1,0,0,0,1),2,3, byrow=TRUE),
+			                                    matrix(c(0,1,0,0,0,1),2,3, byrow=TRUE),
+			                                    matrix(c(0,1,0,0,0,1),2,3, byrow=TRUE),
+			                                    matrix(c(0,1,0,0,0,1),2,3, byrow=TRUE),
+			                                    matrix(c(0,1,0,0,0,1),2,3, byrow=TRUE),
+			                                    matrix(c(0,1,0,0,1,1),2,3, byrow=TRUE),
+			                                    matrix(c(0,1,0,0,0,1),2,3, byrow=TRUE),
+			                                    matrix(c(0,1,0,0,0,0,1,0,0,0,0,1),3,4, byrow=TRUE),
+			                                    matrix(c(0,1,0,0,0,0,1,0,0,0,0,1),3,4, byrow=TRUE),
+			                                    matrix(c(0,1,0,0,0,0,1,0,0,0,0,1),3,4, byrow=TRUE),
+			                                    matrix(c(0,1,0,0,0,0,1,0,0,0,0,1),3,4, byrow=TRUE),
+			                                    matrix(c(0,1,0,0,0,0,1,0,0,0,1,1),3,4, byrow=TRUE),
+			                                    matrix(c(0,1,0,0,0,0,1,0,0,0,0,1),3,4, byrow=TRUE),
+			                                    matrix(c(0,1,0,0,0,0,1,0,0,0,0,1),3,4, byrow=TRUE),
+			                                    matrix(c(0,1,0,0,0,0,1,0,0,0,1,1),3,4, byrow=TRUE),
+			                                    matrix(c(0,1,0,0,0,0,1,0,0,0,0,1),3,4, byrow=TRUE),
+			                                    matrix(c(0,1,0,0,0,0,1,0,0,1,0,1),3,4, byrow=TRUE),
+			                                    matrix(c(0,1,0,0,0, 0,0,1,0,0, 0,0,0,1,0, 0,0,0,0,1),4,5, byrow=TRUE),
+			                                    matrix(c(0,1,0,0,0, 0,0,1,0,0, 0,0,0,1,0, 0,0,0,1,1),4,5, byrow=TRUE),
+			                                    matrix(c(0,1,0,0,0, 0,0,1,0,0, 0,0,0,1,0, 0,0,0,0,1),4,5, byrow=TRUE),		
+			                                    matrix(c(0,1,0,0,0, 0,0,1,0,0, 0,0,0,1,0, 0,0,1,0,1),4,5, byrow=TRUE),
+			                                    matrix(c(0,1,0,0,0, 0,0,1,0,0, 0,0,0,1,0, 0,0,1,0,1),4,5, byrow=TRUE))
+					}
+					else{
+						allConstr <- list(NULL,
+			                                    matrix(c(0,1),1,2),
+			                                    NULL,
+			                                    matrix(c(0,1),1,2),
+			                                    NULL,
+			                                    matrix(c(0,1),1,2),
+			                                    matrix(c(0,1,0),1,3, byrow=TRUE),
+			                                    matrix(c(0,1,0,0,0,1),2,3, byrow=TRUE),
+			                                    matrix(c(0,1,0),1,3, byrow=TRUE),
+			                                    matrix(c(0,1,0,0,0,1),2,3, byrow=TRUE),
+			                                    matrix(c(0,0,1),1,3, byrow=TRUE),
+			                                    NULL,
+			                                    matrix(c(0,0,1),1,3, byrow=TRUE),
+			                                    matrix(c(0,1,0),1,3, byrow=TRUE),
+			                                    matrix(c(0,1,0,0,1,1),2,3, byrow=TRUE),
+			                                    matrix(c(0,0,1),1,3, byrow=TRUE),
+			                                    matrix(c(0,1,0,0, 0,0,0,1),2,4, byrow=TRUE),
+			                                    matrix(c(0,1,0,0),1,4, byrow=TRUE),
+			                                    matrix(c(0,1,0,0, 0,0,0,1),2,4, byrow=TRUE),
+			                                    matrix(c(0,1,0,0,0,0,1,0),2,4, byrow=TRUE),
+			                                    matrix(c(0,1,0,0,0,0,1,0,0,0,1,1),3,4, byrow=TRUE),
+			                                    matrix(c(0,1,0,0,0,0,0,1),2,4, byrow=TRUE),
+			                                    matrix(c(0,0,1,0),1,4, byrow=TRUE),
+			                                    matrix(c(0,0,1,0,0,0,1,1),2,4, byrow=TRUE),
+			                                    matrix(c(0,0,0,1),1,4, byrow=TRUE),
+			                                    matrix(c(0,1,0,0, 0,1,0,1),2,4, byrow=TRUE),
+			                                    matrix(c(0,1,0,0,0, 0,0,0,1,0),2,5, byrow=TRUE),
+			                                    matrix(c(0,1,0,0,0, 0,0,0,1,0, 0,0,0,1,1),3,5, byrow=TRUE),
+			                                    matrix(c(0,1,0,0,0, 0,0,0,0,1),2,5, byrow=TRUE),		
+			                                    matrix(c(0,1,0,0,0, 0,0,1,0,0, 0,0,1,0,1),3,5, byrow=TRUE),
+			                                    matrix(c(0,0,1,0,0, 0,0,1,0,1),2,5, byrow=TRUE))
+					}
+				}
+				else{
+					# normal-gain-amp: 8 models
+					allLabs <- list("theta0.0",
+                                            c("theta0.0", "theta0.1"),
+                                            c("theta0.0", "theta1.1"),
+                                            c("theta0.0", "theta2.1"),
+                                            c("theta0.0", "theta0.1", "theta1.1"),
+                                            c("theta0.0", "theta0.1", "theta2.1"),
+                                            c("theta0.0", "theta1.1", "theta2.1"),
+                                            c("theta0.0", "theta0.1", "theta1.1", "theta2.1"))
+					allTypes <- c("Intercept", "Linear", "Piecewise linear", "Piecewise linear", "Piecewise linear",
+                                            "Piecewise linear", "Piecewise linear", "Piecewise linear")
+					allConstr <- list(NULL,
+		                                    matrix(c(0,1),1,2),
+		                                    matrix(c(0,1),1,2),
+		                                    matrix(c(0,1),1,2),
+		                                    matrix(c(0,1,0,0,0,1),2,3, byrow=TRUE),
+		                                    matrix(c(0,1,0,0,0,1),2,3, byrow=TRUE),
+		                                    matrix(c(0,1,0,0,1,1),2,3, byrow=TRUE),
+		                                    matrix(c(0,1,0,0,0,0,1,0,0,0,1,1),3,4, byrow=TRUE))
+				}
+			}
+		}
+		if(constr.slopes==1){
+			if(!continuous){
+				# 31 models
+				allLabs <- list(c("theta0.0"),
+                                        c("theta0.0", "theta0.1"),
+                                        c("theta0.0", "theta1.0"),
+                                        c("theta0.0", "theta1.1"),
+                                        c("theta0.0", "theta2.0"),
+                                        c("theta0.0", "theta2.1"),
+                                        c("theta0.0", "theta0.1", "theta1.0"),
+                                        c("theta0.0", "theta0.1", "theta1.1"),
+                                        c("theta0.0", "theta0.1", "theta2.0"),
+                                        c("theta0.0", "theta0.1", "theta2.1"),
+                                        c("theta0.0", "theta1.0", "theta1.1"),
+                                        c("theta0.0", "theta1.0", "theta2.0"),
+                                        c("theta0.0", "theta1.0", "theta2.1"),
+                                        c("theta0.0", "theta1.1", "theta2.0"),
+                                        c("theta0.0", "theta1.1", "theta2.1"),
+                                        c("theta0.0", "theta2.0", "theta2.1"),
+                                        c("theta0.0", "theta0.1", "theta1.0", "theta1.1"),
+                                        c("theta0.0", "theta0.1", "theta1.0", "theta2.0"),
+                                        c("theta0.0", "theta0.1", "theta1.0", "theta2.1"),
+                                        c("theta0.0", "theta0.1", "theta1.1", "theta2.0"),
+                                        c("theta0.0", "theta0.1", "theta1.1", "theta2.1"),
+                                        c("theta0.0", "theta0.1", "theta2.0", "theta2.1"),
+                                        c("theta0.0", "theta1.0", "theta1.1", "theta2.0"),
+                                        c("theta0.0", "theta1.0", "theta1.1", "theta2.1"),
+                                        c("theta0.0", "theta1.0", "theta2.0", "theta2.1"),
+                                        c("theta0.0", "theta1.1", "theta2.0", "theta2.1"),
+                                        c("theta0.0", "theta0.1", "theta1.0", "theta1.1", "theta2.0"),
+                                        c("theta0.0", "theta0.1", "theta1.0", "theta1.1", "theta2.1"),
+                                        c("theta0.0", "theta0.1", "theta1.0", "theta2.0", "theta2.1"),
+                                        c("theta0.0", "theta0.1", "theta1.1", "theta2.0", "theta2.1"),
+                                        c("theta0.0", "theta1.0", "theta1.1", "theta2.0", "theta2.1"))
+				allTypes <- c("Intercept", "Linear", "Piecewise level", "Piecewise linear", "Piecewise level",
+                                      rep("Piecewise linear", 6), "Piecewise level", rep("Piecewise linear", 19))
+				if(object@call.arg$constr){
+					if(constr.intercepts){
+						allConstr <- list(NULL,
+		                                          matrix(c(0,1),1,2),
+		                                          matrix(c(0,1),1,2),
+		                                          matrix(c(0,1),1,2),
+      		                                    matrix(c(0,1),1,2),
+      		                                    matrix(c(0,1),1,2),
+      		                                    matrix(c(0,1,0,0,0,1),2,3, byrow=TRUE),
+      		                                    matrix(c(0,1,0,0,1,1),2,3, byrow=TRUE),
+      		                                    matrix(c(0,1,0,0,0,1),2,3, byrow=TRUE),
+      		                                    matrix(c(0,1,0,0,1,1),2,3, byrow=TRUE),
+      		                                    matrix(c(0,1,0,0,0,1),2,3, byrow=TRUE),
+      		                                    matrix(c(0,1,0,0,0,1),2,3, byrow=TRUE),
+      		                                    matrix(c(0,1,0,0,0,1),2,3, byrow=TRUE),
+      		                                    matrix(c(0,1,0,0,0,1),2,3, byrow=TRUE),
+      		                                    matrix(c(0,1,0,0,1,1),2,3, byrow=TRUE),
+      		                                    matrix(c(0,1,0,0,0,1),2,3, byrow=TRUE),
+      		                                    matrix(c(0,1,0,0, 0,1,0,1, 0,0,1,0),3,4, byrow=TRUE),
+      		                                    matrix(c(0,1,0,0, 0,0,1,0, 0,0,0,1),3,4, byrow=TRUE),
+      		                                    matrix(c(0,1,0,0, 0,1,0,1, 0,0,1,0),3,4, byrow=TRUE),
+      		                                    matrix(c(0,1,0,0, 0,1,1,0, 0,0,0,1),3,4, byrow=TRUE),
+      		                                    matrix(c(0,1,0,0, 0,1,1,0, 0,1,1,1),3,4, byrow=TRUE),
+      		                                    matrix(c(0,1,0,0, 0,1,0,1, 0,0,1,0),3,4, byrow=TRUE),
+      		                                    matrix(c(0,1,0,0, 0,0,1,0, 0,0,0,1),3,4, byrow=TRUE),
+      		                                    matrix(c(0,1,0,0, 0,0,1,0, 0,0,1,1),3,4, byrow=TRUE),
+      		                                    matrix(c(0,1,0,0, 0,0,1,0, 0,0,0,1),3,4, byrow=TRUE),
+      		                                    matrix(c(0,1,0,0, 0,1,0,1, 0,0,1,0),3,4, byrow=TRUE),
+      		                                    matrix(c(0,1,0,0,0, 0,1,0,1,0, 0,0,1,0,0, 0,0,0,0,1),4,5, byrow=TRUE),
+      		                                    matrix(c(0,1,0,0,0, 0,1,0,1,0, 0,1,0,1,1, 0,0,1,0,0),4,5, byrow=TRUE),
+      		                                    matrix(c(0,1,0,0,0, 0,1,0,0,1, 0,0,1,0,0, 0,0,0,1,0),4,5, byrow=TRUE),
+      		                                    matrix(c(0,1,0,0,0, 0,1,1,0,0, 0,1,1,0,1, 0,0,0,1,0),4,5, byrow=TRUE),
+      		                                    matrix(c(0,1,0,0,0, 0,0,1,0,0, 0,0,1,0,1, 0,0,0,1,0),4,5, byrow=TRUE))
+					}
+					else{
+						allConstr <- list(NULL,
+      		                                    matrix(c(0,1),1,2),
+      		                                    NULL,
+      		                                    matrix(c(0,1),1,2),
+      		                                    NULL,
+      		                                    matrix(c(0,1),1,2),
+      		                                    matrix(c(0,1,0),1,3, byrow=TRUE),
+      		                                    matrix(c(0,1,0,0,1,1),2,3, byrow=TRUE),
+      		                                    matrix(c(0,1,0),1,3, byrow=TRUE),
+      		                                    matrix(c(0,1,0,0,1,1),2,3, byrow=TRUE),
+      		                                    matrix(c(0,0,1),1,3, byrow=TRUE),
+      		                                    NULL,
+      		                                    matrix(c(0,0,1),1,3, byrow=TRUE),
+      		                                    matrix(c(0,1,0),1,3, byrow=TRUE),
+      		                                    matrix(c(0,1,0,0,1,1),2,3, byrow=TRUE),
+      		                                    matrix(c(0,0,1),1,3, byrow=TRUE),
+      		                                    matrix(c(0,1,0,0, 0,1,0,1),2,4, byrow=TRUE),
+      		                                    matrix(c(0,1,0,0),1,4, byrow=TRUE),
+      		                                    matrix(c(0,1,0,0, 0,1,0,1),2,4, byrow=TRUE),
+      		                                    matrix(c(0,1,0,0, 0,1,1,0),2,4, byrow=TRUE),
+      		                                    matrix(c(0,1,0,0, 0,1,1,0, 0,1,1,1),3,4, byrow=TRUE),
+      		                                    matrix(c(0,1,0,0, 0,1,0,1),2,4, byrow=TRUE),
+      		                                    matrix(c(0,0,1,0),1,4, byrow=TRUE),
+      		                                    matrix(c(0,0,1,0, 0,0,1,1),2,4, byrow=TRUE),
+      		                                    matrix(c(0,0,0,1),1,4, byrow=TRUE),
+      		                                    matrix(c(0,1,0,0, 0,1,0,1),2,4, byrow=TRUE),
+      		                                    matrix(c(0,1,0,0,0, 0,1,0,1,0),2,5, byrow=TRUE),
+      		                                    matrix(c(0,1,0,0,0, 0,1,0,1,0, 0,1,0,1,1),3,5, byrow=TRUE),
+      		                                    matrix(c(0,1,0,0,0, 0,1,0,0,1),2,5, byrow=TRUE),
+      		                                    matrix(c(0,1,0,0,0, 0,1,1,0,0, 0,1,1,0,1),3,5, byrow=TRUE),
+      		                                    matrix(c(0,0,1,0,0, 0,0,1,0,1, 0,0,0,1,0),3,5, byrow=TRUE))
+					}
+				}
+			}
+			else{
+				allLabs <- list(c("theta0.0"),
+                                        c("theta0.0", "theta0.1"),
+                                        c("theta0.0", "theta1.1"),
+                                        c("theta0.0", "theta2.1"),
+                                        c("theta0.0", "theta0.1", "theta1.1"),
+                                        c("theta0.0", "theta0.1", "theta2.1"),
+                                        c("theta0.0", "theta1.1", "theta2.1"),
+                                        c("theta0.0", "theta0.1", "theta1.1", "theta2.1"))
+				if(object@call.arg$constr) allConstr <- list(NULL,
+		                              matrix(c(0,1),1,2),
+		                              matrix(c(0,1),1,2),
+		                              matrix(c(0,1),1,2),
+		                              matrix(c(0,1,0,0,1,1),2,3, byrow=TRUE),
+		                              matrix(c(0,1,0,0,1,1),2,3, byrow=TRUE),
+		                              matrix(c(0,1,0,0,1,1),2,3, byrow=TRUE),
+		                              matrix(c(0,1,0,0, 0,1,1,0, 0,1,1,1),3,4, byrow=TRUE))
+				allTypes <- c("Intercept", "Linear", "Piecewise linear", "Piecewise linear", 
+                                      "Piecewise linear", "Piecewise linear", "Piecewise linear", "Piecewise linear")
+			}
+		}
+	}
+	if(nval==4){
+		if(constr.slopes==2){
+			if(!continuous){
+				# 95 models
+				allLabs <- list(
+					c("theta0.0"),
+					c("theta0.0", "theta0.1"),
+					c("theta0.0", "theta1.0"),
+					c("theta0.0", "theta2.0"),
+					c("theta0.0", "theta2.1"),
+					c("theta0.0", "theta3.0"),
+					c("theta0.0", "theta3.1"),
+					c("theta0.0", "theta0.1", "theta1.0"),
+					c("theta0.0", "theta0.1", "theta1.1"),
+					c("theta0.0", "theta0.1", "theta2.0"),
+					c("theta0.0", "theta0.1", "theta2.1"),
+					c("theta0.0", "theta0.1", "theta3.0"),
+					c("theta0.0", "theta0.1", "theta3.1"),
+					c("theta0.0", "theta1.0", "theta2.0"),
+					c("theta0.0", "theta1.0", "theta2.1"),
+					c("theta0.0", "theta1.0", "theta3.0"),
+					c("theta0.0", "theta1.0", "theta3.1"),
+					c("theta0.0", "theta2.0", "theta2.1"),
+					c("theta0.0", "theta2.0", "theta3.0"),
+					c("theta0.0", "theta2.0", "theta3.1"),
+					c("theta0.0", "theta2.1", "theta3.0"),
+					c("theta0.0", "theta2.1", "theta3.1"),
+					c("theta0.0", "theta3.0", "theta3.1"),
+					c("theta0.0", "theta0.1", "theta1.0", "theta1.1"),
+					c("theta0.0", "theta0.1", "theta1.0", "theta2.0"),
+					c("theta0.0", "theta0.1", "theta1.0", "theta2.1"),
+					c("theta0.0", "theta0.1", "theta1.0", "theta3.0"),
+					c("theta0.0", "theta0.1", "theta1.0", "theta3.1"),
+					c("theta0.0", "theta0.1", "theta1.1", "theta2.0"),
+					c("theta0.0", "theta0.1", "theta1.1", "theta2.1"),
+					c("theta0.0", "theta0.1", "theta1.1", "theta3.0"),
+					c("theta0.0", "theta0.1", "theta1.1", "theta3.1"),
+					c("theta0.0", "theta0.1", "theta2.0", "theta2.1"),
+					c("theta0.0", "theta0.1", "theta2.0", "theta3.0"),
+					c("theta0.0", "theta0.1", "theta2.0", "theta3.1"),
+					c("theta0.0", "theta0.1", "theta2.1", "theta3.0"),
+					c("theta0.0", "theta0.1", "theta2.1", "theta3.1"),
+					c("theta0.0", "theta0.1", "theta3.0", "theta3.1"),
+					c("theta0.0", "theta1.0", "theta2.0", "theta2.1"),
+					c("theta0.0", "theta1.0", "theta2.0", "theta3.0"),
+					c("theta0.0", "theta1.0", "theta2.0", "theta3.1"),
+					c("theta0.0", "theta1.0", "theta2.1", "theta3.0"),
+					c("theta0.0", "theta1.0", "theta2.1", "theta3.1"),
+					c("theta0.0", "theta1.0", "theta3.0", "theta3.1"),
+					c("theta0.0", "theta2.0", "theta2.1", "theta3.0"),
+					c("theta0.0", "theta2.0", "theta2.1", "theta3.1"),
+					c("theta0.0", "theta2.0", "theta3.0", "theta3.1"),
+					c("theta0.0", "theta2.1", "theta3.0", "theta3.1"),
+					c("theta0.0", "theta0.1", "theta1.0", "theta1.1", "theta2.0"),
+					c("theta0.0", "theta0.1", "theta1.0", "theta1.1", "theta2.1"),
+					c("theta0.0", "theta0.1", "theta1.0", "theta1.1", "theta3.0"),
+					c("theta0.0", "theta0.1", "theta1.0", "theta1.1", "theta3.1"),
+					c("theta0.0", "theta0.1", "theta1.0", "theta2.0", "theta2.1"),
+					c("theta0.0", "theta0.1", "theta1.0", "theta2.0", "theta3.0"),
+					c("theta0.0", "theta0.1", "theta1.0", "theta2.0", "theta3.1"),
+					c("theta0.0", "theta0.1", "theta1.0", "theta2.1", "theta3.0"),
+					c("theta0.0", "theta0.1", "theta1.0", "theta2.1", "theta3.1"),
+					c("theta0.0", "theta0.1", "theta1.0", "theta3.0", "theta3.1"),
+					c("theta0.0", "theta0.1", "theta1.1", "theta2.0", "theta2.1"),
+					c("theta0.0", "theta0.1", "theta1.1", "theta2.0", "theta3.0"),
+					c("theta0.0", "theta0.1", "theta1.1", "theta2.0", "theta3.1"),
+					c("theta0.0", "theta0.1", "theta1.1", "theta2.1", "theta3.0"),
+					c("theta0.0", "theta0.1", "theta1.1", "theta2.1", "theta3.1"),
+					c("theta0.0", "theta0.1", "theta1.1", "theta3.0", "theta3.1"),
+					c("theta0.0", "theta0.1", "theta2.0", "theta2.1", "theta3.0"),
+					c("theta0.0", "theta0.1", "theta2.0", "theta2.1", "theta3.1"),
+					c("theta0.0", "theta0.1", "theta2.0", "theta3.0", "theta3.1"),
+					c("theta0.0", "theta0.1", "theta2.1", "theta3.0", "theta3.1"),
+					c("theta0.0", "theta1.0", "theta2.0", "theta2.1", "theta3.0"),
+					c("theta0.0", "theta1.0", "theta2.0", "theta2.1", "theta3.1"),
+					c("theta0.0", "theta1.0", "theta2.0", "theta3.0", "theta3.1"),
+					c("theta0.0", "theta1.0", "theta2.1", "theta3.0", "theta3.1"),
+					c("theta0.0", "theta2.0", "theta2.1", "theta3.0", "theta3.1"),
+					c("theta0.0", "theta0.1", "theta1.0", "theta1.1", "theta2.0", "theta2.1"),
+					c("theta0.0", "theta0.1", "theta1.0", "theta1.1", "theta2.0", "theta3.0"),
+					c("theta0.0", "theta0.1", "theta1.0", "theta1.1", "theta2.0", "theta3.1"),
+					c("theta0.0", "theta0.1", "theta1.0", "theta1.1", "theta2.1", "theta3.0"),
+					c("theta0.0", "theta0.1", "theta1.0", "theta1.1", "theta2.1", "theta3.1"),
+					c("theta0.0", "theta0.1", "theta1.0", "theta1.1", "theta3.0", "theta3.1"),
+					c("theta0.0", "theta0.1", "theta1.0", "theta2.0", "theta2.1", "theta3.0"),
+					c("theta0.0", "theta0.1", "theta1.0", "theta2.0", "theta2.1", "theta3.1"),
+					c("theta0.0", "theta0.1", "theta1.0", "theta2.0", "theta3.0", "theta3.1"),
+					c("theta0.0", "theta0.1", "theta1.0", "theta2.1", "theta3.0", "theta3.1"),
+					c("theta0.0", "theta0.1", "theta1.1", "theta2.0", "theta2.1", "theta3.0"),
+					c("theta0.0", "theta0.1", "theta1.1", "theta2.0", "theta2.1", "theta3.1"),
+					c("theta0.0", "theta0.1", "theta1.1", "theta2.0", "theta3.0", "theta3.1"),
+					c("theta0.0", "theta0.1", "theta1.1", "theta2.1", "theta3.0", "theta3.1"),
+					c("theta0.0", "theta0.1", "theta2.0", "theta2.1", "theta3.0", "theta3.1"),
+					c("theta0.0", "theta1.0", "theta2.0", "theta2.1", "theta3.0", "theta3.1"),
+					c("theta0.0", "theta0.1", "theta1.0", "theta1.1", "theta2.0", "theta2.1", "theta3.0"),
+					c("theta0.0", "theta0.1", "theta1.0", "theta1.1", "theta2.0", "theta2.1", "theta3.1"),
+					c("theta0.0", "theta0.1", "theta1.0", "theta1.1", "theta2.0", "theta3.0", "theta3.1"),
+					c("theta0.0", "theta0.1", "theta1.0", "theta1.1", "theta2.1", "theta3.0", "theta3.1"),
+					c("theta0.0", "theta0.1", "theta1.0", "theta2.0", "theta2.1", "theta3.0", "theta3.1"),
+					c("theta0.0", "theta0.1", "theta1.1", "theta2.0", "theta2.1", "theta3.0", "theta3.1"))
+				allTypes <- c("Intercept", "Linear",
+                                  rep("Piecewise level",2),
+                                  "Piecewise linear",
+                                  "Piecewise level",
+                                  rep("Piecewise linear",7),
+                                  "Piecewise level",
+                                  "Piecewise linear",
+                                  "Piecewise level",
+                                  rep("Piecewise linear",2),
+                                  "Piecewise level",
+                                  rep("Piecewise linear",20),
+                                  "Piecewise level",
+                                  rep("Piecewise linear",55))
+				if(constr.intercepts){
+					allConstr <- list(
+						NULL,
+						matrix(c(0,1),1,2),
+						matrix(c(0,1),1,2),
+						matrix(c(0,1),1,2),
+						matrix(c(0,1),1,2),
+						matrix(c(0,1),1,2),
+						matrix(c(0,1),1,2),
+						matrix(c(0,1,0,0,0,1),2,3, byrow=TRUE),
+						matrix(c(0,1,1,0,0,-1),2,3, byrow=TRUE),
+						matrix(c(0,1,0,0,0,1),2,3, byrow=TRUE),
+						matrix(c(0,1,0,0,0,1),2,3, byrow=TRUE),
+						matrix(c(0,1,0,0,0,1),2,3, byrow=TRUE),
+						matrix(c(0,1,0,0,0,1),2,3, byrow=TRUE),
+						matrix(c(0,1,0,0,0,1),2,3, byrow=TRUE),
+						matrix(c(0,1,0,0,0,1),2,3, byrow=TRUE),
+						matrix(c(0,1,0,0,0,1),2,3, byrow=TRUE),
+						matrix(c(0,1,0,0,0,1),2,3, byrow=TRUE),
+						matrix(c(0,1,0,0,0,1),2,3, byrow=TRUE),
+						matrix(c(0,1,0,0,0,1),2,3, byrow=TRUE),
+						matrix(c(0,1,0,0,0,1),2,3, byrow=TRUE),
+						matrix(c(0,1,0,0,0,1),2,3, byrow=TRUE),
+						matrix(c(0,1,0,0,1,1),2,3, byrow=TRUE),
+						matrix(c(0,1,0,0,0,1),2,3, byrow=TRUE),
+						matrix(c(0,1,0,1, 0,0,1,0, 0,0,0,-1),3,4, byrow=TRUE),
+						matrix(c(0,1,0,0, 0,0,1,0, 0,0,0,1),3,4, byrow=TRUE),
+						matrix(c(0,1,0,0, 0,0,1,0, 0,0,0,1),3,4, byrow=TRUE),
+						matrix(c(0,1,0,0, 0,0,1,0, 0,0,0,1),3,4, byrow=TRUE),
+						matrix(c(0,1,0,0, 0,0,1,0, 0,0,0,1),3,4, byrow=TRUE),
+						matrix(c(0,1,1,0, 0,0,-1,0, 0,0,0,1),3,4, byrow=TRUE),
+						matrix(c(0,1,1,0, 0,0,-1,0, 0,0,0,1),3,4, byrow=TRUE),
+						matrix(c(0,1,1,0, 0,0,-1,0, 0,0,0,1),3,4, byrow=TRUE),
+						matrix(c(0,1,1,0, 0,0,-1,0, 0,0,0,1),3,4, byrow=TRUE),
+						matrix(c(0,1,0,0, 0,0,1,0, 0,0,0,1),3,4, byrow=TRUE),
+						matrix(c(0,1,0,0, 0,0,1,0, 0,0,0,1),3,4, byrow=TRUE),
+						matrix(c(0,1,0,0, 0,0,1,0, 0,0,0,1),3,4, byrow=TRUE),
+						matrix(c(0,1,0,0, 0,0,1,0, 0,0,0,1),3,4, byrow=TRUE),
+						matrix(c(0,1,0,0, 0,0,1,0, 0,0,1,1),3,4, byrow=TRUE),
+						matrix(c(0,1,0,0, 0,0,1,0, 0,0,0,1),3,4, byrow=TRUE),
+						matrix(c(0,1,0,0, 0,0,1,0, 0,0,0,1),3,4, byrow=TRUE),
+						matrix(c(0,1,0,0, 0,0,1,0, 0,0,0,1),3,4, byrow=TRUE),
+						matrix(c(0,1,0,0, 0,0,1,0, 0,0,0,1),3,4, byrow=TRUE),
+						matrix(c(0,1,0,0, 0,0,1,0, 0,0,0,1),3,4, byrow=TRUE),
+						matrix(c(0,1,0,0, 0,0,1,0, 0,0,1,1),3,4, byrow=TRUE),
+						matrix(c(0,1,0,0, 0,0,1,0, 0,0,0,1),3,4, byrow=TRUE),
+						matrix(c(0,1,0,0, 0,0,1,0, 0,0,0,1),3,4, byrow=TRUE),
+						matrix(c(0,1,0,0, 0,0,1,0, 0,0,1,1),3,4, byrow=TRUE),
+						matrix(c(0,1,0,0, 0,0,1,0, 0,0,0,1),3,4, byrow=TRUE),
+						matrix(c(0,1,0,0, 0,0,1,0, 0,0,0,1),3,4, byrow=TRUE),
+						matrix(c(0,1,0,1,0, 0,0,1,0,0, 0,0,0,-1,0, 0,0,0,0,1),4,5, byrow=TRUE),
+						matrix(c(0,1,0,1,0, 0,0,1,0,0, 0,0,0,-1,0, 0,0,0,0,1),4,5, byrow=TRUE),
+						matrix(c(0,1,0,1,0, 0,0,1,0,0, 0,0,0,-1,0, 0,0,0,0,1),4,5, byrow=TRUE),
+						matrix(c(0,1,0,1,0, 0,0,1,0,0, 0,0,0,-1,0, 0,0,0,0,1),4,5, byrow=TRUE),
+						matrix(c(0,1,0,0,0, 0,0,1,0,0, 0,0,0,1,0, 0,0,0,0,1),4,5, byrow=TRUE),
+						matrix(c(0,1,0,0,0, 0,0,1,0,0, 0,0,0,1,0, 0,0,0,0,1),4,5, byrow=TRUE),
+						matrix(c(0,1,0,0,0, 0,0,1,0,0, 0,0,0,1,0, 0,0,0,0,1),4,5, byrow=TRUE),
+						matrix(c(0,1,0,0,0, 0,0,1,0,0, 0,0,0,1,0, 0,0,0,0,1),4,5, byrow=TRUE),
+						matrix(c(0,1,0,0,0, 0,0,1,0,0, 0,0,0,1,0, 0,0,0,1,1),4,5, byrow=TRUE),
+						matrix(c(0,1,0,0,0, 0,0,1,0,0, 0,0,0,1,0, 0,0,0,0,1),4,5, byrow=TRUE),
+						matrix(c(0,1,1,0,0, 0,0,-1,0,0, 0,0,0,1,0, 0,0,0,0,1),4,5, byrow=TRUE),
+						matrix(c(0,1,1,0,0, 0,0,-1,0,0, 0,0,0,1,0, 0,0,0,0,1),4,5, byrow=TRUE),
+						matrix(c(0,1,1,0,0, 0,0,-1,0,0, 0,0,0,1,0, 0,0,0,0,1),4,5, byrow=TRUE),
+						matrix(c(0,1,1,0,0, 0,0,-1,0,0, 0,0,0,1,0, 0,0,0,0,1),4,5, byrow=TRUE),
+						matrix(c(0,1,1,0,0, 0,0,-1,0,0, 0,0,0,1,0, 0,0,0,1,1),4,5, byrow=TRUE),
+						matrix(c(0,1,1,0,0, 0,0,-1,0,0, 0,0,0,1,0, 0,0,0,0,1),4,5, byrow=TRUE),
+						matrix(c(0,1,0,0,0, 0,0,1,0,0, 0,0,0,1,0, 0,0,0,0,1),4,5, byrow=TRUE),
+						matrix(c(0,1,0,0,0, 0,0,1,0,0, 0,0,0,1,0, 0,0,0,1,1),4,5, byrow=TRUE),
+						matrix(c(0,1,0,0,0, 0,0,1,0,0, 0,0,0,1,0, 0,0,0,0,1),4,5, byrow=TRUE),
+						matrix(c(0,1,0,0,0, 0,0,1,0,0, 0,0,0,1,0, 0,0,1,0,1),4,5, byrow=TRUE),
+						matrix(c(0,1,0,0,0, 0,0,1,0,0, 0,0,0,1,0, 0,0,0,0,1),4,5, byrow=TRUE),
+						matrix(c(0,1,0,0,0, 0,0,1,0,0, 0,0,0,1,0, 0,0,0,1,1),4,5, byrow=TRUE),
+						matrix(c(0,1,0,0,0, 0,0,1,0,0, 0,0,0,1,0, 0,0,0,0,1),4,5, byrow=TRUE),
+						matrix(c(0,1,0,0,0, 0,0,1,0,0, 0,0,0,1,0, 0,0,1,0,1),4,5, byrow=TRUE),
+						matrix(c(0,1,0,0,0, 0,0,1,0,0, 0,0,0,1,0, 0,0,1,0,1),4,5, byrow=TRUE),
+						matrix(c(0,1,0,1,0,0, 0,0,1,0,0,0, 0,0,0,-1,0,0, 0,0,0,0,1,0, 0,0,0,0,0,1),5,6, byrow=TRUE),
+						matrix(c(0,1,0,1,0,0, 0,0,1,0,0,0, 0,0,0,-1,0,0, 0,0,0,0,1,0, 0,0,0,0,0,1),5,6, byrow=TRUE),
+						matrix(c(0,1,0,1,0,0, 0,0,1,0,0,0, 0,0,0,-1,0,0, 0,0,0,0,1,0, 0,0,0,0,0,1),5,6, byrow=TRUE),
+						matrix(c(0,1,0,1,0,0, 0,0,1,0,0,0, 0,0,0,-1,0,0, 0,0,0,0,1,0, 0,0,0,0,0,1),5,6, byrow=TRUE),
+						matrix(c(0,1,0,1,0,0, 0,0,1,0,0,0, 0,0,0,-1,0,0, 0,0,0,0,1,0, 0,0,0,0,1,1),5,6, byrow=TRUE),
+						matrix(c(0,1,0,1,0,0, 0,0,1,0,0,0, 0,0,0,-1,0,0, 0,0,0,0,1,0, 0,0,0,0,0,1),5,6, byrow=TRUE),
+						matrix(c(0,1,0,0,0,0, 0,0,1,0,0,0, 0,0,0,1,0,0, 0,0,0,0,1,0, 0,0,0,0,0,1),5,6, byrow=TRUE),
+						matrix(c(0,1,0,0,0,0, 0,0,1,0,0,0, 0,0,0,1,0,0, 0,0,0,0,1,0, 0,0,0,0,1,1),5,6, byrow=TRUE),
+						matrix(c(0,1,0,0,0,0, 0,0,1,0,0,0, 0,0,0,1,0,0, 0,0,0,0,1,0, 0,0,0,0,0,1),5,6, byrow=TRUE),
+						matrix(c(0,1,0,0,0,0, 0,0,1,0,0,0, 0,0,0,1,0,0, 0,0,0,0,1,0, 0,0,0,1,0,1),5,6, byrow=TRUE),
+						matrix(c(0,1,1,0,0,0, 0,0,-1,0,0,0, 0,0,0,1,0,0, 0,0,0,0,1,0, 0,0,0,0,0,1),5,6, byrow=TRUE),
+						matrix(c(0,1,1,0,0,0, 0,0,-1,0,0,0, 0,0,0,1,0,0, 0,0,0,0,1,0, 0,0,0,0,1,1),5,6, byrow=TRUE),
+						matrix(c(0,1,1,0,0,0, 0,0,-1,0,0,0, 0,0,0,1,0,0, 0,0,0,0,1,0, 0,0,0,0,0,1),5,6, byrow=TRUE),
+						matrix(c(0,1,1,0,0,0, 0,0,-1,0,0,0, 0,0,0,1,0,0, 0,0,0,0,1,0, 0,0,0,1,0,1),5,6, byrow=TRUE),
+						matrix(c(0,1,0,0,0,0, 0,0,1,0,0,0, 0,0,0,1,0,0, 0,0,0,0,1,0, 0,0,0,1,0,1),5,6, byrow=TRUE),
+						matrix(c(0,1,0,0,0,0, 0,0,1,0,0,0, 0,0,0,1,0,0, 0,0,0,0,1,0, 0,0,0,1,0,1),5,6, byrow=TRUE),
+						matrix(c(0,1,0,1,0,0,0, 0,0,1,0,0,0,0, 0,0,0,-1,0,0,0, 0,0,0,0,1,0,0, 0,0,0,0,0,1,0, 0,0,0,0,0,0,1),6,7, byrow=TRUE),
+						matrix(c(0,1,0,1,0,0,0, 0,0,1,0,0,0,0, 0,0,0,-1,0,0,0, 0,0,0,0,1,0,0, 0,0,0,0,0,1,0, 0,0,0,0,0,1,1),6,7, byrow=TRUE),
+						matrix(c(0,1,0,1,0,0,0, 0,0,1,0,0,0,0, 0,0,0,-1,0,0,0, 0,0,0,0,1,0,0, 0,0,0,0,0,1,0, 0,0,0,0,0,0,1),6,7, byrow=TRUE),
+						matrix(c(0,1,0,1,0,0,0, 0,0,1,0,0,0,0, 0,0,0,-1,0,0,0, 0,0,0,0,1,0,0, 0,0,0,0,0,1,0, 0,0,0,0,1,0,1),6,7, byrow=TRUE),
+						matrix(c(0,1,0,0,0,0,0, 0,0,1,0,0,0,0, 0,0,0,1,0,0,0, 0,0,0,0,1,0,0, 0,0,0,0,0,1,0, 0,0,0,0,1,0,1),6,7, byrow=TRUE),
+						matrix(c(0,1,1,0,0,0,0, 0,0,-1,0,0,0,0, 0,0,0,1,0,0,0, 0,0,0,0,1,0,0, 0,0,0,0,0,1,0, 0,0,0,0,1,0,1),6,7, byrow=TRUE))
+				}
+				else{
+					allConstr <- list(
+						NULL,
+						matrix(c(0,1),1,2),
+						NULL,
+						NULL,
+						matrix(c(0,1),1,2),
+						NULL,
+						matrix(c(0,1),1,2),
+						matrix(c(0,1,0),1,3, byrow=TRUE),
+						matrix(c(0,1,1,0,0,-1),2,3, byrow=TRUE),
+						matrix(c(0,1,0),1,3, byrow=TRUE),
+						matrix(c(0,1,0,0,0,1),2,3, byrow=TRUE),
+						matrix(c(0,1,0),1,3, byrow=TRUE),
+						matrix(c(0,1,0,0,0,1),2,3, byrow=TRUE),
+						NULL,
+						matrix(c(0,0,1),1,3, byrow=TRUE),
+						NULL,
+						matrix(c(0,0,1),1,3, byrow=TRUE),
+						matrix(c(0,0,1),1,3, byrow=TRUE),
+						NULL,
+						matrix(c(0,0,1),1,3, byrow=TRUE),
+						matrix(c(0,1,0),1,3, byrow=TRUE),
+						matrix(c(0,1,0,0,1,1),2,3, byrow=TRUE),
+						matrix(c(0,0,1),1,3, byrow=TRUE),
+						matrix(c(0,1,0,1, 0,0,0,-1),2,4, byrow=TRUE),
+						matrix(c(0,1,0,0),1,4, byrow=TRUE),
+						matrix(c(0,1,0,0, 0,0,0,1),2,4, byrow=TRUE),
+						matrix(c(0,1,0,0),1,4, byrow=TRUE),
+						matrix(c(0,1,0,0, 0,0,0,1),2,4, byrow=TRUE),
+						matrix(c(0,1,1,0, 0,0,-1,0),2,4, byrow=TRUE),
+						matrix(c(0,1,1,0, 0,0,-1,0, 0,0,0,1),3,4, byrow=TRUE),
+						matrix(c(0,1,1,0, 0,0,-1,0),2,4, byrow=TRUE),
+						matrix(c(0,1,1,0, 0,0,-1,0, 0,0,0,1),3,4, byrow=TRUE),
+						matrix(c(0,1,0,0, 0,0,0,1),2,4, byrow=TRUE),
+						matrix(c(0,1,0,0),1,4, byrow=TRUE),
+						matrix(c(0,1,0,0, 0,0,0,1),2,4, byrow=TRUE),
+						matrix(c(0,1,0,0, 0,0,1,0),2,4, byrow=TRUE),
+						matrix(c(0,1,0,0, 0,0,1,0, 0,0,1,1),3,4, byrow=TRUE),
+						matrix(c(0,1,0,0, 0,0,0,1),2,4, byrow=TRUE),
+						matrix(c(0,0,0,1),1,4, byrow=TRUE),
+						NULL,
+						matrix(c(0,0,0,1),1,4, byrow=TRUE),
+						matrix(c(0,0,1,0),1,4, byrow=TRUE),
+						matrix(c(0,0,1,0, 0,0,1,1),2,4, byrow=TRUE),
+						matrix(c(0,0,0,1),2,4, byrow=TRUE),
+						matrix(c(0,0,1,0),1,4, byrow=TRUE),
+						matrix(c(0,0,1,0, 0,0,1,1),2,4, byrow=TRUE),
+						matrix(c(0,0,0,1),1,4, byrow=TRUE),
+						matrix(c(0,1,0,0, 0,0,0,1),2,4, byrow=TRUE),
+						matrix(c(0,1,0,1,0, 0,0,0,-1,0),2,5, byrow=TRUE),
+						matrix(c(0,1,0,1,0, 0,0,0,-1,0, 0,0,0,0,1),3,5, byrow=TRUE),
+						matrix(c(0,1,0,1,0, 0,0,0,-1,0),2,5, byrow=TRUE),
+						matrix(c(0,1,0,1,0, 0,0,0,-1,0, 0,0,0,0,1),3,5, byrow=TRUE),
+						matrix(c(0,1,0,0,0, 0,0,0,0,1),2,5, byrow=TRUE),
+						matrix(c(0,1,0,0,0),1,5, byrow=TRUE),
+						matrix(c(0,1,0,0,0, 0,0,0,0,1),2,5, byrow=TRUE),
+						matrix(c(0,1,0,0,0, 0,0,0,1,0),2,5, byrow=TRUE),
+						matrix(c(0,1,0,0,0, 0,0,0,1,0, 0,0,0,1,1),3,5, byrow=TRUE),
+						matrix(c(0,1,0,0,0, 0,0,0,0,1),2,5, byrow=TRUE),
+						matrix(c(0,1,1,0,0, 0,0,-1,0,0, 0,0,0,0,1),3,5, byrow=TRUE),
+						matrix(c(0,1,1,0,0, 0,0,-1,0,0),2,5, byrow=TRUE),
+						matrix(c(0,1,1,0,0, 0,0,-1,0,0, 0,0,0,0,1),3,5, byrow=TRUE),
+						matrix(c(0,1,1,0,0, 0,0,-1,0,0, 0,0,0,1,0),3,5, byrow=TRUE),
+						matrix(c(0,1,1,0,0, 0,0,-1,0,0, 0,0,0,1,0, 0,0,0,1,1),4,5, byrow=TRUE),
+						matrix(c(0,1,1,0,0, 0,0,-1,0,0, 0,0,0,0,1),3,5, byrow=TRUE),
+						matrix(c(0,1,0,0,0, 0,0,0,1,0),2,5, byrow=TRUE),
+						matrix(c(0,1,0,0,0, 0,0,0,1,0, 0,0,0,1,1),3,5, byrow=TRUE),
+						matrix(c(0,1,0,0,0, 0,0,0,0,1),2,5, byrow=TRUE),
+						matrix(c(0,1,0,0,0, 0,0,1,0,0, 0,0,1,0,1),3,5, byrow=TRUE),
+						matrix(c(0,0,0,1,0),1,5, byrow=TRUE),
+						matrix(c(0,0,0,1,0, 0,0,0,1,1),2,5, byrow=TRUE),
+						matrix(c(0,0,0,0,1),1,5, byrow=TRUE),
+						matrix(c(0,0,1,0,0, 0,0,1,0,1),2,5, byrow=TRUE),
+						matrix(c(0,0,1,0,0, 0,0,1,0,1),2,5, byrow=TRUE),
+						matrix(c(0,1,0,1,0,0, 0,0,0,-1,0,0, 0,0,0,0,0,1),3,6, byrow=TRUE),
+						matrix(c(0,1,0,1,0,0, 0,0,0,-1,0,0),2,6, byrow=TRUE),
+						matrix(c(0,1,0,1,0,0, 0,0,0,-1,0,0, 0,0,0,0,0,1),3,6, byrow=TRUE),
+						matrix(c(0,1,0,1,0,0, 0,0,0,-1,0,0, 0,0,0,0,1,0),3,6, byrow=TRUE),
+						matrix(c(0,1,0,1,0,0, 0,0,0,-1,0,0, 0,0,0,0,1,0, 0,0,0,0,1,1),4,6, byrow=TRUE),
+						matrix(c(0,1,0,1,0,0, 0,0,0,-1,0,0, 0,0,0,0,0,1),3,6, byrow=TRUE),
+						matrix(c(0,1,0,0,0,0, 0,0,0,0,1,0),2,6, byrow=TRUE),
+						matrix(c(0,1,0,0,0,0, 0,0,0,0,1,0, 0,0,0,0,1,1),3,6, byrow=TRUE),
+						matrix(c(0,1,0,0,0,0, 0,0,0,0,0,1),2,6, byrow=TRUE),
+						matrix(c(0,1,0,0,0,0, 0,0,0,1,0,0, 0,0,0,1,0,1),3,6, byrow=TRUE),
+						matrix(c(0,1,1,0,0,0, 0,0,-1,0,0,0, 0,0,0,0,1,0),3,6, byrow=TRUE),
+						matrix(c(0,1,1,0,0,0, 0,0,-1,0,0,0, 0,0,0,0,1,0, 0,0,0,0,1,1),4,6, byrow=TRUE),
+						matrix(c(0,1,1,0,0,0, 0,0,-1,0,0,0, 0,0,0,0,0,1),3,6, byrow=TRUE),
+						matrix(c(0,1,1,0,0,0, 0,0,-1,0,0,0, 0,0,0,1,0,0, 0,0,0,1,0,1),4,6, byrow=TRUE),
+						matrix(c(0,1,0,0,0,0, 0,0,0,1,0,0, 0,0,0,1,0,1),3,6, byrow=TRUE),
+						matrix(c(0,0,0,1,0,0, 0,0,0,1,0,1),2,6, byrow=TRUE),
+						matrix(c(0,1,0,1,0,0,0, 0,0,0,-1,0,0,0, 0,0,0,0,0,1,0),3,7, byrow=TRUE),
+						matrix(c(0,1,0,1,0,0,0, 0,0,0,-1,0,0,0, 0,0,0,0,0,1,0, 0,0,0,0,0,1,1),4,7, byrow=TRUE),
+						matrix(c(0,1,0,1,0,0,0, 0,0,0,-1,0,0,0, 0,0,0,0,0,0,1),3,7, byrow=TRUE),
+						matrix(c(0,1,0,1,0,0,0, 0,0,0,-1,0,0,0, 0,0,0,0,1,0,0, 0,0,0,0,1,0,1),4,7, byrow=TRUE),
+						matrix(c(0,1,0,0,0,0,0, 0,0,0,0,1,0,0, 0,0,0,0,1,0,1),3,7, byrow=TRUE),
+						matrix(c(0,1,1,0,0,0,0, 0,0,-1,0,0,0,0, 0,0,0,0,1,0,0, 0,0,0,0,1,0,1),4,7, byrow=TRUE))
+				}
+			}
+			else{
+				# 11 models
+				allLabs <- list(
+					"theta0.0",
+					c("theta0.0", "theta0.1"),
+					c("theta0.0", "theta2.1"),
+					c("theta0.0", "theta3.1"),
+					c("theta0.0", "theta0.1", "theta1.1"),
+					c("theta0.0", "theta0.1", "theta2.1"),
+					c("theta0.0", "theta0.1", "theta3.1"),
+					c("theta0.0", "theta2.1", "theta3.1"),
+					c("theta0.0", "theta0.1", "theta1.1", "theta2.1"),
+					c("theta0.0", "theta0.1", "theta1.1", "theta3.1"),
+					c("theta0.0", "theta0.1", "theta2.1", "theta3.1"))
+				allConstr <- list(
+						NULL,
+						matrix(c(0,1),1,2),
+						matrix(c(0,1),1,2),
+						matrix(c(0,1),1,2),
+						matrix(c(0,1,1,0,0,-1),2,3, byrow=TRUE),
+						matrix(c(0,1,0,0,0,1),2,3, byrow=TRUE),
+						matrix(c(0,1,0,0,0,1),2,3, byrow=TRUE),
+						matrix(c(0,1,0,0,1,1),2,3, byrow=TRUE),
+						matrix(c(0,1,1,0, 0,0,-1,0, 0,0,0,1),3,4, byrow=TRUE),
+						matrix(c(0,1,1,0, 0,0,-1,0, 0,0,0,1),3,4, byrow=TRUE),
+						matrix(c(0,1,0,0, 0,0,1,0, 0,0,1,1),3,4, byrow=TRUE))
+				allTypes <- c("Intercept", "Linear", rep("Piecewise linear",9))
+			}
+		}
+		if(constr.slopes==1){
+			constrmat <- t(object@QP$Amat)
+			if(!continuous){
+				mylab <- c("theta0.0", "theta0.1", "theta1.0", "theta1.1", "theta2.0", "theta2.1", "theta3.0", "theta3.1")
+				searchAll <- list(NULL,t(combn(1:7,1))+1,t(combn(1:7,2))+1,t(combn(1:7,3))+1,t(combn(1:7,4))+1,t(combn(1:7,5))+1,t(combn(1:7,6))+1)
+
+				nbmodels <- sum(unlist(lapply(searchAll,nrow)))+1
+				ct <- 1
+				allLabs <- allConstr <- rep(list(NULL),nbmodels)
+				for(i in 1:length(searchAll)){
+					if(i==1){
+						allLabs[[ct]] <- "theta0.0"
+						ct <- ct + 1
+					}
+					else{
+						for(j in 1:nrow(searchAll[[i]])){
+							tokeep <- searchAll[[i]][j,]
+							allLabs[[ct]] <- mylab[c(1,tokeep)]
+							if(object@call.arg$constr){
+								if(length(tokeep)==1) allConstr[[ct]] <- matrix(c(0,1),1,2)
+								else{
+									tp <- constrmat[,c(1,tokeep)]
+									tp <- tp[rowSums(tp)!=0,]
+									if(!is.matrix(tp)) tp <- matrix(tp, 1:length(tp), byrow=TRUE)
+									else{
+										tp <- tp[!duplicated(tp),]
+										if(!is.matrix(tp)) tp <- matrix(tp, 1:length(tp), byrow=TRUE)
+										allConstr[[ct]] <- tp
+									}
+								}
+							}
+							ct <- ct + 1
+						}
+					}
+				}
+				allTypes <- c("Intercept", "Linear",
+						"Piecewise level",
+						"Piecewise linear",
+						"Piecewise level",
+						"Piecewise linear",
+						"Piecewise level",
+						rep("Piecewise linear",8),
+						"Piecewise level",
+						"Piecewise linear",
+						"Piecewise level",
+						rep("Piecewise linear",6),
+						"Piecewise level",
+						rep("Piecewise linear",24),
+						"Piecewise level",
+						rep("Piecewise linear",77))
+			}
+			else{
+					allLabs <- list(
+							"theta0.0",
+							c("theta0.0", "theta0.1"),
+							c("theta0.0", "theta1.1"),
+							c("theta0.0", "theta2.1"),
+							c("theta0.0", "theta3.1"),
+							c("theta0.0", "theta0.1", "theta1.1"),
+							c("theta0.0", "theta0.1", "theta2.1"),
+							c("theta0.0", "theta0.1", "theta3.1"),
+							c("theta0.0", "theta1.1", "theta2.1"),
+							c("theta0.0", "theta1.1", "theta3.1"),
+							c("theta0.0", "theta2.1", "theta3.1"),
+							c("theta0.0", "theta0.1", "theta1.1", "theta2.1"),
+							c("theta0.0", "theta0.1", "theta1.1", "theta3.1"),
+							c("theta0.0", "theta0.1", "theta2.1", "theta3.1"),
+							c("theta0.0", "theta1.1", "theta2.1", "theta3.1"))
+					if(object@call.arg$constr) allConstr <- list(
+							NULL,
+							matrix(c(0,1),1,2, byrow=TRUE),
+							matrix(c(0,1),1,2, byrow=TRUE),
+							matrix(c(0,1),1,2, byrow=TRUE),
+							matrix(c(0,1),1,2, byrow=TRUE),
+							matrix(c(0,1,0, 0,0,1),2,3, byrow=TRUE),
+							matrix(c(0,1,0, 0,1,1),2,3, byrow=TRUE),
+							matrix(c(0,1,0, 0,0,1),2,3, byrow=TRUE),
+							matrix(c(0,1,0, 0,0,1),2,3, byrow=TRUE),
+							matrix(c(0,1,0, 0,0,1),2,3, byrow=TRUE),
+							matrix(c(0,1,0, 0,0,1),2,3, byrow=TRUE),
+							matrix(c(0,1,0,0, 0,0,1,0, 0,1,0,1),3,4, byrow=TRUE),
+							matrix(c(0,1,0,0, 0,0,1,0, 0,0,0,1),3,4, byrow=TRUE),
+							matrix(c(0,1,0,0, 0,1,1,0, 0,0,0,1),3,4, byrow=TRUE),
+							matrix(c(0,1,0,0, 0,0,1,0, 0,0,0,1),3,4, byrow=TRUE))
+					allTypes <- c("Intercept", "Linear", rep("Piecewise linear", 13))
+			}
+		}
+	}
+	f <- function(labs, x){
+		if(length(labs)==1){
+			m <- matrix(1,nrow(x),1)
+			colnames(m) <- "theta0.0"
+		}
+		else{
+			m <- x[,colnames(x)%in%labs]
+		}
+		m
+	}
+	allMats <- lapply(allLabs, f, x=object@X)
+	return(list(allLabs=allLabs, allConstr=allConstr, allTypes=allTypes, allMats=allMats))
+}
+
+
diff --git a/R/plrs.test.r b/R/plrs.test.r
new file mode 100644
index 0000000..55012d8
--- /dev/null
+++ b/R/plrs.test.r
@@ -0,0 +1,103 @@
+# Likelihood ratio test for (constrained) plrs models
+
+# Unconstrained:
+# H0: Intercept model
+# H1: Full model
+
+# Constrained:
+# H0: All constraints are actives (=)
+# H1: At least one constraint is strict (>)
+
+# Authors: Gwenael G.R. Leday based on Ulrike Gromping code from package ic.infer
+
+plrs.test <- function(object, alpha=0.05){
+
+	if(class(object)!="plrs") stop("An object of class \"plrs\" is required")
+	if(object@selected){
+		obj <- plrs(expr=object@mdata$mexpr, cghseg=object@mdata$mcghseg, cghcall=object@mdata$mcghcall,
+				probloss = object@mdata$mprobloss,
+				probnorm = object@mdata$mprobnorm,
+				probgain = object@mdata$mprobgain,
+				probamp = object@mdata$mprobamp,
+				knots = NULL,
+				continuous = object@call.arg$continuous,
+				constr = object@call.arg$constr,
+				constr.slopes = object@call.arg$constr.slopes,
+				constr.intercepts = object@call.arg$constr.intercepts)
+	}else{
+		obj <- object
+	}
+	if(length(coef(obj))!=1){
+		if(obj@call.arg$constr){ # Test whether all inequalities are equalities
+			# Aknowledgement: Taken from package ic.infer
+			if(obj@QP$meq == nrow(t(obj@QP$Amat))){
+				stop("Test not possible: only equality constraints")
+			}
+
+			library(mvtnorm)
+
+			# Unconstrained model
+			unconstr <- .plrs.fit(X=obj@X, matconstr=NULL, mytype=obj@type, modelFull=obj, callArg=obj@call.arg)
+			cov.mat <- unconstr@QP$vcov/unconstr@QP$S2
+			Dmat <- obj@QP$Dmat
+			dvec <- as.vector(obj@QP$dvec)
+			uiw <- t(obj@QP$Amat)
+			bvec <- obj@QP$bvec
+			meq <- obj@QP$meq
+
+			# If intercepts are not constrained set them to 0
+			if(!obj@call.arg$constr.intercepts & !obj@call.arg$continuous){
+			     para <- c(3,5,7)[c(3,5,7) < ncol(model.matrix(obj))]
+			     mat <- matrix(0,length(para),ncol(model.matrix(obj)))
+			     mat[1:length(para),para] <- 1
+			     uiw <- rbind(mat,uiw)
+			     bvec <- c(bvec,rep(0,nrow(mat)))
+			     meq <- meq + length(para)
+			}
+
+			# Weights
+			if(meq==0){
+				wt.bar <- ic.weights(crossprod(t(crossprod(t(uiw),cov.mat)),t(uiw)))
+			}
+			else{
+				wt.bar <- ic.weights(solve(solve(crossprod(t(crossprod(t(uiw),cov.mat)),t(uiw)))[-(1:obj@QP$meq), -(1:obj@QP$meq)]))
+			}
+
+			# Stat, quantile and pval
+			b.eq  <- solve.QP(Dmat, dvec, t(uiw), bvec, meq = nrow(uiw))$sol
+			b.eq <- round(b.eq,8)
+			names(b.eq) <- names(coef(obj))
+			stat <- as.vector(crossprod(coef(obj)-b.eq, solve(cov.mat, coef(obj)-b.eq)))
+			stat <- stat/(unconstr@QP$df.error*unconstr@QP$S2+stat)
+			df.bar <- (nrow(uiw)-meq):0
+			pval <- 1 - pbetabar(stat, df1=df.bar/2, df2=unconstr@QP$df.error/2, wt=wt.bar)
+			qbetabar <- .quantileBetabar(df.bar=df.bar, df.error=unconstr@QP$df.error, wt.bar=wt.bar, alpha=alpha)
+			test <- list("stat" = stat, "pvalue" = pval, "wt.bar" = wt.bar, "df.bar" = df.bar, "unconstr" = unconstr, "qbetabar" = qbetabar, "alpha" = alpha)
+		}
+		else{ # intercept vs given model
+			myX <- matrix(1, nrow(obj@X), 1, dimnames=list(NULL,"theta0.0"))
+			modelNull <- .plrs.fit(X=myX, matconstr=NULL, mytype="Intercept", modelFull=obj)
+
+			rss0 <- sum(residuals(modelNull)^2)
+			rss1 <- sum(residuals(obj)^2)
+
+			# Degree of freedom
+			N <- length(obj@data$expr)
+			df0  <- N - 1
+			df1  <- N - ncol(model.matrix(obj))
+			dfnum <- df0 - df1
+			dfden <- df1
+
+			# P-value
+			stat  <- ((rss0 - rss1)/dfnum)/(rss1/dfden)
+			pval <- 1 - pf(stat, dfnum, dfden)
+			test <- list("stat" = stat, "pvalue" = pval)
+		}
+		object@test <- test
+		return(object)
+	}
+	else{
+		return(object)
+	}
+}
+
diff --git a/R/predict.plrs.r b/R/predict.plrs.r
new file mode 100644
index 0000000..5d4be12
--- /dev/null
+++ b/R/predict.plrs.r
@@ -0,0 +1,18 @@
+# Make predictions from a given "plrs" model
+
+# Author: Gwenael G.R. Leday
+
+predict.plrs <- function(object, newcghseg, ...){
+
+	if(class(object)!="plrs") stop("An object of class \"plrs\" is required")
+	if(!is.vector(newcghseg)) stop("Object \"newcghseg\" has to be a vector")
+
+	X <- .Bmat(newcghseg, knots = object@data$knots,
+		    continuous = object@call.arg$continuous,
+		    general.intercept = T)
+
+	if(object@selected) X <- X[,colnames(X)%in%names(coef(object))]
+
+	fit.val <- as.vector(crossprod(t(X), coef(object)))
+	return(fit.val)
+}
diff --git a/R/quantileBetabar.r b/R/quantileBetabar.r
new file mode 100644
index 0000000..67ff756
--- /dev/null
+++ b/R/quantileBetabar.r
@@ -0,0 +1,25 @@
+# Internal function
+# Determine quantile of mixture of beta distributions
+
+# Author: Gwenael G.R. Leday
+
+.quantileBetabar <- function(df.bar, df.error, wt.bar, alpha){
+	statbis <- 0.5
+	stp <- 0.499999
+	vals <- rep(NA, 100)
+	bool <- T
+	ct <- 0
+	while(bool){
+		ct <- ct + 1
+		vals[ct] <- statbis
+		pvalbis <- 1 - pbetabar(statbis, df1=df.bar/2, df2=df.error/2, wt=wt.bar)
+		if(pvalbis>alpha){
+			statbis <- statbis + stp
+		}else{
+			statbis <- statbis - stp
+		}
+		stp <- stp / 2
+		if(ct>1) if(abs(pvalbis-alpha)<=.Machine$double.eps | abs(vals[ct-1]-vals[ct])<=1e-8) bool <- F
+	}
+	vals[ct]
+}
\ No newline at end of file
diff --git a/R/respCI.r b/R/respCI.r
new file mode 100644
index 0000000..d74a42a
--- /dev/null
+++ b/R/respCI.r
@@ -0,0 +1,88 @@
+# Internal function
+# Confidence bands for a given x, optimization with csdp().
+
+# Author: Gwenael G.R. Leday
+
+.respCI <- function(newx, object, alpha){
+
+	# Test including intercept
+	cov.mat <- object@test$unconstr@QP$vcov/object@test$unconstr@QP$S2
+	x <- rep(0, length(coef(object)))
+	newstat <- as.numeric(t(coef(object)-x) %*% solve(cov.mat, coef(object)-x))
+	newstat <- newstat/(object@test$unconstr@QP$df.error*object@test$unconstr@QP$S2+newstat)
+	qbetabar <- .quantileBetabar(df.bar=object@test$df.bar+1, df.error=object@test$unconstr@QP$df.error, wt.bar=object@test$wt.bar, alpha=alpha)
+	lambda <- (object@test$unconstr@QP$df.error*object@test$unconstr@QP$S2)*(qbetabar/(1-qbetabar))
+
+	R <- solve(cov.mat)
+	M <- chol(R)
+	ML <- coef(object)
+	bb <- -M%*%ML
+	constr <- t(object@QP$Amat)
+	npar <- length(ML)
+
+	# myC
+	myC <- list(
+		-1*rbind(cbind(diag(npar), bb), cbind(t(bb), lambda)),
+		matrix(0, npar-1, npar-1)
+		)
+
+	# myA
+	A11 <- rbind(cbind(matrix(0,npar,npar), matrix(M[,1],npar,1)), cbind(matrix(M[,1],1,npar),0))
+	A12 <- rbind(cbind(matrix(0,npar,npar), matrix(M[,2],npar,1)), cbind(matrix(M[,2],1,npar),0))
+	if(npar>2) A13 <- rbind(cbind(matrix(0,npar,npar), matrix(M[,3],npar,1)), cbind(matrix(M[,3],1,npar),0))
+	if(npar>2) A14 <- rbind(cbind(matrix(0,npar,npar), matrix(M[,4],npar,1)), cbind(matrix(M[,4],1,npar),0))
+	if(npar>4) A15 <- rbind(cbind(matrix(0,npar,npar), matrix(M[,5],npar,1)), cbind(matrix(M[,5],1,npar),0))
+	if(npar>4) A16 <- rbind(cbind(matrix(0,npar,npar), matrix(M[,6],npar,1)), cbind(matrix(M[,6],1,npar),0))
+	if(npar>6) A17 <- rbind(cbind(matrix(0,npar,npar), matrix(M[,7],npar,1)), cbind(matrix(M[,7],1,npar),0))
+	if(npar>6) A18 <- rbind(cbind(matrix(0,npar,npar), matrix(M[,8],npar,1)), cbind(matrix(M[,8],1,npar),0))
+
+	if(npar==2){
+		A21 <- constr[,1]
+		A22 <- constr[,2]
+	}
+	else{
+		A21 <- diag(constr[,1])
+		A22 <- diag(constr[,2])
+	}
+	if(npar>2) A23 <- diag(constr[,3])
+	if(npar>2) A24 <- diag(constr[,4])
+	if(npar>4) A25 <- diag(constr[,5])
+	if(npar>4) A26 <- diag(constr[,6])
+	if(npar>6) A27 <- diag(constr[,7])
+	if(npar>6) A28 <- diag(constr[,8])
+
+	myA <- NULL
+	for(i in 1:npar){
+		myA <- c(myA, list(list(get(paste("A1",i,sep="")), get(paste("A2",i,sep="")))))
+	}
+
+	#myK
+	if(npar!=2){
+		myK <- list(type=c("s","s"), size=c(npar+1,npar-1))
+	}
+	else{
+		myK <- list(type=c("s","l"), size=c(npar+1,npar-1))
+	}
+
+	#myb
+	myb <- as.vector(.Bmat(newx, knots=knots(object), object@call.arg$continuous, TRUE))
+
+	# Semidefinite prog for inf
+	res <- csdp(myC,myA,myb,myK)
+	parinf <- res$y
+	inf <- res$dobj
+
+	#myb
+	myb <- -as.vector(.Bmat(newx, knots=knots(object), object@call.arg$continuous, TRUE))
+
+	# Semidefinite prog for sup
+	res <- csdp(myC,myA,myb,myK)
+	parsup <- res$y
+	sup <- -res$dobj
+
+	out <- c(inf,sup)
+	names(out) <- c("inf","sup")
+
+	return(out)
+}
+
diff --git a/data/neveCN17.RData b/data/neveCN17.RData
new file mode 100644
index 0000000..4a71b03
Binary files /dev/null and b/data/neveCN17.RData differ
diff --git a/data/neveGE17.RData b/data/neveGE17.RData
new file mode 100644
index 0000000..ae23410
Binary files /dev/null and b/data/neveGE17.RData differ
diff --git a/inst/doc/mybib.bib b/inst/doc/mybib.bib
new file mode 100644
index 0000000..0bb980b
--- /dev/null
+++ b/inst/doc/mybib.bib
@@ -0,0 +1,61 @@
+
+@article{gromping2010,
+    title = {Inference with Linear Equality and Inequality Constraints Using {R}: The Package ic.infer},
+    author = {Gr\"{o}mping, Ulrike},
+    year = {2010},
+    journal = {J Stat Softw},
+    volume = {33},
+    number = {i10}
+}
+
+@article{leday2012,
+    author = {Leday, Gwena\"el G. R. and van der Vaart, Aad W. and van Wieringen, Wessel N. and van de Wiel, Mark A.},
+    title = {Modeling association between {DNA} copy number and gene expression with constrained piecewise linear regression splines},
+    journal = {Ann Appl Stat},
+    year = {2012},
+    note = {Accepted for publication}
+}
+
+@article{neve2006,
+    author = {Neve, Richard M. and Chin, Koei and Fridlyand, Jane and Yeh, Jennifer and Baehner, Frederick L. and Fevr, Tea and Clark, Laura and Bayani, Nora and Coppe, Jean-Philippe P. and Tong, Frances and Speed, Terry and Spellman, Paul T. and DeVries, Sandy and Lapuk, Anna and Wang, Nick J. and Kuo, Wen-Lin L. and Stilwell, Jackie L. and Pinkel, Daniel and Albertson, Donna G. and Waldman, Frederic M. and McCormick, Frank and Dickson, Robert B. and Johnson, Michael D. and Lippman, Marc and Ethier, Stephen and Gazdar, Adi and Gray, Joe W.},
+    title = {A collection of breast cancer cell lines for the study of functionally distinct cancer subtypes.},
+    journal = {Cancer cell},
+    volume = {10},
+    year = {2006},
+    month = {December},
+    number = {6},
+    pages = {515--527},
+    issn = {1535-6108}
+}
+
+@article{olshen2004,
+    author = {Olshen, A. B. and Venkatraman, E. S. and Lucito, R. and Wigler, M.},
+    journal = {Biostatistics},
+    month = {October},
+    number = {4},
+    pages = {557--572},
+    title = {Circular binary segmentation for the analysis of array-based {DNA} copy number data.},
+    volume = {5},
+    year = {2004}
+}
+
+@article{vandeWiel2007,
+    author = {van de Wiel, Mark A. and Kim, Kyung I. and Vosse, Sjoerd J. and van Wieringen, Wessel N. and Wilting, Saskia M. and Ylstra, Bauke},
+    journal = {Bioinformatics},
+    month = {April},
+    number = {7},
+    pages = {892--894},
+    title = {{CGHcall}: calling aberrations for array {CGH} tumor profiles},
+    volume = {23},
+    year = {2007}
+}
+
+@article{vanWieringen2012,
+	title = {Matching of array {CGH} and gene expression microarray features for the purpose of integrative genomic analyses},
+	author = {van Wieringen, Wessel N. and Unger, Kristian and Leday, Gwenael and Krijgsman, Oscar and de Menezes, Renee and Ylstra, Bauke and van de Wiel, Mark},
+	journal = {BMC Bioinformatics},
+	year = {2012},
+	volume = {13},
+	number = {1},
+	pages = {80}
+}
diff --git a/inst/doc/plrs_vignette.Rnw b/inst/doc/plrs_vignette.Rnw
new file mode 100644
index 0000000..8369adb
--- /dev/null
+++ b/inst/doc/plrs_vignette.Rnw
@@ -0,0 +1,404 @@
+\documentclass[11pt]{article}
+
+%\VignetteIndexEntry{plrs}
+%\VignetteDepends{}
+%\VignetteKeywords{Piecewise Linear Regression Splines for the association between DNA copy number and gene expression.}
+%\VignettePackage{plrs}
+
+\usepackage{geometry}
+\geometry{left=2.5cm,right=2.5cm,top=2.5cm, bottom=2.5cm}
+
+\newcommand{\plrs}{\texttt{PLRS}}
+\SweaveOpts{keep.source=TRUE}
+
+
+\title{\plrs \\ Piecewise Linear Regression Splines for \\
+the association between DNA copy number and \\ gene expression}
+\author{Gwena\"el G.R. Leday\\
+\\
+\texttt{g.g.r.leday@vu.nl}\\
+\\
+\normalsize{Department of Mathematics, VU University}\\
+\normalsize{De Boelelaan 1081a, 1081 HV Amsterdam, The Netherlands}\\}
+\date{}
+
+
+\begin{document}
+
+\maketitle
+
+\null
+\null
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\section{Introduction}
+\null
+\indent The \plrs~package implements the methodology described by \cite{leday2012} for the joint
+analysis of DNA copy number and mRNA expression.
+The framework employs piecewise linear regression splines (PLRS), a broad class of interpretable
+models, to model \textsl{cis}-relationships between the two molecular levels. \\
+In the present vignette, we provide guidance for:
+\begin{enumerate}
+	\item The analysis of a single DNA-mRNA relationship
+		\begin{itemize}
+			\item model specification and fitting,
+			\item selection of an appropriate model,
+			\item testing the strength of the association and
+			\item obtaining uniform confidence bands.
+		\end{itemize}
+	\item The analysis of multiple DNA-mRNA relationships
+\end{enumerate}
+\null
+We first introduce the breast cancer data used to illustrate this vignette. \\
+
+\newpage
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\section{Breast cancer data}
+\null
+We use the breast cancer data of \cite{neve2006} that are available from the Bioconductor package
+``Neve2006''. aCGH and gene expression for 50 samples (cell lines) were profiled with an OncoBAC
+and Affymetrix HG-U133A arrays. We preprocessed and matched data as follows.
+Data were segmented with the CBS algorithm of \cite{olshen2004} and discretized with CGHcall
+\cite{vandeWiel2007} (also available from Bioconductor). Matching of features was done using
+the Bioconductor package \texttt{sigaR} \cite{vanWieringen2012} and the \textit{distanceAny}
+method with a window of 10,000 bp. The present package includes preprocessed data for
+chromosome 17 only.\\
+
+\noindent Data are loaded in the following way:
+
+<<>>=
+library(plrs)
+data(neveCN17)
+data(neveGE17)
+@
+\null
+
+This loads to the current environment an ``ExpressionSet'' object \texttt{neveGE17} and
+a ``cghCall'' object \texttt{neveCN17}, which contain \textbf{preprocessed} expression and copy number
+data with \textbf{matched} features:
+
+<<>>=
+neveGE17
+neveCN17
+@
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\section{Analysis of a single DNA-mRNA relationship}
+\null
+In this section, we focus on a single \textsl{cis}-effect and show how to specify
+and fit a model. Procedures for selecting an appropriate model, testing
+the association and obtaining confidence intervals for the mean response are presented.\\
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\subsection{Obtain data for a single gene}
+\null
+We choose gene PITPNA for illustration:
+
+<<>>=
+# Index of gene PITPNA
+idx <- which(fData(neveGE17)$Symbol=="PITPNA")[1]
+@
+<<>>=
+# Obtain vectors of gene expression (normalized) and
+# copy number (segmented and called) data
+rna <- exprs(neveGE17)[idx,]
+cghseg <- segmented(neveCN17)[idx,]
+cghcall <- calls(neveCN17)[idx,]
+@
+<<>>=
+# Obtain vectors of posterior membership probabilities
+probloss <- probloss(neveCN17)[idx,]
+probnorm <- probnorm(neveCN17)[idx,]
+probgain <- probgain(neveCN17)[idx,]
+probamp <- probamp(neveCN17)[idx,]
+@
+\null
+Posterior probabilities are used for determining knots (or change points) of the model.
+Their calculation is explained in \cite{leday2012}. \\
+
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\subsection{Number of observations per state}
+\null
+
+For a given gene, PLRS accommodates \textsl{differential} DNA-mRNA relationships across
+the different types of genomic aberrations or \textsl{states}. Presently, we (only) distinguish
+four of them: \textbf{loss} (coded as -1), \textbf{normal} (0), \textbf{gain} (1) and \textbf{amplification} (2).
+To ensure that the model is identifiable, the sample size for each state needs not to be too
+small. A minimum number of three samples is required for estimation, however any higher number
+may be chosen/preferred.
+<<>>=
+# Check: how many observations per state?
+table(cghcall)
+@
+Here, it is possible to fit a 3-state model. If the minimum number of observations is not fulfilled for
+a given state, there are two possibilities: discard data corresponding to the given state or merging it
+to an adjacent one. The function \texttt{modify.conf()} accommodates these two options (see below).
+With this function one can actually control the minimum number of samples per state.
+<<>>=
+# Set the minimum to 4 observations per state
+cghcall2 <- modify.conf(cghcall, min.obs=4, discard=FALSE)
+table(cghcall2)
+
+# Set the minimum to 4 observations per state
+cghcall2 <- modify.conf(cghcall, min.obs=4, discard=TRUE)
+table(cghcall2)
+@
+\null
+In practice, the user does not have to call directly \texttt{modify.conf()} as it is implemented
+internally in function \texttt{plrs()}, which is used to fit a model.
+However, one needs to be aware that the minimum number of observations per state (specified via
+arguments \texttt{min.obs} and \texttt{discard.obs} in \texttt{plrs()}) has a strong
+influence on the resulting model.\\
+\null
+\null
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\subsection{Model fitting}
+\null
+Function \texttt{plrs()} allows the fitting of a PLRS model. Different types of models may be
+specified conveniently. For instance, one may choose to fit a \textit{continuous} PLRS (i.e.
+with no state-specific intercepts or discontinuities) or to change the type of restrictions on
+parameters or simply leave them out.
+Although we recommend users to use default argument values, we below describe how to specify
+alternative types of models.\\
+
+\noindent The followings arguments of function \texttt{plrs()} offer flexibility for modeling:
+\begin{itemize}
+	\item \texttt{continuous = TRUE} or \texttt{FALSE} (default)\\
+	Specify whether state-specific intercepts should be included.
+	\item \texttt{constr = TRUE} (default) or \texttt{FALSE}\\
+	Specify whether inequality constraints on parameter should be applied or not.
+	\item \texttt{constr.slopes = 1} or \texttt{2} (default)\\
+	Specify the type of constraints on slopes; options are:
+	\begin{enumerate}
+		\item \texttt{constr.slopes = 1} forces all slopes to be non-negative.
+		\item \texttt{constr.slopes = 2} forces the slope of the "normal" state (coded 0) to
+		be non-negative while all others are forced to be at least equal to the latter.
+	\end{enumerate}
+	\item \texttt{constr.intercepts = TRUE} (default) or \texttt{FALSE}\\
+	Specify whether state-specific intercepts are to be non-negative.
+\end{itemize}
+\noindent Note that when \texttt{constr = FALSE}, \texttt{constr.slopes} and
+\texttt{constr.intercepts} have no effect.\\
+\\
+\\
+With default settings (recommended):
+\begin{center}
+\setkeys{Gin}{width=0.5\textwidth}
+<<fig=TRUE>>=
+# Fit a model
+model <- plrs(rna, cghseg, cghcall, probloss, probnorm, probgain, probamp)
+model
+plot(model)
+@
+\end{center}
+\null
+The \texttt{plrs()} function returns an S4 object of class ``plrs''. Various generic functions
+have been defined on this class:
+\begin{itemize}
+	\item Functions \texttt{print()} and \texttt{summary()} display information about the model
+	(e.g. estimated coefficients, type of constraints, etc...).
+	\item Function \texttt{plot()} displays the fit of the PLRS model along with data.
+	Segmented copy number is on x-axis while (normalized) gene expression is on y-axis;
+	colors indicate the different aberration states, namely ``loss'' (red), ``normal''
+	(blue) and ``gain'' (green); the black line gives the fit of the PLRS model.
+	Colors and symbols may be changed via the appropriate arguments of function \texttt{plrs()}.
+	Note that the argument \texttt{lin}, if set to \texttt{TRUE}, will additionally display a
+	dashed (default) line giving the fit of the simple linear model.
+	\item Other useful functions include \texttt{coef()}, \texttt{fitted()}, \texttt{residuals()},
+	\texttt{model.matrix()}, \texttt{predict()} and \texttt{knots()}. These are standard functions.
+	Information about these can be found in associated help files.
+\end{itemize}
+\null
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\subsection{Select an appropriate model}
+\null
+Model selection is carried out by \texttt{plrs.select()}, which takes has input an object of
+class ``plrs''. Possible model selection criteria are \texttt{OSAIC}, \texttt{AIC}, \texttt{AICC}
+and \texttt{BIC}. When the model is constrained \texttt{OSAIC} is the only appropriate one.
+If the model has no restrictions on parameters, \texttt{OSAIC} and \texttt{AIC} are equivalent.
+See help file and associated references for more information.
+
+\begin{center}
+\setkeys{Gin}{width=0.5\textwidth}
+<<fig=TRUE>>=
+# Model selection
+modelSelection <- plrs.select(model)
+summary(modelSelection)
+
+# Plot selected model
+plot(modelSelection)
+@
+\end{center}
+\null
+The function \texttt{plrs.select()} returns an S4 object of class ``plrs.select'' (here
+\texttt{modelSelection}), which contains information on model selection. Slot \texttt{model}
+contains the selected model as a ``plrs'' object:
+<<>>=
+selectedModel <- modelSelection@model
+selectedModel
+@
+\null
+
+Although \texttt{model} and \texttt{selectedModel} are both objects of class ``plrs'',
+generic functions defined on this class make implicitly the distinction between
+objects resulting from functions \texttt{plrs()} and \texttt{plrs.select()}
+(see above; ``\texttt{Selected spline coefficients}'').\\
+
+$\Rightarrow$ It is important to note that, for correct inference, subsequent test and confidence
+intervals must be computed on the \textbf{full} model rather than the \textbf{selected} one.
+Therefore, we forced functions \texttt{plrs.test()} and \texttt{plrs.cb()} (used hereafter)
+to operate on the full model, regardless of the input model. This implies that inference
+results are the same with either objects (see below).
+If one wishes to obtain results for the selected model, set \texttt{selectedModel@selected} to \texttt{FALSE}
+and then apply the aforementioned functions.
+\null
+
+%This provides correct inference  the user to continue the analysis with either \texttt{model}
+%or \texttt{selectedModel}. 
+%This distinction is important because subsequent statistical inferences (tests, confidence
+%intervals, ...) should be based on the full model rather than the selected one. Doing
+%otherwise would naively suggests that the form of the relationship is known \textsl{a priori}.
+%Hence, one should continue the analysis with the object \texttt{model}, rather than
+%\texttt{selectedModel}. However, for practical convenience the \texttt{PLRS} package allows
+%to continue the analysis with \texttt{selectedModel} and keep information from the selection
+%procedure.
+%To do so, generic functions associated with the class of objects ``plrs'' make implicitly
+%the distinction between full and selected models, i.e. between objects resulting from
+%\texttt{plrs()} and those from \texttt{plrs.select()}. This can be seen when objects are printed
+%for example (see above; ``\texttt{Selected spline coefficients}''). Thereby, function
+%\texttt{plrs.test}, which will be used for testing model parameters, always consider the full model
+%and tests the association in its generality even though a selected model has been provided.
+%Similarly, function \texttt{plrs.cb} determine confidence intervals on the mean response based on
+%the full model.
+
+\newpage
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\subsection{Testing the strength of the association}
+\null
+
+The function \texttt{plrs.test()} implements a likelihood ratio test to evaluate the effect
+of DNA copy number on expression. It tests the hypothesis that all parameters (except the
+overall intercept) of the PLRS model equal 0. See \cite{gromping2010, leday2012} and associated
+references for more details on the test.
+The function \texttt{plrs.test()} takes as input an object class ``plrs'' and outputs an object
+from the same class. Testing results are contained in slot \texttt{test}.
+
+<<>>=
+# Testing the full model with
+model <- plrs.test(model, alpha=0.05)
+model
+
+# or with
+selectedModel <- plrs.test(selectedModel, alpha=0.05)
+selectedModel
+
+# Testing the selected model
+selectedModel2 <- selectedModel
+selectedModel2@selected <- FALSE
+plrs.test(selectedModel2, alpha=0.05)
+@
+\null
+
+The object \texttt{selectedModel} now contains information on testing, which are
+consequently displayed when printing.\\
+
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\subsection{Uniform confidence bands}
+\null
+Confidence intervals for the spline are computed with function \texttt{plrs.cb}. Again, the
+function requires an object of class ``plrs''. However, the input object must result from
+function \texttt{plrs.test()}). This is because information from the test is required
+(mixture's weights). \texttt{plrs.cb} outputs an object of class ``plrs'' and computed lower and
+upper bounds for the mean response are stored in slot \texttt{cb}.
+
+\begin{center}
+\setkeys{Gin}{width=0.5\textwidth}
+<<fig=TRUE>>=
+# Compute and plot CBs
+selectedModel <- plrs.cb(selectedModel, alpha=0.05)
+
+str(selectedModel@cb)
+
+plot(selectedModel)
+@
+\end{center}
+\null
+Confidence bands are automatically plotted. Color may be change with input arguments \texttt{col.cb}.
+For example:
+\begin{center}
+\setkeys{Gin}{width=0.5\textwidth}
+<<fig=TRUE>>=
+plot(selectedModel, col.pts="black", col.cb="pink")
+@
+\end{center}
+
+\newpage
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\section{Analysis of multiple DNA-mRNA relationships}
+\null
+Function \texttt{plrs.series()} allows to apply the different aforementioned procedures to multiple
+relationships. Input data can be given as \emph{matrix} objects or \emph{ExpressionSet} and
+\emph{cghCall} objects. We now show how to: 1) implement the PLRS screening test and 2) implement the
+model selection procedure. For the purpose of speed, we work on a subset of chromosome 17 (first 200
+features).
+
+<<>>=
+# Testing the full model, no model selection (fast)
+neveSeries <- plrs.series(neveGE17[1:200,], neveCN17[1:200,], control.select=NULL)
+
+# Testing the full model and applying model selection
+neveSeries2 <- plrs.series(neveGE17[1:200,], neveCN17[1:200,])
+@
+
+Function \texttt{plrs.series()} has arguments \texttt{control.model}, \texttt{control.select} and
+\texttt{control.test}, which allows to specify the type of model, model selection and whether
+the test and confidence bands should be computed. Argument \texttt{control.output} allows the user
+to save each ``plrs'' objects and/or associated plot in the working directory.\\
+\texttt{plrs.series()} outputs an object of class ``plrs.series''. Generic functions \texttt{print()}
+and \texttt{summary()} display information on fitted models, selected models and/or the testing procedure.
+
+\newpage
+<<>>=
+# Summary of screening test
+neveSeries
+summary(neveSeries)
+
+# Summary of screening test and model selection
+neveSeries2
+summary(neveSeries2)
+@
+
+\noindent Results of testing and model selection procedures are obtained as follows:
+
+<<>>=
+# Testing results
+head(neveSeries@test)
+head(neveSeries2@test)
+
+# Coefficients of selected models
+head(neveSeries2@coefficients)
+
+@
+
+
+
+\newpage
+\bibliographystyle{plain}
+\bibliography{mybib}
+
+\end{document}
diff --git a/man/criteria.Rd b/man/criteria.Rd
new file mode 100644
index 0000000..6348b77
--- /dev/null
+++ b/man/criteria.Rd
@@ -0,0 +1,46 @@
+\name{criteria}
+\alias{criteria}
+\alias{criteria,plrs-method}
+
+\title{
+  Compute AIC, AICC, BIC and OSAIC for a given \code{plrs} model.
+}
+
+\description{Extract AIC, AICC, BIC and OSAIC from an object of class \code{\link{plrs-class}}.}
+
+\usage{
+criteria(obj, crit = "all")
+}
+
+\arguments{
+  \item{obj}{object of class \code{\link{plrs-class}}}
+  \item{crit}{A \code{character} (vector) among \code{"aic"}, \code{"aicc"}, \code{"bic"}, \code{"osaic"} or \code{"all"}.}
+}
+
+\value{
+  A \code{list} with the following components (if specified):
+  \item{aic}{Akaike's information criterion}
+  \item{aicc}{Small sample correction of AIC}
+  \item{bic}{Bayesian Information Criterion}
+  \item{osaic}{One-Sided AIC. See Hughes and King (2003) for more details.}
+}
+
+\references{
+  Hughes, A. W. and King, M. L. (2003). Model selection using AIC in the presence of one-sided information. \emph{J Stat Plan Infer}, 115(2): 397 411.
+}
+
+\author{Gwenael G.R. Leday \email{g.g.r.leday@vu.nl}}
+
+%\note{}
+%\seealso{}
+\examples{
+
+# Simulate data
+sim <- plrs.sim(n=80, states=4, sigma=0.5)
+
+# Fit
+model <- plrs(expr=sim$expr, cghseg=sim$seg, cghcall=sim$cal)
+
+criteria(model)
+
+}
diff --git a/man/modify.conf.Rd b/man/modify.conf.Rd
new file mode 100644
index 0000000..ee809ae
--- /dev/null
+++ b/man/modify.conf.Rd
@@ -0,0 +1,50 @@
+\name{modify.conf}
+\alias{modify.conf}
+
+\title{Modify the configuration (of calls) of the plrs model}
+
+\description{
+This function changes the discrete copy number values for a given gene in order to 
+force a minimum number of observations per state.
+}
+
+\usage{
+modify.conf(cghcall, min.obs = 3, discard = TRUE)
+}
+
+\arguments{
+  \item{cghcall}{Vector of called values}
+  \item{min.obs}{Minimum number of observations per state}
+  \item{discard}{Logical. Whether discrete states with few observations should be discarded from analysis.}
+}
+
+\details{
+Consider that the number of observations of a given state is lower than \code{min.obs}, then:
+
+- if \code{discard = FALSE}, observations are not discarded and a rearrangement of called values is carried out as follows.
+The "normal" copy number state is taken as a reference. If the minimum number of observations is not obtained,
+"losses" will be merged to "normals", "gains" to "normals" and "amplifications" to "gains".
+Note that this modifies the configuration of the model. Thus, after fitting a model using \code{\link{plrs}},
+original and modified data are stored in the resulting \code{\link{plrs-class}} object, respectively under
+slots \code{data} and \code{mdata}.
+
+- if \code{discard = TRUE}, states for which the number of observations is lower than \code{min.obs}
+are discarded (replaced by NAs).
+
+}
+
+\value{
+  \item{val}{Vector of new called values}
+}
+
+\author{Gwenael G.R. Leday \email{g.g.r.leday@vu.nl}}
+
+\note{
+This function is implemented within function \code{\link{plrs}} and \code{\link{plrs.series}}.
+}
+
+\examples{
+called <- sample(c(rep(-1,5),rep(0,15),rep(1,2),rep(2,1)))
+table(called)
+table(modify.conf(called, min.obs=3))
+}
diff --git a/man/neveCN17.Rd b/man/neveCN17.Rd
new file mode 100644
index 0000000..8d0998e
--- /dev/null
+++ b/man/neveCN17.Rd
@@ -0,0 +1,31 @@
+\name{neveCN17}
+\alias{neveCN17}
+\docType{data}
+
+\title{
+  Copy number for chromosome 17.
+}
+
+\description{Preprocessed copy number data of Neve et al. (2006) for chromosome 17.}
+
+\usage{
+neveCN17
+}
+
+\format{An object of class \code{\link[CGHbase]{cghCall}}}
+
+\source{
+M. Neve et al. in Gray Lab at LBL. Neve2006: expression and CGH data on breast cancer cell lines. R package version 0.1.10.
+}
+
+\references{
+	Neve, R.M. et al. (2006). A collection of breast cancer cell lines for the study of functionally distinct cancer subtypes. \emph{Cancer cell}, 10, 515-527.
+}
+
+
+\examples{
+data(neveCN17)
+dim(neveCN17)
+head(fData(neveCN17))
+
+}
diff --git a/man/neveGE17.Rd b/man/neveGE17.Rd
new file mode 100644
index 0000000..05b82d9
--- /dev/null
+++ b/man/neveGE17.Rd
@@ -0,0 +1,31 @@
+\name{neveGE17}
+\alias{neveGE17}
+\docType{data}
+
+\title{
+  mRNA expression for chromosome 17.
+}
+
+\description{Normalized gene expression data of Neve et al. (2006) for chromosome 17.}
+
+\usage{
+neveGE17
+}
+
+\format{An object of class \code{\link[Biobase]{ExpressionSet}}}
+
+\source{
+M. Neve et al. in Gray Lab at LBL. Neve2006: expression and CGH data on breast cancer cell lines. R package version 0.1.10.
+}
+
+\references{
+	Neve, R.M. et al. (2006). A collection of breast cancer cell lines for the study of functionally distinct cancer subtypes. \emph{Cancer cell}, 10, 515-527.
+}
+
+
+\examples{
+data(neveGE17)
+dim(neveGE17)
+head(fData(neveGE17))
+
+}
diff --git a/man/plot-methods.Rd b/man/plot-methods.Rd
new file mode 100644
index 0000000..f02b906
--- /dev/null
+++ b/man/plot-methods.Rd
@@ -0,0 +1,55 @@
+\name{plot-methods}
+\docType{methods}
+\alias{plot-methods}
+\alias{plot,ANY}
+\alias{plot.plrs}
+\alias{plot,plrs,ANY-method}
+\alias{plot,plrs.select,ANY-method}
+
+\title{Plot functions in package 'plrs'}
+
+\description{
+Methods \code{plot} in package 'plrs'
+}
+
+\usage{
+\method{plot}{plrs}(x, col.line = "black", col.pts = c("red", "blue","green2", "green4"),
+col.cb = "yellow", xlim = c(floor(min(x@data$cghseg)),ceiling(max(x@data$cghseg))),
+ylim = c(floor(min(x@data$expr)),ceiling(max(x@data$expr))),
+pch = 16, lwd=4, cex = 1.2, xlab="", ylab="", main = "",                      
+add = FALSE, lty = 1, lin = FALSE, ...)
+
+}
+
+\arguments{
+  \item{x}{An object of class \code{\link{plrs-class}} or \code{\link{plrs.select-class}}}
+  \item{col.line}{Color of the fitted line}
+  \item{col.pts}{Vector of length 4, for colors associated with each state}
+  \item{col.cb}{Color for the confidence band}
+  \item{xlim}{The x limits of the plot}
+  \item{ylim}{The y limits of the plot}
+  \item{pch}{See \code{\link[graphics]{par}}}
+  \item{lwd}{See \code{\link[graphics]{par}}}
+  \item{cex}{See \code{\link[graphics]{par}}}
+  \item{xlab}{Title of the x-axis}
+  \item{ylab}{Title of the y-axis}
+  \item{main}{Main title for the plot}
+  \item{add}{If the plot should be added to the current device. Default is \code{FALSE}}
+  \item{lty}{See \code{\link[graphics]{par}}}
+  \item{lin}{Logical. Whether the simple linear model should also be plotted}
+  \item{...}{Other arguments, see \code{\link[graphics]{par}}}
+}
+
+\section{Methods}{
+	\describe{
+		\item{\code{signature(x = "plrs")}}{Plot observed points and the fitted line}
+				\item{\code{signature(x = "plrs.select")}}{Plot observed points and the fitted line of the selected model.}
+	}
+}
+
+\details{
+\code{plot.plrs} plots the observed points, the fitted line and potentially the confidence band.
+}
+
+\author{Gwenael G.R. Leday \email{g.g.r.leday@vu.nl}}
+
diff --git a/man/plrs-class.Rd b/man/plrs-class.Rd
new file mode 100644
index 0000000..b816279
--- /dev/null
+++ b/man/plrs-class.Rd
@@ -0,0 +1,58 @@
+\name{plrs-class}
+\Rdversion{1.1}
+\docType{class}
+\alias{plrs-class}
+\alias{coef,plrs-method}
+\alias{effects,plrs-method}
+\alias{fitted,plrs-method}
+\alias{knots,plrs-method}
+\alias{model.matrix,plrs-method}
+
+\alias{predict,plrs-method}
+\alias{print,plrs-method}
+\alias{residuals,plrs-method}
+\alias{show,plrs-method}
+\alias{summary,plrs-method}
+
+\title{Class \code{plrs}}
+\description{
+An S4 class representing the output of the \code{\link{plrs}} function.
+}
+
+\section{Slots}{
+  \describe{
+    \item{\code{coefficients}:}{Object of class \code{numeric} containing spline coefficients}
+    \item{\code{fitted.values}:}{Object of class \code{numeric} containing the fitted values}
+    \item{\code{residuals}:}{Object of class \code{numeric} containing the residuals}
+    \item{\code{X}:}{Object of class \code{matrix} containing the design matrix}
+    \item{\code{data}:}{Object of class \code{list} containing input data}
+    \item{\code{mdata}:}{Object of class \code{list} containing (possibly modified) data used to fit the model (See \code{\link{modify.conf}}).}
+    \item{\code{QP}:}{Object of class \code{list} containing input elements used for quadratic programming. If the model is unconstrained this contains a light version of an \code{lm} object.}
+    \item{\code{test}:}{Object of class \code{list} containing results from testing.}
+    \item{\code{cb}:}{Object of class \code{list} containing lower and upper bounds for predicted values.}   
+    \item{\code{selected}:}{Object of class \code{logical} indicating whether the model results from a selection procedure.}
+    \item{\code{type}:}{Object of class \code{character} giving the type of model}
+    \item{\code{call.arg}:}{Object of class \code{list} containing the input arguments (for reproducibility)}
+  }
+}
+
+\section{Methods}{
+  \describe{
+    \item{coef}{Returns the coefficients}
+    \item{criteria}{See \code{\link{criteria}}}
+    \item{effects}{Returns matrix of effects}
+    \item{fitted}{Returns the fitted values}
+    \item{knots}{Returns the knots}
+    \item{model.matrix}{Returns the design matrix}
+    \item{plot}{See \code{\link{plot.plrs}}}
+    \item{predict}{See \code{\link{predict.plrs}}}
+    \item{print}{Print the object information}
+    \item{residuals}{Returns the residuals}
+    \item{show}{Print the object information}
+    \item{summary}{Print a summary of the object information}
+  }
+}
+
+\author{Gwenael G.R. Leday \email{g.g.r.leday@vu.nl}}
+
+
diff --git a/man/plrs-package.Rd b/man/plrs-package.Rd
new file mode 100644
index 0000000..2a67b0a
--- /dev/null
+++ b/man/plrs-package.Rd
@@ -0,0 +1,41 @@
+\name{plrs-package}
+\alias{plrs-package}
+\docType{package}
+
+\title{Piecewise Linear Regression Splines (PLRS) for the association between DNA copy number and mRNA expression}
+
+\description{
+The present package implements a framework for modeling the relationship between DNA copy number and gene expression data using Piecewise Linear Regression Splines (PLRS). It includes (point and interval) estimation, model selection and testing procedures for such models (possibly under biologically motivated constraints).
+}
+
+\details{
+The use of the present package can be divided into two approaches:\cr
+\cr
+1. Analysis of a single DNA-mRNA relationship\cr
+\cr
+Main functions are:\cr
+\code{\link{plrs}}:	Fit a single plrs model.\cr
+\code{\link{plrs.select}}: Model selection based on AIC, AICC, OSAIC or BIC.\cr
+\code{\link{plrs.test}}:	Likelihood ratio test for a given plrs model.\cr
+\code{\link{plrs.cb}}: Confidence bands for a plrs model.\cr
+\cr
+\cr
+2. Analysis of multiple DNA-mRNA relationships sequentially\cr
+\cr
+Main function is:\cr
+\code{\link{plrs.series}}: point and interval estimation, model selection and testing of DNA-mRNA association for a series of arrays.\cr
+\cr
+Note: This function extend the aforementioned univariate analysis genomewise in the same spirit as some functions of the \pkg{limma} package do. 
+}
+
+\author{
+Gwenael G.R. Leday
+
+Maintainer: Gwenael G.R. Leday \email{g.g.r.leday@vu.nl}
+}
+
+\references{
+Leday GGR, Van der Vaart AW, Van Wieringen WN, Van de Wiel MA. Modeling association between DNA copy number and gene expression with constrained piecewise linear regression splines. Accepted for publication. \emph{Ann Appl Stat}. (2012).
+}
+
+\keyword{copy number, gene expression, regression splines, model selection, constrained inference.}
\ No newline at end of file
diff --git a/man/plrs.Rd b/man/plrs.Rd
new file mode 100644
index 0000000..b6b9ab2
--- /dev/null
+++ b/man/plrs.Rd
@@ -0,0 +1,79 @@
+\name{plrs}
+\alias{plrs}
+\alias{plrs,ANY}
+
+\title{Fit a (constrained) piecewise linear regression spline}
+
+\description{The function fits a piecewise linear regression spline to explain 
+gene expression by the segmented DNA copy number. The called copy number values 
+are used as a template for model building.}
+
+\usage{
+plrs(expr, cghseg, cghcall=NULL, probloss = NULL, probnorm = NULL,
+probgain = NULL, probamp = NULL, knots = NULL, continuous = FALSE,
+constr = TRUE, constr.slopes = 2, constr.intercepts = TRUE,
+min.obs = 3, discard.obs = TRUE)
+}
+
+\arguments{
+  \item{expr}{Vector of gene expression values}
+  \item{cghseg}{Vector of segmented copy number values}
+  \item{cghcall}{Vector of called copy number values. If not provided, we are reduced to a simple linear model.}
+  \item{probloss}{Vector of call probabilities associated with state "loss". Default is \code{NULL}.}
+  \item{probnorm}{Vector of call probabilities associated with state "normal". Default is \code{NULL}.}
+  \item{probgain}{Vector of call probabilities associated with state "gain". Default is \code{NULL}.}
+  \item{probamp}{Vector of call probabilities associated with state "amplification". Default is \code{NULL}.}
+  \item{knots}{knots or change points. If \code{NULL} (default), there are estimated. See details.}
+  \item{continuous}{Logical, whether the model is continuous (no jump) or not.}
+  \item{constr}{Logical, whether the model is constrained or not. (this has been implemented to turn on and off easily the constraints)}
+  \item{constr.slopes}{Type of non-negativity constraints applied on slopes. Either 1 or 2 (default). See details.}
+  \item{constr.intercepts}{If \code{TRUE} (default) jumps from state to state are also constrained to be non-negative}
+  \item{min.obs}{See \code{\link{modify.conf}}}
+  \item{discard.obs}{See \code{\link{modify.conf}}}
+}
+
+\details{
+If \code{cghcall=NULL}, discrete copy number values are omitted, which results in fitting a simple linear model.\cr
+\cr
+If \code{constr.slopes=1}, all slopes are constrained to be non-negative.
+If \code{constr.slopes=2}, the slope associated with state "normal" is constrained to be non-negative and all others are forced to be at least equal to the latter.\cr
+\cr
+Two methods are implemented for the estimation of knots. If call probabilities are provided,
+a knot is determined so that the sum of (the two adjacent) states membership probabilities is maximized.
+Otherwise, this is defined as the midpoint of the interval between the two consecutive states.\cr
+
+The constrained least squares problem is solved using function \code{solve.QP} of package \pkg{quadprog}.
+}
+
+\value{An object of class \code{\link{plrs-class}}}
+
+%\references{}
+
+\author{Gwenael G.R. Leday \email{g.g.r.leday@vu.nl}}
+
+%\note{}
+%\seealso{}
+
+\examples{
+
+# Simulate data
+sim <- plrs.sim(n=80, states=4, sigma=0.5)
+
+# Fit a model 
+model <- plrs(expr=sim$expr, cghseg=sim$seg, cghcall=sim$cal)
+model
+
+# Methods
+coef(model)
+effects(model)
+fitted(model)
+knots(model)
+model.matrix(model)
+plot(model)
+predict(model, newcghseg=seq(0,5, length.out=100))   
+residuals(model)
+summary(model)
+}
+
+
+
diff --git a/man/plrs.cb.Rd b/man/plrs.cb.Rd
new file mode 100644
index 0000000..504d0c3
--- /dev/null
+++ b/man/plrs.cb.Rd
@@ -0,0 +1,54 @@
+\name{plrs.cb}
+\alias{plrs.cb}
+
+\title{Uniform confidence bands (CB) for plrs models}
+
+\description{
+Determine uniform confidence intervals for predicted values of a 'plrs' model.
+}
+
+\usage{
+plrs.cb(object, alpha=0.05, newcgh=NULL)
+}
+
+\arguments{
+  \item{object}{An object of class \code{\link{plrs-class}}.}
+  \item{alpha}{Significance level}
+  \item{newcgh}{Vector of segmented values. Support for building CB.}
+}
+
+\details{
+The input object of class \code{\link{plrs-class}} has to result from function \code{\link{plrs.test}}.\cr
+
+The problem of finding (at a given x) a confidence interval for the mean response is expressed as a
+semi-definite optimization problem and solved using function \code{csdp} of package \pkg{Rcsdp}.
+}
+
+\value{
+An object of class \code{\link{plrs-class} that contains CB information.}
+}
+
+\references{
+Leday GGR, Van der Vaart AW, Van Wieringen WN, Van de Wiel MA. Modeling association between DNA copy number and gene expression with constrained piecewise linear regression splines. Accepted for publication. \emph{Ann Appl Stat}. (2012).
+}
+
+\author{Gwenael G.R. Leday \email{g.g.r.leday@vu.nl}}
+
+\seealso{
+\code{\link{plrs.test}}
+}
+
+\examples{
+
+# Simulate data
+sim <- plrs.sim(n=80, states=4, sigma=0.5)
+
+# Fit a model 
+model <- plrs(expr=sim$expr, cghseg=sim$seg, cghcall=sim$cal)
+
+# Confidence bands
+model <- plrs.test(model)
+model <- plrs.cb(model, alpha=0.05)
+plot(model)
+}
+
diff --git a/man/plrs.select-class.Rd b/man/plrs.select-class.Rd
new file mode 100644
index 0000000..eb71a70
--- /dev/null
+++ b/man/plrs.select-class.Rd
@@ -0,0 +1,33 @@
+\name{plrs.select-class}
+\Rdversion{1.1}
+\docType{class}
+\alias{plrs.select-class}
+\alias{print,plrs.select-method}
+\alias{show,plrs.select-method}
+\alias{summary,plrs.select-method}
+
+\title{Class \code{plrs.select}}
+\description{
+An S4 class representing the output of the \code{\link{plrs.select}} function.
+}
+
+\section{Slots}{
+  \describe{
+    \item{\code{table}:}{Object of class \code{matrix} containing the criterion value for all models}
+    \item{\code{model}:}{Object of class \code{plrs} containing the selected model}
+    \item{\code{crit}:}{Object of class \code{character} containing the criterion used for model selection}
+  }
+}
+
+\section{Methods}{
+  \describe{
+    \item{plot}{See \code{\link{plot.plrs}}}
+    \item{print}{Print the object information}
+    \item{show}{Print the object information}
+    \item{summary}{Print a summary of the object information}
+  }
+}
+
+\author{Gwenael G.R. Leday \email{g.g.r.leday@vu.nl}}
+
+
diff --git a/man/plrs.select.Rd b/man/plrs.select.Rd
new file mode 100644
index 0000000..80600f9
--- /dev/null
+++ b/man/plrs.select.Rd
@@ -0,0 +1,30 @@
+\name{plrs.select}
+\alias{plrs.select}
+\alias{plrs.select,ANY}
+
+\title{Model selection}
+
+\description{Selection of a model based on an information criterion (AIC, AICC, BIC or OSAIC).}
+
+\usage{
+plrs.select(object, crit = ifelse(object@call.arg$constr,"osaic","aic"))
+}
+
+\arguments{
+  \item{object}{An object of class \code{\link{plrs-class}}}
+  \item{crit}{Character corresponding to the criterion to use. See \code{\link{criteria}}.}
+}
+
+%\details{
+
+%}
+
+\value{An object of class \code{\link{plrs.select-class}}
+}
+
+\author{Gwenael G.R. Leday \email{g.g.r.leday@vu.nl}}
+
+%\note{}
+%\seealso{}
+%\examples{}
+%\keyword{}
diff --git a/man/plrs.series-class.Rd b/man/plrs.series-class.Rd
new file mode 100644
index 0000000..8739abd
--- /dev/null
+++ b/man/plrs.series-class.Rd
@@ -0,0 +1,35 @@
+\name{plrs.series-class}
+\Rdversion{1.1}
+\docType{class}
+\alias{plrs.series-class}
+\alias{print,plrs.series-method}
+\alias{show,plrs.series-method}
+\alias{summary,plrs.series-method}
+
+\title{Class \code{plrs.series}}
+\description{
+An S4 class representing the output of the \code{\link{plrs.series}} function.
+}
+
+\section{Slots}{
+  \describe{
+    \item{\code{coefficients}:}{Matrix containing coefficients of models}
+    \item{\code{effects}:}{List containing effects}                                                                
+    \item{\code{test}:}{Matrix containing results from testing.}
+    \item{\code{general}:}{Matrix providing the distribution of the number genes and arrays regarding the copy number states}
+    \item{\code{modelsType}:}{List providing models' type}
+    \item{\code{call.arg}:}{List providing details on the type of models that have been fitted.}
+  }
+}
+
+\section{Methods}{
+  \describe{
+    \item{print}{Print the object information}
+    \item{show}{Print the object information}
+    \item{summary}{Print a summary of the object information}
+  }
+}
+
+\author{Gwenael G.R. Leday \email{g.g.r.leday@vu.nl}}
+
+
diff --git a/man/plrs.series.Rd b/man/plrs.series.Rd
new file mode 100644
index 0000000..4004927
--- /dev/null
+++ b/man/plrs.series.Rd
@@ -0,0 +1,101 @@
+\name{plrs.series}
+\alias{plrs.series}
+
+\title{
+Fit plrs models for a series of arrays.
+}
+
+\description{
+The function fits \code{plrs} models for a series of arrays. Model selection and testing procedures may be applied.
+}
+
+\usage{
+plrs.series(expr, cghseg, cghcall=NULL,
+probloss = NULL, probnorm = NULL, probgain = NULL, probamp = NULL,
+control.model  = list(continuous = FALSE,
+                      constr = TRUE,
+                      constr.slopes = 2,
+                      constr.intercepts = TRUE,
+                      min.obs = 3,
+                      discard.obs = TRUE),
+control.select = list(crit = ifelse(control.model$constr, "osaic","aic")),
+control.test   = list(testing = TRUE,
+                      cb = FALSE,
+                      alpha = 0.05),
+control.output = list(save.models = FALSE,
+                      save.plots = FALSE,
+                      plot.lin = FALSE,
+                      type = "jpeg"))
+}
+
+\arguments{
+  \item{expr}{Either a matrix of expression profiles or an \code{\link[Biobase]{ExpressionSet}} object.}
+  \item{cghseg}{Either a matrix of segmented copy number values or objects of class \code{\link[CGHbase]{cghSeg}} or \code{\link[CGHbase]{cghCall}}}
+  \item{cghcall}{Matrix of called copy number}
+  \item{probloss}{Matrix of call probabilities associated with state "loss". Default is \code{NULL}.}
+  \item{probnorm}{Matrix of call probabilities associated with state "normal". Default is \code{NULL}.}
+  \item{probgain}{Matrix of call probabilities associated with state "gain". Default is \code{NULL}.}
+  \item{probamp}{Matrix of call probabilities associated with state "amplification". Default is \code{NULL}.}
+  \item{control.model}{See details}
+  \item{control.select}{See details}
+  \item{control.test}{See details}
+  \item{control.output}{See details}
+}
+
+\details{
+If DNA and mRNA input data are matrices, rows should correspond to genes and columns to arrays.
+Alternatively, expression data may be provided as an \code{\link[Biobase]{ExpressionSet}} object and aCGH data
+as \code{\link[CGHbase]{cghSeg}} or \code{\link[CGHbase]{cghCall}} objects. A \code{\link[CGHbase]{cghCall}} object
+contain all data from the calling step, thus arguments \code{probloss}, \code{probnorm}, \code{probnorm} and \code{probamp}
+can be omitted. An object of class \code{\link[CGHbase]{cghSeg}} does not contain such data so only simple linear models
+will be fitted.\cr
+\cr
+\code{control.model} allows the user to specify the type of model that has to be fitted. 
+This must be a \code{list} with one or more of the following components:
+\code{constr}, \code{constr.slopes}, \code{constr.intercepts}, \code{min.obs} and \code{discard.obs}.
+See functions \code{\link{plrs}} and \code{\link{modify.conf}} for more details.\cr
+\cr
+\code{control.select} allows the user to specify whether model selection should be done and how.
+This must be a \code{list} with a component named \code{crit}. See function \code{\link{plrs.select}} for more details.
+If \code{control.select = NULL} then no model selection is done. \cr
+\cr
+\code{control.output} allows the user to plot and save each \code{plrs} model. This must be a list with components:\cr
+\code{save.models}, a logical. This will create within the work directory a new directory named "plrsSeriesObjects" that will contain all objects.\cr
+\code{save.plots}, a logical. This will create within the work directory a new directory named "plrsSeriesPlots" that will contains all saved plots.\cr
+\code{plot.lin}, a logical. Whether the simple linear model should aslo be plotted.\cr
+\code{type}, a character. Format of file. To pass through function \code{\link[grDevices]{savePlot}}.\cr
+}
+
+\value{An object of class \code{\link{plrs.series-class}}}
+
+\author{Gwenael G.R. Leday \email{g.g.r.leday@vu.nl}}
+
+%\note{}
+%\seealso{}
+
+\examples{
+
+# Simulate data
+ngenes <- 10
+narray <- 48
+rna <- dnaseg <- dnacal <- matrix(NA, ngenes, narray)
+idx <- sample(1:4, ngenes, replace=TRUE, prob=rep(1/4,4))
+for(i in 1:ngenes){
+	Sim <- plrs.sim(n=narray, states=idx[i], sigma=0.5)
+	rna[i,] <- Sim$expr
+	dnaseg[i,] <- Sim$seg
+	dnacal[i,] <- Sim$cal
+}
+
+
+# Screening procedure with linear model
+series <- plrs.series(expr = rna, cghseg = dnaseg, cghcall = NULL, control.select = NULL)
+
+# Screening procedure with full plrs model
+series <- plrs.series(expr = rna, cghseg = dnaseg, cghcall = dnacal, control.select = NULL)
+
+# Model selection
+series <- plrs.series(expr = rna, cghseg = dnaseg, cghcall = dnacal)
+
+}
+
diff --git a/man/plrs.sim.Rd b/man/plrs.sim.Rd
new file mode 100644
index 0000000..9fd5c54
--- /dev/null
+++ b/man/plrs.sim.Rd
@@ -0,0 +1,55 @@
+\name{plrs.sim}
+\alias{plrs.sim}
+
+\title{
+Simulation of a plrs model
+}
+
+\description{
+Simulation of a piecewise relationship.\cr
+\cr
+The function has been only implemented for convenience of simulations and R examples.\cr
+}
+
+\usage{
+plrs.sim(n = 80, states = 4, sigma = 01, x = NULL)
+}
+
+\arguments{
+  \item{n}{Number of simulated data points}
+  \item{states}{Number of states for the model}
+  \item{sigma}{Noise}
+  \item{x}{Segmented values.}
+}
+
+\details{
+To be written...
+}
+
+%\value{}
+
+\author{Gwenael G.R. Leday \email{g.g.r.leday@vu.nl}}
+
+\examples{
+
+# Simulate 1-state model
+sim <- plrs.sim(n=80, states=1, sigma=0.5)
+model <- plrs(expr=sim$expr, cghseg=sim$seg, cghcall=sim$cal)
+plot(model)
+
+# Simulate 2-state model
+sim <- plrs.sim(n=80, states=2, sigma=0.5)
+model <- plrs(expr=sim$expr, cghseg=sim$seg, cghcall=sim$cal)
+plot(model)
+
+# Simulate 3-state model
+sim <- plrs.sim(n=90, states=3, sigma=0.5)
+model <- plrs(expr=sim$expr, cghseg=sim$seg, cghcall=sim$cal)
+plot(model)
+
+# Simulate 4-state model
+sim <- plrs.sim(n=80, states=4, sigma=0.5)
+model <- plrs(expr=sim$expr, cghseg=sim$seg, cghcall=sim$cal)
+plot(model)
+
+}
diff --git a/man/plrs.test.Rd b/man/plrs.test.Rd
new file mode 100644
index 0000000..ae90982
--- /dev/null
+++ b/man/plrs.test.Rd
@@ -0,0 +1,70 @@
+\name{plrs.test}
+\alias{plrs.test}
+
+\title{Likelihood ratio test for a \code{plrs} model}
+
+\description{Test whether copy number has an effect on mRNA expression.}
+
+\usage{
+plrs.test(object, alpha=0.05)
+}
+
+\arguments{
+  \item{object}{An object of class \code{\link{plrs-class}}}
+  \item{alpha}{Significance level}
+}
+
+\details{
+Two cases present themselves:\cr
+\cr
+1. The model is unconstrained. Thus, the model under the null hypothesis is the intercept
+and an F-test is performed.\cr
+\cr
+2. The model is constrained and the following hypothesis are tested: \cr
+H0: All constraints are actives (=)\cr
+H1: At least one constraint is strict (>)\cr
+Under H0, we always have the intercept model. Indeed, if \code{constr.slopes = 1} (or 2) and
+\code{constr.intercepts = T}, then the only parameter free of inequality constraint is the 
+overall intercept. If \code{constr.intercepts = F}, the local intercepts are additionnaly 
+constrained to be 0 in order to obtain the intercept model under the null.
+The likelihood ratio statistic (unknown variance) is asymptotically distributed as a
+weighted mixture of Beta distribution (cf Gromping (2010)). Calculation of p-values is based on
+functions \code{ic.weights} and \code{pbetabar} of package \pkg{ic.infer}. The package 
+\pkg{mvtnorm} is also involved.\cr
+\cr
+In both cases the input model is taken as the model under the alternative.
+}
+
+\value{A \code{list} object with the following components:
+  \item{stat}{Test statistic}
+  \item{pvalue}{Calculated pvalue}
+  \item{wt.bar}{Weights (if the model is constrained)}
+  \item{df.bar}{Degrees of freedom.}
+  \item{unconstr}{Unconstrained model of class \code{\link{plrs-class}}}
+  \item{qbetabar}{(1-\code{alpha}) quantile of the beta mixture distribution}
+  \item{alpha}{Significance level}
+}
+
+\references{
+Gromping, U. (2010). Inference with linear equality and inequality constraints using
+R: The package ic.infer. \emph{J Stat Softw}, 33(i10).
+}
+
+\author{Gwenael G.R. Leday  \email{g.g.r.leday@vu.nl}}
+
+%\note{}
+%\seealso{}
+\examples{
+
+# Simulate data
+sim <- plrs.sim(n=80, states=2, sigma=0.5)
+
+# Fit a model 
+model <- plrs(expr=sim$expr, cghseg=sim$seg, cghcall=sim$cal)
+
+# Testing
+model <- plrs.test(model)
+model
+
+}
+
diff --git a/man/predict.plrs.Rd b/man/predict.plrs.Rd
new file mode 100644
index 0000000..df05381
--- /dev/null
+++ b/man/predict.plrs.Rd
@@ -0,0 +1,30 @@
+\name{predict.plrs}
+\alias{predict.plrs}
+
+\title{Predict method for \code{plrs} models}
+
+\description{Determine predicted values based on a given plrs model}
+
+\usage{
+\method{predict}{plrs}(object, newcghseg, ...)
+}
+
+\arguments{
+  \item{object}{An object of class \code{\link{plrs-class}}}
+  \item{newcghseg}{A vector of new segmented CGH values}
+  \item{...}{further arguments}
+}
+
+%\details{}
+
+\value{A \code{vector} containing the fitted values.
+}
+
+%\references{}
+
+\author{Gwenael G.R. Leday  \email{g.g.r.leday@vu.nl}}
+
+%\note{}
+%\seealso{}
+%\examples{}
+%\keyword{}