.packageName <- "glmnet"
coef.glmnet=function(object,s=object$lambda,exact=FALSE,...)
  predict(object,s=s,type="coefficients",exact=exact)

cv.glmnet=function(x,y,...,nfolds=10,foldid,type=c("response","class")){
  type=match.arg(type)
###Fit the model once to get dimensions etc of output
  glmnet.object=glmnet(x,y,...)
  lambda=glmnet.object$lambda
  N=nrow(x)
  if(missing(foldid)) foldid=sample(rep(seq(nfolds),length=N)) else nfolds=max(foldid)
  predmat=predict(glmnet.object,x,lambda=lambda)
  nz=sapply(predict(glmnet.object,type="nonzero"),length)

  for(i in seq(nfolds)){
    which=foldid==i
    cvfit=glmnet(x[!which,],y[!which],lambda=lambda,...)
    predmat[which,]=predict(cvfit,x[which,],type=type)
  }
  ###What to do depends on the type and the model fit
 cvstuff= switch(class(glmnet.object)[[2]],
         elnet=cvelnet(y,predmat),
         lognet=cvlognet(y,predmat,type),
         multnet=stop("not yet implemented")
         )
 cvm=cvstuff$cvm
 cvsd=cvstuff$cvsd
  cvname=cvstuff$name
  
out=list(lambda=lambda,cvm=cvm,cvsd=cvsd,cvup=cvm+cvsd,cvlo=cvm-cvsd,nzero=nz,name=cvname)
  lamin=getmin(lambda,cvm,cvsd)
  obj=c(out,as.list(lamin))
  class(obj)="cv.glmnet"
  obj
}
cvelnet=function(y,predmat){
  N=length(y)
   cvraw=(y-predmat)^2
  cvm=apply(cvraw,2,mean)
  cvsd=sqrt(apply(cvraw,2,var)/N)
  list(cvm=cvm,cvsd=cvsd,name="Mean Squared Error")
}
cvlognet=function(y,predmat,type=c("response","class")){
  N=dim(predmat)[[1]]
  ##assume for now that y is binary
  y=as.numeric(as.factor(y))
   cvraw=switch(type,
     response =-2*((y==2)*log(predmat)+(y==1)*log(1-predmat)),
     class=y!=predmat
     )
  cvm=apply(cvraw,2,mean)
  cvsd=sqrt(apply(cvraw,2,var)/N)
  cvname=switch(type,response="Deviance",class="Misclassification Error")
  list(cvm=cvm,cvsd=cvsd,name=cvname)
}
error.bars <-
function(x, upper, lower, width = 0.02, ...)
{
	xlim <- range(x)
	barw <- diff(xlim) * width
	segments(x, upper, x, lower, ...)
	segments(x - barw, upper, x + barw, upper, ...)
	segments(x - barw, lower, x + barw, lower, ...)
	range(upper, lower)
}

fix.lam=function(lam){
  llam=log(lam)
  lam[1]=exp(2*llam[2]-llam[3])
  lam
}
getmin=function(lambda,cvm,cvsd){
  cvmin=min(cvm)
  idmin=cvm<=cvmin
  lambda.min=max(lambda[idmin])
  idmin=match(lambda.min,lambda)
  semin=(cvm+cvsd)[idmin]
  idmin=cvm<semin
  lambda.1se=max(lambda[idmin])
  list(lambda.min=lambda.min,lambda.1se=lambda.1se)
}
 
