R/uni.smooth.const-with-po.r

In scam: Shape Constrained Additive Models

## (c) Natalya Pya (2012-2024). Provided under GPL 2.
## routines for univariate SCOP-spline construction 
## with sum-to-zero identifiability constraints (2023)



#####################################################
### Adding Monotone increasing SCOP-spline construction 
######################################################

smooth.construct.mpi.smooth.spec<- function(object, data, knots)
## construction of the monotone increasing smooth
{ 
   #require(splines)
  m <- object$p.order[1]
  if (is.na(m)) m <- 2 ## default for cubic spline
  if (m < 0) stop("silly m supplied")
  if (object$bs.dim<0) object$bs.dim <- 10 ## default
  q <- object$bs.dim 
  nk <- q+m+2 ## number of knots
  if (nk<=0) stop("k too small for m")
  x <- data[[object$term]]  ## the data
  xk <- knots[[object$term]] ## will be NULL if none supplied
  if (is.null(xk)) # space knots through data
     { n<-length(x)
       xk<-rep(0,q+m+2)
       xk[(m+2):(q+1)]<-seq(min(x),max(x),length=q-m)
       for (i in 1:(m+1)) {xk[i]<-xk[m+2]-(m+2-i)*(xk[m+3]-xk[m+2])}
       for (i in (q+2):(q+m+2)) {xk[i]<-xk[q+1]+(i-q-1)*(xk[m+3]-xk[m+2])}
     }
  if (length(xk)!=nk) # right number of knots?
   stop(paste("there should be ", nk," supplied knots"))
  if (!is.null(object$point.con[[1]])) ## a point constraint is supplied?
    stop(paste("'mpi' smooth does not work with a point constraint; use 'miso' for a start-at-zero constraint, or 'mifo' for a finish-at-zero constraint"))

  #######################################################################
  ## indentifiability constraint by first dropping the first columns of XSigma (setting beta_1=0)
  ## and then applying sum-to-zero (centering) constraint...
  #  get model matrix-------------
  X1 <- splineDesign(xk,x,ord=m+2)
  # get matrix Sigma and remove the first column for the model matrix XSig
 ##  Sig <- matrix(as.numeric(rep(1:q,q)>=rep(1:q,each=q)),q,q) ## coef summation matrix
  Sig <- matrix(1,q,q)  ## coef summation matrix
  Sig[upper.tri(Sig)] <-0
  X <- X1%*%Sig
  X <- X[,-1]
  
  ## applying sum-to-zero (centering) constraint...
  cmx <- colMeans(X)
  X <- sweep(X,2,cmx) ## subtract cmx from columns 
  object$X <- X # the final model matrix
  object$cmX <- cmx

  object$P <- list()
  object$S <- list()
  object$Sigma <- Sig[-1,-1]

  if (!object$fixed) # create the penalty matrix
  { P <- diff(diag(q-1),difference=1)
    object$P[[1]] <- P
    object$S[[1]] <- crossprod(P)
  }
  object$p.ident <- rep(TRUE,q-1) ## p.ident is an indicator of which coefficients must be positive (exponentiated)

  object$rank <- ncol(object$X)-1  # penalty rank
  object$null.space.dim <- 2 ##  ##m+1  # dim. of unpenalized space, 2 as the basis of a straight line is two-dimensional
  object$C <- matrix(0, 0, ncol(X)) # to have no other constraints 
  object$knots <- xk;
  object$m <- m;
  object$df<-ncol(object$X)     # maximum DoF 

  ## get model matrix for 1st and 2nd derivatives of the smooth...
  h <- (max(x)-min(x))/(q-m-1) ## distance between two adjacent knots
  object$Xdf1 <- splineDesign(xk,x,ord=m+1)[,2:(q-1)]/h ## ord is by one less for the 1st derivative
  object$Xdf2 <- if (m==0) matrix(0,nrow(X),ncol(object$Xdf1)-1) ## for piecewise linear splines
                 else splineDesign(xk,x,ord=m)[,2:(q-2)]/h^2 ## ord is by two less for the 2nd derivative 

  class(object)<-"mpi.smooth"  # Give object a class
  object
}


## Prediction matrix for the `mpi` smooth class... 

Predict.matrix.mpi.smooth<-function(object,data)
## prediction method function for the `mpi' smooth class
{ m <- object$m # spline order, m+1=3 default for cubic spline
  q <- object$df +1
 # Sig <- matrix(as.numeric(rep(1:q,q)>=rep(1:q,each=q)),q,q) ## coef summation matrix
  Sig <- matrix(1,q,q)  ## coef summation matrix
  Sig[upper.tri(Sig)] <-0
  ## find spline basis inner knot range...
  ll <- object$knots[m+2];ul <- object$knots[length(object$knots)-m-1]
  x <- data[[object$term]]
  n <- length(x)
  ind <- x<=ul & x>=ll ## data in range
  if (sum(ind)==n) { ## all in range
     X <- spline.des(object$knots,x,m+2)$design
     X <- X%*%Sig 
  } else { ## some extrapolation needed 
     ## matrix mapping coefs to value and slope at end points...
     D <- spline.des(object$knots,c(ll,ll,ul,ul),m+2,c(0,1,0,1))$design
     X <- matrix(0,n,ncol(D)) ## full predict matrix
     if (sum(ind)> 0)  X[ind,] <- spline.des(object$knots,x[ind],m+2)$design ## interior rows
     ## Now add rows for linear extrapolation...
     ind <- x < ll 
     if (sum(ind)>0) X[ind,] <- cbind(1,x[ind]-ll)%*%D[1:2,]
     ind <- x > ul
     if (sum(ind)>0) X[ind,] <- cbind(1,x[ind]-ul)%*%D[3:4,]
     X <- X%*%Sig 
  ## X <- sweep(X,2,c(0,object$cmX))
  }
  X <- sweep(X,2,c(0,object$cmX))
  X
}

####################################################################################################
### Adding monotone increasing SCOP-spline construction without applying identifiability constraints
### to be used with numeric 'by' variable...
####################################################################################################
## when 'by' variable takes more than one value, the smooth terms are identifiable without a
## 'zero intercept' constraint, so they are left unconstrained.  

smooth.construct.mpiBy.smooth.spec<- function(object, data, knots)
## construction of the monotone increasing smooth
{ 
  m <- object$p.order[1]
  if (is.na(m)) m <- 2 ## default for cubic spline
  if (m<0) stop("silly m supplied")
  if (object$bs.dim<0) object$bs.dim <- 10 ## default
  q <- object$bs.dim 
  nk <- q+m+2 ## number of knots
  if (nk<=0) stop("k too small for m")
  x <- data[[object$term]]  ## the data
  xk <- knots[[object$term]] ## will be NULL if none supplied
  if (is.null(xk)) # space knots through data
     { n<-length(x)
       xk<-rep(0,q+m+2)
       xk[(m+2):(q+1)]<-seq(min(x),max(x),length=q-m)
       for (i in 1:(m+1)) {xk[i]<-xk[m+2]-(m+2-i)*(xk[m+3]-xk[m+2])}
       for (i in (q+2):(q+m+2)) {xk[i]<-xk[q+1]+(i-q-1)*(xk[m+3]-xk[m+2])}
     }
  if (length(xk)!=nk) # right number of knots?
   stop(paste("there should be ", nk," supplied knots"))
  if (!is.null(object$point.con[[1]])) ## a point constraint is supplied?
    stop(paste("'mpi' smooth does not work with a point constraint; use 'miso' for a start-at-zero constraint, or 'mifo' for a finish-at-zero constraint"))

  #  get model matrix-------------
  X1 <- splineDesign(xk,x,ord=m+2)
  # get matrix Sigma...
  Sig <- matrix(1,q,q)  ## coef summation matrix
  Sig[upper.tri(Sig)] <-0
  X <- X1%*%Sig 
  object$X <- X # the final model matrix
  object$P <- list()
  object$S <- list()
  object$Sigma <- Sig

  if (!object$fixed) # create the penalty matrix
  { P <- diff(diag(q-1),difference=1)
    P <- rbind(rep(0,q-1),P) ## adding 1st row of zeros
    P <- cbind(rep(0,q-1),P) ## adding first column of zeros
    object$P[[1]] <- P
    object$S[[1]] <- crossprod(P)
  }
  object$p.ident <- c(FALSE,rep(TRUE,q-1)) ## p.ident is an indicator of which coefficients must be positive (exponentiated)

  object$rank <- ncol(object$X)  # penalty rank
  object$null.space.dim <- 2 ## m+1  # dim. of unpenalized space
  object$C <- matrix(0, 0, ncol(X)) # to have no other constraints 
  object$knots <- xk;
  object$m <- m;
  object$df<-ncol(object$X)     # maximum DoF 

  ## get model matrix for 1st and 2nd derivatives of the smooth...
  h <- (max(x)-min(x))/(q-m-1) ## distance between two adjacent knots
  object$Xdf1 <- splineDesign(xk,x,ord=m+1)[,1:(q-1)]/h ## ord is by one less for the 1st derivative
  object$Xdf2 <- if (m==0) matrix(0,nrow(X),ncol(object$Xdf1)-1) ## for piecewise linear splines
                 else splineDesign(xk,x,ord=m)[,2:(q-2)]/h^2 ## ord is by two less for the 2nd derivative 
  class(object)<-"mpiBy.smooth"  # Give object a class
  object
}


## Prediction matrix for the `mpiBy` smooth class... 
Predict.matrix.mpiBy.smooth<-function(object,data)
## prediction method function for the `mpiBy' smooth class
{ m <- object$m # spline order, m+1=3 default for cubic spline
  q <- object$df 
  Sig <- matrix(1,q,q)  ## coef summation matrix
  Sig[upper.tri(Sig)] <-0
  ## find spline basis inner knot range...
  ll <- object$knots[m+2];ul <- object$knots[length(object$knots)-m-1]
  x <- data[[object$term]]
  n <- length(x)
  ind <- x<=ul & x>=ll ## data in range
  if (sum(ind)==n) { ## all in range
     X <- spline.des(object$knots,x,m+2)$design
     X <- X%*%Sig 
  } else { ## some extrapolation needed 
     ## matrix mapping coefs to value and slope at end points...
     D <- spline.des(object$knots,c(ll,ll,ul,ul),m+2,c(0,1,0,1))$design
     X <- matrix(0,n,ncol(D)) ## full predict matrix
     if (sum(ind)> 0)  X[ind,] <- spline.des(object$knots,x[ind],m+2)$design ## interior rows
     ## Now add rows for linear extrapolation...
     ind <- x < ll 
     if (sum(ind)>0) X[ind,] <- cbind(1,x[ind]-ll)%*%D[1:2,]
     ind <- x > ul
     if (sum(ind)>0) X[ind,] <- cbind(1,x[ind]-ul)%*%D[3:4,]
     X <- X%*%Sig 
  }
  X
}


##############################################################################################
### Adding Monotone increasing SCOP-spline construction with a 'start at zero' constraint,
### a SCOP-spline additionally constrained to be zero at a start, at the left-end point of the covariate range...
#############################################################################################

smooth.construct.miso.smooth.spec<- function(object, data, knots)
## construction of the monotone increasing smooth with a 'start-at-zero' constraint;
## achieved simply by setting the first (m+1) spline coefficients to zero...
{ 
  m <- object$p.order[1]
  if (is.na(m)) m <- 2 ## default for cubic spline
  if (m<0) stop("silly m supplied")
  if (object$bs.dim<0) object$bs.dim <- 10 ## default
  q <- object$bs.dim 
  nk <- q+m+2 ## number of knots
  if (nk<=0) stop("k too small for m")
  x <- data[[object$term]]  ## the data
  xk <- knots[[object$term]] ## will be NULL if none supplied
  if (is.null(xk)) { # space knots evenly through data
       n<-length(x)
       xk<-rep(0,q+m+2)
       xk[(m+2):(q+1)]<-seq(min(x),max(x),length=q-m)
       for (i in 1:(m+1)) {xk[i]<-xk[m+2]-(m+2-i)*(xk[m+3]-xk[m+2])}
       for (i in (q+2):(q+m+2)) {xk[i]<-xk[q+1]+(i-q-1)*(xk[m+3]-xk[m+2])}
     }
  ## move m knots, xk[m+3],...,xk[m+3+m-1] (2nd,...,(m+1) inner knots), to the first inner knot xk[m+2],
  ## join these m knots with the constrained xk[m+2] knot to avoid a plateau start...
  xk[(m+3):(m+2+m)] <- xk[m+2]
  if (length(xk)!=nk) # right number of knots?
   stop(paste("there should be ", nk," supplied knots"))
  if (!is.null(object$point.con[[1]])) ## a point constraint is supplied?
    stop(paste("'miso' smooth works only with a 'start at zero' constraint; 'pc' argument does not work here"))
  ## get model matrix...
  X1 <- splineDesign(xk,x,ord=m+2)
  ## get unconstrained matrix Sigma and remove the first m+1 columns and rows...
  Sig <- matrix(1,q,q)  ## coef summation matrix
  Sig[upper.tri(Sig)] <-0
  ind <- 1:(m+1) ## 1:3 if m=2;
  ## Sig[,ind] <- 0; Sig[ind,] <- 0
  Sig <- Sig[-ind,-ind]
  X <- X1[,-ind]%*%Sig # drop (m+1) start terms, model submatrix for the scop-term
  object$X <- X # the final model matrix
  object$P <- list()
  object$S <- list()
  object$Sigma <- Sig
 
  if (!object$fixed) # create the penalty matrix
  { P <- diff(diag(q-length(ind)),difference=1)
    object$P[[1]] <- P
    object$S[[1]] <- crossprod(P)
  }
  object$p.ident <- rep(TRUE,q-length(ind)) ## p.ident is an indicator of which coefficients must be positive (exponentiated)
  object$n.zero <- ind
  object$rank <- ncol(object$X)-1  # penalty rank
  object$null.space.dim <- 2 ## m+1  # dim. of unpenalized space
  object$C <- matrix(0, 0, ncol(X)) # to have no other constraints 
  object$knots <- xk;
  object$m <- m;
  object$df<-ncol(object$X)     # maximum DoF 

  ## get model matrix for 1st and 2nd derivatives of the smooth...
  h <- xk[q]-xk[q-1] ##(max(x)-min(x))/(q-m-1) ## distance between two adjacent knots (both expressions give the same distance)
  ## Xdf1, Xdf2 need to be checked...
  object$Xdf1 <- splineDesign(xk,x,ord=m+1)[,(m+2):(q-1)]/h ## ord is by one less for the 1st derivative
  object$Xdf2 <- if (m==0) matrix(0,nrow(X),ncol(object$Xdf1)-1) ## for piecewise linear splines
                 else splineDesign(xk,x,ord=m)[,2:(q-2)]/h^2 ## ord is by two less for the 2nd derivative 
  class(object)<-"miso.smooth"  # Give object a class
  object
}


