monotonic GAM with more than one term

Hi,

Does anyone know how to fit a GAM where one or more smooth terms are
constrained to be monotonic, in the presence of "by" variables or
other terms? I looked at the example in ?pcls but so far have not been
able to adapt it to the case where there is more than one predictor.
For example,

require(mgcv)

set.seed(0)
n<-100


# Generate data from a monotonic truth.
x<-runif(100)*4-1
x<-sort(x)

m <- structure(rep(1:2,50), .Label=c("one","two"), class="factor")

f<-as.integer(m) * exp(4*x)/(1+exp(4*x))
y<-f+rnorm(100)*0.1
plot(x,y,xlim=c(min(x), max(x)*2))
dat<-data.frame(x=x,m=m,y=y)



# Show regular spline fit (and save fitted object)
f.ug<-gam(y~m+s(x,k=10,by=m,bs="cr"))

bool <- m=="one"
yhat <- fitted(f.ug)

lines(x[bool],yhat[bool])
lines(x[!bool],yhat[!bool])

xx <- seq(max(x), 2*max(x), length=100)
mm <- structure(rep(1:2,50), .Label=c("one","two"), class="factor")
yy <- predict(f.ug, newdata=data.frame(m=mm,x=xx))

bool <- mm=="one"

lines(xx[bool],yy[bool], lty=2) # show effect of extrapolation
lines(xx[!bool],yy[!bool], lty=2)


# this is where I start running into trouble
   sm<-smoothCon(s(x,k=10,by=m,bs="cr"),dat,knots=NULL)[[1]]
   FF<-mono.con(sm$xp);   # get constraints
   G<-list(y=y,
           w=rep(1, n),
           X=sm$X,
           C=matrix(0,0,0),
           S = sm$S,
           off = 0,
           sp=f.ug$sp,
           p=sm$xp,
           Ain = FF$A,
           bin = FF$b
           )


p<-pcls(G)  # fit spline (using s.p. from unconstrained fit)
fv<-Predict.matrix(sm,data.frame(x=x))%*%p

# can we do this without calling smoothCon directly ?
# also having trouble here.
   f.nofit<-gam(y~m+s(x,k=10,by=m,bs="cr"),fit=FALSE)
   FF2 <- mono.con(f.nofit$smooth[[1]]$xp)

   stopifnot(identical(FF, FF2))

   G2 <- list(y = f.nofit$y,
              w = f.nofit$w,
              X = f.nofit$X,
              C = f.nofit$C,
              S = f.nofit$smooth[[1]]$S,
              off = f.nofit$off,
              sp = f.ug$sp,
              p = f.nofit$smooth[[1]]$xp,
              Ain = FF2$A,
              bin = FF2$b
              )

p2 <- pcls(G2)
fv2<-Predict.matrix(f.nofit$smooth[[1]],data.frame(x=x))%*%p2



Many thanks
-Ben