swh:1:snp:16c54c84bc54885e783d4424d714e5cc82f479a1
Tip revision: 61e5d2e598fb5e53345fa9414a615a4d7b524808 authored by Roger Koenker on 18 October 2020, 10:50:03 UTC
version 5.74
version 5.74
Tip revision: 61e5d2e
sfn.R
#WARNING -- needs considerable work to finish implementing sfn.control fix.
#-----------------------------------------------------------------------------
# Storage parameters and other controls for sfn fitting
# nsubmax = upper bound of dimension of lindx
# tmpmax = upper bound for the dimension of tmpvec
# nnzlmax = maximum number of non-zero elements in L
# cachsz = size of the cache memory; it's machine dependent
# small = the tolerance used to check for convergence
# maxiter = allowed limit for sfn iterations
# warn.mesg = flag to control printing of warnings
#
sfn.control <- function( nsubmax = NULL, tmpmax = NULL, nnzlmax = NULL,
cachsz = 64, small = 1e-6, maxiter=100, warn.mesg=TRUE)
list(nsubmax = nsubmax, tmpmax = tmpmax, nnzlmax = nnzlmax, cachsz = cachsz,
small = small, maxiter = maxiter, warn.mesg = warn.mesg)
#################################################################
# Interface for a sparse implementation of LMS interior point method
# a = structure of the design matrix A' stored in csr format
# y = pseudo response vector
# tau = the desired regression quantile
# rhs = the rhs of the dual problem -- specify at your own risk
#################################################################
rq.fit.sfn <- function(a,y,tau=.5, rhs = (1-tau)*c(t(a) %*% rep(1,length(y))), control)
{
y <- -y
n <- length(y)
m <- a@dimension[2]
if(n != a@dimension[1])
stop("Dimensions of design matrix and the response vector not compatible")
u <- rep(1,length=n)
x <- rep((1-tau),length=n)
nnzdmax <- nnza <- a@ia[n+1]-1
iwmax <- 7*m+3
ao <- t(a)
e <- ao %*% a
nnzemax <- e@ia[m+1]-1
ctrl <- sfn.control()
if (!missing(control)) {
control <- as.list(control)
ctrl[names(control)] <- control
}
nsubmax <- ctrl$nsubmax
tmpmax <- ctrl$tmpmax
nnzlmax <- ctrl$nnzlmax
if (is.null(ctrl$nsubmax)) ctrl$nsubmax <- nsubmax <- nnzemax
if (is.null(ctrl$tmpmax)) ctrl$tmpmax <- tmpmax <- 6 * m
if (is.null(ctrl$nnzlmax)) ctrl$nnzlmax <- nnzlmax <- 4 * nnzdmax
wwm <- vector("numeric",3*m)
s <- u - x
b1 <- solve(e, ao %*% y, tmpmax=tmpmax,nnzlmax=nnzlmax,nsubmax=nsubmax)
r <- y - a %*% b1
z <- ifelse(abs(r)<ctrl$small,(r*(r>0)+ctrl$small),r*(r>0))
w <- z - r
wwn <- matrix(0,n,14)
wwn[,1] <- r
wwn[,2] <- z
wwn[,3] <- w
fit <- .Fortran("srqfn",
n = as.integer(n),
m = as.integer(m),
nnza = as.integer(nnza),
a = as.double(a@ra),
ja = as.integer(a@ja),
ia = as.integer(a@ia),
ao = as.double(ao@ra),
jao = as.integer(ao@ja),
iao = as.integer(ao@ia),
nnzdmax = as.integer(nnzdmax),
d = double(nnzdmax),
jd = integer(nnzdmax),
id = integer(m+1),
dsub = double(nnzemax+1),
jdsub = integer(nnzemax+1),
nnzemax = as.integer(nnzemax),
e = as.double(e@ra),
je = as.integer(e@ja),
ie = as.integer(e@ia),
nsubmax = as.integer(nsubmax),
lindx = integer(nsubmax),
xlindx = integer(m+1),
nnzlmax = as.integer(nnzlmax),
lnz = double(nnzlmax),
xlnz = integer(m+1),
iw = integer(m*5),
iwmax = as.integer(iwmax),
iwork = integer(iwmax),
xsuper = integer(m+1),
tmpmax = as.integer(tmpmax),
tmpvec = double(tmpmax),
wwm = as.double(wwm),
wwn = as.double(wwn),
cachsz = as.integer(ctrl$cachsz),
level = as.integer( 8 ),
x = as.double(x),
s = as.double(s),
u = as.double(u),
c = as.double(y),
sol = as.double(b1),
rhs = as.double(rhs),
small = as.double(ctrl$small),
ierr = integer(1),
maxiter = as.integer(ctrl$maxiter),
time = double(7))[c("sol","ierr","maxiter","time")]
ierr <- fit$ierr
if(!(ierr==0) && ctrl$warn.mesg)
warning(sfnMessage(ierr))
coefficients <- -fit$sol
residuals <- -y - a %*% coefficients
list(coefficients = coefficients,
residuals = residuals,
control = ctrl,
ierr = ierr,
it = fit$maxiter)
}
#------------------------------------------------------------------------------
#################################################################
# Interface for a sparse implementation of LMS interior point method
# x = structure of A1' stored in csr format
# y = pseudo response vector
# m = column dimension of A1' in full mode
# R = structure of the constraint matrix A2' stored in csr format
# r = rhs of the inequality constraints
# tau = the desired regression quantile
# rhs = the rhs of the dual problem -- specify at your own risk
#################################################################
rq.fit.sfnc <- function(x, y, R, r, tau = 0.5,
rhs = (1-tau)*c(t(x) %*% rep(1,length(y))),control)
{
y <- -y
r <- -r
n1 <- length(y)
m <- x@dimension[2]
if(n1 != x@dimension[1])
stop("The design matrix A1' and response vector y are not compatible")
n2 <- length(r)
if(n2 != R@dimension[1])
stop("The constraint matrix A2' and constraint rhs are not compatible")
maxn1n2 <- max(n1,n2)
u <- rep(1,length=n1)
x1 <- rep(1-tau,length=n1)
x2 <- rep(1,length=n2)
wwm <- vector("numeric",6*m)
wwm[1:m] <- rhs
nnzx <- x@ia[x@dimension[1]+1]-1
nnzR <- R@ia[R@dimension[1]+1]-1
nnzdmax <- max(nnzx,nnzR)
iwmax <- 7*m+3
ao1 <- t(x)
ao2 <- t(R)
e <- ao1 %*% x
g <- ao2 %*% R
h <- e + g
nnzemax <- e@ia[e@dimension[1]+1]-1
nnzgmax <- g@ia[g@dimension[1]+1]-1
nnzhmax <- h@ia[h@dimension[1]+1]-1
ctrl <- sfn.control()
if (!missing(control)) {
control <- as.list(control)
ctrl[names(control)] <- control
}
nsubmax <- ctrl$nsubmax
tmpmax <- ctrl$tmpmax
nnzlmax <- ctrl$nnzlmax
if (is.null(ctrl$nsubmax)) nsubmax <- nnzhmax
if (is.null(ctrl$tmpmax)) tmpmax <- 6 * m
if (is.null(ctrl$nnzlmax)) nnzlmax <- 4 * nnzdmax
s <- u - x1
chol.o <- chol(e, tmpmax=tmpmax, nsubmax=nsubmax, nnzlmax=nnzlmax)
b <- backsolve(chol.o, ao1 %*% y)
r1 <- y - x %*% b
z1 <- ifelse(abs(r1) < ctrl$small, (r1*(r1>0)+ctrl$small), r1*(r1>0))
w <- z1 - r1
z2 <- rep(1,n2)
wwn1 <- matrix(0,n1,10)
wwn1[,1] <- z1
wwn1[,2] <- w
wwn2 <- matrix(0,n2,7)
wwn2[,2] <- z2
fit <- .Fortran("srqfnc",
n1 = as.integer(n1),
m = as.integer(m),
nnzx = as.integer(nnzx),
x = as.double(x@ra),
jx = as.integer(x@ja),
ix = as.integer(x@ia),
ao1 = as.double(ao1@ra),
jao1 = as.integer(ao1@ja),
iao1 = as.integer(ao1@ia),
n2 = as.integer(n2),
nnzR = as.integer(nnzR),
R = as.double(R@ra),
jR = as.integer(R@ja),
iR = as.integer(R@ia),
ao2 = as.double(ao2@ra),
jao2 = as.integer(ao2@ja),
iao2 = as.integer(ao2@ia),
nnzdmax = as.integer(nnzdmax),
d = double(nnzdmax),
jd = integer(nnzdmax),
id = integer(m+1),
dsub = double(nnzhmax+1),
jdsub = integer(nnzhmax+1),
nnzemax = as.integer(nnzemax),
e = as.double(e@ra),
je = as.integer(e@ja),
ie = as.integer(e@ia),
nnzgmax = as.integer(nnzgmax),
g = double(nnzgmax),
jg = integer(nnzgmax),
ig = integer(m+1),
nnzhmax = as.integer(nnzhmax),
h = double(nnzhmax),
jh = integer(nnzhmax),
ih = integer(m+1),
nsubmax = as.integer(nsubmax),
lindx = integer(nsubmax),
xlindx = integer(m+1),
nnzlmax = as.integer(nnzlmax),
lnz = double(nnzlmax),
xlnz = integer(m+1),
iw = integer(m*5),
iwmax = as.integer(iwmax),
iwork = integer(iwmax),
xsuper = integer(m+1),
tmpmax = as.integer(tmpmax),
tmpvec = double(tmpmax),
maxn1n2 = as.integer(maxn1n2),
ww1 = double(maxn1n2),
wwm = as.double(wwm),
wwn1 = as.double(wwn1),
wwn2 = as.double(wwn2),
cachsz = as.integer(ctrl$cachsz),
level = as.integer( 8 ),
x1 = as.double(x1),
x2 = as.double(x2),
s = as.double(s),
u = as.double(u),
c1 = as.double(y),
c2 = as.double(r),
sol = as.double(b),
small = as.double(ctrl$small),
ierr = integer(1),
maxiter = as.integer(ctrl$maxiter),
time = double(7))[c("sol","ierr","maxiter","time")]
ierr <- fit$ierr
if(ierr == 13)# stop()
stop("Increase nnzh.factor")
if(!(ierr==0) && ctrl$warn.mesg)
warning(sfnMessage(ierr))
coefficients <- -fit$sol
residuals <- -y - x %*% coefficients
list(coefficients = coefficients,
residuals = residuals,
control = ctrl,
ierr = ierr,
it = fit$maxiter)
}
"sfnMessage" <- function(ierr){
switch(ierr,
"insufficient storage (work space) when calling extract\n",
"nnzd > nnzdmax\n",
"insufficient storage in iwork when calling ordmmd\n",
"insufficient storage in iwork when calling sfinit\n",
"nnzl > nnzlmax when calling sfinit\n",
"nsub > nsubmax when calling sfinit\n",
"insufficient work space in iwork when calling symfct\n",
"inconsistancy in input when calling symfct\n",
"tmpsiz > tmpmax when calling bfinit; increase tmpmax\n",
"nonpositive diagonal encountered blkfct() matrix is not positive definite\n",
"insufficient work storage in tmpvec when calling blkfct\n",
"insufficient work storage in iwork when calling blkfct\n",
"impossible error condition",
"impossible error condition",
"impossible error condition",
"impossible error condition",
"tiny diagonals replaced with Inf when calling blkfct\n")
}