## Prediction matrix for the `miso` smooth class... 
Predict.matrix.miso.smooth<-function(object,data)
## prediction method function for the `miso' smooth class
{ m <- object$m # spline order, m+1=3 default for cubic spline
  q <- object$bs.dim ## object$df +1 ## ?????
 # Sig <- matrix(as.numeric(rep(1:q,q)>=rep(1:q,each=q)),q,q) ## coef summation matrix
  Sig <- matrix(1,q,q)  ## coef summation matrix
  Sig[upper.tri(Sig)] <-0
  ind <- object$n.zero ## 1:(m+1)
  Sig[,ind] <- 0; Sig[ind,] <- 0
  ## find spline basis inner knot range...
  ll <- object$knots[m+2];ul <- object$knots[length(object$knots)-m-1]
  x <- data[[object$term]]
  n <- length(x)
  ind <- x<=ul & x>=ll ## data in range
  if (sum(ind)==n) { ## all in range
     X <- spline.des(object$knots,x,m+2)$design
     X <- X%*%Sig 
  } else { ## some extrapolation needed 
     ## matrix mapping coefs to value and slope at end points...
     D <- spline.des(object$knots,c(ll,ll,ul,ul),m+2,c(0,1,0,1))$design
     X <- matrix(0,n,ncol(D)) ## full predict matrix
     if (sum(ind)> 0)  X[ind,] <- spline.des(object$knots,x[ind],m+2)$design ## interior rows
     ## Now add rows for linear extrapolation...
     ind <- x < ll 
     if (sum(ind)>0) X[ind,] <- cbind(1,x[ind]-ll)%*%D[1:2,]
     ind <- x > ul
     if (sum(ind)>0) X[ind,] <- cbind(1,x[ind]-ul)%*%D[3:4,]
     X <- X%*%Sig 
  }
  X
}


##############################################################################################
### Adding Monotone increasing SCOP-spline construction with an 'finish at zero' constraint,
### a SCOP-spline additionally constrained to be zero at the end, at the right-end point of 
### the covariate range...
#############################################################################################

smooth.construct.mifo.smooth.spec<- function(object, data, knots)
## construction of the monotone increasing smooth with a 'finish-at-zero' constraint;
## achieved simply by setting the last (m+1) spline coefficients to zero...
{ 
  if (!is.null(object$point.con[[1]])) ## a point constraint is supplied?
    stop(paste("'mifo' smooth works only with a 'finish at zero' constraint; 'pc' argument does not work here"))
  m <- object$p.order[1]
  if (is.na(m)) m <- 2 ## default for cubic spline
  if (m<0) stop("silly m supplied")
  if (object$bs.dim<0) object$bs.dim <- 10 ## default
  q <- object$bs.dim 
  nk <- q+m+2 ## number of knots
  if (nk<=0) stop("k too small for m")
  x <- data[[object$term]]  ## the data
  xk <- knots[[object$term]] ## will be NULL if none supplied
  if (is.null(xk)) { # space knots evenly through data
       n<-length(x)
       xk<-rep(0,q+m+2)
       xk[(m+2):(q+1)]<-seq(min(x),max(x),length=q-m)
       for (i in 1:(m+1)) {xk[i]<-xk[m+2]-(m+2-i)*(xk[m+3]-xk[m+2])}
       for (i in (q+2):(q+m+2)) {xk[i]<-xk[q+1]+(i-q-1)*(xk[m+3]-xk[m+2])}
     }
  ## move m knots, xk[q-m+1],...,xk[q] (the m pre-last inner knots), to the last inner knot xk[q+1],
  ## join these m knots with the constrained xk[q+1] knot to avoid a plateau end...
  xk[(q-m+1):q] <- xk[q+1]
  if (length(xk)!=nk) # right number of knots?
   stop(paste("there should be ", nk," supplied knots"))
  ##  get model matrix...
  X1 <- splineDesign(xk,x,ord=m+2)
  ## get unconstrained matrix Sigma and remove the last m+1 columns and rows...
  Sig <- matrix(1,q,q)  ## coef summation matrix
  Sig[upper.tri(Sig)] <-0
  ind <- (q-m):q
 ## Sig[,ind] <- 0; Sig[ind,] <- 0
  Sig <- Sig[-ind,-ind]
  X <- X1[,-ind]%*%Sig # drop (m+1) end terms, model submatrix for the scop-term
  object$X <- X # the final model matrix
  object$P <- list()
  object$S <- list()
  object$Sigma <- Sig

  if (!object$fixed) # create the penalty matrix
  { P <- diff(diag(q-length(ind)),difference=1)
    object$P[[1]] <- P
    object$S[[1]] <- crossprod(P)
  }
  object$p.ident <- c(FALSE,rep(TRUE,q-length(ind)-1)) ## p.ident is an indicator of which coefficients must be positive (exponentiated)
  object$n.zero <- ind
  object$rank <- ncol(object$X)-1  # penalty rank
  object$null.space.dim <- 2 ##m+1  # dim. of unpenalized space
  object$C <- matrix(0, 0, ncol(X)) # to have no other constraints 
  object$knots <- xk;
  object$m <- m;
  object$df<-ncol(object$X)     # maximum DoF 

  ## get model matrix for 1st and 2nd derivatives of the smooth...
  h <- (max(x)-min(x))/(q-m-1) ## distance between two adjacent knots
  ## Xdf1, Xdf2 need to be checked...
  object$Xdf1 <- splineDesign(xk,x,ord=m+1)[,1:(q-1-(m+1))]/h ## ord is by one less for the 1st derivative
  object$Xdf2 <- if (m==0) matrix(0,nrow(X),ncol(object$Xdf1)-1) ## for piecewise linear splines
                 else splineDesign(xk,x,ord=m)[,2:(q-2)]/h^2 ## ord is by two less for the 2nd derivative 
  class(object)<-"mifo.smooth"  # Give object a class
  object
}


## Prediction matrix for the `mifo` smooth class... 
Predict.matrix.mifo.smooth<-function(object,data)
## prediction method function for the `mifo' smooth class
{ m <- object$m # spline order, m+1=3 default for cubic spline
  q <- object$bs.dim ## object$df +1 ## ?????
  Sig <- matrix(1,q,q)  ## coef summation matrix
  Sig[upper.tri(Sig)] <-0
  ind <- object$n.zero 
  Sig[,ind] <- 0; Sig[ind,] <- 0
  ## find spline basis inner knot range...
  ll <- object$knots[m+2];ul <- object$knots[length(object$knots)-m-1]
  x <- data[[object$term]]
  n <- length(x)
  ind <- x<=ul & x>=ll ## data in range
  if (sum(ind)==n) { ## all in range
     X <- spline.des(object$knots,x,m+2)$design
     X <- X%*%Sig 
  } else { ## some extrapolation needed 
     ## matrix mapping coefs to value and slope at end points...
     D <- spline.des(object$knots,c(ll,ll,ul,ul),m+2,c(0,1,0,1))$design
     X <- matrix(0,n,ncol(D)) ## full predict matrix
     if (sum(ind)> 0)  X[ind,] <- spline.des(object$knots,x[ind],m+2)$design ## interior rows
     ## Now add rows for linear extrapolation...
     ind <- x < ll 
     if (sum(ind)>0) X[ind,] <- cbind(1,x[ind]-ll)%*%D[1:2,]
     ind <- x > ul
     if (sum(ind)>0) X[ind,] <- cbind(1,x[ind]-ul)%*%D[3:4,]
     X <- X%*%Sig 
  }
  X
}


########################################################
### Adding Monotone decreasing SCOP-spline construction 
########################################################

smooth.construct.mpd.smooth.spec<- function(object, data, knots)
## construction of the monotone decreasing smooth
{ 
  # require(splines)
  m <- object$p.order[1]
  if (is.na(m)) m <- 2 ## default for cubic splines
  if (m<0) stop("silly m supplied")
  if (object$bs.dim<0) object$bs.dim <- 10 ## default
  q <- object$bs.dim 
  nk <- q+m+2 ## number of knots
  if (nk<=0) stop("k too small for m")
  x <- data[[object$term]]  ## the data
  xk <- knots[[object$term]] ## will be NULL if none supplied
  if (is.null(xk)) # space knots through data
     { n <- length(x)
       xk <- rep(0,q+m+2)
       xk[(m+2):(q+1)] <- seq(min(x),max(x),length=q-m)
       for (i in 1:(m+1)) {xk[i] <- xk[m+2]-(m+2-i)*(xk[m+3]-xk[m+2])}
       for (i in (q+2):(q+m+2)) {xk[i] <- xk[q+1]+(i-q-1)*(xk[m+3]-xk[m+2])}
  }
  if (length(xk)!=nk) # right number of knots?
  stop(paste("there should be ",nk," supplied knots"))
  #  get model matrix-------------
  X1 <- splineDesign(xk,x,ord=m+2)
  # get matrix Sigma and remove the first column for the monotone smooth (column and row for Sig)
 # Sig <- matrix(as.numeric(rep(1:q,q)>=rep(1:q,each=q)),q,q) ## coef summation matrix
  Sig <- matrix(-1,q,q)  ## coef summation matrix
  Sig[upper.tri(Sig)] <-0
  Sig[,1] <- -Sig[,1] ## monotone decreasing case
  X <- X1[,-1]%*%Sig[-1,-1] # drop intercept term
  
  ## applying sum-to-zero (centering) constraint...
  cmx <- colMeans(X)
  X <- sweep(X,2,cmx) ## subtract cmx from columns 
  object$X <- X # the final model matrix
  object$cmX <- cmx  

  object$Sigma <- Sig[-1,-1]
  object$P <- list()
  object$S <- list()

  if (!object$fixed) # create the penalty matrix
  { P <- diff(diag(q-1),difference=1)
    object$P[[1]] <- P
    object$S[[1]] <- crossprod(P)
  }
  object$p.ident <- rep(TRUE,q-1)  ## p.ident is an indicator of which coefficients must be positive (exponentiated) 
  object$rank <- ncol(object$X)-1  # penalty rank
  object$null.space.dim <-2 ## m+1  # dim. of unpenalized space
  object$C <- matrix(0, 0, ncol(X)) # to have no other constraints 
  
  ## get model matrix for 1st and 2nd derivatives of the smooth...
  h <- (max(x)-min(x))/(q-m-1) ## distance between two adjacent knots
  object$Xdf1 <- splineDesign(xk,x,ord=m+1)[,2:(q-1)]/h ## ord is by one less for the 1st derivative
  object$Xdf2 <- if (m==0) matrix(0,nrow(X),ncol(object$Xdf1)-1) ## for piecewise linear splines
                 else splineDesign(xk,x,ord=m)[,2:(q-2)]/h^2 ## ord is by two less for the 2nd derivative 

  object$knots <- xk;
  object$m <- m;
  object$df<-ncol(object$X)     # maximum DoF (if unconstrained)
  class(object)<-"mpd.smooth"  # Give object a class
  object
}


## Prediction matrix for the `mpd` smooth class 

Predict.matrix.mpd.smooth<-function(object,data)
## prediction method function for the `mpd' smooth class
{ m <- object$m+1; # spline order, m+1=3 default for cubic spline
  q <- object$df +1
  # elements of matrix Sigma for decreasing smooth...
 # Sig <- matrix(as.numeric(rep(1:q,q)>=rep(1:q,each=q)),q,q) ## coef summation matrix
  Sig <- matrix(-1,q,q)  ## coef summation matrix
  Sig[upper.tri(Sig)] <-0
  Sig[,1] <- -Sig[,1] ## monotone decreasing case
  ## find spline basis inner knot range...
  ll <- object$knots[m+1];ul <- object$knots[length(object$knots)-m]
  m <- m + 1
  x <- data[[object$term]]
  n <- length(x)
  ind <- x<=ul & x>=ll ## data in range
  if (sum(ind)==n) { ## all in range
     X <- spline.des(object$knots,x,m)$design
     X <- X%*%Sig 
  } else { ## some extrapolation needed 
     ## matrix mapping coefs to value and slope at end points...
     D <- spline.des(object$knots,c(ll,ll,ul,ul),m,c(0,1,0,1))$design
     X <- matrix(0,n,ncol(D)) ## full predict matrix
     if (sum(ind)> 0)  X[ind,] <- spline.des(object$knots,x[ind],m)$design ## interior rows
     ## Now add rows for linear extrapolation...
     ind <- x < ll 
     if (sum(ind)>0) X[ind,] <- cbind(1,x[ind]-ll)%*%D[1:2,]
     ind <- x > ul
     if (sum(ind)>0) X[ind,] <- cbind(1,x[ind]-ul)%*%D[3:4,]
     X <- X%*%Sig 
    ## X <- sweep(X,2,c(0,object$cmX)) 
  }
  X <- sweep(X,2,c(0,object$cmX))
  X
}

#######################################################
### Adding Monotone decreasing SCOP-spline construction without applying identifiability constraints
### to be used with numeric 'by' variable...
########################################################

