tsvd <- function(g, X, r){

#%  TSVD  Truncated SVD regularization.
#%
#%  Usage: 
#%         > tmpsvd <- tsvd(g, X, r)
#%         > tmpsvd
#%   or typing 
#%         > tmpsvd$f.r
#%         > tmpsvd$rss
#%         > tmpsvd$f.r.ss   
#%
#%    Given a vector g, a design matrix X, and a truncation 
#%    parameter r, 
#%           tmpsvd <- tsvd(g,X,r)
#%           
#%    returns the truncated SVD estimate (tmpsvd$f.r) of the 
#%    vector f in the linear regression model
#% 
#%            g = X*f + noise.
#%
#%    Also returned are the residual sum of squares (tmpsvd$rss) and 
#%    the sum of squares (tmpsvd$f.r.ss) of the elements of the 
#%    truncated SVD estimate f_r
#%
#%    If r is a vector of truncation parameters, the i-th column
#%    tmpsvd$f.r[,i] is the truncated SVD estimate for the truncation
#%    parameter r[i]; the i-th elements of tmpsvd$rss and tmpsvd$f.r.ss 
#%    are the associated residual sum of squares and estimate sum of squares.
#%
#%    Adapted from the TSVD routine in Per Christian Hansen's
#%    Regularization Toolbox.  
#%
#%    R version by Yajun Mei
  
  #  % Size of inputs
  n <- dim(X)[1]
  p <- dim(X)[2]


  q  <- min(n, p)
  nr <- length(r)

  #% Possible choice of truncation parameter?
  if ( min(r) < 0 | max(r) > q ){
	stop(paste("Impossible truncation parameter r."))
  }
  
  #% Initialize outputs
  f.r 	      <- matrix(0,p,nr)
  rss         <- numeric(nr)
  f.r.ss      <- numeric(nr)

  
  #% Compute SVD of X   
  
  tmp     <- svd(X)
  U          <- tmp$u
  s          <- tmp$d      #% vector of singular values
  V          <- tmp$v         

  
  #% "Fourier" coefficients in expansion of solution in terms of
  #% right singular vectors
  beta        <- t(U[, 1:q]) %*% g
  fc          <- beta / s      
  
  #% Treat each truncation parameter separately.
  for (j in 1:nr){
    k         <- r[j]         #% current truncation parameter
    if (k >= 2)  { f.r[,j]   <- V[,1:k] %*% fc[1:k] }
    if ( k == 1) {   # In this case, use "*" instead of "%*%"
        f.r[,j]  <- V[,1] * fc[1] 
    } 
    f.r.ss[j] <- sum(fc[1:k]^2)
    if (k < q)  { rss[j]    <- sum( beta[(k+1):q]^2 ) }	
  }

  #% In overdetermined case, add rss of least-squares problem
  if (n > p){ 
    rss       <- rss + sum((g - U[, 1:q] %*% beta)^2)
  }

  #%return(f.r,rss,f.r.ss)
  return(list(f.r=f.r,rss=rss,f.r.ss=f.r.ss))
}