glmnet=function(x,y,family=c("gaussian","binomial","multinomial"),weights,alpha=1.0,nlambda=100,lambda.min=ifelse(nobs<nvars,5e-2,1e-4),lambda,standardize=TRUE,thresh=1e-4,dfmax=nvars+1,pmax=min(dfmax*1.2,nvars),exclude,penalty.factor=rep(1,nvars),maxit=100,HessianExact=FALSE,type=c("covariance","naive")){
  family=match.arg(family)
  this.call=match.call()
  nlam=as.integer(nlambda)
  np=dim(x)
  nobs=as.integer(np[1])
  if(missing(weights))weights=rep(1,nobs)
  nvars=as.integer(np[2])
  vnames=colnames(x)
  if(is.null(vnames))vnames=paste("V",seq(nvars),sep="")
  if(family %in% c("binomial","multinomial")){
    nc=dim(y)
    maxit=as.integer(maxit)
    kopt=as.integer(HessianExact)
    if(is.null(nc)){
      ## Need to construct a y matrix, and include the weights
      y=as.factor(y)
      ntab=table(y)
      classnames=names(ntab)
      nc=as.integer(length(ntab))
      y=diag(nc)[as.numeric(y),]
    }
    else{
      noo=nc[1]
      if(noo!=nobs)stop("x and y have different number of rows")
      nc=as.integer(nc[2])
      classnames=colnames(y)
    }
    if(family=="binomial"){
      if(nc>2)stop("More than two classes; use multinomial family instead")
      nc=as.integer(1) # for calling multnet
    }
    if(!missing(weights)){
      ### check if any are zero
      o=weights>0
      if(!all(o)){ #subset the data
        y=y[o,]
        x=x[o,,drop=FALSE]
        weights=weights[o]
        nobs=sum(o)
      }
      y=y*weights
    }
     ### Compute the null deviance
    prior=apply(y,2,sum)
    sumw=sum(y)
    prior=prior/sumw
    nulldev= -2*sum(y*outer(rep(1,nobs),log(prior),"*"))/sumw

  }
else
     {
       weights=as.double(weights)
       ### compute the null deviance
       ybar=weighted.mean(y,weights)
       nulldev=weighted.mean( (y-ybar)^2,weights)
       type=match.arg(type)

       ka=as.integer(switch(type,
         covariance=1,
         naive=2,
         ))
     }
 storage.mode(y)="double"

    ne=as.integer(dfmax)
    nx=as.integer(pmax)
     if(!missing(exclude)){
       jd=match(exclude,seq(nvars),0)
       if(!all(jd>0))stop("Some excluded variables out of range")
       jd=as.integer(c(length(jd),jd))
     }else jd=as.integer(0)
    vp=as.double(penalty.factor)
    isd=as.integer(standardize)
    thresh=as.double(thresh)
    if(missing(lambda)){
      if(lambda.min>=1)stop("lambda.min should be less than 1")
      flmin=as.double(lambda.min)
      ulam=double(1)
    }
    else{
      flmin=as.double(1)    
      if(any(lambda<0))stop("lambdas should be non-negative")
      ulam=as.double(rev(sort(lambda)))
      nlam=as.integer(length(lambda))
    }
    is.sparse=FALSE
    if(class(x)=="dgCMatrix"){##Sparse case
      is.sparse=TRUE
      ix=as.integer(x@p+1)
      jx=as.integer(x@i+1)
      x=as.double(x@x)

      if(family=="gaussian")
              fit=.Fortran("spelnet",
        ka,parm=alpha,nobs,nvars,x,ix,jx,y,weights,jd,vp,ne,nx,nlam,flmin,ulam,thresh,isd,
        lmu=integer(1),
        a0=double(nlam),
        ca=double(nx*nlam),
        ia=integer(nx),
        nin=integer(nlam),
        rsq=double(nlam),
        alm=double(nlam),
        nlp=integer(1),
        jerr=integer(1),PACKAGE="glmnet"
        )

        else
                fit=.Fortran("splognet",
        parm=alpha,nobs,nvars,nc,x,ix,jx,y,jd,vp,ne=ne,nx,nlam,flmin,ulam,thresh,isd,maxit,kopt,
        lmu=integer(1),
        a0=double(nlam*nc),
        ca=double(nx*nlam*nc),
        ia=integer(nx),
        nin=integer(nlam),
        dev=double(nlam),
        alm=double(nlam),
        nlp=integer(1),
        jerr=integer(1),PACKAGE="glmnet"
        )

    }## end of sparse case
  else
     {
       ##regular nonsparse case
       if(family=="gaussian")
         fit=.Fortran("elnet",
          ka,parm=alpha,nobs,nvars,as.double(x),y,weights,jd,vp,ne,nx,nlam,flmin,ulam,thresh,isd,
          lmu=integer(1),
          a0=double(nlam),
          ca=double(nx*nlam),
          ia=integer(nx),
          nin=integer(nlam),
          rsq=double(nlam),
          alm=double(nlam),
          nlp=integer(1),
          jerr=integer(1),PACKAGE="glmnet"
          )
     else
               fit=.Fortran("lognet",
          parm=alpha,nobs,nvars,nc,as.double(x),y,jd,vp,ne,nx,nlam,flmin,ulam,thresh,isd,maxit,kopt,
          lmu=integer(1),
          a0=double(nlam*nc),
          ca=double(nx*nlam*nc),
          ia=integer(nx),
          nin=integer(nlam),
          dev=double(nlam),
          alm=double(nlam),
          nlp=integer(1),
          jerr=integer(1),PACKAGE="glmnet"
          )
     }## end regular nonsparse case
     lmu=fit$lmu
     nin=fit$nin[seq(lmu)]
     ninmax=max(nin)
     lam=fit$alm[seq(lmu)]
     if(missing(lambda))lam=fix.lam(lam)##first lambda is infinity; changed to entry point
  stepnames=paste("s",seq(lmu)-1,sep="")

     errmsg=jerr(fit$jerr,maxit,pmax)### error messages from fortran
     switch(paste(errmsg$n),
            "1"=stop(errmsg$msg,call.=FALSE),
            "-1"=warning(errmsg$msg,call.=FALSE)
            )
if(family=="multinomial"){
      beta.list=as.list(seq(nc))
      names(beta.list)=classnames
      a0=matrix(fit$a0[seq(lmu*nc)],nc,lmu,dimnames=list(classnames,stepnames))
      a0=scale(a0,TRUE,FALSE)
      attr(a0,"scaled:center")=NULL
      dfmat=a0
      dd=c(nvars,lmu)
      if(ninmax>0){
        ca=array(fit$ca[seq(nx*lmu*nc)],c(nx,nc,lmu))[seq(ninmax),,,drop=FALSE]
        ja=fit$ia[seq(ninmax)]#confusing but too hard to change
        oja=order(ja)
        ja=rep(ja[oja],lmu)
        ia=cumsum(c(1,rep(ninmax,lmu)))
        df=apply(abs(ca)>0,c(1,3),any)
        df=apply(df,2,sum)
        for(k in seq(nc)){
          cak=ca[oja,k, ,drop=FALSE]
          dfmat[k,]=apply(abs(cak)>0,2,sum)
          beta=new("dgCMatrix",Dim=dd,Dimnames=list(vnames,stepnames),x=as.vector(cak),p=as.integer(ia-1),i=as.integer(ja-1))
          beta.list[[k]]=beta
        }
   } else{
    for (k in seq(nc)) {
      dfmat[k, ] = rep(0,lmu)
      beta.list[[k]] = zeromat(nvars,lmu,vnames,stepnames)
    }
  }
       
    outlist=list(a0=a0,beta=beta.list,
         dev=fit$dev[seq(lmu)],nulldev=nulldev,dfmat=dfmat,df=df,
         lambda=lam,npasses=fit$nlp,jerr=fit$jerr,dim=dd,call=this.call)
    class(outlist)=c("glmnet","multnet")
    }
     else{
           dd=c(nvars,lmu)
    if(ninmax>0){
           ca=matrix(fit$ca[seq(nx*lmu)],nx,lmu)[seq(ninmax),,drop=FALSE]
           df=apply(abs(ca)>0,2,sum)
           ja=fit$ia[seq(ninmax)]#confusing but too hard to change
           oja=order(ja)
           ja=rep(ja[oja],lmu)
           ia=cumsum(c(1,rep(ninmax,lmu)))
           beta=new("dgCMatrix",Dim=dd,Dimnames=list(vnames,stepnames),x=as.vector(ca[oja,]),p=as.integer(ia-1),i=as.integer(ja-1))
         }else {
           beta = zeromat(nvars,lmu,vnames,stepnames)
           df=rep(0,lmu)
         }

      if(family=="binomial"){
          a0=-fit$a0[seq(lmu)]
          names(a0)=stepnames
           outlist=list(a0=a0,beta=-beta,#sign flips make 2 target class
         dev=fit$dev[seq(lmu)],nulldev=nulldev,df=df,
         lambda=lam,npasses=fit$nlp,jerr=fit$jerr,dim=dd,call=this.call)
    class(outlist)=c("glmnet","lognet")
         }
           else
             {
               a0=fit$a0[seq(lmu)]
               names(a0)=stepnames
                  outlist= list(a0=a0,beta=beta,dev=fit$rsq[seq(lmu)],nulldev=nulldev,df=df,
         lambda=lam,npasses=fit$nlp,jerr=fit$jerr,dim=dd,call=this.call)
    class(outlist)=c("glmnet","elnet")
                }
         }
     outlist
   }