smooth.construct.mpdBy.smooth.spec<- function(object, data, knots)
## construction of the monotone decreasing smooth
{ 
  # require(splines)
  m <- object$p.order[1]
  if (is.na(m)) m <- 2 ## default for cubic splines
  if (m<0) stop("silly m supplied")
  if (object$bs.dim<0) object$bs.dim <- 10 ## default
  q <- object$bs.dim 
  nk <- q+m+2 ## number of knots
  if (nk<=0) stop("k too small for m")
  x <- data[[object$term]]  ## the data
  xk <- knots[[object$term]] ## will be NULL if none supplied
  if (is.null(xk)) # space knots through data
     { n <- length(x)
       xk <- rep(0,q+m+2)
       xk[(m+2):(q+1)] <- seq(min(x),max(x),length=q-m)
       for (i in 1:(m+1)) {xk[i] <- xk[m+2]-(m+2-i)*(xk[m+3]-xk[m+2])}
       for (i in (q+2):(q+m+2)) {xk[i] <- xk[q+1]+(i-q-1)*(xk[m+3]-xk[m+2])}
  }
  if (length(xk)!=nk) # right number of knots?
  stop(paste("there should be ",nk," supplied knots"))
  #  get model matrix-------------
  X1 <- splineDesign(xk,x,ord=m+2)
  # get matrix Sigma and remove the first column for the monotone smooth (column and row for Sig)
 # Sig <- matrix(as.numeric(rep(1:q,q)>=rep(1:q,each=q)),q,q) ## coef summation matrix
  Sig <- matrix(-1,q,q)  ## coef summation matrix
  Sig[upper.tri(Sig)] <-0
  Sig[,1] <- -Sig[,1] ## monotone decrease case

  X <- X1%*%Sig # no drop intercept term
  object$X<-X # the final model matrix
  object$Sigma <- Sig
  object$P <- list()
  object$S <- list()

  if (!object$fixed) # create the penalty matrix
  { P <- diff(diag(q-1),difference=1)
    P <- rbind(rep(0,q-1),P) ## adding 1st row of zeros
    P <- cbind(rep(0,q-1),P) ## adding first column of zeros
    object$P[[1]] <- P
    object$S[[1]] <- crossprod(P)
  }
  object$p.ident <- c(FALSE,rep(TRUE,q-1))  ## p.ident is an indicator of which coefficients must be positive (exponentiated) 
  object$rank <- ncol(object$X)  # penalty rank
  object$null.space.dim <- 2 ## m+1  # dim. of unpenalized space
  object$C <- matrix(0, 0, ncol(X)) # to have no other constraints 
  
  ## get model matrix for 1st and 2nd derivatives of the smooth...
  h <- (max(x)-min(x))/(q-m-1) ## distance between two adjacent knots
  object$Xdf1 <- splineDesign(xk,x,ord=m+1)[,1:(q-1)]/h ## ord is by one less for the 1st derivative
  object$Xdf2 <- if (m==0) matrix(0,nrow(X),ncol(object$Xdf1)-1) ## for piecewise linear splines
                 else splineDesign(xk,x,ord=m)[,2:(q-2)]/h^2 ## ord is by two less for the 2nd derivative 

  object$knots <- xk;
  object$m <- m;
  object$df<-ncol(object$X)     # maximum DoF (if unconstrained)
  class(object)<-"mpdBy.smooth"  # Give object a class
  object
}


## Prediction matrix for the `mpdBy` smooth class 

Predict.matrix.mpdBy.smooth<-function(object,data)
## prediction method function for the `mpdBy' smooth class
{ m <- object$m+1; # spline order, m+1=3 default for cubic spline
  q <- object$df 
  # elements of matrix Sigma for decreasing smooth...
 # Sig <- matrix(as.numeric(rep(1:q,q)>=rep(1:q,each=q)),q,q) ## coef summation matrix
  Sig <- matrix(-1,q,q)  ## coef summation matrix
  Sig[upper.tri(Sig)] <-0
  Sig[,1] <- -Sig[,1] ## monotone decrease case
  ## find spline basis inner knot range...
  ll <- object$knots[m+1];ul <- object$knots[length(object$knots)-m]
  m <- m + 1
  x <- data[[object$term]]
  n <- length(x)
  ind <- x<=ul & x>=ll ## data in range
  if (sum(ind)==n) { ## all in range
     X <- spline.des(object$knots,x,m)$design
     X <- X%*%Sig 
  } else { ## some extrapolation needed 
     ## matrix mapping coefs to value and slope at end points...
     D <- spline.des(object$knots,c(ll,ll,ul,ul),m,c(0,1,0,1))$design
     X <- matrix(0,n,ncol(D)) ## full predict matrix
     if (sum(ind)> 0)  X[ind,] <- spline.des(object$knots,x[ind],m)$design ## interior rows
     ## Now add rows for linear extrapolation...
     ind <- x < ll 
     if (sum(ind)>0) X[ind,] <- cbind(1,x[ind]-ll)%*%D[1:2,]
     ind <- x > ul
     if (sum(ind)>0) X[ind,] <- cbind(1,x[ind]-ul)%*%D[3:4,]
     X <- X%*%Sig 
  }
  X
}


##############################################################
### Smooth constructor for the mixed constrainted smooths ......
##############################################################

###############################################################
### Adding Monotone decreasing & concave P-spline construction 
###############################################################

smooth.construct.mdcv.smooth.spec<- function(object, data, knots)
## construction of the monotone decreasing and concave smooth
{ 
  # require(splines)
  m <- object$p.order[1]
  if (is.na(m)) m <- 2 ## default for cubis spline
  if (m<0) stop("silly m supplied")
  if (object$bs.dim<0) object$bs.dim <- 10 ## default
  q <- object$bs.dim 
  nk <- q+m+2 ## number of knots
  if (nk<=0) stop("k too small for m")
  x <- data[[object$term]]  ## the data
  xk <- knots[[object$term]] ## will be NULL if none supplied
  if (is.null(xk)) # space knots through data
     { n<-length(x)
       xk<-rep(0,q+m+2)
       xk[(m+2):(q+1)]<-seq(min(x),max(x),length=q-m)
       for (i in 1:(m+1)) {xk[i]<-xk[m+2]-(m+2-i)*(xk[m+3]-xk[m+2])}
       for (i in (q+2):(q+m+2)) {xk[i]<-xk[q+1]+(i-q-1)*(xk[m+3]-xk[m+2])}
     }
  if (length(xk)!=nk) # right number of knots?
  stop(paste("there should be ",nk," supplied knots"))
  #  get model matrix-------------
  X1 <- splineDesign(xk,x,ord=m+2)
  # use matrix Sigma and remove the first column for the monotone smooth
  Sig <- matrix(0,(q-1),(q-1))   # Define Matrix Sigma
     # for monotone decreasing & concave smooth
  for (i in 1:(q-1)) Sig[i:(q-1),i]<--c(1:(q-i))
  X <- X1[,2:q]%*%Sig # model submatrix for the constrained term

  ## applying sum-to-zero (centering) constraint...
  cmx <- colMeans(X)
  X <- sweep(X,2,cmx) ## subtract cmx from columns 
  object$X <- X # the final model matrix
  object$cmX <- cmx

  object$Sigma <- Sig
  object$P <- list()
  object$S <- list()

  if (!object$fixed) # create the penalty matrix
  { P <- diff(diag(q-2),difference=1)
    object$P[[1]] <- P
    S <- matrix(0,q-1,q-1)
    S[2:(q-1),2:(q-1)] <- crossprod(P)
    object$S[[1]] <- S
  }
  object$p.ident <- rep(TRUE,q-1)  ## p.ident is an indicator of which coefficients must be positive (exponentiated)  
  object$rank <- ncol(object$X)-1  # penalty rank
  object$null.space.dim <- 2 ##m+1  # dim. of unpenalized space
  object$C <- matrix(0, 0, ncol(X)) # to have no other constraints 
  
  ## get model matrix for 1st and 2nd derivatives of the smooth...
  h <- (max(x)-min(x))/(q-m-1) ## distance between two adjacent knots
  object$Xdf1 <- splineDesign(xk,x,ord=m+1)[,2:(q-1)]/h ## ord is by one less for the 1st derivative
  object$Xdf2 <- if (m==0) matrix(0,nrow(X),ncol(object$Xdf1)-1) ## for piecewise linear splines
                 else splineDesign(xk,x,ord=m)[,2:(q-2)]/h^2 ## ord is by two less for the 2nd derivative 

  object$knots <- xk;
  object$m <- m;
  object$df<-ncol(object$X)     # maximum DoF (if unconstrained)
  class(object)<-"mdcv.smooth"  # Give object a class
  object
}


## Prediction matrix for the `mdcv` smooth class.... 

Predict.matrix.mdcv.smooth<-function(object,data)
## prediction method function for the `mdcv' smooth class
{ m <- object$m+1; # spline order, m+1=3 default for cubic spline
  q <- object$df +1
  Sig <- matrix(0,q,q)   
  Sig[,1] <- rep(1,q)
  for (i in 2:q)  Sig[i:q,i] <- -c(1:(q-i+1))
  ## find spline basis inner knot range...
  ll <- object$knots[m+1];ul <- object$knots[length(object$knots)-m]
  m <- m + 1
  x <- data[[object$term]]
  n <- length(x)
  ind <- x<=ul & x>=ll ## data in range
  if (sum(ind)==n) { ## all in range
     X <- spline.des(object$knots,x,m)$design
     X <- X%*%Sig 
  } else { ## some extrapolation needed 
     ## matrix mapping coefs to value and slope at end points...
     D <- spline.des(object$knots,c(ll,ll,ul,ul),m,c(0,1,0,1))$design
     X <- matrix(0,n,ncol(D)) ## full predict matrix
     if (sum(ind)> 0)  X[ind,] <- spline.des(object$knots,x[ind],m)$design ## interior rows
     ## Now add rows for linear extrapolation...
     ind <- x < ll 
     if (sum(ind)>0) X[ind,] <- cbind(1,x[ind]-ll)%*%D[1:2,]
     ind <- x > ul
     if (sum(ind)>0) X[ind,] <- cbind(1,x[ind]-ul)%*%D[3:4,]
     X <- X%*%Sig 
    ## X <- sweep(X,2,c(0,object$cmX))
  }
  X <- sweep(X,2,c(0,object$cmX))
  X
}


###########################################################
### Adding decreasing & concave SCOP-spline construction without identifiability constraints
### to be used with numeric 'by' variable and linear functional terms...
##########################################################

smooth.construct.mdcvBy.smooth.spec<- function(object, data, knots)
## construction of the monotone decreasing and concave smooth
{ 
  m <- object$p.order[1]
  if (is.na(m)) m <- 2 ## default for cubis spline
  if (m<0) stop("silly m supplied")
  if (object$bs.dim<0) object$bs.dim <- 10 ## default
  q <- object$bs.dim 
  nk <- q+m+2 ## number of knots
  if (nk<=0) stop("k too small for m")
  x <- data[[object$term]]  ## the data
  xk <- knots[[object$term]] ## will be NULL if none supplied
  if (is.null(xk)) # space knots through data
     { n<-length(x)
       xk<-rep(0,q+m+2)
       xk[(m+2):(q+1)]<-seq(min(x),max(x),length=q-m)
       for (i in 1:(m+1)) {xk[i]<-xk[m+2]-(m+2-i)*(xk[m+3]-xk[m+2])}
       for (i in (q+2):(q+m+2)) {xk[i]<-xk[q+1]+(i-q-1)*(xk[m+3]-xk[m+2])}
     }
  if (length(xk)!=nk) # right number of knots?
  stop(paste("there should be ",nk," supplied knots"))
  #  get model matrix-------------
  X <- splineDesign(xk,x,ord=m+2)
  # Define matrix Sigma for monotone decreasing & concave smooth
  Sig <- matrix(0,q,q)   
  Sig[,1] <- rep(1,q)
  for (i in 2:q)  Sig[i:q,i] <- -c(1:(q-i+1))

  X <- X%*%Sig # model matrix for the constrained term
  object$X <- X 
  object$Sigma <- Sig
  object$P <- list()
  object$S <- list()

  if (!object$fixed) # create the penalty matrix
  { P <- diff(diag(q-2),difference=1)
    P <- rbind(matrix(0,2,q-2),P) ## adding first two rows of zeros
    P <- cbind(matrix(0,q-1,2),P) ## adding first two columns of zeros
    object$P[[1]] <- P
    object$S[[1]] <- crossprod(P)
  }
  object$p.ident <- c(FALSE,rep(TRUE,q-1))  ## p.ident is an indicator of which coefficients must be positive (exponentiated)  
  object$rank <- ncol(object$X)  # penalty rank
  object$null.space.dim <- 2 ## m+1  # dim. of unpenalized space
  object$C <- matrix(0, 0, ncol(X)) # to have no other constraints 
  
  ## get model matrix for 1st and 2nd derivatives of the smooth...
  h <- (max(x)-min(x))/(q-m-1) ## distance between two adjacent knots
  object$Xdf1 <- splineDesign(xk,x,ord=m+1)[,1:(q-1)]/h ## ord is by one less for the 1st derivative
  object$Xdf2 <- if (m==0) matrix(0,nrow(X),ncol(object$Xdf1)-1) ## for piecewise linear splines
                 else splineDesign(xk,x,ord=m)[,2:(q-2)]/h^2 ## ord is by two less for the 2nd derivative 

  object$knots <- xk;
  object$m <- m;
  object$df<-ncol(object$X)     # maximum DoF (if unconstrained)
  class(object)<-"mdcvBy.smooth"  # Give object a class
  object
}


## Prediction matrix for the `mdcvBy` smooth class.... 

Predict.matrix.mdcvBy.smooth<-function(object,data)
## prediction method function for the `mdcvBy' smooth class
{ m <- object$m+1; # spline order, m+1=3 default for cubic spline
  q <- object$df 
  Sig <- matrix(0,q,q)   
  Sig[,1] <- rep(1,q)
  for (i in 2:q)  Sig[i:q,i] <- -c(1:(q-i+1))
  ## find spline basis inner knot range...
  ll <- object$knots[m+1];ul <- object$knots[length(object$knots)-m]
  m <- m + 1
  x <- data[[object$term]]
  n <- length(x)
  ind <- x<=ul & x>=ll ## data in range
  if (sum(ind)==n) { ## all in range
     X <- spline.des(object$knots,x,m)$design
     X <- X%*%Sig 
  } else { ## some extrapolation needed 
     ## matrix mapping coefs to value and slope at end points...
     D <- spline.des(object$knots,c(ll,ll,ul,ul),m,c(0,1,0,1))$design
     X <- matrix(0,n,ncol(D)) ## full predict matrix
     if (sum(ind)> 0)  X[ind,] <- spline.des(object$knots,x[ind],m)$design ## interior rows
     ## Now add rows for linear extrapolation...
     ind <- x < ll 
     if (sum(ind)>0) X[ind,] <- cbind(1,x[ind]-ll)%*%D[1:2,]
     ind <- x > ul
     if (sum(ind)>0) X[ind,] <- cbind(1,x[ind]-ul)%*%D[3:4,]
     X <- X%*%Sig 
  }
  X
}


