.packageName <- "rqmcmb2"
scoref<-function(r,tau){
  tt<-signr<-sign(r)
  tt[signr>0]<-tau
  tt[signr<0]<-tau-1
  return(tt)
}

plotDim<-function(p){
  temp<-round(sqrt(p))
  temp1<-(temp+1)
  ifelse(temp*temp1<p, dims<-c(temp1,temp1), dims<-c(temp,temp1))
  if(temp^2==p || p<temp*temp){dims<-c(temp,temp)}
  return(list(dims=dims))
}


#meansd<-function(obj){
#  means <- apply(obj$theta, 2, mean)
#  sds <- sqrt(apply(obj$theta, 2, var))
#  res<-matrix(cbind(means, sds), ncol=2, byrow=FALSE)
#  dimnames(res)<-list(NULL, c("Mean","SD"))
#  return(res)
#}

rqmcmb.ci<-function(obj, alpha=0.1){
  sds <- sqrt(apply(obj$theta, 2, var))
  ciL<-obj$coef+qnorm(alpha/2)* sds
  ciU<-obj$coef-qnorm(alpha/2)* sds
  res<-matrix(cbind(obj$coef, sds, ciL,ciU ), ncol=4, byrow=FALSE)
  dimnames(res)<-list(NULL, c("Coef","SD", "L", "U"))
  return(res)
}


rqmcmb.plot<-function(obj, alpha=.10){
  p <- dim(obj$theta)[2]
  K <- length(obj$theta[,1])
  ci <- rqmcmb.ci(obj)
  lci <- ci[,3]
  uci <- ci[,4]
  #theta <- obj$theta
  thetaOrig <- obj$theta[1,]
  z<-qnorm(alpha/2)
  dev<-dev.cur()
  if(dev==1) {motif()}
  par(mfrow=plotDim(p)$dims)
  plotr <- apply(obj$theta,2,range)
  plotlim <- rbind(pmin(plotr[1,],lci), pmax(plotr[2,],uci))
  for (i in 1:p){
    plot(obj$theta[,i], ylab=paste("theta[",i,"]"),
      ylim=plotlim[,i])
    lines(seq(0, K+1,,10), rep(thetaOrig[i],10),type="b",pch=1,lty=2)
    lines(c(0, K+1), rep(lci[i],2))
    lines(c(0, K+1), rep(uci[i],2))
    title(paste("theta[",i-1,"]"))
  }
  par(mfrow=c(1,1))

}

.First.lib <-
function(lib, pkg) {
   library.dynam("rqmcmb2", pkg, lib)
   print("rqmcmb2 library loaded")}
rqmcmb <- function(x=x, y=y, tau=0.5, K=100, int=TRUE, plotTheta=FALSE) {

if(exists(".Random.seed") == FALSE){ 
  rnorm(1)
  seed <- .Random.seed[2]  
}
else{
  seed <- .Random.seed[2]
}

n<-length(y)
x<-as.matrix(x)
if (nrow(x)!=n) stop("Warning: Dimensions of x and y do not match") 

if(int==TRUE) x <- cbind(rep(1,n), x)

if (n>=200000){
  print("Sample size is too large. Maximum allowed is 200000")}
else {
  if (ncol(x)>=100){
    print("Number of parameters is too large. Maximum allowed is 100")}


#main code


else { 

  if ( n <= 5000) {fit<-rq.fit(x,y, tau=tau, ci=FALSE)} else {
    if ( n <=50000) {fit<-rq.fit.fn(x,y, tau=tau) } else {
      fit<-rq.fit.pfn(x,y, tau=tau) } }


  thetaOrig<-fit$coef
  p<-length(thetaOrig)
  thetaTilda<-thetaOrig


  if ((n*tau < 5*p+1) || (n*(1-tau) < 5*p+1)) {
    print("Warning: May not have enough data to estimate the")
    print("requested quantile reliably")}

  Z <- t(x)%*%x
  cov.svd<-svd(Z)
  A<-cov.svd$u%*%(diag(sqrt(cov.svd$d)^-1))%*%t(cov.svd$v)
  Ainv<-cov.svd$u%*%(diag(sqrt(cov.svd$d)))%*%t(cov.svd$v)
  thetaTilda<-matrix(Ainv%*%thetaOrig, nrow=length(thetaOrig))

  cn<-sqrt(max(cov.svd$d)/min(cov.svd$d))
  if (cn>100){
    print("Warning: Nearly singular design detected;")
    print("the results from rqmcmb may be unreliable")
  }

  psi<-scoref(fit$resid,tau)
  psimat<-matrix(psi,nrow=n,ncol=p,byrow=FALSE)
  ZTilda<-(x%*%A)*psimat

  # re-define x on the theta_tilda scale
  x<-x%*%A
  sumxij<-apply(x,2,sum)
  sumabsxij<-apply(abs(x),2,sum)
  zstar<-.C("rqmcmb",
        array(as.double(t(x))),
        array(as.double(y)),
        array(as.double(tau)),
        array(as.double(thetaTilda)),
        array(as.double(t(A))),
	array(as.double(ZTilda)),
        array(as.double(sumxij)),
        array(as.double(sumabsxij)),
        as.integer(n),
        as.integer(p),
	success=as.integer(1),
        theta=array(as.double(rep(0,K*p+p)),c(p,K+1)),
        as.integer(K),
        as.integer(seed),
        PACKAGE="rqmcmb2"
        )


  if(zstar$success==0){return(list(success=0))}

  thetaTildaMCMB<-t(zstar$theta)
  thetaOrigMCMB<-t(A%*%zstar$theta)

  if(plotTheta==TRUE){
    dev<-dev.cur()
    if(dev==1) {x11()}
    par(mfrow=plotDim(p)$dims)

  for (i in 1:p){
    plot(thetaOrigMCMB[,i], ylab=paste("theta[",i-1,"]"))
    lines(c(0, K+1), rep(thetaOrig[i],2))
    title(paste("theta[",i-1,"]"))
  }

  par(mfrow=c(1,1))
}

return(list(coef=thetaOrig, cov=var(thetaOrigMCMB), theta=thetaOrigMCMB, 
  success=zstar$success, cn=cn))
}

}
}