jerr=function(n,maxit,pmax){
  if(n==0) msg=""
  if(n>0){#fatal error
    if(n<7777)msg="Memory allocation error; contact package maintainer"
    if(n==7777)msg="All used predictors have zero variance"
    if((8000<n) & (n<9000))msg=paste("Null probability for class",n-8000, "< 1.0e-5")
    if((9000<n) & (n<10000))msg=paste("Null probability for class",n-9000, "> 1.0 - 1.0e-5")
    if(n==10000)msg="All penalty factors are <= 0"
    n=1
  msg=paste("in glmnet fortran code -",msg)
  }
  if(n<0){#non fatal error
    if(n>-10000)msg=paste("Convergence for ",-n,"th lambda value not reached after maxit=",maxit," iterations; solutions for larger lambdas returned",sep="")
    if(n < -10000)msg=paste("Number of nonzero coefficients along the path exceeds pmax=",pmax, " at ",-n-10000,"th lambda value; solutions for larger lambdas returned",sep="")
    n=-1
  msg=paste("from glmnet fortran code -",msg)
  }
  list(n=n,msg=msg)
}
                  
                  
     
lambda.interp=function(lambda,s){
### lambda is the index sequence that is produced by the model
### s is the new vector at which evaluations are required.
### the value is a vector of left and right indices, and a vector of fractions.
### the new values are interpolated bewteen the two using the fraction
### Note: lambda decreases. you take:
### sfrac*left+(1-sfrac*right)
  
  if(length(lambda)==1){# degenerate case of only one lambda
    nums=length(s)
    left=rep(1,nums)
    right=left
    sfrac=rep(1,nums)
  }
  else{
      s[s > max(lambda)] = max(lambda)
      s[s < min(lambda)] = min(lambda)
      k=length(lambda)
      sfrac <- (lambda[1]-s)/(lambda[1] - lambda[k])
      lambda <- (lambda[1] - lambda)/(lambda[1] - lambda[k])
      coord <- approx(lambda, seq(lambda), sfrac)$y
      left <- floor(coord)
      right <- ceiling(coord)
      sfrac=(sfrac-lambda[right])/(lambda[left] - lambda[right])
      sfrac[left==right]=1
    }
list(left=left,right=right,frac=sfrac)
}
nonzeroCoef=function(beta,bystep=FALSE){
  ##beta should be in "dgCMatrix" format
  beta=t(beta)
  which=diff(beta@p)
  which=seq(which)[which>0]
  if(bystep){
    nzel=function(x,which)if(any(x))which[x]else NULL
    beta=abs(as.matrix(beta[,which]))>0
    apply(beta,1,nzel,which)
  }
  else which
    
  }