###############################################################
### Adding decreasing & convex SCOP-spline construction 
################################################################

smooth.construct.mdcx.smooth.spec<- function(object, data, knots)
##  the constructor for the monotone decreasing and convex smooth
{ 
  #require(splines)
  m <- object$p.order[1]
  if (is.na(m)) m <- 2 ## default 
  if (m<0) stop("silly m supplied")
  if (object$bs.dim<0) object$bs.dim <- 10 ## default
  q <- object$bs.dim 
  nk <- q+m+2 ## number of knots
  if (nk<=0) stop("k too small for m")
  x <- data[[object$term]]  ## the data
  xk <- knots[[object$term]] ## will be NULL if none supplied
  if (is.null(xk)) # space knots through data
  { n<-length(x)
    xk<-rep(0,q+m+2)
    xk[(m+2):(q+1)]<-seq(min(x),max(x),length=q-m)
    for (i in 1:(m+1)) {xk[i]<-xk[m+2]-(m+2-i)*(xk[m+3]-xk[m+2])}
    for (i in (q+2):(q+m+2)) {xk[i]<-xk[q+1]+(i-q-1)*(xk[m+3]-xk[m+2])}
  }
  if (length(xk)!=nk) # right number of knots?
  stop(paste("there should be ",nk," supplied knots"))
  #  get model matrix-------------
  X1 <- splineDesign(xk,x,ord=m+2)
  # use matrix Sigma and remove the first column for the monotone smooth
  Sig <- matrix(0,(q-1),(q-1))   # Define Matrix Sigma
     # for monotone decreasing & convex smooth
  Sig[1,] <- -rep(1,q-1)
  for (i in 2:(q-1)) {
         Sig[i,1:(q-i)] <- -i;
         Sig[i,(q-i+1):(q-1)] <- -c((i-1):1)
  }
  X <- X1[,2:q]%*%Sig # model submatrix for the constrained term

  ## applying sum-to-zero (centering) constraint...
  cmx <- colMeans(X)
  X <- sweep(X,2,cmx) ## subtract cmx from columns 
  object$X <- X # the final model matrix
  object$cmX <- cmx

  object$Sigma <- Sig
  object$P <- list()
  object$S <- list()

  if (!object$fixed) # create the penalty matrix
  { P <- diff(diag(q-2),difference=1)
    object$P[[1]] <- P
    S <- matrix(0,q-1,q-1)
    S[2:(q-1),2:(q-1)] <- crossprod(P)
    object$S[[1]] <- S
  }
  object$p.ident <- rep(TRUE,q-1)  ## p.ident is an indicator of which coefficients must be positive (exponentiated)  
  object$rank <- ncol(object$X)-1  # penalty rank
  object$null.space.dim <- 2 ## m+1  # dim. of unpenalized space
  object$C <- matrix(0, 0, ncol(X)) # to have no other constraints 
  
  ## get model matrix for 1st and 2nd derivatives of the smooth...
  h <- (max(x)-min(x))/(q-m-1) ## distance between two adjacent knots
  object$Xdf1 <- splineDesign(xk,x,ord=m+1)[,2:(q-1)]/h ## ord is by one less for the 1st derivative
  object$Xdf2 <- if (m==0) matrix(0,nrow(X),ncol(object$Xdf1)-1) ## for piecewise linear splines
                 else splineDesign(xk,x,ord=m)[,2:(q-2)]/h^2 ## ord is by two less for the 2nd derivative 

  object$knots <- xk;
  object$m <- m;
  object$df<-ncol(object$X)     # maximum DoF (if unconstrained)
  class(object)<-"mdcx.smooth"  # Give object a class
  object
}


## Prediction matrix for the `mdcx` smooth class.... 
Predict.matrix.mdcx.smooth<-function(object,data)
## the prediction method for the `mdcx' smooth class
{ x <- data[[object$term]]
  m <- object$m+1; # spline order, m+1=3 default for cubic spline
  q <- object$df +1
  Sig <- matrix(0,q,q)   
  Sig[,1] <- rep(1,q)
  Sig1 <- matrix(0,(q-1),(q-1))   
  Sig1[1,] <- -rep(1,q-1)
  for (i in 2:(q-1)) {
         Sig1[i,1:(q-i)]<--i;
         Sig1[i,(q-i+1):(q-1)]<--c((i-1):1)
  }
  Sig [2:q,2:q] <- Sig1

  ## find spline basis inner knot range...
  ll <- object$knots[m+1];ul <- object$knots[length(object$knots)-m]
  m <- m + 1
  n <- length(x)
  ind <- x<=ul & x>=ll ## data in range
  if (sum(ind)==n) { ## all in range
     X <- spline.des(object$knots,x,m)$design
     X <- X%*%Sig 
  } else { ## some extrapolation needed 
     ## matrix mapping coefs to value and slope at end points...
     D <- spline.des(object$knots,c(ll,ll,ul,ul),m,c(0,1,0,1))$design
     X <- matrix(0,n,ncol(D)) ## full predict matrix
     if (sum(ind)> 0)  X[ind,] <- spline.des(object$knots,x[ind],m)$design ## interior rows
     ## Now add rows for linear extrapolation...
     ind <- x < ll 
     if (sum(ind)>0) X[ind,] <- cbind(1,x[ind]-ll)%*%D[1:2,]
     ind <- x > ul
     if (sum(ind)>0) X[ind,] <- cbind(1,x[ind]-ul)%*%D[3:4,]
     X <- X%*%Sig 
   ##  X <- sweep(X,2,c(0,object$cmX))
  }
  X <- sweep(X,2,c(0,object$cmX))
  X
}

###########################################################
### Adding decreasing & convex SCOP-spline construction without identifiability constraints
### to be used with numeric 'by' variable and linear functional terms...
##########################################################

smooth.construct.mdcxBy.smooth.spec<- function(object, data, knots)
##  the constructor for the monotone decreasing and convex smooth
{ 
  m <- object$p.order[1]
  if (is.na(m)) m <- 2 ## default 
  if (m<0) stop("silly m supplied")
  if (object$bs.dim<0) object$bs.dim <- 10 ## default
  q <- object$bs.dim 
  nk <- q+m+2 ## number of knots
  if (nk<=0) stop("k too small for m")
  x <- data[[object$term]]  ## the data
  xk <- knots[[object$term]] ## will be NULL if none supplied
  if (is.null(xk)) # space knots through data
  { n<-length(x)
    xk<-rep(0,q+m+2)
    xk[(m+2):(q+1)]<-seq(min(x),max(x),length=q-m)
    for (i in 1:(m+1)) {xk[i]<-xk[m+2]-(m+2-i)*(xk[m+3]-xk[m+2])}
    for (i in (q+2):(q+m+2)) {xk[i]<-xk[q+1]+(i-q-1)*(xk[m+3]-xk[m+2])}
  }
  if (length(xk)!=nk) # right number of knots?
  stop(paste("there should be ",nk," supplied knots"))
  #  get model matrix-------------
  X <- splineDesign(xk,x,ord=m+2)
  # define matrix Sigma for monotone decreasing & convex smooth
  Sig <- matrix(0,q,q)   
  Sig[,1] <- rep(1,q)
  Sig1 <- matrix(0,(q-1),(q-1))   
  Sig1[1,] <- -rep(1,q-1)
  for (i in 2:(q-1)) {
         Sig1[i,1:(q-i)]<--i;
         Sig1[i,(q-i+1):(q-1)]<--c((i-1):1)
  }
  Sig [2:q,2:q] <- Sig1

  X <- X%*%Sig # model matrix for the constrained term
  object$X<-X 
  object$Sigma <- Sig
  object$P <- list()
  object$S <- list()

  if (!object$fixed) # create the penalty matrix
  { P <- diff(diag(q-2),difference=1)
    P <- rbind(matrix(0,2,q-2),P) ## adding first two rows of zeros
    P <- cbind(matrix(0,q-1,2),P) ## adding first two columns of zeros
    object$P[[1]] <- P
    object$S[[1]] <- crossprod(P)
  }
  object$p.ident <- c(FALSE,rep(TRUE,q-1))  ## p.ident is an indicator of which coefficients must be positive (exponentiated)  
  object$rank <- ncol(object$X)  # penalty rank
  object$null.space.dim <-2 ## m+1  # dim. of unpenalized space
  object$C <- matrix(0, 0, ncol(X)) # to have no other constraints 
  
  ## get model matrix for 1st and 2nd derivatives of the smooth...
  h <- (max(x)-min(x))/(q-m-1) ## distance between two adjacent knots
  object$Xdf1 <- splineDesign(xk,x,ord=m+1)[,1:(q-1)]/h ## ord is by one less for the 1st derivative
  object$Xdf2 <- if (m==0) matrix(0,nrow(X),ncol(object$Xdf1)-1) ## for piecewise linear splines
                 else splineDesign(xk,x,ord=m)[,2:(q-2)]/h^2 ## ord is by two less for the 2nd derivative 

  object$knots <- xk;
  object$m <- m;
  object$df<-ncol(object$X)     # maximum DoF (if unconstrained)
  class(object)<-"mdcxBy.smooth"  # Give object a class
  object
}


## Prediction matrix for the `mdcxBy` smooth class.... 
Predict.matrix.mdcxBy.smooth<-function(object,data)
## the prediction method for the `mdcxBy' smooth class
{ x <- data[[object$term]]
  m <- object$m+1; # spline order, m+1=3 default for cubic spline
  q <- object$df
  Sig <- matrix(0,q,q)   
  Sig[,1] <- rep(1,q)
  Sig1 <- matrix(0,(q-1),(q-1))   
  Sig1[1,] <- -rep(1,q-1)
  for (i in 2:(q-1)) {
         Sig1[i,1:(q-i)]<--i;
         Sig1[i,(q-i+1):(q-1)]<--c((i-1):1)
  }
  Sig [2:q,2:q] <- Sig1

  ## find spline basis inner knot range...
  ll <- object$knots[m+1];ul <- object$knots[length(object$knots)-m]
  m <- m + 1
  n <- length(x)
  ind <- x<=ul & x>=ll ## data in range
  if (sum(ind)==n) { ## all in range
     X <- spline.des(object$knots,x,m)$design
     X <- X%*%Sig 
  } else { ## some extrapolation needed 
     ## matrix mapping coefs to value and slope at end points...
     D <- spline.des(object$knots,c(ll,ll,ul,ul),m,c(0,1,0,1))$design
     X <- matrix(0,n,ncol(D)) ## full predict matrix
     if (sum(ind)> 0)  X[ind,] <- spline.des(object$knots,x[ind],m)$design ## interior rows
     ## Now add rows for linear extrapolation...
     ind <- x < ll 
     if (sum(ind)>0) X[ind,] <- cbind(1,x[ind]-ll)%*%D[1:2,]
     ind <- x > ul
     if (sum(ind)>0) X[ind,] <- cbind(1,x[ind]-ul)%*%D[3:4,]
     X <- X%*%Sig 
  }
  X
}



################################################################
### Adding monotone increasing & concave SCOP-spline construction ################################################################

smooth.construct.micv.smooth.spec<- function(object, data, knots)
## construction of the monotone increasing and concave smooth
{ 
 # require(splines)
  m <- object$p.order[1]
  if (is.na(m)) m <- 2 ## default 
  if (m<0) stop("silly m supplied")
  if (object$bs.dim<0) object$bs.dim <- 10 ## default
  q <- object$bs.dim 
  nk <- q+m+2 ## number of knots
  if (nk<=0) stop("k too small for m")
  x <- data[[object$term]]  ## the data
  xk <- knots[[object$term]] ## will be NULL if none supplied
  if (is.null(xk)) # space knots through data
  { n<-length(x)
    xk<-rep(0,q+m+2)
    xk[(m+2):(q+1)]<-seq(min(x),max(x),length=q-m)
    for (i in 1:(m+1)) {xk[i]<-xk[m+2]-(m+2-i)*(xk[m+3]-xk[m+2])}
    for (i in (q+2):(q+m+2)) {xk[i]<-xk[q+1]+(i-q-1)*(xk[m+3]-xk[m+2])}
  }
  if (length(xk)!=nk) # right number of knots?
  stop(paste("there should be ",nk," supplied knots"))
  #  get model matrix-------------
  X1 <- splineDesign(xk,x,ord=m+2)
  # use matrix Sigma and remove the first column for the monotone smooth
  Sig <- matrix(0,(q-1),(q-1))   # Define Matrix Sigma
     # for monotone increasing & concave smooth
  Sig[1,]<-rep(1,q-1)
  for (i in 2:(q-1)) {
       Sig[i,1:(q-i)]<-i;
       Sig[i,(q-i+1):(q-1)]<-c((i-1):1)
  }
  X <- X1[,2:q]%*%Sig # model submatrix for the constrained term

  ## applying sum-to-zero (centering) constraint...
  cmx <- colMeans(X)
  X <- sweep(X,2,cmx) ## subtract cmx from columns 
  object$X <- X # the final model matrix
  object$cmX <- cmx

  object$Sigma <- Sig
  object$P <- list()
  object$S <- list()

  if (!object$fixed) # create the penalty matrix
  { P <- diff(diag(q-2),difference=1)
    object$P[[1]] <- P
    S <- matrix(0,q-1,q-1)
    S[2:(q-1),2:(q-1)] <- crossprod(P)
    object$S[[1]] <- S
  }
  object$p.ident <- rep(TRUE,q-1)  ## p.ident is an indicator of which coefficients must be positive (exponentiated) 
  object$rank <- ncol(object$X)-1  # penalty rank
  object$null.space.dim <- 2 ##m+1  # dim. of unpenalized space
  object$C <- matrix(0, 0, ncol(X)) # to have no other constraints 
  
  ## get model matrix for 1st and 2nd derivatives of the smooth...
  h <- (max(x)-min(x))/(q-m-1) ## distance between two adjacent knots
  object$Xdf1 <- splineDesign(xk,x,ord=m+1)[,2:(q-1)]/h ## ord is by one less for the 1st derivative
  object$Xdf2 <- if (m==0) matrix(0,nrow(X),ncol(object$Xdf1)-1) ## for piecewise linear splines
                 else splineDesign(xk,x,ord=m)[,2:(q-2)]/h^2 ## ord is by two less for the 2nd derivative 

  object$knots <- xk;
  object$m <- m;
  object$df<-ncol(object$X)     # maximum DoF (if unconstrained)
  class(object)<-"micv.smooth"  # Give object a class
  object
}


