abcd.lasso<-
function(X,Y,l1,c)
{
  require(quadprog)

  d  <-as.matrix(c(-t(X)%*%Y+l1*c,t(X)%*%Y+l1*c))
  b0 <-rep(0,ncol(X))
  A  <-diag(rep(1,ncol(X)))
  XX <-t(X)%*%X
  D  <-rbind(cbind(XX,-XX),cbind(-XX,XX))
  st <-solve.QP(Dmat=D, dvec=d, Amat=A, bvec=b0, meq=0, 
                factorized=FALSE)$solution

  return(st[1:ncol(X)]-st[(1:ncol(X))+ncol(X)])
}

abcd.augfit<-
function(l1,l2,sigma,abcd.obj)
{
  P<-abcd.obj$P
  Z<-abcd.obj$Z
  R<-abcd.obj$R
  W<-abcd.obj$W
  K<-chol(sigma*abcd.obj$Bi)
  Y<-abcd.obj$Y
  n<-abcd.obj$n

  Y.star <-as.matrix(c(Y,rep(0,nrow(K)+nrow(R))))
  X.star <-rbind(cbind(P,Z,R,W),
                 cbind(matrix(0,ncol=ncol(P),nrow=nrow(K)),K/sqrt(n),matrix(0,ncol=ncol(R)+ncol(W),nrow=nrow(K))),
                 cbind(matrix(0,ncol=ncol(P)+ncol(K),nrow=nrow(R)),R*sqrt(l2),matrix(0,ncol=ncol(R)+ncol(W),nrow=nrow(K)))
                 )
  c.vec  <-rep(c(1,0),times=c(ncol(P),ncol(Z)+ncol(R)+ncol(W)))

  return(abcd.lasso(X=X.star,Y=Y.star,l1=l1,c=c))
}

abcd.crit<-
function(...)
{
  #search through penalty space cleverly.  

  #output the final penalties, model
}



abcd.fit<-
function(...)
{
  require(gss)

  #determine whether to switch to scratch or use an existing fit.
}

abcd.fit.obj<-
function(abcd.obj,l1,l2,sigma)
{

}


#formula for spline part,
#partial is a matrix
#random is a matrix or a factor

abcd.fit.scratch<-
function(formula, data, partial,random,criterion=NULL,l1,l2,sigma)
{
  # formula part
    cl <- match.call()
    mf <- match.call(expand.dots = FALSE)
    m  <- match(c("formula", "data"), names(mf), 0)
    mf <- mf[c(1, m)]
    mf$drop.unused.levels <- TRUE
    mf[[1]] <- as.name("model.frame")
    mf <- eval(mf, parent.frame())
    y  <- model.response(mf, "numeric")
    Y  <- as.matrix(y)
    n  <- nrow(Y)

  # Partial basis part
   if(missing(P)){
     stop("Must specify partial basis")
   }else{
     P    <- as.matrix(partial)
   }

  # random part
   if( is.matrix(random)){
      Z <- random
   }else{ if(is.factor(random)){
      Z <- as.matrix(model.matrix(~random))
   }}
   if(!missing(Z)){ Bi <- solve(Z%*t(Z)) }else{stop("Unrecognized random input")}

   #req for spline part
    rand.parts<-mkran(~1|random,data=data.frame(random=random))

  # Spline parts
   suppressWarnings(
     gss.int <- ssanova(formula=formula,data=data,type="cubic",rand=rand.parts))
   ext.gss <- extract.gss(gss.int)

   R    <- ext.gss$R
   W    <- ext.gss$S

  # criterion
   if(is.null(criterion)){
     crit.fun<-abcd.crit
   }else{
     crit.fun<-match.fun(criterion)
     # functions input here should take the existing forms and find the best
     #  spot based on whatever criterion they evaluate.
   }

   if(is.null(criterion)&any(c(missing(l1),missing(l2),missing(sigma))))
     stop("Must specify either a criterion function or smoothing parameters")
   
  # create an abcd.obj
   abcd.parts <- list(P=P,Z=Z,R=R,W=W,Y=Y,n=n,Bi=Bi)

  # call criterion function
   pen.solved <- crit.fun(l1,l2,sigma,abcd.parts)
}  



extract.gss<-
function(object)
{
  include<-object$terms$labels
  newdata<-object$mf
    nnew <- nrow(newdata)
    nbasis <- length(object$id.basis)
    nnull <- length(object$d)
    nz <- length(object$b)
    nn <- nbasis + nnull + nz
    term <- object$terms
    philist <- rklist <- NULL
    s <- r <- NULL
    nq <- 0
    for (label in include) {
        if (label == "1") {
            philist <- c(philist, term[[label]]$iphi)
            s <- cbind(s, rep(1, len = nnew))
            next
        }
        if (label == "partial") 
            next
        if (label == "offset") 
            next
        xnew <- newdata[, term[[label]]$vlist]
        x <- object$mf[object$id.basis, term[[label]]$vlist]
        nphi <- term[[label]]$nphi
        nrk <- term[[label]]$nrk
        if (nphi) {
            iphi <- term[[label]]$iphi
            phi <- term[[label]]$phi
            for (i in 1:nphi) {
                philist <- c(philist, iphi + (i - 1))
                s <- cbind(s, phi$fun(xnew, nu = i, env = phi$env))
            }
        }
        if (nrk) {
            irk <- term[[label]]$irk
            rk <- term[[label]]$rk
            for (i in 1:nrk) {
                rklist <- c(rklist, irk + (i - 1))
                nq <- nq + 1
                r <- array(c(r, rk$fun(xnew, x, nu = i, env = rk$env, 
                  out = TRUE)), c(nnew, nbasis, nq))
            }
        }
    }
    if (any(include == "partial")) {
        nphi <- term$partial$nphi
        iphi <- term$partial$iphi
        for (i in 1:nphi) philist <- c(philist, iphi + (i - 1))
        s <- cbind(s, newdata$partial)
    }
    r.wk <- matrix(0, nnew, nbasis)
    nq <- 0
    for (i in rklist) {
        nq <- nq + 1
        r.wk <- r.wk + 10^object$theta[i] * r[, , nq]
    }
    R<-r.wk
    S<-s

    return(list(R=R,S=S))
}