plot.cv.glmnet=function(x,sign.lambda=1,...){
  cvobj=x
  xlab="log(Lambda)"
  if(sign.lambda<0)xlab=paste("-",xlab,sep="")
  plot(sign.lambda*log(cvobj$lambda),cvobj$cvm,ylim=range(cvobj$cvup,cvobj$cvlo),xlab=xlab,ylab=cvobj$name,type="n")
error.bars(sign.lambda*log(cvobj$lambda),cvobj$cvup,cvobj$cvlo,width=0.01,col="darkgrey")
  points(sign.lambda*log(cvobj$lambda),cvobj$cvm,pch=20,col="red")
axis(side=3,at=sign.lambda*log(cvobj$lambda),labels=paste(cvobj$nz),tick=FALSE,line=0)
abline(v=sign.lambda*log(cvobj$lambda.min),lty=3)
abline(v=sign.lambda*log(cvobj$lambda.1se),lty=3)
  invisible()
}
plot.glmnet=function(x, xvar=c("norm","lambda","dev"),label=FALSE,...){
  ocl=class(x)[[2]]
  xvar=match.arg(xvar)
  if(ocl=="multnet"){
    beta=x$beta
    if(xvar=="norm"){
      cnorm1=function(beta){
        which=nonzeroCoef(beta)
        beta=as.matrix(beta[which,])
        apply(abs(beta),2,sum)
      }
      norm=apply(sapply(x$beta,cnorm1),1,sum)
    } else norm = NULL
    dfmat=x$dfmat
    ncl=nrow(dfmat)
    clnames=rownames(dfmat)
    for( i in seq(ncl)){
      plotCoef(beta[[i]],norm,x$lambda,dfmat[i,],x$dev,label=label,xvar=xvar,ylab=paste("Coefficients: Class",clnames[i]),...)
    }
  }else plotCoef(x$beta,lambda=x$lambda,df=x$df,dev=x$dev,label=label,xvar=xvar,...)
}
      
      
plotCoef=function(beta,norm,lambda,df,dev,label=FALSE,xvar=c("norm","lambda","dev"),xlab=iname,ylab="Coefficients",...){
  ##beta should be in "dgCMatrix" format
  which=nonzeroCoef(beta)
  beta=as.matrix(beta[which,])
  xvar=match.arg(xvar)
  switch(xvar,
    "norm"={
      index=if(missing(norm))apply(abs(beta),2,sum)else norm
      iname="L1 Norm"
    },
    "lambda"={
      index=log(lambda)
      iname="Log Lambda"
    },
    "dev"= {
      index=dev
      iname="Fraction Deviance Explained"
    }
         )
  dotlist=list(...)
  type=dotlist$type
  if(is.null(type))
    matplot(index,t(beta),lty=1,xlab=xlab,ylab=ylab,type="l",...)
  else matplot(index,t(beta),lty=1,xlab=xlab,ylab=ylab,...)
  atdf=pretty(index)
 prettydf=trunc(approx(x=index,y=df,xout=atdf,rule=2)$y)
 axis(3,at=atdf,label=prettydf,cex.axis=.5,tcl=NA)
 if(label){
   nnz=length(which)
   xpos=max(index)
   pos=4
   if(xvar=="lambda"){
     xpos=min(index)
     pos=2
   }
   xpos=rep(xpos,nnz)
   ypos=beta[,ncol(beta)]
   text(xpos,ypos,paste(which),cex=.5,pos=pos)
 }
         
}
    