## Prediction matrix for the `micv` smooth class....

Predict.matrix.micv.smooth<-function(object,data)
## prediction method function for the `micv' smooth class
{ m <- object$m+1; # spline order, m+1=3 default for cubic spline
  q <- object$df +1
  Sig <- matrix(0,q,q)   
  Sig[,1] <- rep(1,q)
  Sig1 <- matrix(0,(q-1),(q-1))   
  Sig1[1,] <- rep(1,q-1)
  for (i in 2:(q-1)) {
       Sig1[i,1:(q-i)] <- i;
       Sig1[i,(q-i+1):(q-1)] <- c((i-1):1)
  }
  Sig [2:q,2:q] <- Sig1

  ## find spline basis inner knot range...
  ll <- object$knots[m+1];ul <- object$knots[length(object$knots)-m]
  m <- m + 1
  x <- data[[object$term]]
  n <- length(x)
  ind <- x<=ul & x>=ll ## data in range
  if (sum(ind)==n) { ## all in range
     X <- spline.des(object$knots,x,m)$design
     X <- X%*%Sig 
  } else { ## some extrapolation needed 
     ## matrix mapping coefs to value and slope at end points...
     D <- spline.des(object$knots,c(ll,ll,ul,ul),m,c(0,1,0,1))$design
     X <- matrix(0,n,ncol(D)) ## full predict matrix
     if (sum(ind)> 0)  X[ind,] <- spline.des(object$knots,x[ind],m)$design ## interior rows
     ## Now add rows for linear extrapolation...
     ind <- x < ll 
     if (sum(ind)>0) X[ind,] <- cbind(1,x[ind]-ll)%*%D[1:2,]
     ind <- x > ul
     if (sum(ind)>0) X[ind,] <- cbind(1,x[ind]-ul)%*%D[3:4,]
     X <- X%*%Sig 
    ## X <- sweep(X,2,c(0,object$cmX))
  }
  X <- sweep(X,2,c(0,object$cmX))
  X
}


###########################################################
### Adding increasing & concave SCOP-spline construction without identifiability constraints
### to be used with numeric 'by' variable and linear functional terms...
##########################################################

smooth.construct.micvBy.smooth.spec<- function(object, data, knots)
## construction of the monotone increasing and concave smooth
{ 
  m <- object$p.order[1]
  if (is.na(m)) m <- 2 ## default 
  if (m<0) stop("silly m supplied")
  if (object$bs.dim<0) object$bs.dim <- 10 ## default
  q <- object$bs.dim 
  nk <- q+m+2 ## number of knots
  if (nk<=0) stop("k too small for m")
  x <- data[[object$term]]  ## the data
  xk <- knots[[object$term]] ## will be NULL if none supplied
  if (is.null(xk)) # space knots through data
  { n<-length(x)
    xk<-rep(0,q+m+2)
    xk[(m+2):(q+1)]<-seq(min(x),max(x),length=q-m)
    for (i in 1:(m+1)) {xk[i]<-xk[m+2]-(m+2-i)*(xk[m+3]-xk[m+2])}
    for (i in (q+2):(q+m+2)) {xk[i]<-xk[q+1]+(i-q-1)*(xk[m+3]-xk[m+2])}
  }
  if (length(xk)!=nk) # right number of knots?
  stop(paste("there should be ",nk," supplied knots"))
  ##  get model matrix-------------
  X <- splineDesign(xk,x,ord=m+2)
  ## define matrix Sigma for monotone increasing & concave smooth
  Sig <- matrix(0,q,q)   
  Sig[,1] <- rep(1,q)
  Sig1 <- matrix(0,(q-1),(q-1))   
  Sig1[1,] <- rep(1,q-1)
  for (i in 2:(q-1)) {
       Sig1[i,1:(q-i)] <- i;
       Sig1[i,(q-i+1):(q-1)] <- c((i-1):1)
  }
  Sig [2:q,2:q] <- Sig1

  X <- X%*%Sig ## model matrix for the constrained term
  object$X <- X 
  object$Sigma <- Sig
  object$P <- list()
  object$S <- list()

  if (!object$fixed) # create the penalty matrix
  { P <- diff(diag(q-2),difference=1)
    P <- rbind(matrix(0,2,q-2),P) ## adding first two rows of zeros
    P <- cbind(matrix(0,q-1,2),P) ## adding first two columns of zeros
    object$P[[1]] <- P
    object$S[[1]] <- crossprod(P)
  }
  object$p.ident <- c(FALSE,rep(TRUE,q-1))  ## p.ident is an indicator of which coefficients must be positive (exponentiated) 
  object$rank <- ncol(object$X)  # penalty rank
  object$null.space.dim <- 2 ##m+1  # dim. of unpenalized space
  object$C <- matrix(0, 0, ncol(X)) # to have no other constraints 
  
  ## get model matrix for 1st and 2nd derivatives of the smooth...
  h <- (max(x)-min(x))/(q-m-1) ## distance between two adjacent knots
  object$Xdf1 <- splineDesign(xk,x,ord=m+1)[,1:(q-1)]/h ## ord is by one less for the 1st derivative
  object$Xdf2 <- if (m==0) matrix(0,nrow(X),ncol(object$Xdf1)-1) ## for piecewise linear splines
                 else splineDesign(xk,x,ord=m)[,2:(q-2)]/h^2 ## ord is by two less for the 2nd derivative 

  object$knots <- xk;
  object$m <- m;
  object$df<-ncol(object$X)     # maximum DoF (if unconstrained)
  class(object)<-"micvBy.smooth"  # Give object a class
  object
}


## Prediction matrix for the `micvBy` smooth class....
Predict.matrix.micvBy.smooth<-function(object,data)
 ## prediction method function for the `micvBy' smooth class
{ m <- object$m+1; # spline order, m+1=3 default for cubic spline
  q <- object$df
  Sig <- matrix(0,q,q)   
  Sig[,1] <- rep(1,q)
  Sig1 <- matrix(0,(q-1),(q-1))   
  Sig1[1,] <- rep(1,q-1)
  for (i in 2:(q-1)) {
       Sig1[i,1:(q-i)] <- i;
       Sig1[i,(q-i+1):(q-1)] <- c((i-1):1)
  }
  Sig [2:q,2:q] <- Sig1

  ## find spline basis inner knot range...
  ll <- object$knots[m+1];ul <- object$knots[length(object$knots)-m]
  m <- m + 1
  x <- data[[object$term]]
  n <- length(x)
  ind <- x<=ul & x>=ll ## data in range
  if (sum(ind)==n) { ## all in range
     X <- spline.des(object$knots,x,m)$design
     X <- X%*%Sig 
  } else { ## some extrapolation needed 
     ## matrix mapping coefs to value and slope at end points...
     D <- spline.des(object$knots,c(ll,ll,ul,ul),m,c(0,1,0,1))$design
     X <- matrix(0,n,ncol(D)) ## full predict matrix
     if (sum(ind)> 0)  X[ind,] <- spline.des(object$knots,x[ind],m)$design ## interior rows
     ## Now add rows for linear extrapolation...
     ind <- x < ll 
     if (sum(ind)>0) X[ind,] <- cbind(1,x[ind]-ll)%*%D[1:2,]
     ind <- x > ul
     if (sum(ind)>0) X[ind,] <- cbind(1,x[ind]-ul)%*%D[3:4,]
     X <- X%*%Sig 
  }
  X
}


###############################################################
### Adding monotone increasing & convex SCOP-spline construction 
################################################################

smooth.construct.micx.smooth.spec<- function(object, data, knots)
## construction of the monotone increasing and convex smooth
{ 
  # require(splines)
  m <- object$p.order[1]
  if (is.na(m)) m <- 2 ## default 
  if (m < 0) stop("silly m supplied")
  if (object$bs.dim<0) object$bs.dim <- 10 ## default
  q <- object$bs.dim 
  nk <- q+m+2 ## number of knots
  if (nk<=0) stop("k too small for m")
  x <- data[[object$term]]  ## the data
  xk <- knots[[object$term]] ## will be NULL if none supplied
  if (is.null(xk)) # space knots through data
  { n<-length(x)
    xk<-rep(0,q+m+2)
    xk[(m+2):(q+1)]<-seq(min(x),max(x),length=q-m)
    for (i in 1:(m+1)) {xk[i]<-xk[m+2]-(m+2-i)*(xk[m+3]-xk[m+2])}
    for (i in (q+2):(q+m+2)) {xk[i]<-xk[q+1]+(i-q-1)*(xk[m+3]-xk[m+2])}
  }
  if (length(xk)!=nk) # right number of knots?
  stop(paste("there should be ",nk," supplied knots"))
  #  get model matrix-------------
  X1 <- splineDesign(xk,x,ord=m+2)
  # use matrix Sigma and remove the first column for the monotone smooth
  Sig <- matrix(0,(q-1),(q-1))   # Define Matrix Sigma
     # for monotone increasing & convex smooth
  for (i in 1:(q-1)) Sig[i:(q-1),i]<-c(1:(q-i))
  X <- X1[,2:q]%*%Sig # model submatrix for the constrained term

  ## applying sum-to-zero (centering) constraint...
  cmx <- colMeans(X)
  X <- sweep(X,2,cmx) ## subtract cmx from columns 
  object$X <- X # the final model matrix
  object$cmX <- cmx

  object$Sigma <- Sig
  object$P <- list()
  object$S <- list()

  if (!object$fixed) # create the penalty matrix
  { P <- diff(diag(q-2),difference=1)
    object$P[[1]] <- P
    S <- matrix(0,q-1,q-1)
    S[2:(q-1),2:(q-1)] <- crossprod(P)
    object$S[[1]] <- S
  }
  object$p.ident <- rep(TRUE,q-1) ## p.ident is an indicator of which coefficients must be positive (exponentiated) 
  object$rank <- ncol(object$X)-1  # penalty rank
  object$null.space.dim <- 2 ##m+1  # dim. of unpenalized space
  object$C <- matrix(0, 0, ncol(X)) # to have no other constraints 
 
  ## get model matrix for 1st and 2nd derivatives of the smooth...
  h <- (max(x)-min(x))/(q-m-1) ## distance between two adjacent knots
  object$Xdf1 <- splineDesign(xk,x,ord=m+1)[,2:(q-1)]/h ## ord is by one less for the 1st derivative
  object$Xdf2 <- if (m==0) matrix(0,nrow(X),ncol(object$Xdf1)-1) ## for piecewise linear splines
                 else splineDesign(xk,x,ord=m)[,2:(q-2)]/h^2 ## ord is by two less for the 2nd derivative 

  object$knots <- xk;
  object$m <- m;
  object$df<-ncol(object$X)     # maximum DoF (if unconstrained)
  class(object)<-"micx.smooth"  # Give object a class
  object
}

## Prediction matrix for the `micx` smooth class... 

Predict.matrix.micx.smooth<-function(object,data)
## prediction method function for the `micx` smooth class
{ m <- object$m+1; # spline order, m+1=3 default for cubic spline
  q <- object$df +1
  Sig <- matrix(0,q,q)   
  Sig[,1] <- rep(1,q)
  for (i in 2:q) Sig[i:q,i] <- c(1:(q-i+1))
  
  ## find spline basis inner knot range...
  ll <- object$knots[m+1];ul <- object$knots[length(object$knots)-m]
  m <- m + 1
  x <- data[[object$term]]
  n <- length(x)
  ind <- x<=ul & x>=ll ## data in range
  if (sum(ind)==n) { ## all in range
     X <- spline.des(object$knots,x,m)$design
     X <- X%*%Sig 
  } else { ## some extrapolation needed 
     ## matrix mapping coefs to value and slope at end points...
     D <- spline.des(object$knots,c(ll,ll,ul,ul),m,c(0,1,0,1))$design
     X <- matrix(0,n,ncol(D)) ## full predict matrix
     if (sum(ind)> 0)  X[ind,] <- spline.des(object$knots,x[ind],m)$design ## interior rows
     ## Now add rows for linear extrapolation...
     ind <- x < ll 
     if (sum(ind)>0) X[ind,] <- cbind(1,x[ind]-ll)%*%D[1:2,]
     ind <- x > ul
     if (sum(ind)>0) X[ind,] <- cbind(1,x[ind]-ul)%*%D[3:4,]
     X <- X%*%Sig 
    ## X <- sweep(X,2,c(0,object$cmX))
  }
  X <- sweep(X,2,c(0,object$cmX))
  X
}

###########################################################
### Adding increasing & convex SCOP-spline construction without identifiability constraints
### to be used with numeric 'by' variable and linear functional terms...
##########################################################

