# R/l1looDesign.R In cobs99: Constrained B-splines -- outdated 1999 version

#### Documented in l1.designloo.design

```####  Create B-Spline Design matrices for COBS :
####
####  l1.design()  -- L_1         penalty
####  loo.design() -- L_Infinity  penalty

### used to be part of ./cobs.R

l1.design <-
function(x, w, constraint, equal, smaller, greater, gradient,
knots, pw, n.equal, n.smaller, n.greater, n.gradient,
nrq, nl1, neqc, niqc, nvar, lambda)
{
##=########################################################################
##
## Generate the pseudo design matrix for L1 penalty
##
##=########################################################################

## create the pseudo design
##

ks <- 2 ## <==> degree == 1
nk <- length(knots) + 2 ## *(ks - 1)
nkm3 <- nk - 3
nkm4 <- nk - 4
ncoef <- nk - ks
niqc1 <- 0 # per default, e.g., if(constraint == "none") but has 'pointwise'

nobs <- nrq + nl1 + neqc + niqc
X <- matrix(0, nobs, nvar)

ox <- order(x)
sortx <- x[ox]
z1 <- .splBasis(ord = ks, knots, ncoef, xo = sortx)
idx1 <- cbind(rep(ox, rep(ks, nrq)),
c(outer(1:ks, z1\$offsets, "+")))
X[idx1] <- t(t(z1\$design) * rep(w[ox], ks))
##
## formulate inequality constraints for the pseudo X
##
if(lambda != 0 || constraint != "none") {
z2 <- .splBasis(ord = ks, knots, ncoef, xo = knots[1:nkm3],
derivs = rep(1, nkm3))
if(lambda != 0) {
idx2 <- array(rep(1:nkm4, 2), c(nkm4, 2))
idx2[, 1] <- idx2[, 1] + nrq
X[idx2]   <-  - z2\$design[1, 1:nkm4] * lambda
idx2[, 2] <- idx2[, 2] + 2
X[idx2]   <- z2\$design[2, 2:nkm3] * lambda
idx2[, 2] <- idx2[, 2] - 1
X[idx2] <- (z2\$design[1, 2:nkm3] -
z2\$design[2, 1:nkm4]) * lambda
##
## assign different weight to roughness penalty
##
X[nrq + 1:nkm4, ] <- X[nrq + 1:nkm4, ] * pw
}
if(constraint == "increase" || constraint == "decrease") {
niqc1 <- nkm3
idx3 <- cbind(rep(nrq+nl1+neqc + 1:niqc1, rep(ks, niqc1)),
c(outer(1:ks, z2\$offsets, "+")))
X[idx3] <- if(constraint == "increase") z2\$design else -z2\$design
}
else if (constraint == "periodic") {
##
## this portion corresponds to equality constraints
##
neqc3 <- 2
z1.3 <- .splBasis(ord = ks, knots, ncoef,
xo = sortx[c(1,nrq)], derivs = c(1,1))
idx1.3 <- cbind(rep(1:neqc3, rep(ks, neqc3)),
c(outer(1:ks, z1.3\$offsets,"+")))
X.temp <- matrix(0,neqc3,nvar)
X.temp[idx1.3] <- z1.3\$design
X[nrq+nl1+neqc,  ] <- X.temp[2,] - X.temp[1,]
X[nrq+nl1+neqc-1,] <- X[ox[nrq],]- X[ox[1],]
}
else if(constraint == "convex" || constraint == "concave") {
niqc1 <- nkm4
idx3 <- array(rep(1:nkm4, 2), c(nkm4, 2))
idx3[, 1] <- idx3[, 1] + nrq + nl1 + neqc
sgn <- if(constraint == "convex") +1 else -1
X[idx3] <- -sgn* z2\$design[1, 1:nkm4]
idx3[,2] <- idx3[,2] + 2
X[idx3] <- sgn * z2\$design[2, 2:nkm3]
idx3[,2] <- idx3[,2] - 1
X[idx3] <- sgn * (z2\$design[1, 2:nkm3] - z2\$design[2, 1:nkm4])
}
}
niqc2 <- n.smaller
niqc3 <- n.greater
if(n.smaller > 0) {
o.smaller <- order(smaller[,2])
smaller.o <- smaller[,2][o.smaller]
z3.1 <- .splBasis(ord = ks, knots, ncoef, xo = smaller.o)
idx3.1 <- cbind(rep(nrq+nl1+neqc+niqc1 + o.smaller,rep(ks,niqc2)),
c(outer(1:ks,z3.1\$offsets,"+")))
X[idx3.1] <- -z3.1\$design
}
if(n.greater > 0) {
o.greater <- order(greater[,2])
greater.o <- greater[,2][o.greater]
z3.2 <- .splBasis(ord = ks, knots, ncoef, xo = greater.o)
idx3.2 <- cbind(rep(nrq+nl1+neqc+niqc1+niqc2 + o.greater,rep(ks,niqc3)),
c(outer(1:ks,z3.2\$offsets,"+")))
X[idx3.2] <- z3.2\$design
}
neqc1 <- n.equal

if(n.equal > 0) { ## formulate equality constraints for the pseudo X

o.equal <- order(equal[,2])
equal.o <- equal[,2][o.equal]
z1.1 <- .splBasis(ord = ks, knots, ncoef, xo = equal.o)
idx1.1 <- cbind(rep(nrq+nl1+ o.equal, rep(ks, neqc1)),
c(outer(1:ks, z1.1\$offsets,"+")))
X[idx1.1] <- z1.1\$design
}

z1.2 <- .splBasis(ord = ks, knots, ncoef, xo = gradient.o,
derivs = rep(1, neqc2))
idx1.2 <- cbind(rep(nrq+nl1+neqc1+ o.gradient, rep(ks, neqc2)),
c(outer(1:ks, z1.2\$offsets,"+")))
X[idx1.2] <- z1.2\$design
}
return(X)
} ## l1.design

loo.design <- function(x, w, constraint, equal, smaller, greater, gradient,
knots, pw, n.equal, n.smaller, n.greater, n.gradient,
nrq, nl1, neqc, niqc, nvar, lambda)
{
##=########################################################################
##
## Generate the pseudo design matrix for L_oo penalty
##
##=########################################################################

##x <- as.matrix(x, ncol = 1)

ks <- 3 ## <==> degree == 2
nk <- length(knots) + 2*(ks - 1)# = length(new.knots)
ncoef <- nk - ks
nd <- nk - 5
ox <- order(x)
sortx <- x[ox]
nrql1 <- nrq + nl1
nrleq <- nrql1 + neqc
nobs <- nrleq + niqc
X <- matrix(0, nobs, nvar)
z1 <- .splBasis(ord = ks, knots, ncoef, xo = sortx)
idx1 <- cbind(rep(ox, rep(ks, nrq)),
c(outer(1:ks, z1\$offsets, "+")))
X[idx1] <- t(t(z1\$design) * rep(w[ox], ks))
##
## formulate the inequality constraints -- s''()
##
if(lambda != 0) {
niqc1 <- 2*nd
z2 <- .splBasis(ord = ks, knots, ncoef, xo = knots[1:nd],
derivs = rep(2, nd))
X[(nrq+1):nrql1,] <- cbind(c(rep(0,nvar-1), lambda))
idx2 <- cbind(rep(nrleq + 1:nd, rep(ks,nd)),
c(outer(1:ks,z2\$offsets,"+")))
X[idx2] <- z2\$design
idx2[,1] <- idx2[,1]+nd
X[idx2] <- -z2\$design
## assign different weight to roughness penalty
X[nrleq + 1:niqc1, ] <-
X[nrleq + 1:niqc1, ] * rep(pw, 2)
X[nrleq + 1:niqc1, nvar ] <- rep(1,niqc1)
} else niqc1 <- 0
niqc2 <- n.smaller
niqc3 <- n.greater
niqc4 <- niqc - niqc1 -niqc2 -niqc3
if(constraint != "none") {
if(constraint == "convex" || constraint == "concave") {
if(lambda == 0) # z2 not yet above
z2 <- .splBasis(ord = ks, knots, ncoef, xo = knots[1:niqc4],
derivs = rep(2, niqc4))
idx3 <- cbind(rep(nrleq+niqc1 + 1:niqc4, rep(ks,niqc4)),
c(outer(1:ks,z2\$offsets,"+")))
X[idx3] <- if(constraint == "convex") z2\$design else -z2\$design
}
else if(constraint == "periodic") {
neqc3 <- 2
z2 <- .splBasis(ord = ks, knots, ncoef,
xo = sortx[c(1,nrq)], derivs = c(1,1))
idx2 <- cbind(rep(1:neqc3, rep(ks,neqc3)),
c(outer(1:ks,z2\$offsets,"+")))
X.temp <- matrix(0,neqc3,nvar)
X.temp[idx2] <- z2\$design
X[nrleq,  ] <- X.temp[2,]-X.temp[1,]
X[nrleq-1,] <- X[ox[nrq],]-X[ox[1],]
}
else {
z3 <- .splBasis(ord = ks, knots, ncoef, xo = knots[1:niqc4],
derivs = rep(1, niqc4))
idx3 <- cbind(rep(nrleq+niqc1 + 1:niqc4, rep(ks,niqc4)),
c(outer(1:ks,z3\$offsets,"+")))
X[idx3] <- if(constraint == "increase") z3\$design else -z3\$design
}
}
if(n.smaller > 0) {
o.smaller <- order(smaller[,2])
smaller.o <- smaller[,2][o.smaller]
z3.1 <- .splBasis(ord = ks, knots, ncoef, xo = smaller.o)
idx3.1 <- cbind(rep(nrleq+niqc1+niqc4 + o.smaller,rep(ks,niqc2)),
c(outer(1:ks,z3.1\$offsets,"+")))
X[idx3.1] <- -z3.1\$design
}
if(n.greater > 0) {
o.greater <- order(greater[,2])
greater.o <- greater[,2][o.greater]
z3.2 <- .splBasis(ord = ks, knots, ncoef, xo = greater.o)
idx3.2 <- cbind(rep(nrleq+niqc1+niqc4+niqc2 + o.greater, rep(ks,niqc3)),
c(outer(1:ks,z3.2\$offsets,"+")))
X[idx3.2] <- z3.2\$design
}
##
## formulate the equality constraints
##
neqc1 <- n.equal
if(n.equal > 0) {
derivs <- rep(0, neqc1)
o.equal <- order(equal[,2])
equal.o <- equal[,2][o.equal]
z1.1 <- .splBasis(ord = ks, knots, ncoef, xo = equal.o)
idx1.1 <- cbind(rep(nrql1 + o.equal, rep(ks, neqc1)),
c(outer(1:ks, z1.1\$offsets,"+")))
X[idx1.1] <- z1.1\$design
}
z1.2 <- .splBasis(ord = ks, knots, ncoef, xo = gradient.o,
derivs = rep(1, neqc2))
idx1.2 <- cbind(rep(nrql1+neqc1 + o.gradient, rep(ks, neqc2)),
c(outer(1:ks, z1.2\$offsets,"+")))
X[idx1.2] <- z1.2\$design
}
return(X)
} ## loo.design
```

## Try the cobs99 package in your browser

Any scripts or data that you put into this service are public.

cobs99 documentation built on May 31, 2017, 4:39 a.m.