predict.elnet=function(object,newx,s=object$lambda,type=c("link","response","coefficients","nonzero"),exact=FALSE,...){
  type=match.arg(type)
  a0=t(as.matrix(object$a0))
  rownames(a0)="(Intercept)"
  nbeta=rbind2(a0,object$beta)
  if(!missing(s)){
    vnames=dimnames(nbeta)[[1]]
    dimnames(nbeta)=list(NULL,NULL)
    lambda=object$lambda
    lamlist=lambda.interp(lambda,s)
    nbeta=nbeta[,lamlist$left,drop=FALSE]*lamlist$frac +nbeta[,lamlist$right,drop=FALSE]*(1-lamlist$frac)
    dimnames(nbeta)=list(vnames,paste(seq(along=s)))
  }
  switch(type,
         "coefficients"=nbeta,
         "link"=as.matrix(cbind2(1,newx)%*%nbeta),
         "response"=as.matrix(cbind2(1,newx)%*%nbeta),
         "nonzero"=nonzeroCoef(nbeta[-1,,drop=FALSE],bystep=TRUE)
         )
}  
    
predict.glmnet=function(object,newx,s=object$lambda,type=c("link","response","coefficients","class","nonzero"),exact=FALSE,...)
NextMethod("predict")
predict.lognet=function(object,newx,s=object$lambda,type=c("link","response","coefficients","class","nonzero"),exact=FALSE,...){
  type=match.arg(type)
  a0=t(as.matrix(object$a0))
  rownames(a0)="(Intercept)"
  nbeta=rbind2(a0,object$beta)
  if(!missing(s)){
    vnames=dimnames(nbeta)[[1]]
    dimnames(nbeta)=list(NULL,NULL)
   lambda=object$lambda
    lamlist=lambda.interp(lambda,s)
    nbeta=nbeta[,lamlist$left,drop=FALSE]*lamlist$frac +nbeta[,lamlist$right,drop=FALSE]*(1-lamlist$frac)
    dimnames(nbeta)=list(vnames,paste(seq(along=s)))
  }
  ### remember that although the fortran lognet makes predictions
  ### for the first class, we make predictions for the second class
  ### to avoid confusion with 0/1 responses.
  ### glmnet flipped the signs of the coefficients 
  if(type=="coefficients")return(nbeta)
  if(type=="nonzero")return(nonzeroCoef(nbeta[-1,,drop=FALSE],bystep=TRUE))
  nfit=as.matrix(cbind2(1,newx)%*%nbeta)
  switch(type,
         response={
           pp=exp(-nfit)
           1/(1+pp)
         },
         link=nfit,
         class=ifelse(nfit>0,2,1)
         )
}  
    