smooth.construct.micxBy.smooth.spec<- function(object, data, knots)
## construction of the monotone increasing and convex smooth
{ 
  m <- object$p.order[1]
  if (is.na(m)) m <- 2 ## default 
  if (m < 0) stop("silly m supplied")
  if (object$bs.dim<0) object$bs.dim <- 10 ## default
  q <- object$bs.dim 
  nk <- q+m+2 ## number of knots
  if (nk<=0) stop("k too small for m")
  x <- data[[object$term]]  ## the data
  xk <- knots[[object$term]] ## will be NULL if none supplied
  if (is.null(xk)) # space knots through data
  { n<-length(x)
    xk<-rep(0,q+m+2)
    xk[(m+2):(q+1)]<-seq(min(x),max(x),length=q-m)
    for (i in 1:(m+1)) {xk[i]<-xk[m+2]-(m+2-i)*(xk[m+3]-xk[m+2])}
    for (i in (q+2):(q+m+2)) {xk[i]<-xk[q+1]+(i-q-1)*(xk[m+3]-xk[m+2])}
  }
  if (length(xk)!=nk) # right number of knots?
  stop(paste("there should be ",nk," supplied knots"))
  ##  get model matrix...
  X <- splineDesign(xk,x,ord=m+2)
  ## Define matrix Sigma for monotone increasing & convex smooth
  Sig <- matrix(0,q,q)   
  Sig[,1] <- rep(1,q)
  for (i in 2:q) Sig[i:q,i] <- c(1:(q-i+1))

  X <- X%*%Sig # model matrix for the constrained term
  object$X <- X 
  object$Sigma <- Sig
  object$P <- list()
  object$S <- list()

  if (!object$fixed) # create the penalty matrix
  { P <- diff(diag(q-2),difference=1)
    P <- rbind(matrix(0,2,q-2),P) ## adding first two rows of zeros
    P <- cbind(matrix(0,q-1,2),P) ## adding first two columns of zeros
    object$P[[1]] <- P
    object$S[[1]] <- crossprod(P)
  }
  object$p.ident <- c(FALSE,rep(TRUE,q-1)) ## p.ident is an indicator of which coefficients must be positive (exponentiated) 
  object$rank <- ncol(object$X)  # penalty rank
  object$null.space.dim <-2 ## m+1  # dim. of unpenalized space
  object$C <- matrix(0, 0, ncol(X)) # to have no other constraints 
 
  ## get model matrix for 1st and 2nd derivatives of the smooth...
  h <- (max(x)-min(x))/(q-m-1) ## distance between two adjacent knots
  object$Xdf1 <- splineDesign(xk,x,ord=m+1)[,1:(q-1)]/h ## ord is by one less for the 1st derivative
  object$Xdf2 <- if (m==0) matrix(0,nrow(X),ncol(object$Xdf1)-1) ## for piecewise linear splines
                 else splineDesign(xk,x,ord=m)[,2:(q-2)]/h^2 ## ord is by two less for the 2nd derivative 

  object$knots <- xk;
  object$m <- m;
  object$df<-ncol(object$X)     # maximum DoF (if unconstrained)
  class(object)<- "micxBy.smooth"  # Give object a class
  object
}


## Prediction matrix for the `micxBy` smooth class... 

Predict.matrix.micxBy.smooth<-function(object,data)
## prediction method function for the `micxBy` smooth class
{ m <- object$m+1; # spline order, m+1=3 default for cubic spline
  q <- object$df 
  Sig <- matrix(0,q,q)   
  Sig[,1] <- rep(1,q)
  for (i in 2:q) Sig[i:q,i] <- c(1:(q-i+1))
  
  ## find spline basis inner knot range...
  ll <- object$knots[m+1];ul <- object$knots[length(object$knots)-m]
  m <- m + 1
  x <- data[[object$term]]
  n <- length(x)
  ind <- x<=ul & x>=ll ## data in range
  if (sum(ind)==n) { ## all in range
     X <- spline.des(object$knots,x,m)$design
     X <- X%*%Sig 
  } else { ## some extrapolation needed 
     ## matrix mapping coefs to value and slope at end points...
     D <- spline.des(object$knots,c(ll,ll,ul,ul),m,c(0,1,0,1))$design
     X <- matrix(0,n,ncol(D)) ## full predict matrix
     if (sum(ind)> 0)  X[ind,] <- spline.des(object$knots,x[ind],m)$design ## interior rows
     ## Now add rows for linear extrapolation...
     ind <- x < ll 
     if (sum(ind)>0) X[ind,] <- cbind(1,x[ind]-ll)%*%D[1:2,]
     ind <- x > ul
     if (sum(ind)>0) X[ind,] <- cbind(1,x[ind]-ul)%*%D[3:4,]
     X <- X%*%Sig 
  }
  X
}


############################################################
### Smooth construct for the convex/concave smooths ......
###########################################################

###########################################################
### Adding concave SCOP-spline construction *************
###########################################################

smooth.construct.cv.smooth.spec<- function(object, data, knots)
## construction of the concave smooth
{ 
 # require(splines)
  m <- object$p.order[1]
  if (is.na(m)) m <- 2 ## default 
  if (m<0) stop("silly m supplied")
  if (object$bs.dim<0) object$bs.dim <- 10 ## default
  q <- object$bs.dim 
  nk <- q+m+2 ## number of knots
  if (nk<=0) stop("k too small for m")
  x <- data[[object$term]]  ## the data
  xk <- knots[[object$term]] ## will be NULL if none supplied
  if (is.null(xk)) # space knots through data
  { n<-length(x)
    xk<-rep(0,q+m+2)
    xk[(m+2):(q+1)]<-seq(min(x),max(x),length=q-m)
    for (i in 1:(m+1)) {xk[i]<-xk[m+2]-(m+2-i)*(xk[m+3]-xk[m+2])}
    for (i in (q+2):(q+m+2)) {xk[i]<-xk[q+1]+(i-q-1)*(xk[m+3]-xk[m+2])}
  }
  if (length(xk)!=nk) # right number of knots?
  stop(paste("there should be ",nk," supplied knots"))
  #  get model matrix-------------
  X1 <- splineDesign(xk,x,ord=m+2)
  # use matrix Sigma and remove the first column for the monotone smooth
  Sig <- matrix(0,(q-1),(q-1))  # Define Sigma for concave smooth
  Sig[1:(q-1),1]<- c(1:(q-1))
  for (i in 2:(q-1)) Sig[i:(q-1),i]<--c(1:(q-i))
  
  X <- X1[,2:q]%*%Sig # model submatrix for the constrained term

  ## applying sum-to-zero (centering) constraint...
  cmx <- colMeans(X)
  X <- sweep(X,2,cmx) ## subtract cmx from columns 
  object$X <- X # the final model matrix
  object$cmX <- cmx

  object$Sigma <- Sig
  object$P <- list()
  object$S <- list()

  if (!object$fixed) # create the penalty matrix
  { P <- diff(diag(q-2),difference=1)
    object$P[[1]] <- P
    S <- matrix(0,q-1,q-1)
    S[2:(q-1),2:(q-1)] <- crossprod(P)
    object$S[[1]] <- S
  }
  object$p.ident <- rep(TRUE,q-1) ## p.ident is an indicator of which coefficients must be positive (exponentiated) 
  object$rank <- ncol(object$X)-1  # penalty rank
  object$null.space.dim <- 2 ##m+1  # dim. of unpenalized space
  object$C <- matrix(0, 0, ncol(X)) # to have no other constraints 
  
  ## get model matrix for 1st and 2nd derivatives of the smooth...
  h <- (max(x)-min(x))/(q-m-1) ## distance between two adjacent knots
  object$Xdf1 <- splineDesign(xk,x,ord=m+1)[,2:(q-1)]/h ## ord is by one less for the 1st derivative
  object$Xdf2 <- if (m==0) matrix(0,nrow(X),ncol(object$Xdf1)-1) ## for piecewise linear splines
                 else splineDesign(xk,x,ord=m)[,2:(q-2)]/h^2 ## ord is by two less for the 2nd derivative 

  object$knots <- xk;
  object$m <- m;
  object$df<-ncol(object$X)     # maximum DoF (if unconstrained)
  class(object)<-"cv.smooth"  # Give object a class
  object
}


## Prediction matrix for the `cv` smooth class******************

Predict.matrix.cv.smooth<-function(object,data)
## prediction method function for the `cv' smooth class
{ m <- object$m+1; # spline order, m+1=3 default for cubic spline
  q <- object$df +1
  Sig <- matrix(0,q,q)   
  Sig[,1] <- rep(1,q)
  Sig[2:q,2]<- c(1:(q-1))
  for (i in 3:q) Sig[i:q,i] <- -c(1:(q-i+1))

  ## find spline basis inner knot range...
  ll <- object$knots[m+1];ul <- object$knots[length(object$knots)-m]
  m <- m + 1
  x <- data[[object$term]]
  n <- length(x)
  ind <- x<=ul & x>=ll ## data in range
  if (sum(ind)==n) { ## all in range
     X <- spline.des(object$knots,x,m)$design
     X <- X%*%Sig 
  } else { ## some extrapolation needed 
     ## matrix mapping coefs to value and slope at end points...
     D <- spline.des(object$knots,c(ll,ll,ul,ul),m,c(0,1,0,1))$design
     X <- matrix(0,n,ncol(D)) ## full predict matrix
     if (sum(ind)> 0)  X[ind,] <- spline.des(object$knots,x[ind],m)$design ## interior rows
     ## Now add rows for linear extrapolation...
     ind <- x < ll 
     if (sum(ind)>0) X[ind,] <- cbind(1,x[ind]-ll)%*%D[1:2,]
     ind <- x > ul
     if (sum(ind)>0) X[ind,] <- cbind(1,x[ind]-ul)%*%D[3:4,]
     X <- X%*%Sig 
    ## X <- sweep(X,2,c(0,object$cmX))
  }
  X <- sweep(X,2,c(0,object$cmX))
  X
}

####################################################################################################
### Adding concave SCOP-spline construction without applying identifiability constraints
### to be used with numeric 'by' variable...
####################################################################################################
## when 'by' variable takes more than one value, the smooth terms are identifiable without a
## 'zero intercept' constraint, so they are left unconstrained...

smooth.construct.cvBy.smooth.spec<- function(object, data, knots)
## construction of the concave smooth
{ 
 # require(splines)
  m <- object$p.order[1]
  if (is.na(m)) m <- 2 ## default 
  if (m<0) stop("silly m supplied")
  if (object$bs.dim<0) object$bs.dim <- 10 ## default
  q <- object$bs.dim 
  nk <- q+m+2 ## number of knots
  if (nk<=0) stop("k too small for m")
  x <- data[[object$term]]  ## the data
  xk <- knots[[object$term]] ## will be NULL if none supplied
  if (is.null(xk)) # space knots through data
  { n<-length(x)
    xk<-rep(0,q+m+2)
    xk[(m+2):(q+1)]<-seq(min(x),max(x),length=q-m)
    for (i in 1:(m+1)) {xk[i]<-xk[m+2]-(m+2-i)*(xk[m+3]-xk[m+2])}
    for (i in (q+2):(q+m+2)) {xk[i]<-xk[q+1]+(i-q-1)*(xk[m+3]-xk[m+2])}
  }
  if (length(xk)!=nk) # right number of knots?
  stop(paste("there should be ",nk," supplied knots"))
  #  get model matrix-------------
  X <- splineDesign(xk,x,ord=m+2)
  # get matrix Sigma for concave smooth...
  Sig <- matrix(0,q,q)   
  Sig[,1] <- rep(1,q)
  Sig[2:q,2]<- c(1:(q-1))
  for (i in 3:q) Sig[i:q,i] <- -c(1:(q-i+1))
  
  X <- X%*%Sig # model matrix for the constrained term
  object$X<-X 
  object$Sigma <- Sig
  object$P <- list()
  object$S <- list()

  if (!object$fixed) # create the penalty matrix
  { P <- diff(diag(q-2),difference=1)
    P <- rbind(matrix(0,2,q-2),P) ## adding first two rows of zeros
    P <- cbind(matrix(0,q-1,2),P) ## adding first two columns of zeros
    object$P[[1]] <- P
    object$S[[1]] <- crossprod(P)
  }
  object$p.ident <- c(FALSE,rep(TRUE,q-1)) ## p.ident is an indicator of which coefficients must be positive (exponentiated) 
  object$rank <- ncol(object$X)  # penalty rank
  object$null.space.dim <- 2 ##m+1  # dim. of unpenalized space
  object$C <- matrix(0, 0, ncol(X)) # to have no other constraints 
  
  ## get model matrix for 1st and 2nd derivatives of the smooth...
  h <- (max(x)-min(x))/(q-m-1) ## distance between two adjacent knots
  object$Xdf1 <- splineDesign(xk,x,ord=m+1)[,1:(q-1)]/h ## ord is by one less for the 1st derivative
  object$Xdf2 <- if (m==0) matrix(0,nrow(X),ncol(object$Xdf1)-1) ## for piecewise linear splines
                 else splineDesign(xk,x,ord=m)[,2:(q-2)]/h^2 ## ord is by two less for the 2nd derivative 

  object$knots <- xk;
  object$m <- m;
  object$df<-ncol(object$X)     # maximum DoF (if unconstrained)
  class(object)<-"cvBy.smooth"  # Give object a class
  object
}


## Prediction matrix for the `cvBy` smooth class******************

Predict.matrix.cvBy.smooth<-function(object,data)
## prediction method function for the `cvBy' smooth class
{ m <- object$m+1; # spline order, m+1=3 default for cubic spline
  q <- object$df 
  Sig <- matrix(0,q,q)   
  Sig[,1] <- rep(1,q)
  Sig[2:q,2]<- c(1:(q-1))
  for (i in 3:q) Sig[i:q,i] <- -c(1:(q-i+1))

  ## find spline basis inner knot range...
  ll <- object$knots[m+1];ul <- object$knots[length(object$knots)-m]
  m <- m + 1
  x <- data[[object$term]]
  n <- length(x)
  ind <- x<=ul & x>=ll ## data in range
  if (sum(ind)==n) { ## all in range
     X <- spline.des(object$knots,x,m)$design
     X <- X%*%Sig 
  } else { ## some extrapolation needed 
     ## matrix mapping coefs to value and slope at end points...
     D <- spline.des(object$knots,c(ll,ll,ul,ul),m,c(0,1,0,1))$design
     X <- matrix(0,n,ncol(D)) ## full predict matrix
     if (sum(ind)> 0)  X[ind,] <- spline.des(object$knots,x[ind],m)$design ## interior rows
     ## Now add rows for linear extrapolation...
     ind <- x < ll 
     if (sum(ind)>0) X[ind,] <- cbind(1,x[ind]-ll)%*%D[1:2,]
     ind <- x > ul
     if (sum(ind)>0) X[ind,] <- cbind(1,x[ind]-ul)%*%D[3:4,]
     X <- X%*%Sig 
  }
  X
}



###########################################################
### Adding convex SCOP-spline construction *************
##########################################################