predict.multnet=function(object,newx,s=object$lambda,type=c("link","response","coefficients","class","nonzero"),exact=FALSE,...){
  type=match.arg(type)
  a0=object$a0
  rownames(a0)=rep("(Intercept)",nrow(a0))
  nbeta=object$beta
  klam=dim(a0)
  nclass=klam[[1]]
  nlambda=length(s)
  if(!missing(s)){
    lambda=object$lambda
    lamlist=lambda.interp(lambda,s)
    for(i in seq(nclass)){
      kbeta=rbind2(a0[i,,drop=FALSE],nbeta[[i]])
      vnames=dimnames(kbeta)[[1]]
    dimnames(kbeta)=list(NULL,NULL)
      kbeta=kbeta[,lamlist$left,drop=FALSE]*lamlist$frac +kbeta[,lamlist$right,drop=FALSE]*(1-lamlist$frac)
     dimnames(kbeta)=list(vnames,paste(seq(along=s)))
     nbeta[[i]]=kbeta
    }
  }
  else
    {
      for(i in seq(nclass))nbeta[[i]]=rbind2(a0[i,],nbeta[[i]])
    }
if(type=="coefficients")return(nbeta)
if(type=="nonzero")return(
     lapply(nbeta,function(x)nonzeroCoef(x[-1,,drop=FALSE],bystep=TRUE)))

  dd=dim(newx)
  npred=dd[[1]]
  dn=list(names(nbeta),dimnames(nbeta[[1]])[[2]],dimnames(newx)[[1]])
  dp=array(0,c(nclass,nlambda,npred),dimnames=dn)
  for(i in seq(nclass)){
    fitk=cbind2(1,newx)%*%(nbeta[[i]])
    dp[i,,]=dp[i,,]+as.matrix(t(fitk))
  }
  switch(type,
         response={
           pp=exp(dp)
           psum=apply(pp,c(2,3),sum)
           aperm(pp/rep(psum,rep(nclass,nlambda*npred)),c(3,1,2))
         },
         link=aperm(dp,c(3,1,2)),
         class={
           dp=aperm(dp,c(3,1,2))
           apply(dp,3,softmax)
         }
         )
}  
    
print.glmnet=function(x,digits = max(3, getOption("digits") - 3),...){
      cat("\nCall: ", deparse(x$call), "\n\n")

  print(cbind(Df=x$df,"%Dev"=signif(x$dev,digits),Lambda=signif(x$lambda,digits)))
    }
`softmax` <-
function (x, gap = FALSE) 
{
    d <- dim(x)
    maxdist <- x[, 1]
    pclass <- rep(1, d[1])
    for (i in seq(2, d[2])) {
        l <- x[, i] > maxdist
        pclass[l] <- i
        maxdist[l] <- x[l, i]
    }
    dd <- dimnames(x)[[2]]
    if (gap) {
        x <- abs(maxdist - x)
        x[cbind(seq(d[1]), pclass)] <- drop(x %*% rep(1, d[2]))
        gaps <- do.call("pmin", data.frame(x))
    }
    pclass <- if (is.null(dd) || !length(dd)) 
        pclass
    else factor(pclass, levels = seq(d[2]), labels = dd)
    if (gap) 
        list(class = pclass, gaps = gaps)
    else pclass
}
zeromat=function(nvars,lmu,vnames,stepnames){
  ca=rep(0,lmu)
  ia=seq(lmu+1)
  ja=rep(1,lmu)
  dd=c(nvars,lmu)
  new("dgCMatrix", Dim = dd,
      Dimnames = list(vnames,stepnames),
      x = as.vector(ca),
      p = as.integer(ia - 1), i = as.integer(ja - 1))
}