smooth.construct.cx.smooth.spec<- function(object, data, knots)
## construction of the convex smooth
{ 
  # require(splines)
  m <- object$p.order[1]
  if (is.na(m)) m <- 2 ## default 
  if (m<0) stop("silly m supplied")
  if (object$bs.dim<0) object$bs.dim <- 10 ## default
  q <- object$bs.dim 
  nk <- q+m+2 ## number of knots
  if (nk<=0) stop("k too small for m")
  x <- data[[object$term]]  ## the data
  xk <- knots[[object$term]] ## will be NULL if none supplied
  if (is.null(xk)) # space knots through data
  { n<-length(x)
    xk<-rep(0,q+m+2)
    xk[(m+2):(q+1)]<-seq(min(x),max(x),length=q-m)
    for (i in 1:(m+1)) {xk[i]<-xk[m+2]-(m+2-i)*(xk[m+3]-xk[m+2])}
    for (i in (q+2):(q+m+2)) {xk[i]<-xk[q+1]+(i-q-1)*(xk[m+3]-xk[m+2])}
  }
  if (length(xk)!=nk) # right number of knots?
  stop(paste("there should be ",nk," supplied knots"))
  #  get model matrix-------------
  X1 <- splineDesign(xk,x,ord=m+2)
  # use matrix Sigma and remove the first column for the monotone smooth
  Sig <- matrix(0,(q-1),(q-1))   # Define Sigma for convex smooth
  Sig[1:(q-1),1]<- -c(1:(q-1))
  for (i in 2:(q-1)) Sig[i:(q-1),i]<- c(1:(q-i))
  
  X <- X1[,2:q]%*%Sig # model submatrix for the constrained term

  ## applying sum-to-zero (centering) constraint...
  cmx <- colMeans(X)
  X <- sweep(X,2,cmx) ## subtract cmx from columns 
  object$X <- X # the final model matrix
  object$cmX <- cmx

  object$Sigma <- Sig
  object$P <- list()
  object$S <- list()

  if (!object$fixed) # create the penalty matrix
  { P <- diff(diag(q-2),difference=1)
    object$P[[1]] <- P
    S <- matrix(0,q-1,q-1)
    S[2:(q-1),2:(q-1)] <- crossprod(P)
    object$S[[1]] <- S
  }
  object$p.ident <- rep(TRUE,q-1) ## p.ident is an indicator of which coefficients must be positive (exponentiated) 
  object$rank <- ncol(object$X)-1  # penalty rank
  object$null.space.dim <- 2 ##m+1  # dim. of unpenalized space
  object$C <- matrix(0, 0, ncol(X)) # to have no other constraints 
  
  ## get model matrix for 1st and 2nd derivatives of the smooth...
  h <- (max(x)-min(x))/(q-m-1) ## distance between two adjacent knots
  object$Xdf1 <- splineDesign(xk,x,ord=m+1)[,2:(q-1)]/h ## ord is by one less for the 1st derivative
  object$Xdf2 <- if (m==0) matrix(0,nrow(X),ncol(object$Xdf1)-1) ## for piecewise linear splines
                 else splineDesign(xk,x,ord=m)[,2:(q-2)]/h^2 ## ord is by two less for the 2nd derivative 

  object$knots <- xk;
  object$m <- m;
  object$df<-ncol(object$X)     # maximum DoF (if unconstrained)
  class(object)<-"cx.smooth"  # Give object a class
  object
}


## Prediction matrix for the `cx` smooth class *************************

Predict.matrix.cx.smooth<-function(object,data)
## prediction method function for the `cx' smooth class
{ m <- object$m+1; # spline order, m+1=3 default for cubic spline
  q <- object$df +1
  Sig <- matrix(0,q,q)   
  Sig[,1] <- rep(1,q)
  Sig[2:q,2]<- -c(1:(q-1))
  for (i in 3:q) Sig[i:q,i] <- c(1:(q-i+1))
  
  ## find spline basis inner knot range...
  ll <- object$knots[m+1];ul <- object$knots[length(object$knots)-m]
  m <- m + 1
  x <- data[[object$term]]
  n <- length(x)
  ind <- x<=ul & x>=ll ## data in range
  if (sum(ind)==n) { ## all in range
     X <- spline.des(object$knots,x,m)$design
     X <- X%*%Sig 
  } else { ## some extrapolation needed 
     ## matrix mapping coefs to value and slope at end points...
     D <- spline.des(object$knots,c(ll,ll,ul,ul),m,c(0,1,0,1))$design
     X <- matrix(0,n,ncol(D)) ## full predict matrix
     if (sum(ind)> 0)  X[ind,] <- spline.des(object$knots,x[ind],m)$design ## interior rows
     ## Now add rows for linear extrapolation...
     ind <- x < ll 
     if (sum(ind)>0) X[ind,] <- cbind(1,x[ind]-ll)%*%D[1:2,]
     ind <- x > ul
     if (sum(ind)>0) X[ind,] <- cbind(1,x[ind]-ul)%*%D[3:4,]
     X <- X%*%Sig 
    ## X <- sweep(X,2,c(0,object$cmX))
  }
  X <- sweep(X,2,c(0,object$cmX))
  X
}

###########################################################
### Adding convex SCOP-spline construction without identifiability constraints
### to be used with numeric 'by' variable and linear functional terms...
##########################################################

smooth.construct.cxBy.smooth.spec<- function(object, data, knots)
## construction of the convex smooth
{ 
  m <- object$p.order[1]
  if (is.na(m)) m <- 2 ## default 
  if (m<0) stop("silly m supplied")
  if (object$bs.dim<0) object$bs.dim <- 10 ## default
  q <- object$bs.dim 
  nk <- q+m+2 ## number of knots
  if (nk<=0) stop("k too small for m")
  x <- data[[object$term]]  ## the data
  xk <- knots[[object$term]] ## will be NULL if none supplied
  if (is.null(xk)) # space knots through data
  { n<-length(x)
    xk<-rep(0,q+m+2)
    xk[(m+2):(q+1)]<-seq(min(x),max(x),length=q-m)
    for (i in 1:(m+1)) {xk[i]<-xk[m+2]-(m+2-i)*(xk[m+3]-xk[m+2])}
    for (i in (q+2):(q+m+2)) {xk[i]<-xk[q+1]+(i-q-1)*(xk[m+3]-xk[m+2])}
  }
  if (length(xk)!=nk) # right number of knots?
  stop(paste("there should be ",nk," supplied knots"))
  #  get model matrix-------------
  X <- splineDesign(xk,x,ord=m+2)
  # get matrix Sigma for convex  smooth...
  Sig <- matrix(0,q,q)   
  Sig[,1] <- 1
  Sig[2:q,2]<- -c(1:(q-1))
  for (i in 3:q) Sig[i:q,i] <- c(1:(q-i+1))
  
  X <- X%*%Sig # model submatrix for the constrained term
  object$X<-X # the final model matrix
  object$Sigma <- Sig
  object$P <- list()
  object$S <- list()

  if (!object$fixed) # create the penalty matrix
  { P <- diff(diag(q-2),difference=1)
    P <- rbind(matrix(0,2,q-2),P) ## adding first two rows of zeros
    P <- cbind(matrix(0,q-1,2),P) ## adding first two columns of zeros
    object$P[[1]] <- P
    object$S[[1]] <- crossprod(P)
  }
  object$p.ident <- c(FALSE,rep(TRUE,q-1)) ## p.ident is an indicator of which coefficients must be positive (exponentiated) 
  object$rank <- ncol(object$X)  # penalty rank
  object$null.space.dim <- 2 ##m+1  # dim. of unpenalized space
  object$C <- matrix(0, 0, ncol(X)) # to have no other constraints 
  
  ## get model matrix for 1st and 2nd derivatives of the smooth...
  h <- (max(x)-min(x))/(q-m-1) ## distance between two adjacent knots
  object$Xdf1 <- splineDesign(xk,x,ord=m+1)[,1:(q-1)]/h ## ord is by one less for the 1st derivative
  object$Xdf2 <- if (m==0) matrix(0,nrow(X),ncol(object$Xdf1)-1) ## for piecewise linear splines
                 else splineDesign(xk,x,ord=m)[,2:(q-2)]/h^2 ## ord is by two less for the 2nd derivative 

  object$knots <- xk;
  object$m <- m;
  object$df<-ncol(object$X)     # maximum DoF (if unconstrained)
  class(object)<-"cxBy.smooth"  # Give object a class
  object
}


## Prediction matrix for the `cxBy` smooth class *************************

Predict.matrix.cxBy.smooth<-function(object,data)
## prediction method function for the `cxBy' smooth class
{ m <- object$m+1; # spline order, m+1=3 default for cubic spline
  q <- object$df 
  Sig <- matrix(0,q,q)   
  Sig[,1] <- rep(1,q)
  Sig[2:q,2]<- -c(1:(q-1))
  for (i in 3:q) Sig[i:q,i] <- c(1:(q-i+1))
  
  ## find spline basis inner knot range...
  ll <- object$knots[m+1];ul <- object$knots[length(object$knots)-m]
  m <- m + 1
  x <- data[[object$term]]
  n <- length(x)
  ind <- x<=ul & x>=ll ## data in range
  if (sum(ind)==n) { ## all in range
     X <- spline.des(object$knots,x,m)$design
     X <- X%*%Sig 
  } else { ## some extrapolation needed 
     ## matrix mapping coefs to value and slope at end points...
     D <- spline.des(object$knots,c(ll,ll,ul,ul),m,c(0,1,0,1))$design
     X <- matrix(0,n,ncol(D)) ## full predict matrix
     if (sum(ind)> 0)  X[ind,] <- spline.des(object$knots,x[ind],m)$design ## interior rows
     ## Now add rows for linear extrapolation...
     ind <- x < ll 
     if (sum(ind)>0) X[ind,] <- cbind(1,x[ind]-ll)%*%D[1:2,]
     ind <- x > ul
     if (sum(ind)>0) X[ind,] <- cbind(1,x[ind]-ul)%*%D[3:4,]
     X <- X%*%Sig 
  }
  X
}


#####################################################
### Adding positive SCOP-spline construction 
######################################################


smooth.construct.po.smooth.spec<- function(object, data, knots)
## construction of the positively costrained smooth
{ 
  m <- object$p.order[1]
  if (is.na(m)) m <- 2 ## default for cubic spline
  if (m<0) stop("silly m supplied")
  if (object$bs.dim<0) object$bs.dim <- 10 ## default
  q <- object$bs.dim 
  nk <- q+m+2 ## number of knots
  if (nk<=0) stop("k too small for m")
  x <- data[[object$term]]  ## the data
  xk <- knots[[object$term]] ## will be NULL if none supplied
  if (is.null(xk)) # space knots through data
     { n<-length(x)
       xk<-rep(0,q+m+2)
       xk[(m+2):(q+1)]<-seq(min(x),max(x),length=q-m)
       for (i in 1:(m+1)) {xk[i]<-xk[m+2]-(m+2-i)*(xk[m+3]-xk[m+2])}
       for (i in (q+2):(q+m+2)) {xk[i]<-xk[q+1]+(i-q-1)*(xk[m+3]-xk[m+2])}
     }
  if (length(xk)!=nk) # right number of knots?
  stop(paste("there should be ", nk," supplied knots"))
  #  get model matrix-------------
  X <- splineDesign(xk,x,ord=m+2)
 
  Sig <- diag(1,(q-1))  ## identity matrix 
 ## Sig <- diag(1,q)  ## identity matrix
     
  X <- X[,-1] ## X[,1:(q-1)] 
  object$X <- X # the final model matrix
  object$P <- list()
  object$S <- list()
  object$Sigma <- Sig

  if (!object$fixed) # create the penalty matrix
  { P <- diff(diag(q-1),difference=1)
   ## P <- diff(diag(q),difference=1)
    object$P[[1]] <- P
    object$S[[1]] <- crossprod(P)
  }
  object$p.ident <- rep(TRUE,q-1) ## p.ident is an indicator of which coefficients must be positive (exponentiated)
 ## object$p.ident <- rep(TRUE,q)
  object$rank <- ncol(object$X)-1  # penalty rank
 ## object$rank <- ncol(object$X)  # penalty rank
  object$null.space.dim <-2 ## m+1  # dim. of unpenalized space
  object$C <- matrix(0, 0, ncol(X)) # to have no other constraints 
  object$knots <- xk;
  object$m <- m;
  object$df<-ncol(object$X)     # maximum DoF (if unconstrained)

  ## get model matrix for 1st and 2nd derivatives of the smooth...
  h <- (max(x)-min(x))/(q-m-1) ## distance between two adjacent knots
  object$Xdf1 <- splineDesign(xk,x,ord=m+1)[,2:(q-1)]/h ## ord is by one less for the 1st derivative
  object$Xdf2 <-  if (m>0) splineDesign(xk,x,ord=m)[,2:(q-2)]/h^2 ## ord is by two less for the 2nd deriv
                  else matrix(0,nrow(X),ncol(object$Xdf1)-1) ## for piecewise linear splines    

  class(object)<-"po.smooth"  # Give object a class
  object
}

## NOte: maybe no need to set the first coefficient to zero, as no intercept is required in the model with 'po' smooth
## Prediction matrix for the `po` smooth class... 

Predict.matrix.po.smooth<-function(object,data)
## prediction method function for the `po' smooth class
{ m <- object$m+1; # spline order, m+1=3 default for cubic spline
 ## q <- object$df +1
 ## Sig <- diag(1,q)   # Define Matrix Sigma
  ## find spline basis inner knot range...
  ll <- object$knots[m+1];ul <- object$knots[length(object$knots)-m]
  m <- m + 1
  x <- data[[object$term]]
  n <- length(x)
  ind <- x<=ul & x>=ll ## data in range
  if (sum(ind)==n) { ## all in range
     X <- spline.des(object$knots,x,m)$design
  } else { ## some extrapolation needed 
     ## matrix mapping coefs to value and slope at end points...
     D <- spline.des(object$knots,c(ll,ll,ul,ul),m,c(0,1,0,1))$design
     X <- matrix(0,n,ncol(D)) ## full predict matrix
     if (sum(ind)> 0)  X[ind,] <- spline.des(object$knots,x[ind],m)$design ## interior rows
     ## Now add rows for linear extrapolation...
     ind <- x < ll 
     if (sum(ind)>0) X[ind,] <- cbind(1,x[ind]-ll)%*%D[1:2,]
     ind <- x > ul
     if (sum(ind)>0) X[ind,] <- cbind(1,x[ind]-ul)%*%D[3:4,]
  }
  X
}

##################################################################
### Adding SCOP-spline with decreasing and positivity constraints
##################################################################

smooth.construct.dpo.smooth.spec<- function(object, data, knots)
## construction of the positively costrained smooth
{ 
  m <- object$p.order[1]
  if (is.na(m)) m <- 2 ## default for cubic spline
  if (m<0) stop("silly m supplied")
  if (object$bs.dim<0) object$bs.dim <- 10 ## default
  q <- object$bs.dim 
  nk <- q+m+2 ## number of knots
  if (nk<=0) stop("k too small for m")
  x <- data[[object$term]]  ## the data
  xk <- knots[[object$term]] ## will be NULL if none supplied
  if (is.null(xk)) # space knots through data
     { n<-length(x)
       xk<-rep(0,q+m+2)
       xk[(m+2):(q+1)]<-seq(min(x),max(x),length=q-m)
       for (i in 1:(m+1)) {xk[i]<-xk[m+2]-(m+2-i)*(xk[m+3]-xk[m+2])}
       for (i in (q+2):(q+m+2)) {xk[i]<-xk[q+1]+(i-q-1)*(xk[m+3]-xk[m+2])}
     }
  if (length(xk)!=nk) # right number of knots?
  stop(paste("there should be ", nk," supplied knots"))
  #  get model matrix-------------
  X1 <- splineDesign(xk,x,ord=m+2)
   
  # get matrix Sigma and remove the last column for the decreasing and positive smooth 
  Sig <- matrix(1,q,q)  ## coef summation matrix
  Sig[lower.tri(Sig)] <-0
  X <- X1%*%Sig
  X <- X[,-q]
  
  object$X <- X # the final model matrix
 
  object$Sigma <- Sig[-q,-q]
  object$P <- list()
  object$S <- list()

  if (!object$fixed) # create the penalty matrix
  { P <- diff(diag(q-1),difference=1)
    object$P[[1]] <- P
    object$S[[1]] <- crossprod(P)
  }   
   
  object$p.ident <- rep(TRUE,q-1) ## p.ident is an indicator of which coefficients must be positive (exponentiated)
 ## object$p.ident <- rep(TRUE,q)
  object$rank <- ncol(object$X)-1  # penalty rank
 ## object$rank <- ncol(object$X)  # penalty rank
  object$null.space.dim <-2 ## m+1  # dim. of unpenalized space
  object$C <- matrix(0, 0, ncol(X)) # to have no other constraints 
  object$knots <- xk;
  object$m <- m;
  object$df<-ncol(object$X)     # maximum DoF (if unconstrained)

  ## get model matrix for 1st and 2nd derivatives of the smooth...
  h <- (max(x)-min(x))/(q-m-1) ## distance between two adjacent knots
  object$Xdf1 <- splineDesign(xk,x,ord=m+1)[,2:(q-1)]/h ## ord is by one less for the 1st derivative
  object$Xdf2 <-  if (m>0) splineDesign(xk,x,ord=m)[,2:(q-2)]/h^2 ## ord is by two less for the 2nd deriv
                  else matrix(0,nrow(X),ncol(object$Xdf1)-1) ## for piecewise linear splines    

  class(object)<-"dpo.smooth"  # Give object a class
  object
}

## Prediction matrix for the `dpo` smooth class... 
Predict.matrix.dpo.smooth<-function(object,data)
## prediction method function for the `dpo' smooth class
{ m <- object$m+1; # spline order, m+1=3 default for cubic spline
  q <- object$df +1
  # elements of matrix Sigma for decreasing smooth...
 # Sig <- matrix(as.numeric(rep(1:q,q)>=rep(1:q,each=q)),q,q) ## coef summation matrix
  Sig <- matrix(1,q,q)  ## coef summation matrix
  Sig[lower.tri(Sig)] <-0
  
  ## find spline basis inner knot range...
  ll <- object$knots[m+1];ul <- object$knots[length(object$knots)-m]
  m <- m + 1
  x <- data[[object$term]]
  n <- length(x)
  ind <- x<=ul & x>=ll ## data in range
  if (sum(ind)==n) { ## all in range
     X <- spline.des(object$knots,x,m)$design
     X <- X%*%Sig 
  } else { ## some extrapolation needed 
     ## matrix mapping coefs to value and slope at end points...
     D <- spline.des(object$knots,c(ll,ll,ul,ul),m,c(0,1,0,1))$design
     X <- matrix(0,n,ncol(D)) ## full predict matrix
     if (sum(ind)> 0)  X[ind,] <- spline.des(object$knots,x[ind],m)$design ## interior rows
     ## Now add rows for linear extrapolation...
     ind <- x < ll 
     if (sum(ind)>0) X[ind,] <- cbind(1,x[ind]-ll)%*%D[1:2,]
     ind <- x > ul
     if (sum(ind)>0) X[ind,] <- cbind(1,x[ind]-ul)%*%D[3:4,]
     X <- X%*%Sig     
  }
  X
}

##################################################################
### Adding SCOP-spline with inreasing and positivity constraints
##################################################################

smooth.construct.ipo.smooth.spec<- function(object, data, knots)
## construction of the increasing and positively costrained smooth
{ 
  m <- object$p.order[1]
  if (is.na(m)) m <- 2 ## default for cubic spline
  if (m<0) stop("silly m supplied")
  if (object$bs.dim<0) object$bs.dim <- 10 ## default
  q <- object$bs.dim 
  nk <- q+m+2 ## number of knots
  if (nk<=0) stop("k too small for m")
  x <- data[[object$term]]  ## the data
  xk <- knots[[object$term]] ## will be NULL if none supplied
  if (is.null(xk)) # space knots through data
     { n<-length(x)
       xk<-rep(0,q+m+2)
       xk[(m+2):(q+1)]<-seq(min(x),max(x),length=q-m)
       for (i in 1:(m+1)) {xk[i]<-xk[m+2]-(m+2-i)*(xk[m+3]-xk[m+2])}
       for (i in (q+2):(q+m+2)) {xk[i]<-xk[q+1]+(i-q-1)*(xk[m+3]-xk[m+2])}
     }
  if (length(xk)!=nk) # right number of knots?
  stop(paste("there should be ", nk," supplied knots"))
  #  get model matrix-------------
  X1 <- splineDesign(xk,x,ord=m+2)
   
  # get matrix Sigma and remove the first column for the increasing and positive smooth 
  Sig <- matrix(1,q,q)  ## coef summation matrix
  Sig[upper.tri(Sig)] <-0
  X <- X1%*%Sig
  X <- X[,-1]
  object$X <- X # the final model matrix
  object$Sigma <- Sig[-1,-1]
  object$P <- list()
  object$S <- list()

  if (!object$fixed) # create the penalty matrix
  { P <- diff(diag(q-1),difference=1)
    object$P[[1]] <- P
    object$S[[1]] <- crossprod(P)
  }   
   
  object$p.ident <- rep(TRUE,q-1) ## p.ident is an indicator of which coefficients must be positive (exponentiated)
 ## object$p.ident <- rep(TRUE,q)
  object$rank <- ncol(object$X)-1  # penalty rank
 ## object$rank <- ncol(object$X)  # penalty rank
  object$null.space.dim <-2 ## m+1  # dim. of unpenalized space
  object$C <- matrix(0, 0, ncol(X)) # to have no other constraints 
  object$knots <- xk;
  object$m <- m;
  object$df<-ncol(object$X)     # maximum DoF (if unconstrained)

  ## get model matrix for 1st and 2nd derivatives of the smooth...
  h <- (max(x)-min(x))/(q-m-1) ## distance between two adjacent knots
  object$Xdf1 <- splineDesign(xk,x,ord=m+1)[,2:(q-1)]/h ## ord is by one less for the 1st derivative
  object$Xdf2 <-  if (m>0) splineDesign(xk,x,ord=m)[,2:(q-2)]/h^2 ## ord is by two less for the 2nd deriv
                  else matrix(0,nrow(X),ncol(object$Xdf1)-1) ## for piecewise linear splines    

  class(object)<-"ipo.smooth"  # Give object a class
  object
}

## Prediction matrix for the `ipo` smooth class... 
Predict.matrix.ipo.smooth<-function(object,data)
## prediction method function for the `ipo' smooth class
{ m <- object$m+1; # spline order, m+1=3 default for cubic spline
  q <- object$df +1
  # elements of matrix Sigma for increasing smooth...
 # Sig <- matrix(as.numeric(rep(1:q,q)>=rep(1:q,each=q)),q,q) ## coef summation matrix
  Sig <- matrix(1,q,q)  ## coef summation matrix
  Sig[upper.tri(Sig)] <-0
  
  ## find spline basis inner knot range...
  ll <- object$knots[m+1];ul <- object$knots[length(object$knots)-m]
  m <- m + 1
  x <- data[[object$term]]
  n <- length(x)
  ind <- x<=ul & x>=ll ## data in range
  if (sum(ind)==n) { ## all in range
     X <- spline.des(object$knots,x,m)$design
     X <- X%*%Sig 
  } else { ## some extrapolation needed 
     ## matrix mapping coefs to value and slope at end points...
     D <- spline.des(object$knots,c(ll,ll,ul,ul),m,c(0,1,0,1))$design
     X <- matrix(0,n,ncol(D)) ## full predict matrix
     if (sum(ind)> 0)  X[ind,] <- spline.des(object$knots,x[ind],m)$design ## interior rows
     ## Now add rows for linear extrapolation...
     ind <- x < ll 
     if (sum(ind)>0) X[ind,] <- cbind(1,x[ind]-ll)%*%D[1:2,]
     ind <- x > ul
     if (sum(ind)>0) X[ind,] <- cbind(1,x[ind]-ul)%*%D[3:4,]
     X <- X%*%Sig     
  }
  X
}

##################################################################
### Adding cyclic P-spline with positivity constraint
## based on  R routines for the package mgcv (c) Simon Wood 2000-2019
##################################################################
cSplineDes <- function (x, knots, ord = 4,derivs=0)
{ ## cyclic version of spline design... the package mgcv (c) Simon Wood 2000-2019
  ##require(splines)
  nk <- length(knots)
  if (ord<2) stop("order too low")
  if (nk<ord) stop("too few knots")
  knots <- sort(knots)
  k1 <- knots[1]
  if (min(x)<k1||max(x)>knots[nk]) stop("x out of range")
  xc <- knots[nk-ord+1] ## wrapping involved above this point
  ## copy end intervals to start, for wrapping purposes...
  knots <- c(k1-(knots[nk]-knots[(nk-ord+1):(nk-1)]),knots)
  ind <- x>xc ## index for x values where wrapping is needed
  X1 <- splineDesign(knots,x,ord,outer.ok=TRUE,derivs=derivs)
  x[ind] <- x[ind] - max(knots) + k1
  if (sum(ind)) { ## wrapping part...
    X2 <- splineDesign(knots,x[ind],ord,outer.ok=TRUE,derivs=derivs) 
    X1[ind,] <- X1[ind,] + X2
  }
  X1 ## final model matrix
} ## cSplineDes


smooth.construct.cpop.smooth.spec<- function(object, data, knots)
## construction of the cyclic and positively costrained smooth
##`s(x,bs="cpop",m=c(2,1))' would be a cubic B-spline basis with a 1st order difference 
## penalty with positivity constraint. m==c(0,0) would be linear splines with a ridge penalty).
{ 
 if (length(object$p.order)==1) m <- rep(object$p.order,2) 
  else m <- object$p.order  ## m[1] - basis order, m[2] - penalty order
  m[is.na(m)] <- 2 ## default
  object$p.order <- m
  if (object$bs.dim<0) object$bs.dim <- max(10,m[1]) ## default
  nk <- object$bs.dim +1  ## number of interior knots
  if (nk<=m[1]) stop("basis dimension too small for b-spline order")
  if (length(object$term)!=1) stop("Basis only handles 1D smooths")
  x <- data[[object$term]]    # find the data
  k <- knots[[object$term]]
  
  if (is.null(k)) { x0 <- min(x);x1 <- max(x) } else
  if (length(k)==2) { 
    x0 <- min(k);x1 <- max(k);
    if (x0>min(x)||x1<max(x)) stop("knot range does not include data")
  } 
  if (is.null(k)||length(k)==2) {
     k <- seq(x0,x1,length=nk)  
  } else {
    if (length(k)!=nk) 
    stop(paste("there should be ",nk," supplied knots"))
  }

  if (length(k)!=nk) stop(paste("there should be",nk,"knots supplied"))

  X <- cSplineDes(x,k,ord=m[1]+2)  ## model matrix
  
  if (!is.null(k)) {
    if (sum(colSums(X)==0)>0) warning("knot range is so wide that there is *no* information about some basis coefficients")
  }  

  p.ord <- m[2]
  np <- ncol(X)
  if (p.ord>np-1) stop("penalty order too high for basis dimension")
  
  Sig <- diag(1,(np-1))  ## identity matrix 
  object$X <- X[,-1] 
  object$Sigma <- Sig

  object$p.ident <- rep(TRUE,np-1) ## p.ident is an indicator of which coefficients must be positive (exponentiated)
  object$C <- matrix(0, 0, ncol(object$X)) # to have no other constraints 


  ## now construct penalty...  
  np <- np-1
  
  De <- diag(np + p.ord)
  if (p.ord>0) { 
    for (i in 1:p.ord) De <- diff(De)
    D <- De[,-(1:p.ord)]
    D[,(np-p.ord+1):np] <-  D[,(np-p.ord+1):np] + De[,1:p.ord]
  } else D <- De
  object$S <- list(t(D)%*%D)  # get penalty

  ## other stuff...
  object$rank <- ncol(object$X)-1 # penalty rank
  object$null.space.dim <- 1    # dimension of unpenalized space
  object$knots <- k; 
  object$m <- m      # store p-spline specific info.   
  object$df<-ncol(object$X)     # maximum DoF (if unconstrained)
   
  class(object)<-"cpopspline.smooth"  # Give object a class
  object
}


Predict.matrix.cpopspline.smooth<-function(object,data)
## prediction method function for the cpopspline smooth class
{ x <- data[[object$term]] 
  k0 <- min(object$knots);k1 <- max(object$knots) 
  if (min(x)<k0||max(x)>k1) x <- cwrap(k0,k1,x)
  X <- cSplineDes(x,object$knots,object$m[1]+2)
  X
}