Nothing
R/qDEA.R
In qDEA: Quantile Data Envelopment Analysis
Defines functions
LP_highs
LPSOLVER
sort_list
merge_lists
rbindSM
cbindSM
SM2SM
SM2A
A2SM
iter_delete
Write.LP
my_seconds
lagMat
my_lag
qDEAbuild
DEAbuild
nCm_mpick
qDEA_mlist
qDEA_solve
qDEA
Documented in
A2SM
cbindSM
DEAbuild
iter_delete
lagMat
LP_highs
LPSOLVER
merge_lists
my_lag
my_seconds
nCm_mpick
qDEA
qDEAbuild
qDEA_mlist
qDEA_solve
rbindSM
SM2A
SM2SM
sort_list
Write.LP
###############################################################################
#' qDEA: Calling function for set of DEA and qDEA processes
#'
#' Note: Optional arguments in function call have default values of 'NULL'
#'
#' @param X Reference dmu's = ndmu x number of inputs input matrix.
#' @param Y Reference dmu's = ndmu x number of outputs output matrix.
#' @param qout Maximal proportion of dmu's allowed external to DEA hull.
#' Default = 1/ndmu
#' @param qoutS Proportion of external points to identify using qDEA slicing
#' (Atwood and Shaik 2020 EJOR)with qoutS <= qout. Default = qout
#' @param X0 Inputs for set of ndmu0 dmu's to be processed. Default = X
#' @param Y0 Outputs for set of ndmu0 dmu's to be processed. Default = Y
#' @param DX0 Input directions for ndmu0 dmu's in X0 and Y0.
#' (Must be provided if orient = 'ddea'). Default = NULL
#' @param DY0 Output directions for ndmu0 dmu's in X0 and Y0.
#' (Must be provided if orient = 'ddea') Default = NULL
#' @param orient Model orientation ('in','out','inout','oneone','ddea').
#' (!!If orient='ddea',DX0,and DY0 must be provided)
#' Default='out'
#' @param RTS Returns to scale.('CRS,'VRS,'DRS','IRS) Default 'CRS'.
#' @param dmulist Index vector of dmus in (X0,Y0) to process. NULL = process all.
#' @param nqiter Maximal number of qDEA iterations. Default = 1.
#' @param nboot Number of bootstrap replications. Default = 0.
#' @param transform Transform DDEA distance to traditional efficiency metrics
#' (input or output orientations only) Default=TRUE
#' @param mcells Number of subsample sizes for bootstapping. Default = 5
#' @param mlist Optional list of user chosen subsample sizes. Default = NULL
#' @param seedval Seed value for random number generator - used in bootstapping.
#' Default = 1001.
#' @param qtol q search tolerance with iterative qDEA. Default = 1E-6
#' @param BIGM Default Big M in RHS of qDEA stage 2 process. Default = 1E9
#' @param eps Search tolerance in qDEA improvement tests. Default = 1E-6
#' @param skipzprob Skip qDEA if qout=0. Default=TRUE.
#' @param replaceA1 Put dmu0's data in first row of reference sets. Default=FALSE
#' @param baseqDEA Use basic qDEA model from EJOR article. Default=FALSE
#' @param unbounded qDEA reported as unbounded if obj<=unbounded. Default = (-1E3)
#' @param obj2test Convergance tol for objective in iterative qDEA.Default = 1E-4
#' @param replaceM Subsample with replacement in subsample bootstrap. Default=FALSE
#' @param alpha Alpha level for Simar-Wilson(2011)subsample size selection
#' procedure. Default = 0.05
#' @param betaq qDEA convergence rate. 0.5 => "root n" convergence.
#' @param siglist Vector of user's desired confidence interval widths
#' Default = c(0.10, 0.05, 0.01)
#' @param CILag CILag level for Simar-Wilson (2011) subsample size
#' selection procedure. Default = 1
#' @param printlog Progress of DDEA dmu's solved when (X0,Y0) >1 dmus.
#' @param prntmod Print progress every prntmod dmus. Default=100.
#' @param printtxt Additional text to print with progress printlog.
#' @param getproject Compute projected values for dmu's in dmulist. Default=FALSE
#' @param getbootpeers Pull and store peers for each bootstrapped solution.
#' Default=FALSE
#' @param solver LP solver Default='highs'
#'
#' @return '###################################################################'
#' @return A list of individual items and lists of addiitional/optional output
#'
#' @return effvals = vector of DDEA distances
#' (efficiencies if transform=TRUE and input or output orientation).
#' @return effvalsq = vector of qDDEA distances
#' (efficiencies if transform=TRUE and input or output orientation).
#' @return distvals = vector of DDEA distances
#' @return distvalsq = vector of qDEA distances
#' @return distMAT = ndmu0 by 3 matrix with DDEA,qDDEA-S1,qDDEA-S2 distances.
#' '###########################################################################'
#' @return INPUT_DATA = A list containing:
#'
#' @return X = reference dmu's = ndmu x number of inputs input matrix.
#' @return Y = reference dmu's = ndmu x number of outputs output matrix.
#' @return X0 = inputs for set of ndmu0 dmu's processed.
#' @return Y0 = outputs for set of ndmu0 dmu's processed.
#' @return DX0 = input directions for ndmu0 dmu's processed.
#' @return DY0 = output directions for ndmu0 dmu's processed.
#' @return qout = maximal proportion of dmu's allowed external to DEA hull.
#' @return qoutS = proportion of external points to identify using qDEA slicing.
#' @return RTS = returns to scale.
#' @return orient = model orientation. ('ddea','in','out','inout','oneone')
#' @return baseqDEA = use basic qDEA model from EJOR article
#' @return dmulist0 = DMU0 index in originally inputs X0,Y0,DX0,and DY0
#' @return '###################################################################'
#' @return BOOT_DATA = A list containing:
#'
#' @return effvals.bc = vector of DDEA bias corrected distance or efficiency metrics
#' (efficiencies if transform=TRUE and input or output orientation)
#' @return effvalsq.bc = vector of qDEA bias corrected distance or efficiency metrics
#' (efficiencies if transform=TRUE and input or output orientation)
#' @return distvals.bc = vector of bias corrected DDEA distances
#' @return distvalsq.bc = vector of bias corrected qDEA distances
#' @return mcells = number of subsample sizes for bootstapping
#' @return mlist = list of subsample sizes used in bootstraps
#' @return BOOTdmus = matrix containing indexes of reference dmus chosen for each bootstrap
#' @return BOOT = nboot by mcells by ndmu0 array of bootstrapped DDEA distances
#' (efficiencies if transform=TRUE and input or output orientation)
#' @return BOOTS = nboot by mcells by ndmu0 array of DDEA bootstrapped s-statistics
#' (statn = (m(n)/n)^beta (statm - statn)
#' @return BOOTs = nboot by ndmu0 matrix of DDEA bootstrapped s-statistics
#' (at sample size chosen by Simar and Wilson (2011) suggested process)
#' @return BOOTq = nboot by mcells by ndmu0 array of bootstrapped qDEA distances
#' (efficiencies if transform=TRUE and input or output orientation)
#' @return BOOTSq = nboot by mcells by ndmu0 array of qDEA bootstrapped s-statistics
#' (statn = (m(n)/n)^beta (statm - statn)
#' @return BOOTsq = nboot by ndmu0 matrix of qDEA bootstrapped s-statistics
#' (at sample size chosen by Simar and Wilson (2011) suggested process
#' @return mpick = length(ndmu0) index vector of DDEA bootstrap sample size chosen from mlist
#' using Simar and Wilson (2011) suggested sample size selection process
#' @return mpickq = length(ndmu0) index vector of qDDEA bootstrap sample size chosen from mlist
#' using Simar and Wilson (2011) suggested sample size selection process
#' @return beta = DDEA convergence rate indicated Simar-Wilson (2011)
#' @return betaq = qDEA convergence rate indicated Atwood and Shaik(2020) = 0.5
#' @return siglist = vector of user's desired confidence interval widths
#' @return CI = ndmu0 by 2 by length(siglist) array of DDEA confidence intervals
#' (estimated using quantiles on bias corrcted "s" statistics)
#' @return CIq = ndmu0 by 2 by length(siglist) array of qDEA confidence intervals
#' (estimated using quantiles on bias corrcted "s" statistics)
#' @return CIq_norm = ndmu0 by 2 by length(siglist) array of qDEA confidence intervals
#' (estimated using sample mean and standard error of bootstrapped
#' "s" statistics and assuming normality) (Atwood and Shaik 2020)
#' @return '###################################################################'
#' @return PEER_DATA = A list containing:
#'
#' @return PEERS = dataframe of DDEA peers and projection weights for each dmu0
#' @return PEERSq = dataframe of qDDEA-S2 peers and projection weights each dmu0
#' @return DOUTq = data frame indicating external dmus for each dmu's qDEA solution
#' @return BPEERS = ndmu0 by mcells by (nIN+nOUT) by 2(dmu-z,zweight) array
#' of DDEA subsampled peers (dmu-z) and projection weights (z)
#' if(nboot>0&getbootpeers==TRUE)
#' @return BPEERSq = ndmu0 by mcells by (nIN+nOUT) by 2(dmu-z,zweight) array
#' of qDEA subsampled peers (dmu-z) and projection weights (z)
#' if(nboot>0&getbootpeers==TRUE)
#' @return BPEERS1 = data.frame of DDEA subsampled peers and projection weights (z)
#' if(nboot>0&getbootpeers==TRUE) for chosen subsample size
#' data.frame: nb = bootstrap rep, dmu0 = given dmu,
#' dmuz = reference dmu, z = reference dmu projection weight
#' @return BPEERS1q = data.frame of qDEA subsampled peers and projection weights (z)
#' if(nboot>0&getbootpeers==TRUE) for chosen subsample size
#' data.frame: nb = bootstrap rep, dmu0 = given dmu,
#' dmuz = reference dmu, z = reference dmu projection weight
#' @return '###################################################################'
#' @return PROJ_DATA = Projected Values. A list containing:
#'
#' @return X0HAT = DDEA projected input levels (if getproject=TRUE)
#' @return Y0HAT = DDEA projected output levels (if getproject=TRUE)
#' @return X0HAT.bc = bias corrected DDEA projected input levels
#' (if nboot>0 and getproject=TRUE)
# @return Y0HAT.bc = bias corrected DDEA projected output levels
#' (if nboot>0 and getproject=TRUE) )
#' @return X0HATq = qDEA projected input levels (if getproject=TRUE)
#' @return Y0HATq = qDEA projected output levels (if getproject=TRUE)
#' @return X0HATq.bc = bias corrected qDEA projected input levels
#' (if nboot>0 and getproject=TRUE)
#' @return Y0HATq.bc = bias corrected qDEA projected output levels
#' (if nboot>0 and getproject=TRUE) )
#' @return '###################################################################'
#' @return LP_DATA = LP Models, data, and results. A list containing:
#'
#' @return status = ndmu0 by 3 matrix with LP status of DDEA,qDDEA-S1,qDDEA-S2
#' @return qhat = proportion of dmu's external to hull in qDEA solution
#' (NA indicates qDEA for given DMU was unbounded)
#' @return qiter = number of qDEA iterations completed.
#' @return PSOL = ndmu0 by ? matrix with DDEA LP solutions for each dmu0
#' @return PSOLq1 = ndmu0 by ? matrix with qDDEA-S1 LP solutions for each dmu0
#' @return PSOLq2 = ndmu0 by ? matrix with qDDEA-S2 Lp solutions for each dmu0
#' @return RCOST = ndmu0 by ? matrix of LP reduced costs for DDEA solutions
#' @return RCOSTq1 = ndmu0 by ? matrix of LP reduced costs for qDDEA-S1 solutions
#' @return RCOSTq2 = ndmu0 by ? matrix of LP reduced costs for qDDEA-S2 solutions
#' @return LPModels = list of LP0, LP1, and LP2 LP objects from qDEA_solve function
#' @examples
#' # Examples from CST(2006): Cooper, W., Seiford, L., and Tone, K., 2006.
#' # Introduction to Data Envelopment Analysis and Its Uses. Springer, New York.
#'
#' # CST One Input - One Output Example - Table 1.1
#' data(CST11)
#' X = as.matrix(CST11$EMPLOYEES)
#' Y = as.matrix(CST11$SALES_EJOR)
#' qout = 1/nrow(X)
#' sol = qDEA(X, Y, RTS = 'crs', orient = 'in', qout = qout)
#' # DEA efficiency scores
#' sol$effvals
#' # qDEA input efficiency scores allowing one external point
#' sol$effvalsq
#'
#' # CST Two Input - One Output Example - Table 1.3
#' data(CST21)
#' X = as.matrix(cbind(CST21$EMPLOYEES, CST21$FLOOR_AREA))
#' Y = as.matrix(CST21$SALES)
#' qout = 1/nrow(X)
#' sol = qDEA(X, Y, RTS = 'crs', orient = 'in', qout = qout)
#' sol$effvals
#' sol$effvalsq
#'
#' # CST One Input - Two Output Example - Table 1.4
#' data(CST12)
#' X = as.matrix(CST12$EMPLOYEES)
#' Y = as.matrix(cbind(CST12$CUSTOMERS, CST12$SALES))
#' qout = 1/nrow(X)
#' sol = qDEA(X, Y, RTS = 'crs', orient = 'out', qout = qout)
#' sol$effvals
#' sol$effvalsq
#'
#' # CST Two Input - Two Output Example - Table 1.5
#' data(CST22)
#' X = as.matrix(cbind(CST22$DOCTORS, CST22$NURSES))
#' Y = as.matrix(cbind(CST22$OUT_PATIENTS, CST22$IN_PATIENTS))
#' qout = 1/nrow(X)
#' sol = qDEA(X, Y, RTS = 'crs', orient = 'in', qout = qout)
#' # DEA efficiency scores - see table 1.6 CCR estimates
#' round(sol$effvals, 2)
#' # qDEA efficiency scores allowing one external point
#' round(sol$effvalsq, 2)
#'
#' \donttest{
#' # Atwood-Shaik EJOR qDEA Examples (longer running examples)
#' data(CST11)
#' X = as.matrix(CST11$EMPLOYEES)
#' Y = as.matrix(CST11$SALES_EJOR)
#' ###############################################################################
#' # EJOR Efficiency Results Table 1 Input Orientation
#' tmpC1=qDEA(X,Y,RTS='crs',orient='in',qout=1/8,getproject=TRUE)
#' tmpC2=qDEA(X,Y,RTS='crs',orient='in',qout=2/8,getproject=TRUE)
#'
#' # Table 1 Input Orientation
#' round(cbind(tmpC1$distvals,tmpC1$PROJ_DATA$X0HAT,tmpC1$PROJ_DATA$Y0HAT,
#' tmpC1$distvalsq,tmpC1$PROJ_DATA$X0HATq,tmpC1$PROJ_DATA$Y0HATq,
#' tmpC2$distvalsq,tmpC2$PROJ_DATA$X0HATq,tmpC2$PROJ_DATA$Y0HATq),3)
#' ###############################################################################
#' # EJOR Efficiency Results Table 1 Output Orientation
#' tmpC3=qDEA(X,Y,RTS='crs',orient='out',qout=1/8,getproject=TRUE)
#' tmpC4=qDEA(X,Y,RTS='crs',orient='out',qout=2/8,getproject=TRUE)
#'
#' # Table 1 Output Orientation
#' round(cbind(tmpC3$distvals,tmpC3$PROJ_DATA$X0HAT,tmpC3$PROJ_DATA$Y0HAT,
#' tmpC3$distvalsq,tmpC3$PROJ_DATA$X0HATq,tmpC3$PROJ_DATA$Y0HATq,
#' tmpC4$distvalsq,tmpC4$PROJ_DATA$X0HATq,tmpC4$PROJ_DATA$Y0HATq),3)
#' ################################################################################
#' # EJOR Efficiency Results Table 1 one-one Orientation
#' tmpC5=qDEA(X,Y,RTS='crs',orient='oneone',qout=1/8,getproject=TRUE)
#' tmpC6=qDEA(X,Y,RTS='crs',orient='oneone',qout=2/8,getproject=TRUE)
#'
#' # Table 1 one-one Orientation
#' round(cbind(tmpC5$distvals,tmpC5$PROJ_DATA$X0HAT,tmpC5$PROJ_DATA$Y0HAT,
#' tmpC5$distvalsq,tmpC5$PROJ_DATA$X0HATq,tmpC5$PROJ_DATA$Y0HATq,
#' tmpC6$distvalsq,tmpC6$PROJ_DATA$X0HATq,tmpC6$PROJ_DATA$Y0HATq),3)
#' ############################################################################
#' }
###############################################################################
#
###############################################################################
# debugging function qDEA
###############################################################################
###############################################################################
## default settings
# X=X;Y=Y;X0=NULL;Y0=NULL;DX0=NULL;DY0=NULL
###############################################################################
# qout=1/nrow(X);qoutS=qout
###############################################################################
# orient='out';RTS='CRS';dmulist=NULL;nqiter=1;nboot=0;transform=TRUE;mcells=5
# mlist=NULL;seedval=1001;qtol=1E-6;BIGM=1E9;eps=1E-6;skipzprob=TRUE;replaceA1=FALSE
# baseqDEA=FALSE;unbounded=(-1E3);obj2test=1E-4;replaceM=FALSE;alpha=0.05;betaq=0.5
# siglist=c(0.10,0.05,0.01);CILag=1;printlog=TRUE;prntmod=100;printtxt=''
# getproject=FALSE;getbootpeers=FALSE
###############################################################################
# MRE plots
# tmpC=qDEA(X,Y,qout=qout,RTS=RTS,nqiter=nqiter,orient='ddea',
# X0=X0,Y0=Y0,DX0=DX0,DY0=DY0,baseqDEA=TRUE,getproject=TRUE)
#
# RTS='VRS';nboot=0; mcells=1
###############################################################################
# John-Illex-Data
#tmp2=qDEA(X,Y,X0=X,Y0=Y,orient='ddea',DX0=DX0,DY0=DY0,
# RTS='VRS',qout=qout, nboot=nboot,mcells=1,
# getbootpeers=TRUE, replaceA1=TRUE, getproject=TRUE)
###############################################################################
# nqiter=4; RTS='VRS'; nboot=50; mcells=1
# #############################################################################
# qoutS=qout;orient='ddea';dmulist=NULL;transform=TRUE
# mlist=NULL;seedval=1001;qtol=1E-6;BIGM=1E9;eps=1E-6;skipzprob=TRUE;replaceA1=TRUE
# baseqDEA=TRUE;unbounded=(-1E3);obj2test=1E-4;replaceM=FALSE;alpha=0.05;betaq=0.5
# siglist=c(0.10,0.05,0.01);CILag=1;printlog=TRUE;prntmod=100;printtxt=''
# getproject=TRUE;getbootpeers=TRUE;solver='highs'
###############################################################################
#' @export
###############################################################################
qDEA=function(X,Y,qout=1/nrow(X),qoutS=qout,X0=NULL,Y0=NULL,DX0=NULL,DY0=NULL,
orient='out',RTS='CRS',dmulist=NULL,nqiter=1,nboot=0,transform=TRUE,mcells=5,
mlist=NULL,seedval=1001,qtol=1E-6,BIGM=1E9,eps=1E-6,skipzprob=TRUE,replaceA1=FALSE,
baseqDEA=FALSE,unbounded=(-1E3),obj2test=1E-4,replaceM=FALSE,alpha=0.05,betaq=0.5,
siglist=c(0.10,0.05,0.01),CILag=1,printlog=TRUE,prntmod=100,printtxt='',
getproject=FALSE,getbootpeers=FALSE,solver='highs'){
###############################################################################
# print(paste('printtxt = ',printtxt))
# print(paste('solver = ',solver))
###############################################################################
###############################################################################
#setDTthreads(1)
if(baseqDEA==TRUE) qoutS=qout
# truncate qout
if(qout<1E-9) qout = 1E-9
if(qoutS<1E-9) qoutS = 1E-9
if(qoutS>qout) qoutS = qout
###############################################################################
# data check and construction
#######################################
if(!is.matrix(X)) X=as.matrix(X)
if(!is.matrix(Y)) Y=as.matrix(Y)
if(nrow(X)!=nrow(Y)) stop('STOP: nrow(X)!=nrow(Y)')
(nrX=nrow(X))
(ncX=ncol(X));(ncY=ncol(Y))
(ncYX=ncY+ncX)
#######################################
if(length(which(is.na(as.vector(X)))) >0) {
stop('STOP: The X matrix cannot contain NA elements')
}
if(length(which(is.na(as.vector(Y)))) >0) {
stop('STOP: The Y matrix cannot contain NA elements')
}
###############################################################################
###############################################################################
dmulist0=dmulist
orient
if(orient=='ddea'|orient=='DDEA'){
if(is.null(DX0)|is.null(DY0)|is.null(DX0)|is.null(DY0)) {
stop('STOP: with orient = ddea, X0, Y0, DX0 and DY0 must be non-null')
}
if(length(which(is.na(c(as.vector(X0),as.vector(Y0),
as.vector(DX0),as.vector(DY0)) ))) >0) {
stop('STOP: X0,Y0,DX0, and DY0 cannot contain NA elelemts')
}
if(is.vector(X0) & is.vector(Y0)&is.vector(DX0) & is.vector(DY0)){
ndmu0=1
if(length(X0)!=ncol(X) | length(Y0)!=ncol(Y)|length(DX0)!=ncol(X)
| length(DY0)!=ncol(Y)) {
stop('STOP: X0, Y0, DX0, DY0 vectors must of appropriate length')
}
X0=matrix(X0,nrow=ndmu0)
Y0=matrix(Y0,nrow=ndmu0)
DX0=matrix(DX0,nrow=ndmu0)
DY0=matrix(DY0,nrow=ndmu0)
} # end if(is.vector(X0,Y0,DX0 & DY0))
(ndmu0=nrow(X0))
if(ncol(DX0)!=ncol(X)) stop('STOP:ncol(DX0)!=ncol(X)')
if(ncol(DY0)!=ncol(Y)) stop('STOP:ncol(DY0)!=ncol(Y)')
if(!is.null(dmulist)){
dmulist0=dmulist
if(max(dmulist)>ndmu0) stop('STOP: max dmulist>nrow(x0)')
ndmu0=length(dmulist)
X0=matrix(X0[dmulist,],nrow=ndmu0)
Y0=matrix(Y0[dmulist,],nrow=ndmu0)
DX0=matrix(DX0[dmulist,],nrow=ndmu0)
DY0=matrix(DY0[dmulist,],nrow=ndmu0)
dmulist=1:nrow(X0)
}
if(is.null(dmulist)) dmulist = 1:nrow(X0)
} # end if(orient=='ddea')
################################################################################
if(orient!='ddea'&orient!='DDEA'){
DX0=NULL; DY0=NULL
}
################################################################################
################################################################################
if(!is.null(X0)|!is.null(Y0)){
if(is.vector(X0)&is.vector(Y0)){
X0=matrix(X0,nrow=1)
Y0=matrix(Y0,nrow=1)
} # end if(is.vector(X0)&is.vector(Y0)){
if(nrow(X0)!=nrow(Y0)) stop('STOP: nrow(X0)!=nrow(Y0)')
if(ncol(X0)!=ncol(X)) stop('STOP: ncol(X0)!=ncol(X)')
if(ncol(Y0)!=ncol(Y)) stop('STOP: ncol(Y0)!=ncol(Y)')
if(length(which(is.na(as.vector(X0)))) >0) {
stop('STOP: The X0 matrix cannot contain NA elements')
}
if(length(which(is.na(as.vector(Y0)))) >0) {
stop('STOP: The Y0 matrix cannot contain NA elements')
}
} # end if(!is.null(X0)|!is.null(Y0)){
###############################################################################
if(is.null(dmulist)&(is.null(X0))){
dmulist = 1:nrow(X)
X0=X; Y0=Y
}
###############################################################################
if(!is.null(dmulist)&(!is.null(X0))){
################################
if(length(dmulist)>nrow(X0)){
stop('STOP: length(dmulist)>nrow(X0)')
}
if(max(dmulist)>nrow(X0)){
stop('STOP: max(dmulist)>nrow(X0)')
}
if(length(dmulist)<nrow(X0)){
dmulist0=dmulist
ndmu0=length(dmulist)
X0=matrix(X0[dmulist,],nrow=ndmu0)
Y0=matrix(Y0[dmulist,],nrow=ndmu0)
dmulist=1:nrow(X0)
}
###################################
} # end if(!is.null(dmulist)&(!is.null(X0)|!is.null(Y0)))
#######################################
if(!is.null(dmulist)&is.null(X0)){
ndmu0=length(dmulist)
dmulist0=dmulist
X0=matrix(X[dmulist,],nrow=ndmu0)
Y0=matrix(Y[dmulist,],nrow=ndmu0)
} else{
if(is.null(dmulist)&!is.null(X0)){
dmulist=1:nrow(X0)
dmulist0=dmulist
}
}
###############################################################################
###############################################################################
(ndmu0=nrow(X0))
#######################################
if(is.null(DX0)|is.null(DY0)) {
(DX0=matrix(0,ndmu0,ncX))
(DY0=matrix(0,ndmu0,ncY))
if(orient=='in'|orient=='IN'){
(DX0=matrix(X0,nrow=ndmu0))
}
if(orient=='out'|orient=='OUT'){
(DY0=matrix(Y0,nrow=ndmu0))
}
if(orient=='inout'|orient=='INOUT'){
(DX0=matrix(X0,nrow=ndmu0))
(DY0=matrix(Y0,nrow=ndmu0))
}
if(orient=='oneone'|orient=='ONEONE'){
DX0=X0^0
DY0=Y0^0
}
} # end if is.null(DX0 or DY0)
(nDMU0=nrow(X0))
###############################################################################
###############################################################################
## highs blows up if abs of any sparse matrix ra values are <1e-9
#X[abs(X)< 2e-9] = 2e-9; Y[abs(Y)< 2e-9] = 2e-9
#X0[abs(X0)< 2e-9] = 2e-9; Y0[abs(Y0)< 2e-9] = 2e-9
#DX0[abs(DX0)< 2e-9] = 2e-9; DY0[abs(DY0)< 2e-9] = 2e-9
###############################################################################
time1=NULL;time2=NULL;time3=NULL
###############################################################################
time10=my_seconds()
tmp0=qDEA_solve(X,Y,qout=qout,qoutS=qoutS,X0=X0,Y0=Y0,DX0=DX0,DY0=DY0,RTS=RTS,
nqiter=nqiter,qtol=qtol,BIGM=BIGM,eps=eps,skipzprob=skipzprob,
unbounded=unbounded,obj2test=obj2test,baseqDEA=baseqDEA,replaceA1=replaceA1,
printlog=printlog,prntmod=prntmod,printtxt=printtxt,dmulist0=dmulist0,
solver=solver)
time11=my_seconds()
(time1=time11-time10)
###############################################################################
###############################################################################
BOOT='NOT RUN'; BOOTq='NOT RUN'; BOOT_M='NOT RUN'; BOOT_dmus='NOT RUN'
BPEERS=list('NOT RUN'); BPEERSq=list('NOT RUN')
###############################################################################
#debugoff
if(nboot>0){
###############################################################################
##
(Bdrop = which(is.na(tmp0$effvals[,1])))
if(length(Bdrop)==0) Bdrop=NULL
(Bdropq = which(is.na(tmp0$effvals[,3])))
##
(n=nrX)
if(!is.null(mlist)) mcells=length(mlist)
##
if(is.null(mlist[1])){
# compute mlist values
if(mcells==1){
if(qout<1e-6) mlist=floor(2.5*sqrt(nrX))
if(qout>=1e-6) mlist=floor(max(2.5*sqrt(nrX),1/qout))+1
} # end if(mcells==1){
if(mcells>=2){
(n1=floor(1.0*sqrt(nrX)))
(if(qout>0) n1=floor(max(1/qout,sqrt(nrX)))+1)
(n2=floor(min(5*n1,nrX/2))+1)
if(n2>0.5*nrX) (n2=floor(0.5*nrX)+1)
(mlist=round(exp(seq(log(n1),log(n2),length.out=mcells)),0))
if(max(table(mlist)>1)) (mlist=floor(seq(n1,n2,length.out=mcells)))
} # end if(mcells>=2)
mlist
} # end if(is.null(mlist[1]))
###############################################################################
mlist
(mmax=max(mlist))
if(mmax>nrX){
stop('STOP: User must provide subsample size vector with largest element < nrow(X)')
}
BOOT_dmus = matrix(NA,nboot,mmax)
BOOT=array(NA,dim=c(nboot,mcells,nDMU0))
BOOTq=BOOT
time20=my_seconds()
###############################################################################
nb=1
# debugoff
for(nb in 1:nboot){
if(printlog==TRUE&nb%%prntmod==0){
print(paste(printtxt,'bootstrap',nb,'of',nboot))
}
set.seed(seedval+nb)
mpick=sample(1:nrX,mmax,replace=replaceM)
BOOT_dmus[nb,]=mpick
YM=as.matrix(Y[mpick,]); XM=as.matrix(X[mpick,])
tmp2=qDEA_mlist(XM=XM,YM=YM,mpick=mpick,Bdrop=Bdrop,Bdropq=Bdropq,qout=qout,
qoutS=qoutS,X0=X0,Y0=Y0,DX0=DX0,DY0=DY0,mlist=mlist,RTS=RTS,nqiter=nqiter,
unbounded=unbounded,replaceA1=replaceA1,baseqDEA=baseqDEA,
getbootpeers=getbootpeers,dmulist0=dmulist0,solver=solver)
BOOT[nb,,]=t(tmp2$EFFmlist)
BOOTq[nb,,]=t(tmp2$EFFmlistq)
peers=NULL; peersq=NULL
if(getbootpeers==TRUE){
peers = tmp2$peers
peers$nb = nb
peersq = tmp2$peersq
peersq$nb = nb
BPEERS[[nb]]=peers
BPEERSq[[nb]]=peersq
} # end if(getbootpeers==TRUE)
###############################################################################
} # end loop on nb
###############################################################################
if(getbootpeers==TRUE){
BPEERS = as.data.frame(dplyr::bind_rows(BPEERS))
BPEERSq = as.data.frame(dplyr::bind_rows(BPEERSq))
BPEERS$dmuz=ifelse(BPEERS$dmuz==0,BPEERS$dmu0,BPEERS$dmuz)
BPEERSq$dmuz=ifelse(BPEERSq$dmuz==0,BPEERSq$dmu0,BPEERSq$dmuz)
} # end if (getbootpeers==TRUE)
if(getbootpeers!=TRUE) {BPEERS=list('NOT RUN'); BPEERSq=list('NOT RUN')}
time21=my_seconds()
(time2=time21-time20)
###############################################################################
###############################################################################
} # end if(nboot>0)
###############################################################################
###############################################################################
# process results
###############################################################################
# DEA and qDEA results
###############################################################################
distMAT=tmp0$effvals;qhat=tmp0$qhat;qiter=tmp0$qiter; status=tmp0$status
PSOL=tmp0$PSOL;PSOLq1=tmp0$PSOLq1;PSOLq2=tmp0$PSOLq2
PEERS=tmp0$PEERS;PEERSq=tmp0$PEERSq
RCOST=tmp0$RCOST;RCOSTq1=tmp0$RCOSTq1;RCOSTq2=tmp0$RCOSTq2
DOUTq=tmp0$D2OUT;devstart=tmp0$devstart;devend=tmp0$devend
LP0=tmp0$LP0;LP1=tmp0$LP1;LP2=tmp0$LP2
################################################################################
PEERS$dmuz=ifelse(PEERS$dmuz==0,PEERS$dmu0,PEERS$dmuz)
if(!is.character(PEERSq)) {
PEERSq$dmuz=ifelse(PEERSq$dmuz==0,PEERSq$dmu0,PEERSq$dmuz)
}
################################################################################
distvals=distMAT[,1]
distvalsq=distMAT[,3]
effvals=distvals
effvalsq=distvalsq
######################################
if(transform==TRUE){
if(orient=='in'|orient=='IN') {
effvals=1-distvals
effvalsq=1-distvalsq
}
if(orient=='out'|orient=='OUT') {
effvals=1+distvals
effvalsq=1+distvalsq
}
} # end if(transform==TRUE){
##############################################################################
##############################################################################
# bootstrap and bias correction results
###############################################################################
distvals.bc=rep(NA,length(distvals));distvalsq.bc=rep(NA,length(distvalsq))
effvals.bc=distvals.bc;effvalsq.bc=distvalsq.bc
mpick=NA; mpickq=NA; time_boot=NA
BOOTS='NOT RUN'; BOOTs='NOT RUN'
BOOTSq='NOT RUN'; BOOTsq='NOT RUN'
CI='NOT RUN'; CIq='NOT RUN'; CIq_norm='NOT RUN'
###############################################################################
if(nboot>1){
(nsigs = length(siglist))
###############################################################################
# bias correction and confidence intervals for DEA efficiency estimates
###############################################################################
BOOTS=array(NA,dim=c(nboot,mcells,nDMU0)) #set of S matrices m in mlist by DMU
BOOTs=matrix(NA,nboot,nDMU0) #set of chosen S vectors by DMU
CI=array(NA,dim=c(nDMU0,2,nsigs))
#############################
# beta values for DDEA estimates - see Simar-Wilson ?? 2009 ??
if(RTS!='CRS'&RTS!='crs') (beta = 2/(ncX + ncY + 1))
if(RTS=='CRS'|RTS=='crs') (beta = 2/(ncX + ncY))
#############################
i=1 # for debugging
for(i in 1:nDMU0){
distvals[i]
if(!is.na(distvals[i])){
tmp=nCm_mpick(stat=distvals[i],boot=BOOT[,,i],n,mlist,beta=beta,
alpha=alpha,CILag=CILag)
(mpick[i]=tmp$mpick)
S=tmp$S
s=tmp$s
BOOTS[,,i]=S
BOOTs[,i]=s
(distvals.bc[i]=mean(s,na.rm=TRUE))
j=1 # for debugging
for(j in 1:nsigs){
(sig=siglist[j])
(CI[i,1,j]=stats::quantile(s,probs=sig/2,na.rm=TRUE))
(CI[i,2,j]=stats::quantile(s,probs=1-sig/2,na.rm=TRUE))
} #end loop for(j in 1:length(siglist))
} # end if(is.na(distvals[i]))
} # end loop for(i in 1:nDMU0){
effvals.bc=distvals.bc
###############################################################################
# bias correction and confidence intervals for qDEA estimates
###############################################################################
BOOTSq=array(NA,dim=c(nboot,mcells,nDMU0)) #set of S matrices by DMU
BOOTsq=matrix(NA,nboot,nDMU0) #set of chosen S vectors by DMU
CIq=array(NA,dim=c(nDMU0,2,nsigs))
CIq_norm=array(NA,dim=c(nDMU0,2,nsigs))
i=1 # for debugging
for(i in 1:nDMU0){
distvalsq[i]
if(!is.na(distvalsq[i])){
tmpq=nCm_mpick(stat=distvalsq[i],boot=BOOTq[,,i],n,mlist,beta=betaq,
alpha=alpha,CILag=CILag)
(mpickq[i]=tmpq$mpick)
Sq=tmpq$S
sq=tmpq$s
BOOTSq[,,i]=Sq
BOOTsq[,i]=sq
(distvalsq.bc[i]=mean(sq,na.rm=TRUE))
#j=1 # for debugging
for(j in 1:nsigs){
(sig=siglist[j])
(CIq[i,1,j]=stats::quantile(sq,probs=sig/2,na.rm=TRUE))
(CIq[i,2,j]=stats::quantile(sq,probs=1-sig/2,na.rm=TRUE))
(CIq_norm[i,1,j]=stats::qnorm(sig/2, mean(sq),stats::sd(sq,na.rm=TRUE)))
(CIq_norm[i,2,j]=stats::qnorm(1-sig/2,mean(sq),stats::sd(sq,na.rm=TRUE)))
} #end loop for(j in 1:length(siglist))
} # end if(!is.na(distvalssq[i])){
} # end loop for(i in 1:nDMU0){
effvalsq.bc=distvalsq.bc
#summary(as.vector(CIq))
#(pick=which(mpickq==1))
#CIq[pick,,]
#effq
#effq.bc
###############################################################################
# transform to efficiencies
###############################################################################
if(transform==TRUE){
if(orient=='in'|orient=='IN') {
effvals.bc=1-distvals.bc
BOOT=1-BOOT
BOOTS=1-BOOTS
BOOTs=1-BOOTs
CI=1-CI
tmp=CI[,1,]
CI[,1,]=CI[,2,]
CI[,2,]=tmp
effvalsq.bc=1-distvalsq.bc
BOOTq=1-BOOTq
BOOTSq=1-BOOTSq
BOOTsq=1-BOOTsq
CIq=1-CIq
tmp=CIq[,1,]
CIq[,1,]=CIq[,2,]
CIq[,2,]=tmp
tmp=CIq_norm[,1,]
CIq_norm[,1,]=CIq_norm[,2,]
CIq_norm[,2,]=tmp
} # end if(orient=='in'|orient=='IN')
if(orient=='out'|orient=='OUT') {
effvals.bc=1+distvals.bc
BOOT = 1 + BOOT
BOOTS = 1 + BOOTS
BOOTs = 1 + BOOTs
CI = 1 + CI
effvalsq.bc=1+distvalsq.bc
BOOTq = 1 + BOOTq
BOOTSq = 1 + BOOTSq
BOOTsq = 1 + BOOTsq
CIq = 1 + CIq
CIq_norm = 1 + CIq_norm
} # end if(orient=='out'|orient=='OUT')
} # end if(transform==TRUE)
###############################################################################
} # end if(nboot>0)
###############################################################################
###############################################################################
Y0HAT='NOT RUN'; X0HAT='NOT RUN'; Y0HATq='NOT RUN'; X0HATq='NOT RUN'
Y0HAT.bc='NOT RUN'; X0HAT.bc='NOT RUN'; Y0HATq.bc='NOT RUN'; X0HATq.bc='NOT RUN'
if(getproject==TRUE) {
Y0HAT=Y0+distvals*DY0
X0HAT=X0-distvals*DX0
Y0HATq=Y0+distvalsq*DY0
X0HATq=X0-distvalsq*DX0
if(nboot>0){
Y0HAT.bc=Y0+distvals.bc*DY0
X0HAT.bc=X0-distvals.bc*DX0
Y0HATq.bc=Y0+distvalsq.bc*DY0
X0HATq.bc=X0-distvalsq.bc*DX0
} # end if(nboot>0)
} # end if(getproject==TRUE)
###############################################################################
###############################################################################
BPEERS1='NOT RUN'
BPEERSq1='NOT RUN'
if(nboot>0&getbootpeers==TRUE){
dftmp=data.frame(cbind(m=mlist[mpick],dmu0=1:nDMU0))
BPEERS1=merge(dftmp,BPEERS)
BPEERS1=doBy::orderBy(~nb+dmu0+dmuz,data=BPEERS1)
dftmp=data.frame(cbind(m=mlist[mpickq],dmu0=1:nDMU0))
BPEERSq1=merge(dftmp,BPEERSq)
BPEERSq1=doBy::orderBy(~nb+dmu0+dmuz,data=BPEERSq1)
} # end if(getbootpeers==TRUE)
##############################################################################
BOOT_DATA = list(effvals.bc=effvals.bc,effvalsq.bc=effvalsq.bc,
distvals.bc=distvals.bc,distvalsq.bc=distvalsq.bc,mcells=mcells,
mlist=mlist,BOOT_dmus=BOOT_dmus,BOOT=BOOT,BOOTS=BOOTS,BOOTs=BOOTs,
BOOTq=BOOTq,BOOTSq=BOOTSq,BOOTsq=BOOTsq,BPEERS=BPEERS,
BPEERSq=BPEERSq,BPEERS1=BPEERS1,BPEERSq1=BPEERSq1,mpick=mpick,
mpickq=mpickq,beta=beta,betaq=betaq,siglist=siglist,CI=CI,CIq=CIq,
CIq_norm=CIq_norm) # end list BOOTDATA
##
INPUT_DATA = list(X=X,Y=Y,X0=X0,Y0=Y0,DX0=DX0,DY0=DY0,qout=qout,qoutS=qoutS,
RTS=RTS,orient=orient,baseqDEA=baseqDEA,dmulist=dmulist,
dmulist0=dmulist0) # end list INPUT_DATA
##
LP_DATA = list(status=status,qhat=qhat,qiter=qiter,PSOL=PSOL,PSOLq1=PSOLq1,
PSOLq2=PSOLq2,RCOST=RCOST,RCOSTq1=RCOSTq1,RCOSTq2=RCOSTq2,
LPModels=list(LP0=tmp0$LP0,LP1=tmp0$LP1,LP2=tmp0$LP2))# endlist LPDATA
##
PEER_DATA = list(PEERS=PEERS,PEERSq=PEERSq,DOUTq=DOUTq,BPEERS=BPEERS,
BPEERSq=BPEERSq,BPEERS1=BPEERS1,BPEERSq1=BPEERSq1) # end list PEER_DATA
##
PROJ_DATA = list(X0HAT=X0HAT,Y0HAT=Y0HAT,X0HAT.bc=X0HAT.bc,Y0HAT.bc=Y0HAT.bc,
X0HATq=X0HATq,Y0HATq=Y0HATq,X0HATq.bc=X0HATq.bc,
Y0HATq.bc=Y0HATq.bc) # end list PROJ_DATA
###############################################################################
###############################################################################
list(effvals=effvals,effvalsq=effvalsq,distvals=distvals,
distvalsq=distvalsq,distMAT=distMAT,INPUT_DATA=INPUT_DATA,BOOT_DATA=BOOT_DATA,
PEER_DATA=PEER_DATA,PROJ_DATA=PROJ_DATA,LP_DATA=LP_DATA)
###############################################################################
} # end function qDEA
###############################################################################
###############################################################################
#' qDEA_solve: Sets up LP objects and solves DDEA dual DEA and qDEA using highs
#' !!!! Intended to be called from qDEA function.!!!!
#' Note: ndmu=number of dmus in reference set, ndmu0 = number dmus to process.
#' @param X Reference dmu's = ndmu x number of inputs input matrix.
#' @param Y Reference dmu's = ndmu x number of outputs output matrix.
#' @param qout Maximal proportion of dmu's allowed external to DEA hull.
#' @param qoutS Proportion of external points to identify using EJOR qDEA slicing.
#' @param X0 Inputs for set of ndmu0 dmu's to be processed.
#' @param Y0 Outputs for set of ndmu0 dmu's to be processed.
#' @param DX0 Input directions for ndmu0 dmu's in X0 and Y0.
#' @param DY0 Output directions for ndmu0 dmu's in X0 and Y0.
#' @param RTS Returns to scale default='CRS' options are 'CRS,'VRS,'DRS','IRS.
#' @param nqiter Maximal number of qDEA iterations.
#' @param qtol q search tolerance with iterative qDEA.
#' @param BIGM Default Big M in RHS of qDEA stage 2 process.
#' @param eps Zero-test criterion. abs(x) > eps implies x!=0
#' @param skipzprob Skip qDEA if qout=0.
#' @param unbounded qDEA reported as unbounded if obj<=unbounded.
#' @param obj2test Convergence tol for objective in iterative qDEA.
#' @param replaceA1 Put dmu0's data in first row of reference sets.
#' @param baseqDEA Use basic qDEA model from EJOR article
#' @param printlog Progress of dmu's solved when (X0,Y0) >1 dmus.
#' @param prntmod Print progress every prntmod dmus. default=100
#' @param printtxt Additional text to print with progress printlog.
#' @param dmulist0 DMU0 index in originally inputs X0,Y0,DX0,and DY0
#' @param solver LP solver Default='highs'
#' @return A list containing the following components:
#' @return effvals = ndmu0 by 3 matrix with DDEA,qDDEA-S1,qDDEA-S2 distances
#' (NA indicates qDEA for given DMU was unbounded)
#' @return qhat = proportion of dmu's external to hull in qDEA solution
#' (NA indicates qDEA for given DMU was unbounded)
#' @return qiter = number of qDEA iterations completed.
#' @return D2OUT = data frame indicating qDEA external reference dmus
#' @return status = ndmu0 by 3 matrix with LP status of DDEA,qDDEA-S1,qDDEA-S2
#' @return PSOL = ndmu0 by ? matrix with DDEA (dual-side) solutions for each dmu0
#' @return PSOLq1 = ndmu0 by ? matrix with qDDEA-S1 solutions for each dmu0
#' @return PSOLq2 = ndmu0 by ? matrix with qDDEA-S2 solutions for each dmu0
#' @return PEERS = dataframe of non-zero DDEA dmuz's and projection weights for each dmu0
#' @return PEERSq = dataframe of non-zero qDDEA-S2 dmuz's projection weights for each dmu0
#' @return RCOST = ndmu0 by ? matrix of LP reduced costs for DDEA solutions
#' @return RCOSTq1 = ndmu0 by ? matrix of LP reduced costs for qDDEA-S1 solutions
#' @return RCOSTq2 = ndmu0 by ? matrix of LP reduced costs for qDDEA-S2 solutions
#' @return devstart = LP index for column of first qDDEA-S1 LPM "deviation" value
#' @return devend = LP index for column of last qDDEA-S1 LPM "deviation" value
#' @return LP0 = DDEA LP object for last dmu0 processed
#' @return LP1 = qDDEA-S1 LP object for last dmu0 processed
#' @return LP2 = qDDEA-S2 LP object for last dmu0 processed
##############################################################################
# debugging qDEA_solve
##############################################################################
# X0=X00;Y0=Y00;DX0=DX00;DY0=DY00
##############################################################################
# (qout = 2/8); (qoutS = 1/8)
# (qout = 1/8); (qoutS = 0/8)
##############################################################################
##############################################################################
# X=Xm;Y=Ym
##############################################################################
# (qoutS = qout);(RTS='VRS');(nqiter=4);(qtol=1e-6);(BIGM=1E9);(eps=1e-6)
# (skipzprob=TRUE);(unbounded = (-1E3));(obj2test=1e-4);(printlog=TRUE)
# (prntmod=1); (printtxt=''); (replaceA1=TRUE)
# (baseqDEA=FALSE)
##############################################################################
#' @export
###############################################################################
qDEA_solve=function(X,Y,qout=0.10,qoutS=qout,X0,Y0,DX0,DY0,RTS='CRS',
nqiter=1,qtol=1E-6,BIGM=1E9,eps=1E-6,skipzprob=TRUE,unbounded=(-1E3),
obj2test=1E-4,replaceA1=FALSE,baseqDEA=FALSE,printlog=TRUE,prntmod=100,printtxt='',
dmulist0=1:nrow(X0),solver='highs'){
###############################################################################
if(nrow(X)!=nrow(Y)) stop('Error: nrow(X)!= nrow(Y)')
if(nrow(X0)!=nrow(Y0)) stop('Error: nrow(X0)!= nrow(Y0)')
if(ncol(X)!=ncol(X0)) stop('Error: ncol(X)!=ncol(X0)')
if(ncol(Y)!=ncol(Y0)) stop('Error: ncol(Y)!=ncol(Y0)')
if(nrow(X0)!=nrow(DX0)) stop('Error: nrow(X0)!=nrow(DX0)')
if(nrow(Y0)!=nrow(DY0)) stop('Error: nrow(Y0)!=nrow(DY0)')
if(ncol(X0)!=ncol(DX0)) stop('Error: ncol(X0)!=ncol(DX0)')
if(ncol(Y0)!=ncol(DY0)) stop('Error: ncol(Y0)!=ncol(DY0)')
###############################################################################
(ncY=ncol(Y));(ncX=ncol(X))
(ncYX=ncY+ncX)
(nobs=nrow(Y))
(devstart=ncYX+2)
(devend=devstart+nobs-1)
(nDMU0=nrow(X0))
###############################################################################
###############################################################################
# LP output data
D2OUT=matrix(NA,nDMU0,nobs)
effvals=matrix(NA,nDMU0,3)
status=matrix(NA,nDMU0,3)
qhat=rep(NA,nDMU0)
qiter=rep(0,nDMU0)
#####################################################################
LP0='NOT RUN'; LP1='NOT RUN'; LP2='NOT RUN'
PSOL='NOT RUN'; PSOLq1='NOT RUN';PSOLq2='NOT RUN'
PEERS='NOT RUN'; DSOL1='NOT RUN';PEERSq='NOT RUN'
RCOST='NOT RUN';RCOSTq1='NOT RUN';RCOSTq2='NOT RUN'
#####################################################################
###################################################################
# SETUP SECTIONS
###################################################################
# debug off
# set up DDEA dual model
time0=my_seconds()
LP0=DEAbuild(X,Y,X0,Y0,DX0,DY0,dmu0=1,RTS=RTS,
unbounded=unbounded,solver=solver)
# Matrices to store primal and dual solutions
PSOL=matrix(NA,nDMU0,LP0$nc)
PEERS=matrix(NA,nDMU0,LP0$nr)
RCOST=matrix(NA,nDMU0,LP0$nc)
#####################################################################
#####################################################################
# set up for qDEA
#debug off
if(qout>=qtol|skipzprob==FALSE) {
# set up qDEA LP models with SMM CMO
LP1=qDEAbuild(X,Y,X0,Y0,DX0,DY0,dmu0=1,qout=qoutS,RTS=RTS,
unbounded=unbounded,solver=solver)
LP2=DEAbuild(X,Y,X0,Y0,DX0,DY0,dmu0=1,RTS=RTS,
unbounded=unbounded,solver=solver)
# Matrices to store primal and dual solutions
PSOLq1=matrix(NA,nDMU0,LP1$nc)
DSOL1=matrix(NA,nDMU0,LP1$nr)
RCOSTq1=matrix(NA,nDMU0,LP1$nc)
PSOLq2=matrix(NA,nDMU0,LP2$nc)
PEERSq=matrix(NA,nDMU0,LP2$nr)
RCOSTq2=matrix(NA,nDMU0,LP2$nc)
tau0=LP1$tau0
(tau0inv=LP1$ra[LP1$qchng[1]])
# save initial values for row boundaries
rhs1=LP1$rhs
rhs2=LP2$rhs
rlower1=LP1$rlower
rlower2=LP2$rlower
rupper1=LP1$rupper
rupper2=LP2$rupper
###############################################################################
} # end setup #if(qoutS>=qtol|skipzprob==FALSE) {
###############################################################################
#####################################################################
# run dual DDEA and qDEA models
#####################################################################
i=1
# debug off
for(i in 1:nDMU0){
if(printlog==TRUE){
if(i==1|i%%prntmod==0|i==nDMU0) print(paste(printtxt,'on dmu',i,'of',nDMU0))
} # end if(printlog==TRUE)
# DDEA DUAL section
# modify lp objective and matrix coefficients
(LP0$obj[1:ncYX]=LP0$DOBJ[i,])
LP0$ra[LP0$dyxchngC]= LP0$DYX[i,]
LP0$ra[LP0$yxchngR] = -LP0$DOBJ[i,]
if(replaceA1==TRUE){
yx = -LP0$DOBJ[i,]
LP0$ra[LP0$yxchngC]=yx
}
LP0$dmu0=i
LPSOL0=LPSOLVER(LP0,solver=solver)
names(LPSOL0)
# [1] "status" "objval" "solution" "dual" "rcost" "lhs"
# [7] "TIME" "SMtime" "proctime" "LP"
(status[i,1]=LPSOL0$status)
(obj0=LPSOL0$objval)
(effvals[i,1]=obj0)
(PSOL[i,]=LPSOL0$solution)
(PEERS[i,]=LPSOL0$dual)
(RCOST[i,]=LPSOL0$rcost)
if(obj0<=unbounded+obj2test){
effvals[i,1]=NA
PSOL[i,] = NA
PEERS[i,] = NA
RCOST[i,]= NA
} # end if(obj0<=unbounded+obj2test)
LP2 = LP0
#############################################
# end dual DDEA LP section
#############################################
#############################################
# start qDEA section
#############################################
###############################################
# debugoff
if(qoutS>=qtol|skipzprob==FALSE) {
# modify lp objectives
(LP1$obj[1:ncYX]=LP1$DOBJ[i,])
(LP2$obj[1:ncYX]=LP2$DOBJ[i,])
# modify/restore lp matrix coefficients
# restore original tau0inv
LP1$ra[LP1$qchng[1]]=tau0inv
# modify with DMU i's directions
dyx = LP1$DYX[i,]
dobj = -LP1$DOBJ[i,]
LP1$ra[LP1$dyxchngC] = dyx
LP1$ra[LP1$yxchngR] = dobj
LP2$ra[LP2$dyxchngC] = dyx
LP2$ra[LP2$yxchngR] = dobj
if(replaceA1==TRUE){
# yx = c(Y0[i,],-X0[i,])
yx = -LP1$DOBJ[i,]
LP1$ra[LP1$yxchngC]=yx
LP2$ra[LP2$yxchngC]=yx
}
(LP1$dmu0=i)
(LP2$dmu0=i)
# restore original lp rhs
LP1$rupper=rupper1
LP2$rupper=rupper2
# qDEA stage 1
(pointsout0=floor(qoutS*nobs))
LPSOL1=LPSOLVER(LP1,solver=solver)
names(LPSOL1)
# [1] "status" "objval" "solution" "dual" "rcost" "lhs"
# [7] "TIME" "SMtime" "proctime" "LP"
(status[i,2]=LPSOL1$status)
(obj1=LPSOL1$objval)
effvals[i,2]=obj1
PSOLq1[i,]=LPSOL1$solution
DSOL1[i,]=LPSOL1$dual
RCOSTq1[i,]=LPSOL1$rcost
xvals1=LPSOL1$solution
devs1=xvals1[devstart:devend]
D1out=ifelse(devs1>eps,1,0)
sum(D1out)
# qDEA stage 2
LP2$rhs[1:nobs]=D1out*BIGM;
LP2$rupper[1:nobs]=D1out*BIGM;
LPSOL2=LPSOLVER(LP2,solver=solver)
names(LPSOL2)
# [1] "status" "objval" "solution" "dual" "rcost" "lhs"
# [7] "TIME" "SMtime" "proctime" "LP"
(status[i,3]=LPSOL2$status)
(obj2=LPSOL2$objval)
effvals[i,3]=obj2
PSOLq2[i,]=LPSOL2$solution
(PEERSq[i,]=LPSOL2$dual)
RCOSTq2[i,]=LPSOL2$rcost
(xvals2=LPSOL2$solution)
LHS2=LPSOL2$lhs
D2out=ifelse(LHS2[1:nobs]>eps,1,0)
sum(D2out)
(qhat[i]=sum(D2out)/nobs)
qiter[i]=1
(qleft= qoutS-qhat[i])
(pointsout=as.integer(qleft*nobs))
obj2diff = 1E3
# debug off
if(nqiter>1&abs(qoutS-qhat[i])>qtol) {
# debugoff
while(qiter[i]<nqiter & pointsout>0&qleft>eps) {
(qleft= qoutS-qhat[i])
floor(qleft*nobs)
########################################################################
(pointsout=floor(qleft*nobs))
########################################################################
(q0 = qleft-eps)
(q1 = (floor(qleft*nobs)+1)/nobs - eps)
floor(q1*nobs)
tau2 = q1
(tau2inv=1/tau2)
(LP1$ra[LP1$qchng[1]]=tau2inv)
##############################################
LP1$rhs[1:nobs]=D2out*BIGM;
LP1$rupper[1:nobs]=D2out*BIGM;
LPSOL1=LPSOLVER(LP1,solver=solver)
names(LPSOL1)
# [1] "status" "objval" "solution" "dual" "rcost" "lhs"
# [7] "TIME" "SMtime" "proctime" "LP"
(status2=LPSOL1$status)
(obj1=LPSOL1$objval)
xvals1=LPSOL1$solution
devs1=xvals1[devstart:devend]
D1out=ifelse(devs1>eps,1,0)
sum(D1out)
# qDEA-stage2
tmp=D2out+D1out
tmp2=ifelse(tmp>0,1,0)
sum(tmp2)
LP2$rhs[1:nobs]=tmp*BIGM;
LP2$rupper[1:nobs]=tmp*BIGM;
LPSOL2=LPSOLVER(LP2,solver=solver)
names(LPSOL2)
# [1] "status" "objval" "solution" "dual" "rcost" "lhs"
# [7] "TIME" "SMtime" "proctime" "LP"
(status3=LPSOL2$status)
(obj2=LPSOL2$objval)
# ###########################################################################
# debugoff
if(baseqDEA!=TRUE) {
if(abs(obj2-effvals[i,3])<=eps & pointsout>=1){
xvals2=LPSOL2$solution
LHS2=LPSOL2$lhs
psol2=LPSOL2$solution
dsol2=LPSOL2$dual
d2out=ifelse(LHS2[1:nobs]>eps,1,0)
LPlist=list(NULL)
(dualList=which(abs(dsol2[1:nobs])>eps))
(nduals=length(dualList))
# debugoff
if(nduals>=1){
qefftest=100
nd=1
#debugoff
for(nd in 1:nduals){
LP22=LP2
(ndpick=dualList[nd])
LP22$rhs[ndpick]=BIGM
LP22$rupper[ndpick]=BIGM
LPSOL=LPSOLVER(LP22,solver=solver)
(qefftest[nd]=LPSOL$objval)
LPlist[[nd]]=list(LP=LP22)
} # end loop for(nd in 1:nduals)
(ndpick=which(qefftest==min(qefftest))[1])
(obj22=qefftest[ndpick])
#debugoff
if(obj22<(obj2-eps)){
LP2=LPlist[[ndpick]]$LP
LPSOL2=LPSOLVER(LP2,solver=solver)
status3=LPSOL2$status
obj2=LPSOL2$objval
} # end if(obj22<(obj2-eps)){
} # end loop if(length(nduals)>=1)
} # end if(abs(obj2-effvals[i,3])<=eps)
} # end if(baseqDEA!=TRUE)
############################################################################
(obj2diff=abs(obj2-effvals[i,3]))
xvals2=LPSOL2$solution
LHS2=LPSOL2$lhs
(qiter[i]=qiter[i]+1)
if(obj2diff>obj2test){
(status[i,2]=status2)
(status[i,3]=status3)
(effvals[i,2]=obj1)
(effvals[i,3]=obj2)
(PSOLq2[i,]=LPSOL2$solution)
(PEERSq[i,]=LPSOL2$dual)
RCOSTq2[i,]=LPSOL2$rc
D2out=ifelse(LHS2[1:nobs]>eps,1,0)
(qhat[i]=sum(D2out)/nobs)
} # end if(obj2diff>obj2test)
(qleft= qoutS-qhat[i])
(floor(qleft*nobs))
(pointsout=floor(qleft*nobs))
} # end #while(qiter[i]<nqiter & abs(qoutS-qhat[i])>qtol & pointsout>0 & obj2diff>eps) {
}# end if(nqiter>1&abs(qhat[i]-qoutS)>qtol)
D2OUT[i,]=D2out
} # end if(qoutS>(1.0/nobs)|skipzprob==FALSE)
###############################################################################
# if(qout>(qoutS+eps)) {
if(baseqDEA==FALSE&qout>(qoutS+eps)) {
(ndmus=nrow(X))
(nout0=floor(qout*ndmus))
tmp2=iter_delete(LP2,ndmus,nout0,iterlim=max(100,nout0),
BIGM=BIGM,unbounded=unbounded,eps=eps,solver=solver)
LP2=tmp2$LP2
LP2sol=tmp2$LP2sol
(status[i,3]=LP2sol$status)
(obj2=tmp2$obj2)
effvals[i,3]=obj2
if(!is.na(obj2)){
PSOLq2[i,]=LP2sol$solution
PEERSq[i,]=LP2sol$dual
RCOSTq2[i,]=LP2sol$rcost
(xvals2=LP2sol$solution)
(outlist=tmp2$outlist)
D2out=rep(0,ndmus)
D2out[outlist]=1
sum(D2out)
(qhat[i]=sum(D2out)/nobs)
D2OUT[i,]=D2out
} # end if(qout>0)
(qiter[i]=qiter[i]+tmp2$iter)
} # end if(qout>(qoutS+eps))
################################################################################
if(qout>=qtol|skipzprob==FALSE) {
if(is.na(obj2)|obj2<=unbounded+obj2test){
PSOLq2[i,] = NA
PEERSq[i,] = NA
RCOSTq2[i,]= NA
# D2OUT[i,] = NA
} # end if(obj2<=unbounded)
} # end if(qout>=qtol|skipzprob==FALSE)
} # end loop on i
###############################################################################
(status[,1]=ifelse(effvals[,1]<=unbounded+obj2test,NA,status[,1]))
(effvals[,1]=ifelse(effvals[,1]<=unbounded+obj2test,NA,effvals[,1]))
(status[,2]=ifelse(effvals[,3]<=unbounded+obj2test,NA,status[,2]))
(status[,3]=ifelse(effvals[,3]<=unbounded+obj2test,NA,status[,3]))
(effvals[,2]=ifelse(effvals[,3]<=unbounded+obj2test,NA,effvals[,2]))
(effvals[,3]=ifelse(effvals[,3]<=unbounded+obj2test,NA,effvals[,3]))
(napick= which(is.na(effvals[,3])))
(qhat[napick]=NA)
#####
(tmp=PEERS)
(PEERS=A2SM(matrix(PEERS[,1:nrow(X)],nrow=nrow(PEERS)),SMM='RCI',NA_flag=-1E12))
(PEERS=as.data.frame(cbind(PEERS$ia,PEERS$ja,-1*PEERS$ra)))
names(PEERS)=c('dmu0','dmuz','z')
if(replaceA1==TRUE) PEERS$dmuz=ifelse(PEERS$dmuz==1,0,PEERS$dmuz)
PEERS$dmuz = ifelse(is.na(PEERS$z),NA,PEERS$dmuz)
(PEERS = unique(PEERS))
if(!is.null(dmulist0)) PEERS$dmu0=dmulist0[PEERS$dmu0]
PEERS
#######################
if(qout>=qtol|skipzprob==FALSE) {
## process peer data
(tmpq=PEERSq)
PEERSq[is.na(PEERSq)]=0
PEERSq=A2SM(matrix(PEERSq[,1:nrow(X)],nrow=nrow(PEERSq)),SMM='RCI',NA_flag=-1E12)
(PEERSq=as.data.frame(cbind(PEERSq$ia,PEERSq$ja,-1*PEERSq$ra)))
names(PEERSq)=c('dmu0','dmuz','z')
if(replaceA1==TRUE) PEERSq$dmuz=ifelse(PEERSq$dmuz==1,0,PEERSq$dmuz)
PEERSq$dmuz = ifelse(is.na(PEERSq$z),NA,PEERSq$dmuz)
(PEERSq = unique(PEERSq))
if(!is.null(dmulist0)) PEERSq$dmu0=dmulist0[PEERSq$dmu0]
PEERSq
## process D2OUT data
tmpD2=D2OUT
D2OUT[is.na(D2OUT)]=0
D2OUT=A2SM(as.matrix(D2OUT),SMM='RCI',NA_flag=-1E12)
D2OUT=as.data.frame(cbind(D2OUT$ia,D2OUT$ja,D2OUT$ra))
names(D2OUT)=c('dmu0','dmuout','dout')
if(replaceA1==TRUE) D2OUT$dmuout=ifelse(D2OUT$dmuout==1,0,D2OUT$dmuout)
D2OUT$dmuout = ifelse(is.na(D2OUT$dout),NA,D2OUT$dmuout)
D2OUT = unique(D2OUT)
D2OUT = D2OUT[,c('dmu0','dmuout')]
} # end if(qout>=qtol|skipzprob==FALSE)
#######################
###############################################################################
list(effvals=effvals,qhat=qhat,qiter=qiter,D2OUT=D2OUT,status=status,
PSOL=PSOL,PSOLq1=PSOLq1,PSOLq2=PSOLq2,PEERS=PEERS,PEERSq=PEERSq,
RCOST=RCOST,RCOSTq1=RCOSTq1,RCOSTq2=RCOSTq2,devstart=devstart,devend=devend,
LP0=LP0,LP1=LP1,LP2=LP2)
###############################################################################
} # end function qDEA_solve
###############################################################################
###############################################################################
#' qDEA_mlist: Obtain subsample qDEA results for a set (vector) of subsample sizes
#' !!!Intended to be called from qDEA function!!!!
#'
#' mcells = number of subsample sizes
#' @param XM Input levels for subsample of reference set X.
#' @param YM Output levels for subsample of reference set Y.
#' @param mpick Index of subsample dmus in original (X,Y) reference set.
#' @param qout Maximal proportion of dmu's allowed external to DEA hull.
#' @param qoutS Proportion of external points to identify using qDEA slicing.
#' @param Bdrop Index list of DMU's whose initial DDEA solutions were unbounded.
#' @param Bdropq Index list of DMU's whose initial qDEA solutions were unbounded.
#' @param X0 Inputs for set of ndmu0 dmu's to be processed.
#' @param Y0 Outputs for set of ndmu0 dmu's to be processed.
#' @param DX0 Input directions for ndmu0 dmu's in X0 and Y0.
#' @param DY0 Output directions for ndmu0 dmu's in X0 and Y0.
#' @param mlist Vector of subsample sizes
#' @param RTS Returns to scale default='CRS' options are 'CRS,'VRS,'DRS','IRS.
#' @param nqiter Maximal number of qDEA iterations.
#' @param qtol qout search tolerance with iterative qDEA.
#' @param BIGM Default Big M in RHS of qDEA stage 2 process.
#' @param eps Search tolerance in qDEA improvement tests.
#' @param skipzprob Skip qDEA if qout=0.
#' @param unbounded qDEA reported as unbounded if obj<=unbounded.
#' @param obj2test Converge tol for objective in iterative qDEA.
#' @param replaceA1 Put dmu0's data in first row of reference sets. .
#' @param baseqDEA Use basic qDEA model from EJOR article
#' @param printlog Progress of dmu's solved when (X0,Y0) >1 dmus.
#' @param prntmod Print progress every prntmod dmus.
#' @param getbootpeers Return dmu projection weights for each dmu and each bootstrap.
#' @param dmulist0 DMU0 index in originally inputs X0,Y0,DX0,and DY0
#' @param solver LP solver Default='highs'
#' @return A list containing the following components:
#' @return EFFmlist = ndmu0 by mcell matrix of subsample DDEA distances.
#' @return EFFmlistq = ndmu0 by mcell matrix of subsample qDEA distances.
#' @return XM = Input levels for subsample of reference set X.
#' @return YM = Output levels for subsample of reference set X.
#' @return qout = Maximal proportion of dmu's allowed external to DEA hull.
#' @return qoutS = Proportion of external points to identify using qDEA slicing.
#' @return X0 = Inputs for set of ndmu0 dmu's to be processed.
#' @return Y0 = Outputs for set of ndmu0 dmu's to be processed.
#' @return DX0 = Input directions for ndmu0 dmu's in X0 and Y0.
#' @return DY0 = Output directions for ndmu0 dmu's in X0 and Y0.
#' @return mlist = Vector of subsample sizes
#' @return mcells = Number of subsample sizes
#' @return peers = A data frame containing DDEA peer dmus and projection weights
#' @return peersq = A data frame containing qDEA peer dmus and projection weights
###############################################################################
# debug qDEA_mlist
###############################################################################
# XM=XM;YM=YM;X0=X0;Y0=Y0;DX0=DX0;DY0=DY0
# (qout=qout); (qoutS=qoutS)
# Bdrop=Bdrop;Bdropq=Bdropq
##############################################################################
# RTS=RTS;nqiter=nqiter;qtol=qtol;BIGM=BIGM;eps=eps;skipzprob=TRUE;unbounded=unbounded
# obj2test=obj2test;printlog=printlog;prntmod=prntmod;printtxt=printtxt;replaceA1=FALSE
# getbootpeers=getbootpeers
###############################################################################
#' @export
###############################################################################
qDEA_mlist=function(XM,YM,mpick,qout=0.10,qoutS=qout,Bdrop=NULL,Bdropq=NULL,
X0=XM,Y0=YM,DX0=XM,DY0=YM,mlist,RTS='CRS',nqiter=1,qtol=1e-6,BIGM=1E9,eps=1e-6,
skipzprob=TRUE,unbounded=(-1E3),obj2test=1e-4,replaceA1=replaceA1,baseqDEA=baseqDEA,
printlog=FALSE,prntmod=100,getbootpeers=FALSE,dmulist0=1:nrow(X0),solver='highs') {
###############################################################################
#setDTthreads(1)
###############################################################################
if(nrow(XM)!=nrow(YM)) stop('Error: nrow(X)!= nrow(Y)')
if(nrow(X0)!=nrow(Y0)) stop('Error: nrow(X0)!= nrow(Y0)')
if(ncol(XM)!=ncol(X0)) stop('Error: ncol(X)!=ncol(X0)')
if(ncol(YM)!=ncol(Y0)) stop('Error: ncol(Y)!=ncol(Y0)')
if(nrow(X0)!=nrow(DX0)) stop('Error: nrow(X0)!=nrow(DX0)')
if(nrow(Y0)!=nrow(DY0)) stop('Error: nrow(Y0)!=nrow(DY0)')
if(ncol(X0)!=ncol(DX0)) stop('Error: ncol(X0)!=ncol(DX0)')
if(ncol(Y0)!=ncol(DY0)) stop('Error: ncol(Y0)!=ncol(DY0)')
###############################################################################
(mcells=length(mlist))
(ncY=ncol(YM));(ncX=ncol(XM))
(ncYX=ncY+ncX)
(nDMU0=nrow(X0))
EFFmlist=matrix(NA,nDMU0,mcells)
EFFmlistq=matrix(NA,nDMU0,mcells)
# data frame for matching dmu index from XM with dmu index from original X
dmuset = data.frame(cbind(dmuz=mpick,dmuzm=1:length(mpick)))
#head(dmuset)
#data frame for putting unbounded dmus back in PEERSq object with NA entries
dmu0_ub = data.frame(cbind(dmu0=1:nDMU0,ub=rep(NA,nDMU0)))
#head(dmu_ub)
peers = 'NOT RUN'
peersq = 'NOT RUN'
if(getbootpeers==TRUE){
peers = list('NOT RUN')
peersq =list('NOT RUN')
}
#######################################
j=1
# debugoff
for(j in 1:mcells){
(m=mlist[j])
Ym=as.matrix(YM[1:m,]);Xm=as.matrix(XM[1:m,])
tmp=qDEA_solve(X=Xm,Y=Ym,qout=qout,qoutS=qoutS,X0=X0,Y0=Y0,DX0=DX0,DY0=DY0,
RTS=RTS,nqiter=nqiter,qtol=qtol,BIGM=BIGM,eps=eps,skipzprob=skipzprob,
unbounded=unbounded,obj2test=obj2test,printlog=printlog,
prntmod=prntmod,replaceA1=replaceA1,baseqDEA=baseqDEA,dmulist0=NULL,solver=solver)
EFFmlist[,j]=tmp$effvals[,1]
EFFmlistq[,j]=tmp$effvals[,3]
# debugoff
if(getbootpeers==TRUE){
PEERS=tmp$PEERS
names(PEERS) = c('dmu0','dmuzm','z')
PEERS=merge(PEERS,dmuset,all.x=TRUE)
PEERS$dmuz=ifelse(is.na(PEERS$dmuz),PEERS$dmu0,PEERS$dmuz)
PEERS$m = m
PEERS=PEERS[,c('m','dmu0','dmuz','z')]
PEERS=doBy::orderBy(~dmu0+dmuz,data=PEERS)
########################
PEERSq=tmp$PEERSq
names(PEERSq) = c('dmu0','dmuzm','z')
PEERSq=merge(PEERSq,dmuset,all.x=TRUE)
PEERSq$dmuz=ifelse(is.na(PEERSq$dmuz),PEERSq$dmu0,PEERSq$dmuz)
PEERSq$m = m
PEERSq=PEERSq[,c('m','dmu0','dmuz','z')]
PEERSq=doBy::orderBy(~dmu0+dmuz,data=PEERSq)
########################
peers[[j]] = PEERS
peersq[[j]] = PEERSq
#########################
} # end if(getbootpeers==TRUE)
} # end loop on j
if(getbootpeers==TRUE){
peers = dplyr::bind_rows(peers)
peersq = dplyr::bind_rows(peersq)
if(length(Bdrop)>=1){
EFFmlist[Bdrop,] = NA
pick1=match(peers$dmu0,Bdrop)
pick2=which(!is.na(pick1))
peers$dmuz[pick2] = NA
peers$z[pick2] = NA
peers=unique(peers)
} # end if(length(Bdrop)>=1)
if(length(Bdropq)>=1) {
EFFmlistq[Bdropq,] = NA
pick1=match(peersq$dmu0,Bdropq)
pick2=which(!is.na(pick1))
peersq$dmuz[pick2] = NA
peersq$z[pick2] = NA
peersq=unique(peersq)
} # end if(length(Bdropq)>=1)
} #end if(getbootpeers==TRUE)
peers=stats::na.omit(peers)
peersq=stats::na.omit(peersq)
#################################################################################
list(EFFmlist=EFFmlist,EFFmlistq=EFFmlistq,XM=XM,YM=YM,qout=qout,qoutS=qoutS,
X0=X0,Y0=Y0,DX0=DX0,DY0=DY0,mlist=mlist,mcells=mcells,peers=peers,peersq=peersq)
###############################################################################
} # end function qDEA_mlist
###############################################################################
###############################################################################
#' nCM_mpick: Select subsample size m from mlist
#' !!!!Intended to be called from qDEA function!!!!
#'
#' Pulls associated set of nboot bootstrapped values and computes bias corected
#' "s(m)" values.
#' ** Uses Simar and Wilson (2011) suggested procedure for selecting m **
#' @param stat DDEA or qDEA distance point estimate
#' @param boot nboot by mcells matrix of bootstrapped DDEA and qDEA distance
#' @param n Number of dmus in full reference set (X,Y)
#' @param mlist Vector of subsample sizes
#' @param beta Convergence rate. 0.5 => "root n" convergence
#' @param alpha Alpha level for Simar-Wilson (2011) subsample size selection procedure
#' @param CILag CILag level for Simar-Wilson (2011) subsample size selection procedure
#' @return mpick = selected subsample size,
#' @return S = nboot by mcell matrix of
#' s = stat(n) - (m/n)^beta(boot(m)-stat(n)
#' bias corrected values for each sample size m
#' @return s = vector (of length nboot) of 'bias.corected' s-values for
#' sample size 'mpick'
#' @return stat.bc = bias corrected point estimate of stat = mean(s)
###############################################################################
#' @export
###############################################################################
nCm_mpick=function(stat,boot,n,mlist,beta=0.5,alpha=0.05,CILag=1){
###############################################################################
boot=as.matrix(boot)
(mcells=length(mlist))
if(mcells!=ncol(boot)) stop('length of mlist != ncol(boot)')
###############################################################################
# Estimate confidence intervals
(m = mlist)
(nM = length(m))
S = 0 * boot
for (j in 1:nM) {
S[, j] = stat- ((m[j]/n)^beta) * (boot[, j] - stat)
}
UL = apply(S, 2, stats::quantile, probs = (1 - alpha/2),na.rm=TRUE)
LL = apply(S, 2, stats::quantile, probs = (alpha/2),na.rm=TRUE)
summary(UL)
summary(LL)
CL = stat - UL
CU = stat - LL
summary(CL)
summary(CU)
if(nM>2){
MCL = lagMat(CL, -CILag:CILag)
MCU = lagMat(CU, -CILag:CILag)
CLsd = apply(MCL, 1, stats::sd)
CUsd = apply(MCU, 1, stats::sd)
Csd = CLsd + CUsd
(mpick0 = which(Csd == suppressWarnings(min(Csd, na.rm=TRUE))))
} # end if(nM>2)
if(nM==1){
mpick0=1
} # end if(nM==1)
mpick=mpick0[1]
s=S[,mpick]
###############################
list(mpick=mpick,S=S,s=s,stat.bc=mean(s))
###############################################################################
} # end function nCm_mpick
###############################################################################
###############################################################################
#' DEAbuild: Builds DDEA LP object for use in qDEA_solve function
#'
#' @param X Reference dmu's = ndmu x number of inputs input matrix.
#' @param Y Reference dmu's = ndmu x number of outputs output matrix.
#' @param X0 Inputs for set of ndmu0 dmu's to be processed.
#' @param Y0 Outputs for set of ndmu0 dmu's to be processed.
#' @param DX0 Input directions for ndmu0 dmu's in X0 and Y0.
#' @param DY0 Output directions for ndmu0 dmu's in X0 and Y0.
#' @param dmu0 Row in (X0,Y0,DX0,DY0) to use as given dmu
#' @param RTS Returns to scale: 'CRS,'VRS,'DRS','IRS.
#' @param unbounded DEA obj restricted >= to unbounded. Default = -1E3
#' @param solver LP solver Default='highs'
#' @return Returns an LP list object containing the following elements:
#' @return LPsense = 'max' or 'min'
#' @return nnz = number of nonzero elements in 'A' matrix
#' @return nr = number or rows in 'A' matrix
#' @return nc = number of columns in 'A' matrix
#' @return obj = nc length vector of objective coefficients
#' @return ra = nonzero coefficients in 'A' matrix
#' @return ia = row indexes for non zero elements in 'A' matrix
#' @return ja = column indexes for non zero elements in 'A' matrix
#' @return SMM = Sparse matrix method: 'CMO','RMO','CRI',or 'RCI'
#' @return ZINDEX = 'zero indexing' (T) or 'one indexing' (F)
#' @return dirs = nr length vector of constraint signs ('<=','>=', or '=')
#' @return rhs = nr length vector of RHS coefficients
#' @return xlower = nc vector of lower bounds on 'x' choice variables
#' @return xupper = nc vector of upper bounds on 'x' choice variables
#' @return rlower = nr vector of Ax lower bounds i.e. rlower <= Ax
#' @return rupper = nr vector of Ax upper bounds i.e. Ax <= rupper
#' @return vartypes = nc vector of variable types: ('C','B','I') -- qDEA:rep('C',nc)
#' @return yxchngC = index of given dmu's output-input values locations in ra vector
#' (used to edit LP problem as we loop through DMU's in (X0,Y0,DX0,DY0)
#' @return dyxchngC = index of given dmu's direction values locations in ra vector
#' (used to edit LP problem as we loop through DMU's in (X0,Y0,DX0,DY0)
#' @return yxchngR = index of given dmu's obj restriction output-input values locations in ra vector
#' (used to edit LP problem as we loop through DMU's in (X0,Y0,DX0,DY0)
#' @return iyxchngC = row indexes associated with yxchngC cells
#' @return idyxchngC = row indexes associated with dyxchngC cells
#' @return iyxchngR = row indexes associated with yxchngR cells
#' @return DOBJ = matrix of obj values to change by dmu with DOBJ=cbind(-Y0,X0)
#' @return DYX = matrix of (dy,dx) values to change by dmu with DYX=cbind(DY0,DX0)
#' @return RTS = Returns to scale: 'CRS,'VRS,'DRS','IRS.
#' @return dmu0 = Row in (X0,Y0,DX0,DY0) to use as given dmu
###############################################################################
# debugging DEAbuild
###############################################################################
# debug on
# (dmu0=1); (RTS='VRS'); (SMM='CMO'); (ZINDEX=TRUE)
###############################################################################
#' @export
###############################################################################
DEAbuild=function(X,Y,X0,Y0,DX0,DY0,dmu0=1,RTS='CRS',
unbounded=-1E3,solver='highs'){
###############################################################################
# obj restriction
SMM='CMO'; ZINDEX=TRUE
(nrY=nrow(Y));(nrX=nrow(X))
(ncY=ncol(Y));(ncX=ncol(X))
if(nrY!=nrX) stop('nrY!=nrX')
DOBJ=cbind(-Y0,X0)
DYX=cbind(DY0,DX0)
y0=Y0[dmu0,];x0=X0[dmu0,];dy=DY0[dmu0,];dx=DX0[dmu0,]
(ncYX=ncY+ncX)
(nc=ncY+ncX)
(nr=nrY+2)
# Note: At this stage ja contains "columns counts"
# Include obj restriction
(M1=as.matrix(as.matrix(rbind(Y,dy,y0))))
(ra=as.vector(M1))
(ia=rep(1:(nrY+2),ncY))
# Note: At this stage ja contains "columns counts"
(ja=rep(nrY+2,ncY))
(M2=as.matrix(as.matrix(rbind(-X,dx,-x0))))
(ra=c(ra,as.vector(M2)))
(ia=c(ia,rep(1:(nrX+2),ncX)))
(ja=c(ja,rep((nrX+2),ncX)))
if(RTS!='CRS'&RTS!='crs'){
(ra=c(ra,rep(-1,nrY),-1))
(ia=c(ia,1:nrY,nrY+2))
# Note: At this stage ja contains "columns counts"
(ja=c(ja,nrY+1))
}
# compute CMO ja vector from counts
(ja=cumsum(c(1,ja)))
# change locations in constraint matrix
(yxchngC = cumsum(c(1,rep((nrY+2),(ncYX-1)))))
(dyxchngC = cumsum(c(nrY+1,rep((nrY+2),(ncYX-1)))))
(yxchngR = cumsum(c(nrY+2,rep((nrY+2),(ncYX-1)))))
###########################################################
LPsense='min'
(obj=c(-y0,x0))
(xlower=rep(0,nc))
(xupper=rep(1E8,nc))
(dirs=c(rep('<=',nrY),'=','<='))
(rhs=c(rep(0,nrY),1,-unbounded))
if(RTS!='CRS'&RTS!='crs'){
nc=nc+1
obj=c(obj,1)
if(RTS=='VRS'|RTS=='vrs'){
xlower=c(xlower,-1E8)
xupper=c(xupper,1E8)
}
if(RTS=='DRS'|RTS=='drs'){
xlower=c(xlower,0)
xupper=c(xupper,1E8)
}
if(RTS=='IRS'|RTS=='irs'){
xlower=c(xlower,0)
xupper=c(xupper,1E8)
}
} # end if(RTS!='CRS'&RTS!='crs')
rlower=rhs
rupper=rhs
rlower=ifelse(dirs=='<='|dirs=='<',-1E8,rlower)
rupper=ifelse(dirs=='<='|dirs=='<',rhs,rupper)
rlower=ifelse(dirs=='>='|dirs=='>',rhs,rlower)
rupper=ifelse(dirs=='>='|dirs=='>',1E8,rupper)
rlower=ifelse(dirs=='=='|dirs=='=',rhs,rlower)
rupper=ifelse(dirs=='=='|dirs=='=',rhs+1E-8,rupper)
rlower;rupper
###############################################################################
iyxchngC=ia[yxchngC]
idyxchngC=ia[dyxchngC]
iyxchngR=ia[yxchngR]
###############################################################################
if(ZINDEX==TRUE) {ia=ia-1;ja=ja-1}
vartypes=rep("C",nc)
###############################################################################
###############################################################################
LP=list(LPsense=LPsense,nnz=length(ra),nr=nr,nc=nc,obj=obj,ra=ra,ia=ia,ja=ja,
SMM=SMM,ZINDEX=ZINDEX,dirs=dirs,rhs=rhs,xlower=xlower,xupper=xupper,
rlower=rlower,rupper=rupper,vartypes=vartypes,
yxchngC=yxchngC,dyxchngC=dyxchngC,yxchngR=yxchngR,iyxchngC=iyxchngC,
idyxchngC=idyxchngC,iyxchngR=iyxchngR,DOBJ=DOBJ,DYX=DYX,RTS=RTS,dmu0=dmu0)
###############################################################################
LP=SM2SM(LP,SMM2='CRI',ZINDEX2=FALSE)
###############################################################################
LP
###############################################################################
} # end function DEAbuild
###############################################################################
###############################################################################
#' qDEAbuild: Builds qDEA LP object for use in qDEA_solve function
#'
#' @param X Reference dmu's = ndmu x number of inputs input matrix.
#' @param Y Reference dmu's = ndmu x number of outputs output matrix.
#' @param X0 Inputs for set of ndmu0 dmu's to be processed.
#' @param Y0 Outputs for set of ndmu0 dmu's to be processed.
#' @param DX0 Input directions for ndmu0 dmu's in X0 and Y0.
#' @param DY0 Output directions for ndmu0 dmu's in X0 and Y0.
#' @param qout Maximal proportion of dmu's allowed external to DEA hull.
#' @param dmu0 Row in (X0,Y0,DX0,DY0) to use as given dmu
#' @param RTS Returns to scale: 'CRS,'VRS,'DRS','IRS.
#' @param unbounded DEA obj restricted >= to unbounded. Default = -1E3
#' @param solver LP solver Default='highs'
#' @return Returns an LP list object containing the following elements:
#' @return LPsense = 'max' or 'min'
#' @return nnz = number of nonzero elements in 'A' matrix
#' @return nr = number or rows in 'A' matrix
#' @return nc = number of columns in 'A' matrix
#' @return obj = nc length vector of objective coefficients
#' @return ra = nonzero coefficients in 'A' matrix
#' @return ia = row indexes for non zero elements in 'A' matrix
#' @return ja = column indexes for non zero elements in 'A' matrix
#' @return SMM = Sparse matrix method: 'CMO','RMO','CRI',or 'RCI'
#' @return ZINDEX = 'zero indexing' (T) or 'one indexing' (F)
#' @return dirs = nr length vector of constaint signs ('<=','>=', or '=')
#' @return rhs = nr length vector of RHS coefficients
#' @return xlower = nc vector of lower bounds on 'x' choice variables
#' @return xupper = nc vector of upper bounds on 'x' choice variables
#' @return rlower = nr vector of Ax lower bounds i.e. rlower <= Ax
#' @return rupper = nr vector of Ax upper bounds i.e. Ax <= rupper
#' @return vartypes = nc vector of variable types: ('C','B','I') -- qDEA:rep('C',nc)
#' @return yxchngC = index of given dmu's output-input values locations in ra vector
#' (used to edit LP problem as we loop through DMU's in (X0,Y0,DX0,DY0)
#' @return dyxchngC = index of given dmu's direction values locations in ra vector
#' (used to edit LP problem as we loop through DMU's in (X0,Y0,DX0,DY0)
#' @return yxchngR = index of given dmu's output-input values locations in ra vector
#' (used to edit LP problem as we loop through DMU's in (X0,Y0,DX0,DY0)
#' @return iyxchngC = row indexes associated with yxchngC cells
#' @return idyxchngC = row indexes associated with dyxchngC cells
#' @return iyxchngR = row indexes associated with yxchngR cells
#' @return DOBJ = matrix of obj values to change by dmu with DOBJ=cbind(-Y0,X0)
#' @return DYX = matrix of (dy,dx) values to change by dmu with DYX=cbind(DY0,DX0)
#' @return qchng = index of 1/q value location in ra vector
#' (used to edit LP problem as we iterate qDEA
#' @return RTS = Returns to scale: 'CRS,'VRS,'DRS','IRS.
#' @return dmu0 = Row in (X0,Y0,DX0,DY0) to use as given dmu.
#' @return devstart = LP index for column of first qDDEA-S1 LPM "deviation" value
#' @return devend = LP index for column of last qDDEA-S1 LPM "deviation" value
#' @return tau0 = parameter used in iterative qDEA process
###############################################################################
# debug qDEAbuild
# (qout=0.10); (dmu0=1); (RTS='VRS'); (SMM='CMO'); (ZINDEX=TRUE)
###############################################################################
#' @export
###############################################################################
qDEAbuild=function(X,Y,X0,Y0,DX0,DY0,qout=0.10,dmu0=1,RTS='CRS',
unbounded=-1E3,solver='highs'){
###############################################################################
SMM='CMO'; ZINDEX=TRUE
eps = (1E-6)
DOBJ=cbind(-Y0,X0)
DYX=cbind(DY0,DX0)
(nrY=nrow(Y));(nrX=nrow(X))
(ncY=ncol(Y));(ncX=ncol(X))
if(nrY!=nrX) stop('nrY!=nrX')
(ncYX=ncY+ncX)
(nobs=nrow(Y))
(devstart=ncYX+2)
(devend=devstart+nobs-1)
y0=Y0[dmu0,];x0=X0[dmu0,];dy=DY0[dmu0,];dx=DX0[dmu0,]
q0 = qout-eps
q1 = (floor(qout*nrY)+0.95)/nrY - eps
tau0 = q1
(ncYX=ncY+ncX)
(nc=ncY+ncX+1+nrY+1)
(nr=nrY+4)
# Include obj restriction
M1=as.matrix(as.matrix(rbind(Y,dy,y0)))
ra=as.vector(M1)
ia=rep(1:(nrY+2),ncY)
# initially ja is just a count of nz elements in each column
ja=rep(nrY+2,ncY)
M2=as.matrix(rbind(-X,dx,-x0))
ra=c(ra,as.vector(M2))
ia=c(ia,rep(1:(nrX+2),ncX))
ja=c(ja,rep((nrY+2),ncX))
ra=c(ra,rep(1,nrY),-1)
ia=c(ia,1:nrY,nrY+4)
ja=c(ja,nrY+1)
ra=c(ra,rep(c(-1,1/nrY),nrY))
tmp=rep(1:nrY,each=2)
pick=seq(2,2*nrY,2)
tmp[pick]=nrY+3
ia=c(ia,tmp)
ja=c(ja,rep(2,nrY))
ra=c(ra,-1,1/tau0)
ia=c(ia,nrY+3,nrY+4)
ja=c(ja,2)
qchng=length(ra)
if(RTS!='CRS'& RTS!='crs'){
ra=c(ra,rep(-1,nrY),-1)
ia=c(ia,1:nrY,nrY+2)
ja=c(ja,nrY+1)
}
# compute CMO ja vector from counts
ja=cumsum(c(1,ja))
# change locations in constraint matrix
(yxchngC = cumsum(c(1,rep((nrY+2),(ncYX-1)))))
(dyxchngC = cumsum(c(nrY+1,rep((nrY+2),(ncYX-1)))))
(yxchngR = cumsum(c(nrY+2,rep((nrY+2),(ncYX-1)))))
###########################################################
LPsense='min'
(obj=c(-y0,x0,0,rep(0,nrY),0))
(xlower=rep(0,nc))
(xupper=rep(1E8,nc))
(dirs=c(rep('<=',nrY),'=','<=','<=','<='))
(rhs=c(rep(0,nrY),1,-unbounded,0,0))
if(RTS!='CRS'&RTS!='crs'){
nc=nc+1
obj=c(obj,1)
if(RTS=='VRS'|RTS=='vrs'){
xlower=c(xlower,-1E8)
xupper=c(xupper,1E8)
}
if(RTS=='DRS'|RTS=='drs'){
xlower=c(xlower,0)
xupper=c(xupper,1E8)
}
if(RTS=='IRS'|RTS=='irs'){
xlower=c(xlower,0)
xupper=c(xupper,1E8)
}
} # end if(RTS!='CRS'&RTS!='crs')
rlower=rhs
rupper=rhs
rlower=ifelse(dirs=='<='|dirs=='<',-1E8,rlower)
rupper=ifelse(dirs=='<='|dirs=='<',rhs,rupper)
rlower=ifelse(dirs=='>='|dirs=='>',rhs,rlower)
rupper=ifelse(dirs=='>='|dirs=='>',1E8,rupper)
rlower=ifelse(dirs=='=='|dirs=='=',rhs,rlower)
rupper=ifelse(dirs=='=='|dirs=='=',rhs+1E-8,rupper)
rlower;rupper
###############################################################################
iyxchngC=ia[yxchngC][1]
idyxchngC=ia[dyxchngC][1]
iyxchngR=ia[yxchngR][1]
###############################################################################
if(ZINDEX==TRUE) {ia=ia-1;ja=ja-1}
vartypes=rep("C",nc)
###############################################################################
###############################################################################
LP=list(LPsense=LPsense,nnz=length(ra),nr=nr,nc=nc,obj=obj,ra=ra,ia=ia,ja=ja,
SMM=SMM,ZINDEX=ZINDEX,dirs=dirs,rhs=rhs,xlower=xlower,xupper=xupper,
rlower=rlower,rupper=rupper,vartypes=vartypes,yxchngC=yxchngC,
dyxchngC=dyxchngC,yxchngR=yxchngR,iyxchngC=iyxchngC,idyxchngC=idyxchngC,
iyxchngR=iyxchngR,DOBJ=DOBJ,DYX=DYX,qchng=qchng,RTS=RTS,qout=qout,dmu0=dmu0,
devstart=devstart,devend=devend,qout=qout,tau0=tau0)
###############################################################################
LP=SM2SM(LP,SMM2='CRI',ZINDEX2=FALSE)
###############################################################################
LP
###############################################################################
} # end function qDEAbuild
###############################################################################
###############################################################################
#' my_lag: Lag a vector of values
#'
#' @param x A vector of values to be lagged
#' @param nlag Length of lag: negative value indicates an 'uplag'
#' @return = a lagged vector
###############################################################################
#' @export
###############################################################################
my_lag=function(x,nlag=1) {
nobx=length(x)
if(nlag<0) x=rev(x)
x0=rep(NA,abs(nlag))
x1=x[1:(nobx-abs(nlag))]
lagx=c(x0,x1)
if(nlag<0){
x=rev(x)
lagx=rev(lagx)
}
lagx
} # end function my_lag
###############################################################################
###############################################################################
#' lagMat: Create a matrix of lags
#'
#' @param x Vector to be 'lagged'
#' @param lags Number of lagged columns in resulting matrix
#' @param Lzero (F= do not include x[i] in row i) (T = include x[i] in row i)
#' @return = lagged matrix
###############################################################################
lagMat=function(x,lags=2,Lzero=FALSE) {
if(length(lags)==1){
if(lags[1]>0&Lzero==FALSE) lags=1:lags
if(lags[1]<0&Lzero==FALSE) lags=-1:lags
if(Lzero!=FALSE) lags=0:lags
} # end if(length(lags)==1)
nobx=length(x)
nlags=length(lags)
X=matrix(0,nobx,nlags)
for(j in 1:nlags){
X[,j]=my_lag(x,lags[j])
}
X
} # end function lagMat
###############################################################################
###############################################################################
#' my_seconds: Function to pull seconds from proc.time function
#'
#'@return = seconds from proc.time function
###############################################################################
#' @export
###############################################################################
my_seconds=function(){
time1=proc.time()
as.numeric(time1[3])
} # end function my_seconds
###############################################################################
###############################################################################
#' Write.LP
#'
#' Write a linear programming model to a csv file.
#'
#' This function converts a linear programming model to a dense
#' representation and writes it to a csv file that can be opened or solved
#' using spreadsheet based solvers.
#'
#' @param LP A linear programming object.
#' @param path Directory to write LP .csv file to. REQUIRED to use function
#' @param filename Character string specifying the csv file name.
#'
#' @return A csv file written to disk.
#'
#' @export
###############################################################################
Write.LP=function(LP,path=NULL,filename='LP.csv'){
###############################################################################
if (is.null(path)) {
stop("Please specify a directory path/folder to save the file in.")
}
###############################################################################
(fname=paste(path,'/',filename,sep=''))
###############################################################################
A = SM2A(LP)
(nc=ncol(A))
(nr=nrow(A))
if('cnames' %in% names(LP) == FALSE) LP$cnames=paste('C_',1:nc,sep='')
if('rnames' %in% names(LP) == FALSE) LP$rnames=paste('R_',1:nr,sep='')
if('rhs' %in% names(LP)==FALSE&'rupper' %in% names(LP)==FALSE){
stop('rhs and (rlower,rupper )cannot both be missing in LP object')
}
(LPsense=LP$LPsense)
obj=LP$obj
if('vartypes' %in% names(LP) == FALSE) LP$vartypes = rep('C',length(obj))
if('xlower' %in% names(LP) == FALSE) LP$xlower = rep(0,length(obj))
if('xupper' %in% names(LP) == FALSE) LP$xupper = rep(1E10,length(obj))
if('rlower' %in% names(LP) ==TRUE){if(is.na(LP$rlower[1])) LP$rlower=NULL}
if('rupper' %in% names(LP) ==TRUE){if(is.na(LP$rupper[1])) LP$rupper=NULL}
if('rlower' %in% names(LP) == FALSE) {
rlower = LP$rhs
rlower = ifelse(LP$dirs=='<='|LP$dirs=='<',-1E10,rlower)
rlower = ifelse(LP$dirs=='>='|LP$dirs=='>',LP$rhs,rlower)
LP$rlower = rlower
}
if('rupper' %in% names(LP) == FALSE) {
rupper = LP$rhs
rupper = ifelse(LP$dirs=='<='|LP$dirs=='<',LP$rhs,rupper)
rupper = ifelse(LP$dirs=='>='|LP$dirs=='>',1E10,rupper)
LP$rupper = rupper
}
solution=0*obj
if('solution' %in% names(LP) == TRUE) solution=LP$solution
cnames=LP$cnames
rnames=LP$rnames
vartypes=LP$vartypes
xlower=LP$xlower
xupper=LP$xupper
xlower[xlower<=-Inf] = -1E10
xupper[xupper>=Inf] = 1E10
dirs=LP$dirs
rhs=NA
if('rhs' %in% names(LP) == TRUE) rhs = LP$rhs
rlower=LP$rlower
rupper=LP$rupper
M=A
M=cbind(M,'|','lhsA','|')
(tmp=c(cnames,'|','LHS','|'))
if(exists('rhs')){
M=cbind(M,dirs,'|',rhs,'|')
tmp=c(tmp,'dirs','|','rhs','|')
}
if(exists('rlower')){
M=cbind(M,rlower,'|',rupper,'|')
tmp=c(tmp,'rlower','|','rupper','|')
}
M=cbind(rnames,M)
tmp=c('',tmp)
M=rbind(rep('========',ncol(M)),M)
M=rbind(tmp,M)
M=rbind(rep('========',ncol(M)),M)
M=rbind(M,rep('========',ncol(M)))
(tmp=c('OBJ',obj,'|','lhsOBJ','=',paste(LPsense)))
(tmp=c(tmp,rep('',ncol(M)-length(tmp))))
M=rbind(M,tmp)
M=rbind(M,rep('========',ncol(M)))
(tmp=c('nrows =',nr,'ncols =',nc))
M=rbind(c(tmp,rep('',ncol(M)-length(tmp))),M)
(tmp=c('xlower',xlower,'|','xlower'))
M=rbind(M,c(tmp,rep('',ncol(M)-length(tmp))))
(tmp=c('xupper',xupper,'|','xupper'))
M=rbind(M,c(tmp,rep('',ncol(M)-length(tmp))))
M=rbind(M,rep('========',ncol(M)))
(tmp=c('xvals',solution,'|','xvalues'))
M=rbind(M,c(tmp,rep('',ncol(M)-length(tmp))))
M=rbind(M,rep('========',ncol(M)))
M=rbind(rep('========',ncol(M)),M)
(tmp=c('VARTYPE',vartypes,'|','vartypes'))
M=rbind(M,c(tmp,rep('',ncol(M)-length(tmp))))
M=rbind(M,rep('========',ncol(M)))
rm(A,nr,nc)
utils::write.table(M,file=fname,row.names=FALSE,col.names=FALSE,sep=',')
###############################################################################
} # end function Write.LP
###############################################################################
#################################################################################
#' iter_delete: Function used for 'peeling' that supplements or supplants qDEA
#' 'slicing' procedure. Intended to be called from qDEA function
#'
#' @param LP0 qDEA Stage-2 LP object.
#' @param ndmus Number of dmus in reference data (X,Y).
#' @param nout0 Number of external points in current LP object's solution.
#' @param iterlim Maximal number of 'peeling' iterations.
#' @param BIGM Default Big M in RHS of qDEA stage 2 process.
#' @param unbounded LP reported as unbounded if obj<=unbounded.
#' @param eps iter_delete convergence test parameter
#' @param solver LP solver Default='highs'
#' @return obj2 = best objective level found
#' @return nout = number of external points
#' @return outlist = list of external dmus
#' @return iter = number of iterations completed
#' @return LP2 = LP object for best solution found
#' @return LP2sol = optimal LP solution
#################################################################################
# (ndmus=nrow(X))
# nout0=1
# iterlim=100
# BIGM=1E6
# unbounded=1E-3
# eps=1e-6
#################################################################################
#' @export
###############################################################################
iter_delete=function(LP0,ndmus,nout0,iterlim=max(250,nout0),BIGM=1E6,
unbounded=1E-3,eps=1e-6,solver='highs'){
###############################################################################
LP2=LP0
LP2sol=LPSOLVER(LP2,solver=solver)
(obj2=LP2sol$objval)
nobs=ndmus
lhs2=LP2sol$lhs[1:nobs]
dual2=LP2sol$dual[1:nobs]
(pick20=which(lhs2 > eps))
(nout=length(pick20))
(nleft=nout0-nout)
(pick21=which(dual2 < (-eps)))
pick22=sort(c(pick20,pick21))
(npick22=length(pick22))
go=1
iter=1
while(iter <= iterlim & npick22<=nout0){
LP2=LP0
LP2$rhs[pick22]=BIGM
LP2$rupper[pick22]=BIGM
LP2sol=LPSOLVER(LP2,solver=solver)
lhs2=LP2sol$lhs[1:nobs]
dual2=LP2sol$dual[1:nobs]
(obj2=LP2sol$objval)
(pick20=which(lhs2 > eps))
(nout=length(pick20))
(nleft=nout0-nout)
(pick21=which(dual2 < (-eps)))
(npick21=length(pick21))
(pick22=sort(c(pick20,pick21)))
(npick22=length(pick22))
(iter=iter+1)
} # end while(iter <= iterlim & npick22<=nleft){
##############################################################
iter
nout
nout0
(npick21=length(pick21))
while(iter <= iterlim & nout<nout0){
LPlist=list('NOT RUN')
objhat=0
pick20
j=1
for(j in 1:npick21){
(pick22=sort(c(pick20,pick21[j])))
LP2=LP0
LP2$rhs[pick22]=BIGM
LP2$rupper[pick22]=BIGM
LP2sol=LPSOLVER(LP2,solver=solver)
(objhat[j]=LP2sol$objval)
LPlist[[j]]=list(LP2=LP2,LP2sol=LP2sol)
} # end loop for(jon in 1:npick21)
objhat
(jpick=which(objhat==min(objhat))[1])
LP2=LPlist[[jpick]]$LP2
LP2sol=LPlist[[jpick]]$LP2sol
lhs2=LP2sol$lhs[1:nobs]
dual2=LP2sol$dual[1:nobs]
(obj2=LP2sol$objval)
(pick20=which(lhs2 > eps))
(nout=length(pick20))
(nleft=nout0-nout)
(pick21=which(dual2 < (-eps)))
(npick21=length(pick21))
(pick22=sort(c(pick20,pick21)))
(npick22=length(pick22))
(iter=iter+1)
} #end while(iter <= iterlim & nout<nout0){
list(obj2=obj2,nout=nout,outlist=pick20,iter=iter,LP2=LP2,LP2sol=LP2sol)
} # end function iter_delete
###############################################################################
###############################################################################
#require(data.table);setDTthreads(1)
###############################################################################
#' A2SM: Convert a matrix A to sparse matrix object ASM
#'
#' @param A The matrix A.
#' @param SMM Sparse matrix method with SMM='A',CRI','RCI','CMO', or 'RMO'
#' @param ZINDEX T = Use zero indexing or F = Use one indexing.
#' @param eps Value to use in non-zero test.
#' @param NA_flag Number to uses as NA flag
#'
#' @return A list object containing the matrix components:
#' @return nnz = the number of non-zero elements in ra.
#' @return nr = the number of rows in matrix A.
#' @return nc = the number of columns in matrix A.
#' @return ia = the row index.
#' @return ja = the column index.
#' @return ra = the non-zero coefficients in A
#' @return rnames = matrix row names -- may be ""
#' @return cnames = matrix column names -- may be ""
#' @return SMM = the sparse matrix type.
#' @return ZINDEX with T = Use zero indexing or F = Use one indexing.
#' @examples
#' (A = matrix(c(1,0,0,2,0,3,1,0,0,5,0,6),3,4,byrow=TRUE))
#' (ASM1 = A2SM(A,SMM='CMO',ZINDEX=TRUE))
#' (ASM2=SM2SM(ASM1,SMM2='CRI',ZINDEX2=FALSE))
#' SM2A(ASM1)
#' SM2A(ASM2)
#'
#' (A = matrix(c(1,0,0,2,0,3,1,0,0,5,0,6),3,4,byrow=TRUE)); A[2,3]=NA; A
#' (ASM1 = A2SM(A,SMM='CMO',ZINDEX=TRUE)); A; SM2A(ASM1)
#' (ASM2 = SM2SM(ASM1,SMM2='CRI',ZINDEX2=TRUE)); A ; SM2A(ASM2)
#'
#' #CMO documentation example
#' nr=5
#' nc=8
#' ra=c(3.0,5.6,1.0,2.0,1.1,1.0,-2.0,2.8,-1.0,1.0,1.0,-1.2,-1.0,1.9)
#' ia=c(0,4,0,1,1,2,0,3,0,4,2,3,0,4)
#' ja=c(0,2,4,6,8,10,11,12,14)
#' SMM='CMO'
#' ZINDEX=TRUE
#' ASM=list(nr=nr,nc=nc,ra=ra,ia=ia,ja=ja,SMM=SMM,ZINDEX=ZINDEX)
#' (A=SM2A(ASM))
###############################################################################
#' @export
###############################################################################
A2SM = function(A, SMM = 'CRI',ZINDEX=FALSE,
eps = .Machine$double.eps,NA_flag=(-1E12)){
##############################################################################
A[is.na(A)] = NA_flag
##############################################################################
(nr = nrow(A))
(nc = ncol(A))
if(is.null(rownames(A))) rownames(A)=rep('',nr)
if(is.null(colnames(A))) colnames(A)=rep('',nc)
(rnames=rownames(A))
(cnames=colnames(A))
##############################################################################
# initially uses SMM RCI
tmp=list(NULL); eps=1E-9
j=1
for(j in 1:nc){
a=as.vector(A[,j])
ipick=which(abs(a)>eps)
ra=a[ipick]
ia=ipick
ja=rep(j,length(ipick))
as.data.frame(cbind(ia,ja,ra))
tmp[[j]]=as.data.frame(cbind(ia,ja,ra))
} # end loop for(j in 1:nc)
tmp=dplyr::bind_rows(tmp)
(nnz=nrow(tmp))
ASM=list(nnz=nnz,nr=nr,nc=nc,ia=tmp$ia,ja=tmp$ja,ra=tmp$ra,
rnames=rnames,cnames=cnames,SMM='CRI',ZINDEX=FALSE)
rm(tmp)
# convert to user specified sparse matrix format
SM2SM(ASM,SMM2=SMM,ZINDEX2=ZINDEX)
} # end function A2SM
###############################################################################
###############################################################################
#' SM2A: Convert a sparse matrix object ASM into a matrix A
#'
#' ASM is a sparse matrix object containing:
#' nr = number of rows, nc = number of columns, ra = nonzero coeff,
#' ia = row indices, ja = col indices,
#' rnames = row names,cnames = column names # Note these may be missing or ""
#' SMM = space matrix method, and ZINDEX = use 0-indexing
#' @param ASM A sparse matrix object
#' @param NA_flag Number to uses as NA flag
#'
#' @return The matrix A.
#' @examples
#' (A = matrix(c(1,0,0,2,0,3,1,0,0,5,0,6),3,4,byrow=TRUE))
#' (ASM1 = A2SM(A,SMM='CMO',ZINDEX=TRUE))
#' (ASM2=SM2SM(ASM1,SMM2='CRI',ZINDEX2=FALSE))
#' SM2A(ASM1)
#' SM2A(ASM2)
#'
#' (A = matrix(c(1,0,0,2,0,3,1,0,0,5,0,6),3,4,byrow=TRUE)); A[2,3]=NA; A
#' (ASM1 = A2SM(A,SMM='CMO',ZINDEX=TRUE)); A; SM2A(ASM1)
#' (ASM2 = SM2SM(ASM1,SMM2='CRI',ZINDEX2=TRUE)); A ; SM2A(ASM2)
#'
#' #CMO documentation example
#' nr=5
#' nc=8
#' ra=c(3.0,5.6,1.0,2.0,1.1,1.0,-2.0,2.8,-1.0,1.0,1.0,-1.2,-1.0,1.9)
#' ia=c(0,4,0,1,1,2,0,3,0,4,2,3,0,4)
#' ja=c(0,2,4,6,8,10,11,12,14)
#' SMM='CMO'
#' ZINDEX=TRUE
#' ASM=list(nr=nr,nc=nc,ra=ra,ia=ia,ja=ja,SMM=SMM,ZINDEX=ZINDEX)
#' (A=SM2A(ASM))
###############################################################################
#' @export
###############################################################################
SM2A = function(ASM,NA_flag=(-1E12)){
###############################################################################
ASM$ra[is.na(ASM$ra)] = NA_flag
##
nr=ASM$nr; nc=ASM$nc
if(!exists('rnames',where=ASM)) ASM$rnames=rep('',nr)
if(!exists('cnames',where=ASM)) ASM$cnames=rep('',nc)
##
(SMM0=ASM$SMM)
if(SMM0=='A') A = ASM$A
if(SMM0!='A'){
ASM2=SM2SM(ASM,SMM2='CRI',ZINDEX2=FALSE)
nr=ASM2$nr;nc=ASM2$nc;
ia=ASM2$ia;ja=ASM2$ja;ra=ASM2$ra
(ra=as.numeric(ra))
(ia=as.integer(ia))
(ja=as.integer(ja))
(nr=as.integer(nr))
(nc=as.integer(nc))
A=matrix(0,nr,nc)
A[cbind(ia,ja)]=ra
} # end if(SMM0!='A')
##
A[A <= NA_flag] = NA
##
if(!is.null(ASM$rnames)&length(ASM$rnames)==nrow(A)) rownames(A)=ASM$rnames
if(!is.null(ASM$cnames)&length(ASM$cnames)==ncol(A)) colnames(A)=ASM$cnames
###############################################################################
A
###############################################################################
} # end function SM2A
###############################################################################
###############################################################################
#' SM2SM: Convert a sparse matrix object ASM into a different sparse matrix form.
#'
#' ASM is a sparse matrix object containing:
#' nr = number of rows, nc = number of columns, ra = nonzero coeff,
#' ia = row indices, ja = col indices, SMM = space matrix method,
#' ZINDEX=use 0-indexing
#' @param ASM A sparse matrix object
#' @param SMM2 sparse matrix method with SMM2 ='A','CRI','RCI','CMO',or'RMO'
#' @param ZINDEX2 T = Use zero indexing or F = Use one indexing.
#' @param NA_flag Number to uses as NA flag
#'
#' @return A list object containing the following sparse matrix components:
#' @return nnz = the number of non-zero elements in ra.
#' @return nr = the number of rows in matrix.
#' @return nc = the number of columns in matrix.
#' @return ia = the revised row index.
#' @return ja = the revised column index.
#' @return rnames = row names -- may be missing or ""
#' @return cnames = column names -- may be missing or ""
#' @return ra = the revised non-zero coefficients in A
#' @return SMM = the revised sparse matrix type.
#' @return ZINDEX = TRUE or F for the revised sparse form.
#' @examples
#' (A = matrix(c(1,0,0,2,0,3,1,0,0,5,0,6),3,4,byrow=TRUE))
#' (ASM1 = A2SM(A,SMM='CMO',ZINDEX=TRUE))
#' (ASM2 = SM2SM(ASM1,SMM2='CRI',ZINDEX2=FALSE))
#' SM2A(ASM1)
#' SM2A(ASM2)
#'
#' (A = matrix(c(1,0,0,2,0,3,1,0,0,5,0,6),3,4,byrow=TRUE)); A[2,3]=NA; A
#' (ASM1 = A2SM(A,SMM='CMO',ZINDEX=TRUE)); SM2A(ASM1)
#' (ASM2 = SM2SM(ASM1,SMM2='CRI',ZINDEX2=TRUE)); A ; SM2A(ASM2)
#'
#' #CMO documentation example
#' nr=5
#' nc=8
#' ra=c(3.0,5.6,1.0,2.0,1.1,1.0,-2.0,2.8,-1.0,1.0,1.0,-1.2,-1.0,1.9)
#' ia=c(0,4,0,1,1,2,0,3,0,4,2,3,0,4)
#' ja=c(0,2,4,6,8,10,11,12,14)
#' SMM='CMO'
#' ZINDEX=TRUE
#' ASM=list(nr=nr,nc=nc,ra=ra,ia=ia,ja=ja,SMM=SMM,ZINDEX=ZINDEX)
#' SM2A(ASM)
###############################################################################
#' @export
###############################################################################
SM2SM = function(ASM, SMM2='CRI', ZINDEX2=FALSE, NA_flag=(-1E12)){
###############################################################################
ASM$ra[is.na(ASM$ra)] = NA_flag
##
nnz=ASM$nnz;nr=ASM$nr; nc=ASM$nc
if(!exists('rnames',where=ASM)) ASM$rnames=rep('',nr)
if(!exists('cnames',where=ASM)) ASM$cnames=rep('',nc)
##
(rnames=ASM$rnames)
(cnames=ASM$cnames)
## original sparse form
(SMM0=ASM$SMM)
## if SMM0 indicates a matrix convert ASM to SMM='CRI' sparse form
if(SMM0=='A') {
(tmp=A2SM(ASM$A, SMM = 'CRI', ZINDEX=FALSE))
ASM=merge_lists(tmp,ASM)
} # end if(SMM0=='A')
## Pull sparse matrix elements from ASM
(nr=ASM$nr);(nc=ASM$nc);(ia=ASM$ia);(ja=ASM$ja);(ra=ASM$ra)
(SMM=ASM$SMM);(ZINDEX=ASM$ZINDEX)
## convert "zero" to "one" indexing if needed
if(ZINDEX==TRUE){
ia=ia+1
ja=ja+1
}# end if (ZINDEX == TRUE)
ia;ja
## convert CMO ja vector or RMO ia vector to complete index vectors
if(SMM=='CMO') ja=rep(1:nc,diff(ja))
if(SMM=='RMO') ia=rep(1:nr,diff(ia))
tmp=as.data.frame(cbind(ia,ja,ra))
## convert ASM to ASM2 form
if(SMM2=='CRI'|SMM2=='CMO'){
# resort to be sure in "complete and contiguous CRI form
(tmp=tmp[with(tmp,order(ja,ia)),])
# pull resorted data
ra=tmp$ra; ia=tmp$ia; ja=tmp$ja
if(SMM2=='CMO'){
(ja=c(match(1:nc,ja),nnz+1))
} # end if(SMM2=='CMO')
} # end if(SMM2=='CRI'|SMM2=='CMO')
if(SMM2=='RCI'|SMM2=='RMO'){
# resort to be sure in "complete and contiguous RCI form
(tmp=tmp[with(tmp,order(ia,ja)),])
# pull resorted data
ra=tmp$ra; ia=tmp$ia; ja=tmp$ja
if(SMM2=='RMO'){
ia=c(match(1:nr,ia),nnz+1)
} # end if(SMM2=='RMO')
} # end if(SMM2=='RCI'|SMM2=='RMO)
(ra=as.numeric(ra))
(ia=as.integer(ia))
(ja=as.integer(ja))
(nr=as.integer(nr))
(nc=as.integer(nc))
if(ZINDEX2 == TRUE) {ia = ia-1; ja = ja-1}
##
##
ASM2=ASM # inherits all objects in ASM
# put revised data into ASM2
ASM2$nnz=nnz;ASM2$nr=nr;ASM2$nc=nc;ASM2$ia=ia;ASM2$ja=ja;ASM2$ra=ra
ASM2$SMM=SMM2;ASM2$ZINDEX=ZINDEX2
if(SMM2=='A'){
ASM2=ASM # inherits all objects in original ASM
ASM2$A = SM2A(ASM)
ASM2$SMM = SMM2
}
#(ASM2 = list(nnz=nnz,nr = nr,nc = nc,ia = ia2,ja = ja2,ra = ra2,
# rnames=rnames,cnames=cnames,SMM = SMM2,ZINDEX = ZINDEX2))
###############################################################################
ASM2
###############################################################################
} # end function SM2SM
###############################################################################
# ###############################################################################
# # Examples:
# ###############################################################################
# # create matrix with missing data
# (A = matrix(c(1,0,0,2,0,3,1,0,0,5,0,6),3,4,byrow=TRUE)); A[2,3]=NA; A
# # no names for in matrix A
# c(rownames(A),colnames(A))
# # add names to A matrix
# rownames(A)=paste('r',1:nrow(A),sep='');colnames(A)=paste('c',1:ncol(A),sep='')
# A
# # Convert to sparse matrix with NA's replaced by NA_flag
# ASM1=A2SM(A,NA_flag=-1000)
# # Note below: SM2A puts " NA's" back in matrix if original matrix had them
# # or if user indicates levels in matrix that designate missing data
# SM2A(ASM1);SM2A(ASM1,NA_flag=-1E4);SM2A(ASM1,NA_flag=-1E3)
# ##############################################################################
#
# ##############################################################################
# A=t(matrix(1:15,5,3))
# A[cbind(c(3,1,2,3,1,3),c(1,2,3,3,4,5))]=0
# A
# #
# # (A=matrix(c(1,0,2,0,0,3,0,4,0,5,6,0),3,4,byrow=TRUE))
# # (A=matrix(c(1,0,0,2,0,3,4,0,0,5,0,6),3,4,byrow=TRUE))
# # (A=matrix(c(1,0,0,2,0,3,0,0,0,5,0,6),3,4,byrow=TRUE))
# # (A=matrix(c(2,3,1,4,1,2,3,4,2),3,3,byrow=TRUE))
# #
# A2SM(A,SMM='CMO',ZINDEX=TRUE)
#
# SMMlist=c('RCI','CRI','RMO','CMO')
# ZINDEXlist=c(FALSE,TRUE)
#
# SMM1 = 'RCI' # 'RCI', 'CRI', 'RMO', 'CMO'
# ZINDEX1=FALSE # T or F
# SMM2='CRI' # 'RCI', 'CRI', 'RMO', 'CMO'
# ZINDEX2=FALSE # T or F
#
# for(SMM1 in SMMlist){
# for(ZINDEX1 in ZINDEXlist) {
# for(SMM2 in SMMlist) {
# for(ZINDEX2 in ZINDEXlist) {
# print(paste(SMM1,ZINDEX1,SMM2,ZINDEX2))
# (ASM=A2SM(A,SMM=SMM1,ZINDEX=ZINDEX1))
# (ASM1=SM2SM(ASM,SMM2=SMM2,ZINDEX2=ZINDEX2))
# (ASM2=A2SM(A,SMM=SMM2,ZINDEX=ZINDEX2))
# print(summary(unlist(ASM1)==unlist(ASM2)))
# } # end for(ZINDEX2 in ZINDEXlist)
# } # end for(SMM2 in SMMlist)
# } # end for(ZINDEX1 in ZINDEXlist)
# } # end for(SMM1 in SMMlist)
#
#
# summary(SM2A(ASM1)==SM2A(ASM2))
# #############################################################
# (A=matrix(c(1,0,2,0,0,3,0,4,0,5,6,0),3,4,byrow=TRUE))
#
# (ASM=A2SM(A,SMM='A'))
# ASM2=SM2SM(ASM,SMM2=SMM2,ZINDEX2=ZINDEX2)
# (A2=SM2A(ASM2))
# summary(as.vector(A)==as.vector(A2))
#
# for(SMM2 in SMMlist) {
# for(ZINDEX2 in ZINDEXlist) {
# print(paste(SMM2,ZINDEX2))
# (ASM=A2SM(A=A,SMM='A'))
# (ASM2=SM2SM(ASM,SMM2=SMM2,ZINDEX2=ZINDEX2))
# (A2=SM2A(ASM2))
# print(summary(as.vector(A)==as.vector(A2)))
# } # end for(ZINDEX2 in ZINDEXlist)
# } # end for(SMM2 in SMMlist)
# ###########################################################
# (A=matrix(c(1,0,2,0,0,3,0,4,0,5,6,0),3,4,byrow=TRUE))
# SMM2='CRI' # 'RCI', 'CRI', 'RMO', 'CMO'
# ZINDEX2=FALSE # T or F
#
# (ASM=A2SM(A,SMM='A'))
# ASM2=SM2SM(ASM,SMM2=SMM2,ZINDEX2=ZINDEX2)
# (A2=SM2A(ASM2))
# summary(as.vector(A)==as.vector(A2))
# (A3=SM2A(ASM))
# summary(as.vector(A)==as.vector(A3))
# ##########################################################
# # CMO documentation example
# nr=5
# nc=8
# ra=c(3.0,5.6,1.0,2.0,1.1,1.0,-2.0,2.8,-1.0,1.0,1.0,-1.2,-1.0,1.9)
# ia=c(0,4,0,1,1,2,0,3,0,4,2,3,0,4)
# ja=c(0,2,4,6,8,10,11,12,14)
# SMM='CMO'
# ZINDEX=TRUE
# ASM=list(nr=nr,nc=nc,ra=ra,ia=ia,ja=ja,SMM=SMM,ZINDEX=ZINDEX)
# (A=SM2A(ASM))
# #############################################################################
###############################################################################
#' cbindSM: "column bind" two sparse matrices
#'
#' @param SM1 First sparse matrix object
#' @param SM2 Second sparse matrix object
#' @param SMM Sparse matrix method with 'CRI','RCI','CMO',or'RMO'
#' @param ZINDEX T = Use zero indexing or F = Use one indexing.
#'
#' @return A 'column bound" sparse matrix object with components"
#' @return nnz = the number of non-zero elements in ra.
#' @return nr = the number of rows in matrix A.
#' @return nc = the number of columns in matrix A.
#' @return ia = the row index.
#' @return ja = the column index.
#' @return ra = the non-zero coefficients in A
#' @return rnames = matrix row names -- may be ''
#' @return cnames = matrix column names -- may be ''
#' @return SMM = the sparse matrix type.
#' @return ZINDEX with T = Use zero indexing or F = Use one indexing.
###############################################################################
#' @export
###############################################################################
cbindSM=function(SM1,SM2,SMM='CRI',ZINDEX=FALSE){
nr1=SM1$nr; nc1=SM1$nc; nr2=SM2$nr; nc2=SM2$nc
if(nr1!=nr2) stop('STOP: row numbers not compatible')
tmp1=SM2SM(SM1,SMM2='CRI',ZINDEX2=FALSE)
tmp2=SM2SM(SM2,SMM2='CRI',ZINDEX2=FALSE)
nr=nr1; nc = nc1+nc2
ia=c(tmp1$ia,tmp2$ia)
ja=c(tmp1$ja,tmp2$ja+nc1)
ra=c(tmp1$ra,tmp2$ra)
tmp1$nnz=length(ra);tmp1$nr=nr;tmp1$nc=nc
tmp1$ia=ia;tmp1$ja=ja;tmp1$ra=ra;tmp1$SMM='CRI';tmp1$ZINDEX=FALSE
if(length(tmp1$cnames)>1) {
if(length(tmp2$cnames)>1) tmp1$cnames=c(tmp1$cnames,tmp2$cnames)
if(length(tmp2$cnames)<=1) tmp1$cnames=c(tmp1$cnames,rep('',nc2))
} # end
SM2SM(tmp1,SMM2=SMM,ZINDEX2=ZINDEX)
} # end function cbindSPM
###############################################################################
###############################################################################
#' rbindSM: "row bind" two sparse matrices
#'
#' @param SM1 First sparse matrix object
#' @param SM2 Second sparse matrix object
#' @param SMM Sparse matrix method with 'CRI','RCI','CMO',or'RMO'
#' @param ZINDEX T = Use zero indexing or F = Use one indexing.
#'
#' @return A 'row bound" sparse matrix object with components"
#' @return nnz = the number of non-zero elements in ra.
#' @return nr = the number of rows in matrix A.
#' @return nc = the number of columns in matrix A.
#' @return ia = the row index.
#' @return ja = the column index.
#' @return ra = the non-zero coefficients in A
#' @return rnames = matrix row names -- may be ''
#' @return cnames = matrix column names -- may be ''
#' @return SMM = the sparse matrix type.
#' @return ZINDEX with T = Use zero indexing or F = Use one indexing.
###############################################################################
#' @export
###############################################################################
rbindSM=function(SM1,SM2,SMM='CRI',ZINDEX=FALSE){
nr1=SM1$nr; nc1=SM1$nc; nr2=SM2$nr; nc2=SM2$nc
if(nc1!=nc2) stop('STOP: column number not compatible')
tmp1=SM2SM(SM1,SMM2='RCI',ZINDEX2=FALSE)
tmp2=SM2SM(SM2,SMM2='RCI',ZINDEX2=FALSE)
nr=nr1+nr2; nc = nc1
ia=c(tmp1$ia,tmp2$ia+nr1)
ja=c(tmp1$ja,tmp2$ja)
ra=c(tmp1$ra,tmp2$ra)
tmp1$nnz=length(ra);tmp1$nr=nr;tmp1$nc=nc
tmp1$ia=ia;tmp1$ja=ja;tmp1$ra=ra;tmp1$SMM='RCI';tmp1$ZINDEX=FALSE
if(length(tmp1$rnames)>1) {
if(length(tmp2$rnames)>1) tmp1$rnames=c(tmp1$rnames,tmp2$rnames)
if(length(tmp2$rnames)<=1) tmp1$rnames=c(tmp1$rnames,rep('',nr2))
} # end
SM2SM(tmp1,SMM2=SMM,ZINDEX2=ZINDEX)
} # end function rbindSPM
###############################################################################
#
#
#
################################################################################
# # Examples of cbindSM and rbindSM
# set.seed(1001)
# nr=7;nc=7
# (A=matrix(round(runif(nr*nc),3),nr,nc))
# A[A<0.5]=0
# colnames(A)=paste('c',1:nc,sep='')
# rownames(A)=paste('r',1:nr,sep='')
# A
#
# (A11=A[1:5,1:5]); (A12=A[1:5,6:7])
# (A21=A[6:7,1:5]); (A22=A[6:7,6:7])
#
#
# (SM11=A2SM(A11))
# (SM12=A2SM(A12))
# (SM21=A2SM(A21))
# (SM22=A2SM(A22))
#
# cbind(A11,A12)
# SM2A(cbindSM(SM11,SM12))
#
#
# rbind(A11,A21)
# SM2A(rbindSM(SM11,SM21))
#
# A
# SM2A(cbindSM(rbindSM(SM11,SM21),rbindSM(SM12,SM22)))
###############################################################################
# A2SM(A,SMM='CMO')
# cbindSM(rbindSM(SM11,SM21),rbindSM(SM12,SM22),SMM='CMO')
###############################################################################
###############################################################################
#' merge_lists: Merges two list objects.
#' This function appends any objects in list2 and not in list 1
#' to list1 with priority given to list 1 components.
#'
#' @param list1 A list object
#' @param list2 A list object
#' @return list3 A merged object
#' @examples
#' L1=list(a=1:3,b=4:6,c=7:9,e=10:12)
#' L2=list(b=13:15,d=16:18,e=19:21,f=22:24)
#' (L3=merge_lists(L1,L2))
#' (L4=merge_lists(L2,L1))
###############################################################################
#' @export
###############################################################################
merge_lists = function(list1, list2){
#####################################################################
# This function appends any objects in list2 and not in list 1
# to list1 with priority given to list 1 components.
#####################################################################
(names1 = names(list1))
(names2 = names(list2))
(list3 = list1)
(n1=length(list1))
(pick2 = which(is.na(match(names2, names1))))
if (length(pick2) > 0) {
for(j in 1:length(pick2)) list3 = c(list3,list2[pick2[j]])
}
sort_list(list3)
} # end function merge_lists
###############################################################################
###############################################################################
#' sort_list: Sorts objects in list by object name
#'
#' @param list1 A list object
#' @return list2 A resorted list object
###############################################################################
#' @export
###############################################################################
sort_list=function(list1){
(n=length(list1))
df1=data.frame(Lnames=names(list1))
df1$cell=1:n; df1$cell=as.numeric(df1$cell)
df1=df1[order(df1$Lnames),]
list2=list1[df1$cell[1]]
for(j in 2:n) list2=c(list2,list1[df1$cell[j]])
list2
} # end function sort_list
################################################################################
###############################################################################
#' LPSOLVER Function to call specified solver to sparse LP problem
#'
#' @param LP An LP object
#' @param solver LP solver Default='highs' Note: Versions for 'clp','highs','glpk' or 'rglpk' are available from author
#' @return A list containing following components:
#' @return status = status of solver
#' @return objval = objective value
#' @return solution = lp primal solution
#' @return dual = lp dual solution
#' @return rcost = lp reduced cost
#' @return lhs = lp lhs associated with primal solution
#' @return TIME = total seconds to execute LP_* function
#' @return SMtime = seconds to transform form of sparse LP (if needed)
#' @return proctime = seconds required by specified solver to solve LP
#' @return LP = LP object
###############################################################################
#' @export
###############################################################################
LPSOLVER=function(LP,solver='highs'){
LP_highs(LP)
} # end function LPSOLVER
###############################################################################
###############################################################################
#' LP_highs Function to solve sparse LP problem using highs package
#'
#' @param LP An LP object
#' @return A list containing following components:
#' @return status = status of highs solver
#' @return objval = objective value
#' @return solution = lp primal solution
#' @return dual = lp dual solution
#' @return rcost = lp reduced cost
#' @return lhs = lp lhs associated with primal solution
#' @return TIME = total seconds to execute LP_highs function
#' @return SMtime = seconds to transform form of sparse LP (if needed)
#' @return proctime = seconds required by highs to solve LP
#' @return LP = LP object (transformed to SMM='CRI', ZINDEX=FALSE if needed)
###############################################################################
#' @export
###############################################################################
LP_highs=function(LP){
###############################################################################
#require(highs);require(Matrix)
###############################################################################
if(!requireNamespace("highs", quietly=TRUE)|!requireNamespace("Matrix", quietly=TRUE)) {
stop('The packages highs and Matrix must be installed to use the option solver = highs')
}
###############################################################################
# Note: This code is set up to allow
# the user to measure the time used in
# different procedures. The following
#lines pull initial times from the system
###############################################################################
time000=my_seconds(); time00=my_seconds()
###############################################################################
#Convert LP matrix format (if necessary)
###############################################################################
# print(paste('SMM',LP$SMM))
###############################################################################
if(LP$SMM!='CRI'|LP$ZINDEX!=FALSE) LP=SM2SM(LP,SMM2='CRI',ZINDEX2=FALSE)
# Compute sparse matrix conversion time used
SMtime=my_seconds()-time00
#######################################
if('vartypes' %in% names(LP) == FALSE) LP$vartypes = rep('C',length(LP$obj))
#######################################
# Construct lower and upper bounds on x
# If xbounds are not specified in LP
# we assume 0 <= x <= Inf
#######################################
if('xlower' %in% names(LP) == FALSE) LP$xlower = rep(0,length(LP$obj))
if('xupper' %in% names(LP) == FALSE) LP$xupper = rep(Inf,length(LP$obj))
#######################################
# Construct lower and upper bounds on Ax
# If rbounds are not specified in LP
# we construct them from the LP list's
# 'rhs' and 'dirs' vectors
#######################################
if('rlower' %in% names(LP) == FALSE) {
rlower = LP$rhs
rlower = ifelse(LP$dirs=='<='|LP$dirs=='<',-Inf,rlower)
rlower = ifelse(LP$dirs=='>='|LP$dirs=='>',LP$rhs,rlower)
LP$rlower = rlower
}
if('rupper' %in% names(LP) == FALSE) {
rupper = LP$rhs
rupper = ifelse(LP$dirs=='<='|LP$dirs=='<',LP$rhs,rupper)
rupper = ifelse(LP$dirs=='>='|LP$dirs=='>',Inf,rupper)
LP$rupper = rupper
}
(LPsense = ifelse(LP$LPsense!='max'|LP$LPsense!='MAX',FALSE,TRUE))
# # highs blows up if abs of any sparse matrix ra values are <1e-9
# LP$ra[abs(LP$ra)<2e-9] = 2e-9
ASM = Matrix::sparseMatrix(LP$ia,LP$ja,x=LP$ra)
time0=my_seconds()
tmp = highs::highs_solve(L = LP$obj, lower = LP$xlower, upper = LP$xupper,
A = ASM, lhs = LP$rlower, rhs = LP$rupper,types=LP$vartypes,
maximum=LPsense,control=list(threads=1))
proctime=my_seconds()-time0
# retrieve the results
(status=tmp$status)
(objval=tmp$objective_value)
(solution=tmp$primal_solution)
rcost=tmp$solver_msg$col_dual
lhs=tmp$solver_msg$row_value
dual=tmp$solver_msg$row_dual
TIME=my_seconds()-time000
###############################################################################
list(status=status,objval=objval,solution=solution,dual=dual,rcost=rcost,
lhs=lhs,TIME=TIME,SMtime=SMtime,proctime=proctime,LP=LP)
###############################################################################
} # end function LP_highs
###############################################################################
Try the qDEA package in your browser
Any scripts or data that you put into this service are public.
qDEA documentation built on April 13, 2026, 5:07 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions LP_highs LPSOLVER sort_list merge_lists rbindSM cbindSM SM2SM SM2A A2SM iter_delete Write.LP my_seconds lagMat my_lag qDEAbuild DEAbuild nCm_mpick qDEA_mlist qDEA_solve qDEA
Documented in A2SM cbindSM DEAbuild iter_delete lagMat LP_highs LPSOLVER merge_lists my_lag my_seconds nCm_mpick qDEA qDEAbuild qDEA_mlist qDEA_solve rbindSM SM2A SM2SM sort_list Write.LP
###############################################################################
#' qDEA: Calling function for set of DEA and qDEA processes
#'
#' Note: Optional arguments in function call have default values of 'NULL'
#'
#' @param X Reference dmu's = ndmu x number of inputs input matrix.
#' @param Y Reference dmu's = ndmu x number of outputs output matrix.
#' @param qout Maximal proportion of dmu's allowed external to DEA hull.
#' Default = 1/ndmu
#' @param qoutS Proportion of external points to identify using qDEA slicing
#' (Atwood and Shaik 2020 EJOR)with qoutS <= qout. Default = qout
#' @param X0 Inputs for set of ndmu0 dmu's to be processed. Default = X
#' @param Y0 Outputs for set of ndmu0 dmu's to be processed. Default = Y
#' @param DX0 Input directions for ndmu0 dmu's in X0 and Y0.
#' (Must be provided if orient = 'ddea'). Default = NULL
#' @param DY0 Output directions for ndmu0 dmu's in X0 and Y0.
#' (Must be provided if orient = 'ddea') Default = NULL
#' @param orient Model orientation ('in','out','inout','oneone','ddea').
#' (!!If orient='ddea',DX0,and DY0 must be provided)
#' Default='out'
#' @param RTS Returns to scale.('CRS,'VRS,'DRS','IRS) Default 'CRS'.
#' @param dmulist Index vector of dmus in (X0,Y0) to process. NULL = process all.
#' @param nqiter Maximal number of qDEA iterations. Default = 1.
#' @param nboot Number of bootstrap replications. Default = 0.
#' @param transform Transform DDEA distance to traditional efficiency metrics
#' (input or output orientations only) Default=TRUE
#' @param mcells Number of subsample sizes for bootstapping. Default = 5
#' @param mlist Optional list of user chosen subsample sizes. Default = NULL
#' @param seedval Seed value for random number generator - used in bootstapping.
#' Default = 1001.
#' @param qtol q search tolerance with iterative qDEA. Default = 1E-6
#' @param BIGM Default Big M in RHS of qDEA stage 2 process. Default = 1E9
#' @param eps Search tolerance in qDEA improvement tests. Default = 1E-6
#' @param skipzprob Skip qDEA if qout=0. Default=TRUE.
#' @param replaceA1 Put dmu0's data in first row of reference sets. Default=FALSE
#' @param baseqDEA Use basic qDEA model from EJOR article. Default=FALSE
#' @param unbounded qDEA reported as unbounded if obj<=unbounded. Default = (-1E3)
#' @param obj2test Convergance tol for objective in iterative qDEA.Default = 1E-4
#' @param replaceM Subsample with replacement in subsample bootstrap. Default=FALSE
#' @param alpha Alpha level for Simar-Wilson(2011)subsample size selection
#' procedure. Default = 0.05
#' @param betaq qDEA convergence rate. 0.5 => "root n" convergence.
#' @param siglist Vector of user's desired confidence interval widths
#' Default = c(0.10, 0.05, 0.01)
#' @param CILag CILag level for Simar-Wilson (2011) subsample size
#' selection procedure. Default = 1
#' @param printlog Progress of DDEA dmu's solved when (X0,Y0) >1 dmus.
#' @param prntmod Print progress every prntmod dmus. Default=100.
#' @param printtxt Additional text to print with progress printlog.
#' @param getproject Compute projected values for dmu's in dmulist. Default=FALSE
#' @param getbootpeers Pull and store peers for each bootstrapped solution.
#' Default=FALSE
#' @param solver LP solver Default='highs'
#'
#' @return '###################################################################'
#' @return A list of individual items and lists of addiitional/optional output
#'
#' @return effvals = vector of DDEA distances
#' (efficiencies if transform=TRUE and input or output orientation).
#' @return effvalsq = vector of qDDEA distances
#' (efficiencies if transform=TRUE and input or output orientation).
#' @return distvals = vector of DDEA distances
#' @return distvalsq = vector of qDEA distances
#' @return distMAT = ndmu0 by 3 matrix with DDEA,qDDEA-S1,qDDEA-S2 distances.
#' '###########################################################################'
#' @return INPUT_DATA = A list containing:
#'
#' @return X = reference dmu's = ndmu x number of inputs input matrix.
#' @return Y = reference dmu's = ndmu x number of outputs output matrix.
#' @return X0 = inputs for set of ndmu0 dmu's processed.
#' @return Y0 = outputs for set of ndmu0 dmu's processed.
#' @return DX0 = input directions for ndmu0 dmu's processed.
#' @return DY0 = output directions for ndmu0 dmu's processed.
#' @return qout = maximal proportion of dmu's allowed external to DEA hull.
#' @return qoutS = proportion of external points to identify using qDEA slicing.
#' @return RTS = returns to scale.
#' @return orient = model orientation. ('ddea','in','out','inout','oneone')
#' @return baseqDEA = use basic qDEA model from EJOR article
#' @return dmulist0 = DMU0 index in originally inputs X0,Y0,DX0,and DY0
#' @return '###################################################################'
#' @return BOOT_DATA = A list containing:
#'
#' @return effvals.bc = vector of DDEA bias corrected distance or efficiency metrics
#' (efficiencies if transform=TRUE and input or output orientation)
#' @return effvalsq.bc = vector of qDEA bias corrected distance or efficiency metrics
#' (efficiencies if transform=TRUE and input or output orientation)
#' @return distvals.bc = vector of bias corrected DDEA distances
#' @return distvalsq.bc = vector of bias corrected qDEA distances
#' @return mcells = number of subsample sizes for bootstapping
#' @return mlist = list of subsample sizes used in bootstraps
#' @return BOOTdmus = matrix containing indexes of reference dmus chosen for each bootstrap
#' @return BOOT = nboot by mcells by ndmu0 array of bootstrapped DDEA distances
#' (efficiencies if transform=TRUE and input or output orientation)
#' @return BOOTS = nboot by mcells by ndmu0 array of DDEA bootstrapped s-statistics
#' (statn = (m(n)/n)^beta (statm - statn)
#' @return BOOTs = nboot by ndmu0 matrix of DDEA bootstrapped s-statistics
#' (at sample size chosen by Simar and Wilson (2011) suggested process)
#' @return BOOTq = nboot by mcells by ndmu0 array of bootstrapped qDEA distances
#' (efficiencies if transform=TRUE and input or output orientation)
#' @return BOOTSq = nboot by mcells by ndmu0 array of qDEA bootstrapped s-statistics
#' (statn = (m(n)/n)^beta (statm - statn)
#' @return BOOTsq = nboot by ndmu0 matrix of qDEA bootstrapped s-statistics
#' (at sample size chosen by Simar and Wilson (2011) suggested process
#' @return mpick = length(ndmu0) index vector of DDEA bootstrap sample size chosen from mlist
#' using Simar and Wilson (2011) suggested sample size selection process
#' @return mpickq = length(ndmu0) index vector of qDDEA bootstrap sample size chosen from mlist
#' using Simar and Wilson (2011) suggested sample size selection process
#' @return beta = DDEA convergence rate indicated Simar-Wilson (2011)
#' @return betaq = qDEA convergence rate indicated Atwood and Shaik(2020) = 0.5
#' @return siglist = vector of user's desired confidence interval widths
#' @return CI = ndmu0 by 2 by length(siglist) array of DDEA confidence intervals
#' (estimated using quantiles on bias corrcted "s" statistics)
#' @return CIq = ndmu0 by 2 by length(siglist) array of qDEA confidence intervals
#' (estimated using quantiles on bias corrcted "s" statistics)
#' @return CIq_norm = ndmu0 by 2 by length(siglist) array of qDEA confidence intervals
#' (estimated using sample mean and standard error of bootstrapped
#' "s" statistics and assuming normality) (Atwood and Shaik 2020)
#' @return '###################################################################'
#' @return PEER_DATA = A list containing:
#'
#' @return PEERS = dataframe of DDEA peers and projection weights for each dmu0
#' @return PEERSq = dataframe of qDDEA-S2 peers and projection weights each dmu0
#' @return DOUTq = data frame indicating external dmus for each dmu's qDEA solution
#' @return BPEERS = ndmu0 by mcells by (nIN+nOUT) by 2(dmu-z,zweight) array
#' of DDEA subsampled peers (dmu-z) and projection weights (z)
#' if(nboot>0&getbootpeers==TRUE)
#' @return BPEERSq = ndmu0 by mcells by (nIN+nOUT) by 2(dmu-z,zweight) array
#' of qDEA subsampled peers (dmu-z) and projection weights (z)
#' if(nboot>0&getbootpeers==TRUE)
#' @return BPEERS1 = data.frame of DDEA subsampled peers and projection weights (z)
#' if(nboot>0&getbootpeers==TRUE) for chosen subsample size
#' data.frame: nb = bootstrap rep, dmu0 = given dmu,
#' dmuz = reference dmu, z = reference dmu projection weight
#' @return BPEERS1q = data.frame of qDEA subsampled peers and projection weights (z)
#' if(nboot>0&getbootpeers==TRUE) for chosen subsample size
#' data.frame: nb = bootstrap rep, dmu0 = given dmu,
#' dmuz = reference dmu, z = reference dmu projection weight
#' @return '###################################################################'
#' @return PROJ_DATA = Projected Values. A list containing:
#'
#' @return X0HAT = DDEA projected input levels (if getproject=TRUE)
#' @return Y0HAT = DDEA projected output levels (if getproject=TRUE)
#' @return X0HAT.bc = bias corrected DDEA projected input levels
#' (if nboot>0 and getproject=TRUE)
# @return Y0HAT.bc = bias corrected DDEA projected output levels
#' (if nboot>0 and getproject=TRUE) )
#' @return X0HATq = qDEA projected input levels (if getproject=TRUE)
#' @return Y0HATq = qDEA projected output levels (if getproject=TRUE)
#' @return X0HATq.bc = bias corrected qDEA projected input levels
#' (if nboot>0 and getproject=TRUE)
#' @return Y0HATq.bc = bias corrected qDEA projected output levels
#' (if nboot>0 and getproject=TRUE) )
#' @return '###################################################################'
#' @return LP_DATA = LP Models, data, and results. A list containing:
#'
#' @return status = ndmu0 by 3 matrix with LP status of DDEA,qDDEA-S1,qDDEA-S2
#' @return qhat = proportion of dmu's external to hull in qDEA solution
#' (NA indicates qDEA for given DMU was unbounded)
#' @return qiter = number of qDEA iterations completed.
#' @return PSOL = ndmu0 by ? matrix with DDEA LP solutions for each dmu0
#' @return PSOLq1 = ndmu0 by ? matrix with qDDEA-S1 LP solutions for each dmu0
#' @return PSOLq2 = ndmu0 by ? matrix with qDDEA-S2 Lp solutions for each dmu0
#' @return RCOST = ndmu0 by ? matrix of LP reduced costs for DDEA solutions
#' @return RCOSTq1 = ndmu0 by ? matrix of LP reduced costs for qDDEA-S1 solutions
#' @return RCOSTq2 = ndmu0 by ? matrix of LP reduced costs for qDDEA-S2 solutions
#' @return LPModels = list of LP0, LP1, and LP2 LP objects from qDEA_solve function
#' @examples
#' # Examples from CST(2006): Cooper, W., Seiford, L., and Tone, K., 2006.
#' # Introduction to Data Envelopment Analysis and Its Uses. Springer, New York.
#'
#' # CST One Input - One Output Example - Table 1.1
#' data(CST11)
#' X = as.matrix(CST11$EMPLOYEES)
#' Y = as.matrix(CST11$SALES_EJOR)
#' qout = 1/nrow(X)
#' sol = qDEA(X, Y, RTS = 'crs', orient = 'in', qout = qout)
#' # DEA efficiency scores
#' sol$effvals
#' # qDEA input efficiency scores allowing one external point
#' sol$effvalsq
#'
#' # CST Two Input - One Output Example - Table 1.3
#' data(CST21)
#' X = as.matrix(cbind(CST21$EMPLOYEES, CST21$FLOOR_AREA))
#' Y = as.matrix(CST21$SALES)
#' qout = 1/nrow(X)
#' sol = qDEA(X, Y, RTS = 'crs', orient = 'in', qout = qout)
#' sol$effvals
#' sol$effvalsq
#'
#' # CST One Input - Two Output Example - Table 1.4
#' data(CST12)
#' X = as.matrix(CST12$EMPLOYEES)
#' Y = as.matrix(cbind(CST12$CUSTOMERS, CST12$SALES))
#' qout = 1/nrow(X)
#' sol = qDEA(X, Y, RTS = 'crs', orient = 'out', qout = qout)
#' sol$effvals
#' sol$effvalsq
#'
#' # CST Two Input - Two Output Example - Table 1.5
#' data(CST22)
#' X = as.matrix(cbind(CST22$DOCTORS, CST22$NURSES))
#' Y = as.matrix(cbind(CST22$OUT_PATIENTS, CST22$IN_PATIENTS))
#' qout = 1/nrow(X)
#' sol = qDEA(X, Y, RTS = 'crs', orient = 'in', qout = qout)
#' # DEA efficiency scores - see table 1.6 CCR estimates
#' round(sol$effvals, 2)
#' # qDEA efficiency scores allowing one external point
#' round(sol$effvalsq, 2)
#'
#' \donttest{
#' # Atwood-Shaik EJOR qDEA Examples (longer running examples)
#' data(CST11)
#' X = as.matrix(CST11$EMPLOYEES)
#' Y = as.matrix(CST11$SALES_EJOR)
#' ###############################################################################
#' # EJOR Efficiency Results Table 1 Input Orientation
#' tmpC1=qDEA(X,Y,RTS='crs',orient='in',qout=1/8,getproject=TRUE)
#' tmpC2=qDEA(X,Y,RTS='crs',orient='in',qout=2/8,getproject=TRUE)
#'
#' # Table 1 Input Orientation
#' round(cbind(tmpC1$distvals,tmpC1$PROJ_DATA$X0HAT,tmpC1$PROJ_DATA$Y0HAT,
#' tmpC1$distvalsq,tmpC1$PROJ_DATA$X0HATq,tmpC1$PROJ_DATA$Y0HATq,
#' tmpC2$distvalsq,tmpC2$PROJ_DATA$X0HATq,tmpC2$PROJ_DATA$Y0HATq),3)
#' ###############################################################################
#' # EJOR Efficiency Results Table 1 Output Orientation
#' tmpC3=qDEA(X,Y,RTS='crs',orient='out',qout=1/8,getproject=TRUE)
#' tmpC4=qDEA(X,Y,RTS='crs',orient='out',qout=2/8,getproject=TRUE)
#'
#' # Table 1 Output Orientation
#' round(cbind(tmpC3$distvals,tmpC3$PROJ_DATA$X0HAT,tmpC3$PROJ_DATA$Y0HAT,
#' tmpC3$distvalsq,tmpC3$PROJ_DATA$X0HATq,tmpC3$PROJ_DATA$Y0HATq,
#' tmpC4$distvalsq,tmpC4$PROJ_DATA$X0HATq,tmpC4$PROJ_DATA$Y0HATq),3)
#' ################################################################################
#' # EJOR Efficiency Results Table 1 one-one Orientation
#' tmpC5=qDEA(X,Y,RTS='crs',orient='oneone',qout=1/8,getproject=TRUE)
#' tmpC6=qDEA(X,Y,RTS='crs',orient='oneone',qout=2/8,getproject=TRUE)
#'
#' # Table 1 one-one Orientation
#' round(cbind(tmpC5$distvals,tmpC5$PROJ_DATA$X0HAT,tmpC5$PROJ_DATA$Y0HAT,
#' tmpC5$distvalsq,tmpC5$PROJ_DATA$X0HATq,tmpC5$PROJ_DATA$Y0HATq,
#' tmpC6$distvalsq,tmpC6$PROJ_DATA$X0HATq,tmpC6$PROJ_DATA$Y0HATq),3)
#' ############################################################################
#' }
###############################################################################
#
###############################################################################
# debugging function qDEA
###############################################################################
###############################################################################
## default settings
# X=X;Y=Y;X0=NULL;Y0=NULL;DX0=NULL;DY0=NULL
###############################################################################
# qout=1/nrow(X);qoutS=qout
###############################################################################
# orient='out';RTS='CRS';dmulist=NULL;nqiter=1;nboot=0;transform=TRUE;mcells=5
# mlist=NULL;seedval=1001;qtol=1E-6;BIGM=1E9;eps=1E-6;skipzprob=TRUE;replaceA1=FALSE
# baseqDEA=FALSE;unbounded=(-1E3);obj2test=1E-4;replaceM=FALSE;alpha=0.05;betaq=0.5
# siglist=c(0.10,0.05,0.01);CILag=1;printlog=TRUE;prntmod=100;printtxt=''
# getproject=FALSE;getbootpeers=FALSE
###############################################################################
# MRE plots
# tmpC=qDEA(X,Y,qout=qout,RTS=RTS,nqiter=nqiter,orient='ddea',
# X0=X0,Y0=Y0,DX0=DX0,DY0=DY0,baseqDEA=TRUE,getproject=TRUE)
#
# RTS='VRS';nboot=0; mcells=1
###############################################################################
# John-Illex-Data
#tmp2=qDEA(X,Y,X0=X,Y0=Y,orient='ddea',DX0=DX0,DY0=DY0,
# RTS='VRS',qout=qout, nboot=nboot,mcells=1,
# getbootpeers=TRUE, replaceA1=TRUE, getproject=TRUE)
###############################################################################
# nqiter=4; RTS='VRS'; nboot=50; mcells=1
# #############################################################################
# qoutS=qout;orient='ddea';dmulist=NULL;transform=TRUE
# mlist=NULL;seedval=1001;qtol=1E-6;BIGM=1E9;eps=1E-6;skipzprob=TRUE;replaceA1=TRUE
# baseqDEA=TRUE;unbounded=(-1E3);obj2test=1E-4;replaceM=FALSE;alpha=0.05;betaq=0.5
# siglist=c(0.10,0.05,0.01);CILag=1;printlog=TRUE;prntmod=100;printtxt=''
# getproject=TRUE;getbootpeers=TRUE;solver='highs'
###############################################################################
#' @export
###############################################################################
qDEA=function(X,Y,qout=1/nrow(X),qoutS=qout,X0=NULL,Y0=NULL,DX0=NULL,DY0=NULL,
orient='out',RTS='CRS',dmulist=NULL,nqiter=1,nboot=0,transform=TRUE,mcells=5,
mlist=NULL,seedval=1001,qtol=1E-6,BIGM=1E9,eps=1E-6,skipzprob=TRUE,replaceA1=FALSE,
baseqDEA=FALSE,unbounded=(-1E3),obj2test=1E-4,replaceM=FALSE,alpha=0.05,betaq=0.5,
siglist=c(0.10,0.05,0.01),CILag=1,printlog=TRUE,prntmod=100,printtxt='',
getproject=FALSE,getbootpeers=FALSE,solver='highs'){
###############################################################################
# print(paste('printtxt = ',printtxt))
# print(paste('solver = ',solver))
###############################################################################
###############################################################################
#setDTthreads(1)
if(baseqDEA==TRUE) qoutS=qout
# truncate qout
if(qout<1E-9) qout = 1E-9
if(qoutS<1E-9) qoutS = 1E-9
if(qoutS>qout) qoutS = qout
###############################################################################
# data check and construction
#######################################
if(!is.matrix(X)) X=as.matrix(X)
if(!is.matrix(Y)) Y=as.matrix(Y)
if(nrow(X)!=nrow(Y)) stop('STOP: nrow(X)!=nrow(Y)')
(nrX=nrow(X))
(ncX=ncol(X));(ncY=ncol(Y))
(ncYX=ncY+ncX)
#######################################
if(length(which(is.na(as.vector(X)))) >0) {
stop('STOP: The X matrix cannot contain NA elements')
}
if(length(which(is.na(as.vector(Y)))) >0) {
stop('STOP: The Y matrix cannot contain NA elements')
}
###############################################################################
###############################################################################
dmulist0=dmulist
orient
if(orient=='ddea'|orient=='DDEA'){
if(is.null(DX0)|is.null(DY0)|is.null(DX0)|is.null(DY0)) {
stop('STOP: with orient = ddea, X0, Y0, DX0 and DY0 must be non-null')
}
if(length(which(is.na(c(as.vector(X0),as.vector(Y0),
as.vector(DX0),as.vector(DY0)) ))) >0) {
stop('STOP: X0,Y0,DX0, and DY0 cannot contain NA elelemts')
}
if(is.vector(X0) & is.vector(Y0)&is.vector(DX0) & is.vector(DY0)){
ndmu0=1
if(length(X0)!=ncol(X) | length(Y0)!=ncol(Y)|length(DX0)!=ncol(X)
| length(DY0)!=ncol(Y)) {
stop('STOP: X0, Y0, DX0, DY0 vectors must of appropriate length')
}
X0=matrix(X0,nrow=ndmu0)
Y0=matrix(Y0,nrow=ndmu0)
DX0=matrix(DX0,nrow=ndmu0)
DY0=matrix(DY0,nrow=ndmu0)
} # end if(is.vector(X0,Y0,DX0 & DY0))
(ndmu0=nrow(X0))
if(ncol(DX0)!=ncol(X)) stop('STOP:ncol(DX0)!=ncol(X)')
if(ncol(DY0)!=ncol(Y)) stop('STOP:ncol(DY0)!=ncol(Y)')
if(!is.null(dmulist)){
dmulist0=dmulist
if(max(dmulist)>ndmu0) stop('STOP: max dmulist>nrow(x0)')
ndmu0=length(dmulist)
X0=matrix(X0[dmulist,],nrow=ndmu0)
Y0=matrix(Y0[dmulist,],nrow=ndmu0)
DX0=matrix(DX0[dmulist,],nrow=ndmu0)
DY0=matrix(DY0[dmulist,],nrow=ndmu0)
dmulist=1:nrow(X0)
}
if(is.null(dmulist)) dmulist = 1:nrow(X0)
} # end if(orient=='ddea')
################################################################################
if(orient!='ddea'&orient!='DDEA'){
DX0=NULL; DY0=NULL
}
################################################################################
################################################################################
if(!is.null(X0)|!is.null(Y0)){
if(is.vector(X0)&is.vector(Y0)){
X0=matrix(X0,nrow=1)
Y0=matrix(Y0,nrow=1)
} # end if(is.vector(X0)&is.vector(Y0)){
if(nrow(X0)!=nrow(Y0)) stop('STOP: nrow(X0)!=nrow(Y0)')
if(ncol(X0)!=ncol(X)) stop('STOP: ncol(X0)!=ncol(X)')
if(ncol(Y0)!=ncol(Y)) stop('STOP: ncol(Y0)!=ncol(Y)')
if(length(which(is.na(as.vector(X0)))) >0) {
stop('STOP: The X0 matrix cannot contain NA elements')
}
if(length(which(is.na(as.vector(Y0)))) >0) {
stop('STOP: The Y0 matrix cannot contain NA elements')
}
} # end if(!is.null(X0)|!is.null(Y0)){
###############################################################################
if(is.null(dmulist)&(is.null(X0))){
dmulist = 1:nrow(X)
X0=X; Y0=Y
}
###############################################################################
if(!is.null(dmulist)&(!is.null(X0))){
################################
if(length(dmulist)>nrow(X0)){
stop('STOP: length(dmulist)>nrow(X0)')
}
if(max(dmulist)>nrow(X0)){
stop('STOP: max(dmulist)>nrow(X0)')
}
if(length(dmulist)<nrow(X0)){
dmulist0=dmulist
ndmu0=length(dmulist)
X0=matrix(X0[dmulist,],nrow=ndmu0)
Y0=matrix(Y0[dmulist,],nrow=ndmu0)
dmulist=1:nrow(X0)
}
###################################
} # end if(!is.null(dmulist)&(!is.null(X0)|!is.null(Y0)))
#######################################
if(!is.null(dmulist)&is.null(X0)){
ndmu0=length(dmulist)
dmulist0=dmulist
X0=matrix(X[dmulist,],nrow=ndmu0)
Y0=matrix(Y[dmulist,],nrow=ndmu0)
} else{
if(is.null(dmulist)&!is.null(X0)){
dmulist=1:nrow(X0)
dmulist0=dmulist
}
}
###############################################################################
###############################################################################
(ndmu0=nrow(X0))
#######################################
if(is.null(DX0)|is.null(DY0)) {
(DX0=matrix(0,ndmu0,ncX))
(DY0=matrix(0,ndmu0,ncY))
if(orient=='in'|orient=='IN'){
(DX0=matrix(X0,nrow=ndmu0))
}
if(orient=='out'|orient=='OUT'){
(DY0=matrix(Y0,nrow=ndmu0))
}
if(orient=='inout'|orient=='INOUT'){
(DX0=matrix(X0,nrow=ndmu0))
(DY0=matrix(Y0,nrow=ndmu0))
}
if(orient=='oneone'|orient=='ONEONE'){
DX0=X0^0
DY0=Y0^0
}
} # end if is.null(DX0 or DY0)
(nDMU0=nrow(X0))
###############################################################################
###############################################################################
## highs blows up if abs of any sparse matrix ra values are <1e-9
#X[abs(X)< 2e-9] = 2e-9; Y[abs(Y)< 2e-9] = 2e-9
#X0[abs(X0)< 2e-9] = 2e-9; Y0[abs(Y0)< 2e-9] = 2e-9
#DX0[abs(DX0)< 2e-9] = 2e-9; DY0[abs(DY0)< 2e-9] = 2e-9
###############################################################################
time1=NULL;time2=NULL;time3=NULL
###############################################################################
time10=my_seconds()
tmp0=qDEA_solve(X,Y,qout=qout,qoutS=qoutS,X0=X0,Y0=Y0,DX0=DX0,DY0=DY0,RTS=RTS,
nqiter=nqiter,qtol=qtol,BIGM=BIGM,eps=eps,skipzprob=skipzprob,
unbounded=unbounded,obj2test=obj2test,baseqDEA=baseqDEA,replaceA1=replaceA1,
printlog=printlog,prntmod=prntmod,printtxt=printtxt,dmulist0=dmulist0,
solver=solver)
time11=my_seconds()
(time1=time11-time10)
###############################################################################
###############################################################################
BOOT='NOT RUN'; BOOTq='NOT RUN'; BOOT_M='NOT RUN'; BOOT_dmus='NOT RUN'
BPEERS=list('NOT RUN'); BPEERSq=list('NOT RUN')
###############################################################################
#debugoff
if(nboot>0){
###############################################################################
##
(Bdrop = which(is.na(tmp0$effvals[,1])))
if(length(Bdrop)==0) Bdrop=NULL
(Bdropq = which(is.na(tmp0$effvals[,3])))
##
(n=nrX)
if(!is.null(mlist)) mcells=length(mlist)
##
if(is.null(mlist[1])){
# compute mlist values
if(mcells==1){
if(qout<1e-6) mlist=floor(2.5*sqrt(nrX))
if(qout>=1e-6) mlist=floor(max(2.5*sqrt(nrX),1/qout))+1
} # end if(mcells==1){
if(mcells>=2){
(n1=floor(1.0*sqrt(nrX)))
(if(qout>0) n1=floor(max(1/qout,sqrt(nrX)))+1)
(n2=floor(min(5*n1,nrX/2))+1)
if(n2>0.5*nrX) (n2=floor(0.5*nrX)+1)
(mlist=round(exp(seq(log(n1),log(n2),length.out=mcells)),0))
if(max(table(mlist)>1)) (mlist=floor(seq(n1,n2,length.out=mcells)))
} # end if(mcells>=2)
mlist
} # end if(is.null(mlist[1]))
###############################################################################
mlist
(mmax=max(mlist))
if(mmax>nrX){
stop('STOP: User must provide subsample size vector with largest element < nrow(X)')
}
BOOT_dmus = matrix(NA,nboot,mmax)
BOOT=array(NA,dim=c(nboot,mcells,nDMU0))
BOOTq=BOOT
time20=my_seconds()
###############################################################################
nb=1
# debugoff
for(nb in 1:nboot){
if(printlog==TRUE&nb%%prntmod==0){
print(paste(printtxt,'bootstrap',nb,'of',nboot))
}
set.seed(seedval+nb)
mpick=sample(1:nrX,mmax,replace=replaceM)
BOOT_dmus[nb,]=mpick
YM=as.matrix(Y[mpick,]); XM=as.matrix(X[mpick,])
tmp2=qDEA_mlist(XM=XM,YM=YM,mpick=mpick,Bdrop=Bdrop,Bdropq=Bdropq,qout=qout,
qoutS=qoutS,X0=X0,Y0=Y0,DX0=DX0,DY0=DY0,mlist=mlist,RTS=RTS,nqiter=nqiter,
unbounded=unbounded,replaceA1=replaceA1,baseqDEA=baseqDEA,
getbootpeers=getbootpeers,dmulist0=dmulist0,solver=solver)
BOOT[nb,,]=t(tmp2$EFFmlist)
BOOTq[nb,,]=t(tmp2$EFFmlistq)
peers=NULL; peersq=NULL
if(getbootpeers==TRUE){
peers = tmp2$peers
peers$nb = nb
peersq = tmp2$peersq
peersq$nb = nb
BPEERS[[nb]]=peers
BPEERSq[[nb]]=peersq
} # end if(getbootpeers==TRUE)
###############################################################################
} # end loop on nb
###############################################################################
if(getbootpeers==TRUE){
BPEERS = as.data.frame(dplyr::bind_rows(BPEERS))
BPEERSq = as.data.frame(dplyr::bind_rows(BPEERSq))
BPEERS$dmuz=ifelse(BPEERS$dmuz==0,BPEERS$dmu0,BPEERS$dmuz)
BPEERSq$dmuz=ifelse(BPEERSq$dmuz==0,BPEERSq$dmu0,BPEERSq$dmuz)
} # end if (getbootpeers==TRUE)
if(getbootpeers!=TRUE) {BPEERS=list('NOT RUN'); BPEERSq=list('NOT RUN')}
time21=my_seconds()
(time2=time21-time20)
###############################################################################
###############################################################################
} # end if(nboot>0)
###############################################################################
###############################################################################
# process results
###############################################################################
# DEA and qDEA results
###############################################################################
distMAT=tmp0$effvals;qhat=tmp0$qhat;qiter=tmp0$qiter; status=tmp0$status
PSOL=tmp0$PSOL;PSOLq1=tmp0$PSOLq1;PSOLq2=tmp0$PSOLq2
PEERS=tmp0$PEERS;PEERSq=tmp0$PEERSq
RCOST=tmp0$RCOST;RCOSTq1=tmp0$RCOSTq1;RCOSTq2=tmp0$RCOSTq2
DOUTq=tmp0$D2OUT;devstart=tmp0$devstart;devend=tmp0$devend
LP0=tmp0$LP0;LP1=tmp0$LP1;LP2=tmp0$LP2
################################################################################
PEERS$dmuz=ifelse(PEERS$dmuz==0,PEERS$dmu0,PEERS$dmuz)
if(!is.character(PEERSq)) {
PEERSq$dmuz=ifelse(PEERSq$dmuz==0,PEERSq$dmu0,PEERSq$dmuz)
}
################################################################################
distvals=distMAT[,1]
distvalsq=distMAT[,3]
effvals=distvals
effvalsq=distvalsq
######################################
if(transform==TRUE){
if(orient=='in'|orient=='IN') {
effvals=1-distvals
effvalsq=1-distvalsq
}
if(orient=='out'|orient=='OUT') {
effvals=1+distvals
effvalsq=1+distvalsq
}
} # end if(transform==TRUE){
##############################################################################
##############################################################################
# bootstrap and bias correction results
###############################################################################
distvals.bc=rep(NA,length(distvals));distvalsq.bc=rep(NA,length(distvalsq))
effvals.bc=distvals.bc;effvalsq.bc=distvalsq.bc
mpick=NA; mpickq=NA; time_boot=NA
BOOTS='NOT RUN'; BOOTs='NOT RUN'
BOOTSq='NOT RUN'; BOOTsq='NOT RUN'
CI='NOT RUN'; CIq='NOT RUN'; CIq_norm='NOT RUN'
###############################################################################
if(nboot>1){
(nsigs = length(siglist))
###############################################################################
# bias correction and confidence intervals for DEA efficiency estimates
###############################################################################
BOOTS=array(NA,dim=c(nboot,mcells,nDMU0)) #set of S matrices m in mlist by DMU
BOOTs=matrix(NA,nboot,nDMU0) #set of chosen S vectors by DMU
CI=array(NA,dim=c(nDMU0,2,nsigs))
#############################
# beta values for DDEA estimates - see Simar-Wilson ?? 2009 ??
if(RTS!='CRS'&RTS!='crs') (beta = 2/(ncX + ncY + 1))
if(RTS=='CRS'|RTS=='crs') (beta = 2/(ncX + ncY))
#############################
i=1 # for debugging
for(i in 1:nDMU0){
distvals[i]
if(!is.na(distvals[i])){
tmp=nCm_mpick(stat=distvals[i],boot=BOOT[,,i],n,mlist,beta=beta,
alpha=alpha,CILag=CILag)
(mpick[i]=tmp$mpick)
S=tmp$S
s=tmp$s
BOOTS[,,i]=S
BOOTs[,i]=s
(distvals.bc[i]=mean(s,na.rm=TRUE))
j=1 # for debugging
for(j in 1:nsigs){
(sig=siglist[j])
(CI[i,1,j]=stats::quantile(s,probs=sig/2,na.rm=TRUE))
(CI[i,2,j]=stats::quantile(s,probs=1-sig/2,na.rm=TRUE))
} #end loop for(j in 1:length(siglist))
} # end if(is.na(distvals[i]))
} # end loop for(i in 1:nDMU0){
effvals.bc=distvals.bc
###############################################################################
# bias correction and confidence intervals for qDEA estimates
###############################################################################
BOOTSq=array(NA,dim=c(nboot,mcells,nDMU0)) #set of S matrices by DMU
BOOTsq=matrix(NA,nboot,nDMU0) #set of chosen S vectors by DMU
CIq=array(NA,dim=c(nDMU0,2,nsigs))
CIq_norm=array(NA,dim=c(nDMU0,2,nsigs))
i=1 # for debugging
for(i in 1:nDMU0){
distvalsq[i]
if(!is.na(distvalsq[i])){
tmpq=nCm_mpick(stat=distvalsq[i],boot=BOOTq[,,i],n,mlist,beta=betaq,
alpha=alpha,CILag=CILag)
(mpickq[i]=tmpq$mpick)
Sq=tmpq$S
sq=tmpq$s
BOOTSq[,,i]=Sq
BOOTsq[,i]=sq
(distvalsq.bc[i]=mean(sq,na.rm=TRUE))
#j=1 # for debugging
for(j in 1:nsigs){
(sig=siglist[j])
(CIq[i,1,j]=stats::quantile(sq,probs=sig/2,na.rm=TRUE))
(CIq[i,2,j]=stats::quantile(sq,probs=1-sig/2,na.rm=TRUE))
(CIq_norm[i,1,j]=stats::qnorm(sig/2, mean(sq),stats::sd(sq,na.rm=TRUE)))
(CIq_norm[i,2,j]=stats::qnorm(1-sig/2,mean(sq),stats::sd(sq,na.rm=TRUE)))
} #end loop for(j in 1:length(siglist))
} # end if(!is.na(distvalssq[i])){
} # end loop for(i in 1:nDMU0){
effvalsq.bc=distvalsq.bc
#summary(as.vector(CIq))
#(pick=which(mpickq==1))
#CIq[pick,,]
#effq
#effq.bc
###############################################################################
# transform to efficiencies
###############################################################################
if(transform==TRUE){
if(orient=='in'|orient=='IN') {
effvals.bc=1-distvals.bc
BOOT=1-BOOT
BOOTS=1-BOOTS
BOOTs=1-BOOTs
CI=1-CI
tmp=CI[,1,]
CI[,1,]=CI[,2,]
CI[,2,]=tmp
effvalsq.bc=1-distvalsq.bc
BOOTq=1-BOOTq
BOOTSq=1-BOOTSq
BOOTsq=1-BOOTsq
CIq=1-CIq
tmp=CIq[,1,]
CIq[,1,]=CIq[,2,]
CIq[,2,]=tmp
tmp=CIq_norm[,1,]
CIq_norm[,1,]=CIq_norm[,2,]
CIq_norm[,2,]=tmp
} # end if(orient=='in'|orient=='IN')
if(orient=='out'|orient=='OUT') {
effvals.bc=1+distvals.bc
BOOT = 1 + BOOT
BOOTS = 1 + BOOTS
BOOTs = 1 + BOOTs
CI = 1 + CI
effvalsq.bc=1+distvalsq.bc
BOOTq = 1 + BOOTq
BOOTSq = 1 + BOOTSq
BOOTsq = 1 + BOOTsq
CIq = 1 + CIq
CIq_norm = 1 + CIq_norm
} # end if(orient=='out'|orient=='OUT')
} # end if(transform==TRUE)
###############################################################################
} # end if(nboot>0)
###############################################################################
###############################################################################
Y0HAT='NOT RUN'; X0HAT='NOT RUN'; Y0HATq='NOT RUN'; X0HATq='NOT RUN'
Y0HAT.bc='NOT RUN'; X0HAT.bc='NOT RUN'; Y0HATq.bc='NOT RUN'; X0HATq.bc='NOT RUN'
if(getproject==TRUE) {
Y0HAT=Y0+distvals*DY0
X0HAT=X0-distvals*DX0
Y0HATq=Y0+distvalsq*DY0
X0HATq=X0-distvalsq*DX0
if(nboot>0){
Y0HAT.bc=Y0+distvals.bc*DY0
X0HAT.bc=X0-distvals.bc*DX0
Y0HATq.bc=Y0+distvalsq.bc*DY0
X0HATq.bc=X0-distvalsq.bc*DX0
} # end if(nboot>0)
} # end if(getproject==TRUE)
###############################################################################
###############################################################################
BPEERS1='NOT RUN'
BPEERSq1='NOT RUN'
if(nboot>0&getbootpeers==TRUE){
dftmp=data.frame(cbind(m=mlist[mpick],dmu0=1:nDMU0))
BPEERS1=merge(dftmp,BPEERS)
BPEERS1=doBy::orderBy(~nb+dmu0+dmuz,data=BPEERS1)
dftmp=data.frame(cbind(m=mlist[mpickq],dmu0=1:nDMU0))
BPEERSq1=merge(dftmp,BPEERSq)
BPEERSq1=doBy::orderBy(~nb+dmu0+dmuz,data=BPEERSq1)
} # end if(getbootpeers==TRUE)
##############################################################################
BOOT_DATA = list(effvals.bc=effvals.bc,effvalsq.bc=effvalsq.bc,
distvals.bc=distvals.bc,distvalsq.bc=distvalsq.bc,mcells=mcells,
mlist=mlist,BOOT_dmus=BOOT_dmus,BOOT=BOOT,BOOTS=BOOTS,BOOTs=BOOTs,
BOOTq=BOOTq,BOOTSq=BOOTSq,BOOTsq=BOOTsq,BPEERS=BPEERS,
BPEERSq=BPEERSq,BPEERS1=BPEERS1,BPEERSq1=BPEERSq1,mpick=mpick,
mpickq=mpickq,beta=beta,betaq=betaq,siglist=siglist,CI=CI,CIq=CIq,
CIq_norm=CIq_norm) # end list BOOTDATA
##
INPUT_DATA = list(X=X,Y=Y,X0=X0,Y0=Y0,DX0=DX0,DY0=DY0,qout=qout,qoutS=qoutS,
RTS=RTS,orient=orient,baseqDEA=baseqDEA,dmulist=dmulist,
dmulist0=dmulist0) # end list INPUT_DATA
##
LP_DATA = list(status=status,qhat=qhat,qiter=qiter,PSOL=PSOL,PSOLq1=PSOLq1,
PSOLq2=PSOLq2,RCOST=RCOST,RCOSTq1=RCOSTq1,RCOSTq2=RCOSTq2,
LPModels=list(LP0=tmp0$LP0,LP1=tmp0$LP1,LP2=tmp0$LP2))# endlist LPDATA
##
PEER_DATA = list(PEERS=PEERS,PEERSq=PEERSq,DOUTq=DOUTq,BPEERS=BPEERS,
BPEERSq=BPEERSq,BPEERS1=BPEERS1,BPEERSq1=BPEERSq1) # end list PEER_DATA
##
PROJ_DATA = list(X0HAT=X0HAT,Y0HAT=Y0HAT,X0HAT.bc=X0HAT.bc,Y0HAT.bc=Y0HAT.bc,
X0HATq=X0HATq,Y0HATq=Y0HATq,X0HATq.bc=X0HATq.bc,
Y0HATq.bc=Y0HATq.bc) # end list PROJ_DATA
###############################################################################
###############################################################################
list(effvals=effvals,effvalsq=effvalsq,distvals=distvals,
distvalsq=distvalsq,distMAT=distMAT,INPUT_DATA=INPUT_DATA,BOOT_DATA=BOOT_DATA,
PEER_DATA=PEER_DATA,PROJ_DATA=PROJ_DATA,LP_DATA=LP_DATA)
###############################################################################
} # end function qDEA
###############################################################################
###############################################################################
#' qDEA_solve: Sets up LP objects and solves DDEA dual DEA and qDEA using highs
#' !!!! Intended to be called from qDEA function.!!!!
#' Note: ndmu=number of dmus in reference set, ndmu0 = number dmus to process.
#' @param X Reference dmu's = ndmu x number of inputs input matrix.
#' @param Y Reference dmu's = ndmu x number of outputs output matrix.
#' @param qout Maximal proportion of dmu's allowed external to DEA hull.
#' @param qoutS Proportion of external points to identify using EJOR qDEA slicing.
#' @param X0 Inputs for set of ndmu0 dmu's to be processed.
#' @param Y0 Outputs for set of ndmu0 dmu's to be processed.
#' @param DX0 Input directions for ndmu0 dmu's in X0 and Y0.
#' @param DY0 Output directions for ndmu0 dmu's in X0 and Y0.
#' @param RTS Returns to scale default='CRS' options are 'CRS,'VRS,'DRS','IRS.
#' @param nqiter Maximal number of qDEA iterations.
#' @param qtol q search tolerance with iterative qDEA.
#' @param BIGM Default Big M in RHS of qDEA stage 2 process.
#' @param eps Zero-test criterion. abs(x) > eps implies x!=0
#' @param skipzprob Skip qDEA if qout=0.
#' @param unbounded qDEA reported as unbounded if obj<=unbounded.
#' @param obj2test Convergence tol for objective in iterative qDEA.
#' @param replaceA1 Put dmu0's data in first row of reference sets.
#' @param baseqDEA Use basic qDEA model from EJOR article
#' @param printlog Progress of dmu's solved when (X0,Y0) >1 dmus.
#' @param prntmod Print progress every prntmod dmus. default=100
#' @param printtxt Additional text to print with progress printlog.
#' @param dmulist0 DMU0 index in originally inputs X0,Y0,DX0,and DY0
#' @param solver LP solver Default='highs'
#' @return A list containing the following components:
#' @return effvals = ndmu0 by 3 matrix with DDEA,qDDEA-S1,qDDEA-S2 distances
#' (NA indicates qDEA for given DMU was unbounded)
#' @return qhat = proportion of dmu's external to hull in qDEA solution
#' (NA indicates qDEA for given DMU was unbounded)
#' @return qiter = number of qDEA iterations completed.
#' @return D2OUT = data frame indicating qDEA external reference dmus
#' @return status = ndmu0 by 3 matrix with LP status of DDEA,qDDEA-S1,qDDEA-S2
#' @return PSOL = ndmu0 by ? matrix with DDEA (dual-side) solutions for each dmu0
#' @return PSOLq1 = ndmu0 by ? matrix with qDDEA-S1 solutions for each dmu0
#' @return PSOLq2 = ndmu0 by ? matrix with qDDEA-S2 solutions for each dmu0
#' @return PEERS = dataframe of non-zero DDEA dmuz's and projection weights for each dmu0
#' @return PEERSq = dataframe of non-zero qDDEA-S2 dmuz's projection weights for each dmu0
#' @return RCOST = ndmu0 by ? matrix of LP reduced costs for DDEA solutions
#' @return RCOSTq1 = ndmu0 by ? matrix of LP reduced costs for qDDEA-S1 solutions
#' @return RCOSTq2 = ndmu0 by ? matrix of LP reduced costs for qDDEA-S2 solutions
#' @return devstart = LP index for column of first qDDEA-S1 LPM "deviation" value
#' @return devend = LP index for column of last qDDEA-S1 LPM "deviation" value
#' @return LP0 = DDEA LP object for last dmu0 processed
#' @return LP1 = qDDEA-S1 LP object for last dmu0 processed
#' @return LP2 = qDDEA-S2 LP object for last dmu0 processed
##############################################################################
# debugging qDEA_solve
##############################################################################
# X0=X00;Y0=Y00;DX0=DX00;DY0=DY00
##############################################################################
# (qout = 2/8); (qoutS = 1/8)
# (qout = 1/8); (qoutS = 0/8)
##############################################################################
##############################################################################
# X=Xm;Y=Ym
##############################################################################
# (qoutS = qout);(RTS='VRS');(nqiter=4);(qtol=1e-6);(BIGM=1E9);(eps=1e-6)
# (skipzprob=TRUE);(unbounded = (-1E3));(obj2test=1e-4);(printlog=TRUE)
# (prntmod=1); (printtxt=''); (replaceA1=TRUE)
# (baseqDEA=FALSE)
##############################################################################
#' @export
###############################################################################
qDEA_solve=function(X,Y,qout=0.10,qoutS=qout,X0,Y0,DX0,DY0,RTS='CRS',
nqiter=1,qtol=1E-6,BIGM=1E9,eps=1E-6,skipzprob=TRUE,unbounded=(-1E3),
obj2test=1E-4,replaceA1=FALSE,baseqDEA=FALSE,printlog=TRUE,prntmod=100,printtxt='',
dmulist0=1:nrow(X0),solver='highs'){
###############################################################################
if(nrow(X)!=nrow(Y)) stop('Error: nrow(X)!= nrow(Y)')
if(nrow(X0)!=nrow(Y0)) stop('Error: nrow(X0)!= nrow(Y0)')
if(ncol(X)!=ncol(X0)) stop('Error: ncol(X)!=ncol(X0)')
if(ncol(Y)!=ncol(Y0)) stop('Error: ncol(Y)!=ncol(Y0)')
if(nrow(X0)!=nrow(DX0)) stop('Error: nrow(X0)!=nrow(DX0)')
if(nrow(Y0)!=nrow(DY0)) stop('Error: nrow(Y0)!=nrow(DY0)')
if(ncol(X0)!=ncol(DX0)) stop('Error: ncol(X0)!=ncol(DX0)')
if(ncol(Y0)!=ncol(DY0)) stop('Error: ncol(Y0)!=ncol(DY0)')
###############################################################################
(ncY=ncol(Y));(ncX=ncol(X))
(ncYX=ncY+ncX)
(nobs=nrow(Y))
(devstart=ncYX+2)
(devend=devstart+nobs-1)
(nDMU0=nrow(X0))
###############################################################################
###############################################################################
# LP output data
D2OUT=matrix(NA,nDMU0,nobs)
effvals=matrix(NA,nDMU0,3)
status=matrix(NA,nDMU0,3)
qhat=rep(NA,nDMU0)
qiter=rep(0,nDMU0)
#####################################################################
LP0='NOT RUN'; LP1='NOT RUN'; LP2='NOT RUN'
PSOL='NOT RUN'; PSOLq1='NOT RUN';PSOLq2='NOT RUN'
PEERS='NOT RUN'; DSOL1='NOT RUN';PEERSq='NOT RUN'
RCOST='NOT RUN';RCOSTq1='NOT RUN';RCOSTq2='NOT RUN'
#####################################################################
###################################################################
# SETUP SECTIONS
###################################################################
# debug off
# set up DDEA dual model
time0=my_seconds()
LP0=DEAbuild(X,Y,X0,Y0,DX0,DY0,dmu0=1,RTS=RTS,
unbounded=unbounded,solver=solver)
# Matrices to store primal and dual solutions
PSOL=matrix(NA,nDMU0,LP0$nc)
PEERS=matrix(NA,nDMU0,LP0$nr)
RCOST=matrix(NA,nDMU0,LP0$nc)
#####################################################################
#####################################################################
# set up for qDEA
#debug off
if(qout>=qtol|skipzprob==FALSE) {
# set up qDEA LP models with SMM CMO
LP1=qDEAbuild(X,Y,X0,Y0,DX0,DY0,dmu0=1,qout=qoutS,RTS=RTS,
unbounded=unbounded,solver=solver)
LP2=DEAbuild(X,Y,X0,Y0,DX0,DY0,dmu0=1,RTS=RTS,
unbounded=unbounded,solver=solver)
# Matrices to store primal and dual solutions
PSOLq1=matrix(NA,nDMU0,LP1$nc)
DSOL1=matrix(NA,nDMU0,LP1$nr)
RCOSTq1=matrix(NA,nDMU0,LP1$nc)
PSOLq2=matrix(NA,nDMU0,LP2$nc)
PEERSq=matrix(NA,nDMU0,LP2$nr)
RCOSTq2=matrix(NA,nDMU0,LP2$nc)
tau0=LP1$tau0
(tau0inv=LP1$ra[LP1$qchng[1]])
# save initial values for row boundaries
rhs1=LP1$rhs
rhs2=LP2$rhs
rlower1=LP1$rlower
rlower2=LP2$rlower
rupper1=LP1$rupper
rupper2=LP2$rupper
###############################################################################
} # end setup #if(qoutS>=qtol|skipzprob==FALSE) {
###############################################################################
#####################################################################
# run dual DDEA and qDEA models
#####################################################################
i=1
# debug off
for(i in 1:nDMU0){
if(printlog==TRUE){
if(i==1|i%%prntmod==0|i==nDMU0) print(paste(printtxt,'on dmu',i,'of',nDMU0))
} # end if(printlog==TRUE)
# DDEA DUAL section
# modify lp objective and matrix coefficients
(LP0$obj[1:ncYX]=LP0$DOBJ[i,])
LP0$ra[LP0$dyxchngC]= LP0$DYX[i,]
LP0$ra[LP0$yxchngR] = -LP0$DOBJ[i,]
if(replaceA1==TRUE){
yx = -LP0$DOBJ[i,]
LP0$ra[LP0$yxchngC]=yx
}
LP0$dmu0=i
LPSOL0=LPSOLVER(LP0,solver=solver)
names(LPSOL0)
# [1] "status" "objval" "solution" "dual" "rcost" "lhs"
# [7] "TIME" "SMtime" "proctime" "LP"
(status[i,1]=LPSOL0$status)
(obj0=LPSOL0$objval)
(effvals[i,1]=obj0)
(PSOL[i,]=LPSOL0$solution)
(PEERS[i,]=LPSOL0$dual)
(RCOST[i,]=LPSOL0$rcost)
if(obj0<=unbounded+obj2test){
effvals[i,1]=NA
PSOL[i,] = NA
PEERS[i,] = NA
RCOST[i,]= NA
} # end if(obj0<=unbounded+obj2test)
LP2 = LP0
#############################################
# end dual DDEA LP section
#############################################
#############################################
# start qDEA section
#############################################
###############################################
# debugoff
if(qoutS>=qtol|skipzprob==FALSE) {
# modify lp objectives
(LP1$obj[1:ncYX]=LP1$DOBJ[i,])
(LP2$obj[1:ncYX]=LP2$DOBJ[i,])
# modify/restore lp matrix coefficients
# restore original tau0inv
LP1$ra[LP1$qchng[1]]=tau0inv
# modify with DMU i's directions
dyx = LP1$DYX[i,]
dobj = -LP1$DOBJ[i,]
LP1$ra[LP1$dyxchngC] = dyx
LP1$ra[LP1$yxchngR] = dobj
LP2$ra[LP2$dyxchngC] = dyx
LP2$ra[LP2$yxchngR] = dobj
if(replaceA1==TRUE){
# yx = c(Y0[i,],-X0[i,])
yx = -LP1$DOBJ[i,]
LP1$ra[LP1$yxchngC]=yx
LP2$ra[LP2$yxchngC]=yx
}
(LP1$dmu0=i)
(LP2$dmu0=i)
# restore original lp rhs
LP1$rupper=rupper1
LP2$rupper=rupper2
# qDEA stage 1
(pointsout0=floor(qoutS*nobs))
LPSOL1=LPSOLVER(LP1,solver=solver)
names(LPSOL1)
# [1] "status" "objval" "solution" "dual" "rcost" "lhs"
# [7] "TIME" "SMtime" "proctime" "LP"
(status[i,2]=LPSOL1$status)
(obj1=LPSOL1$objval)
effvals[i,2]=obj1
PSOLq1[i,]=LPSOL1$solution
DSOL1[i,]=LPSOL1$dual
RCOSTq1[i,]=LPSOL1$rcost
xvals1=LPSOL1$solution
devs1=xvals1[devstart:devend]
D1out=ifelse(devs1>eps,1,0)
sum(D1out)
# qDEA stage 2
LP2$rhs[1:nobs]=D1out*BIGM;
LP2$rupper[1:nobs]=D1out*BIGM;
LPSOL2=LPSOLVER(LP2,solver=solver)
names(LPSOL2)
# [1] "status" "objval" "solution" "dual" "rcost" "lhs"
# [7] "TIME" "SMtime" "proctime" "LP"
(status[i,3]=LPSOL2$status)
(obj2=LPSOL2$objval)
effvals[i,3]=obj2
PSOLq2[i,]=LPSOL2$solution
(PEERSq[i,]=LPSOL2$dual)
RCOSTq2[i,]=LPSOL2$rcost
(xvals2=LPSOL2$solution)
LHS2=LPSOL2$lhs
D2out=ifelse(LHS2[1:nobs]>eps,1,0)
sum(D2out)
(qhat[i]=sum(D2out)/nobs)
qiter[i]=1
(qleft= qoutS-qhat[i])
(pointsout=as.integer(qleft*nobs))
obj2diff = 1E3
# debug off
if(nqiter>1&abs(qoutS-qhat[i])>qtol) {
# debugoff
while(qiter[i]<nqiter & pointsout>0&qleft>eps) {
(qleft= qoutS-qhat[i])
floor(qleft*nobs)
########################################################################
(pointsout=floor(qleft*nobs))
########################################################################
(q0 = qleft-eps)
(q1 = (floor(qleft*nobs)+1)/nobs - eps)
floor(q1*nobs)
tau2 = q1
(tau2inv=1/tau2)
(LP1$ra[LP1$qchng[1]]=tau2inv)
##############################################
LP1$rhs[1:nobs]=D2out*BIGM;
LP1$rupper[1:nobs]=D2out*BIGM;
LPSOL1=LPSOLVER(LP1,solver=solver)
names(LPSOL1)
# [1] "status" "objval" "solution" "dual" "rcost" "lhs"
# [7] "TIME" "SMtime" "proctime" "LP"
(status2=LPSOL1$status)
(obj1=LPSOL1$objval)
xvals1=LPSOL1$solution
devs1=xvals1[devstart:devend]
D1out=ifelse(devs1>eps,1,0)
sum(D1out)
# qDEA-stage2
tmp=D2out+D1out
tmp2=ifelse(tmp>0,1,0)
sum(tmp2)
LP2$rhs[1:nobs]=tmp*BIGM;
LP2$rupper[1:nobs]=tmp*BIGM;
LPSOL2=LPSOLVER(LP2,solver=solver)
names(LPSOL2)
# [1] "status" "objval" "solution" "dual" "rcost" "lhs"
# [7] "TIME" "SMtime" "proctime" "LP"
(status3=LPSOL2$status)
(obj2=LPSOL2$objval)
# ###########################################################################
# debugoff
if(baseqDEA!=TRUE) {
if(abs(obj2-effvals[i,3])<=eps & pointsout>=1){
xvals2=LPSOL2$solution
LHS2=LPSOL2$lhs
psol2=LPSOL2$solution
dsol2=LPSOL2$dual
d2out=ifelse(LHS2[1:nobs]>eps,1,0)
LPlist=list(NULL)
(dualList=which(abs(dsol2[1:nobs])>eps))
(nduals=length(dualList))
# debugoff
if(nduals>=1){
qefftest=100
nd=1
#debugoff
for(nd in 1:nduals){
LP22=LP2
(ndpick=dualList[nd])
LP22$rhs[ndpick]=BIGM
LP22$rupper[ndpick]=BIGM
LPSOL=LPSOLVER(LP22,solver=solver)
(qefftest[nd]=LPSOL$objval)
LPlist[[nd]]=list(LP=LP22)
} # end loop for(nd in 1:nduals)
(ndpick=which(qefftest==min(qefftest))[1])
(obj22=qefftest[ndpick])
#debugoff
if(obj22<(obj2-eps)){
LP2=LPlist[[ndpick]]$LP
LPSOL2=LPSOLVER(LP2,solver=solver)
status3=LPSOL2$status
obj2=LPSOL2$objval
} # end if(obj22<(obj2-eps)){
} # end loop if(length(nduals)>=1)
} # end if(abs(obj2-effvals[i,3])<=eps)
} # end if(baseqDEA!=TRUE)
############################################################################
(obj2diff=abs(obj2-effvals[i,3]))
xvals2=LPSOL2$solution
LHS2=LPSOL2$lhs
(qiter[i]=qiter[i]+1)
if(obj2diff>obj2test){
(status[i,2]=status2)
(status[i,3]=status3)
(effvals[i,2]=obj1)
(effvals[i,3]=obj2)
(PSOLq2[i,]=LPSOL2$solution)
(PEERSq[i,]=LPSOL2$dual)
RCOSTq2[i,]=LPSOL2$rc
D2out=ifelse(LHS2[1:nobs]>eps,1,0)
(qhat[i]=sum(D2out)/nobs)
} # end if(obj2diff>obj2test)
(qleft= qoutS-qhat[i])
(floor(qleft*nobs))
(pointsout=floor(qleft*nobs))
} # end #while(qiter[i]<nqiter & abs(qoutS-qhat[i])>qtol & pointsout>0 & obj2diff>eps) {
}# end if(nqiter>1&abs(qhat[i]-qoutS)>qtol)
D2OUT[i,]=D2out
} # end if(qoutS>(1.0/nobs)|skipzprob==FALSE)
###############################################################################
# if(qout>(qoutS+eps)) {
if(baseqDEA==FALSE&qout>(qoutS+eps)) {
(ndmus=nrow(X))
(nout0=floor(qout*ndmus))
tmp2=iter_delete(LP2,ndmus,nout0,iterlim=max(100,nout0),
BIGM=BIGM,unbounded=unbounded,eps=eps,solver=solver)
LP2=tmp2$LP2
LP2sol=tmp2$LP2sol
(status[i,3]=LP2sol$status)
(obj2=tmp2$obj2)
effvals[i,3]=obj2
if(!is.na(obj2)){
PSOLq2[i,]=LP2sol$solution
PEERSq[i,]=LP2sol$dual
RCOSTq2[i,]=LP2sol$rcost
(xvals2=LP2sol$solution)
(outlist=tmp2$outlist)
D2out=rep(0,ndmus)
D2out[outlist]=1
sum(D2out)
(qhat[i]=sum(D2out)/nobs)
D2OUT[i,]=D2out
} # end if(qout>0)
(qiter[i]=qiter[i]+tmp2$iter)
} # end if(qout>(qoutS+eps))
################################################################################
if(qout>=qtol|skipzprob==FALSE) {
if(is.na(obj2)|obj2<=unbounded+obj2test){
PSOLq2[i,] = NA
PEERSq[i,] = NA
RCOSTq2[i,]= NA
# D2OUT[i,] = NA
} # end if(obj2<=unbounded)
} # end if(qout>=qtol|skipzprob==FALSE)
} # end loop on i
###############################################################################
(status[,1]=ifelse(effvals[,1]<=unbounded+obj2test,NA,status[,1]))
(effvals[,1]=ifelse(effvals[,1]<=unbounded+obj2test,NA,effvals[,1]))
(status[,2]=ifelse(effvals[,3]<=unbounded+obj2test,NA,status[,2]))
(status[,3]=ifelse(effvals[,3]<=unbounded+obj2test,NA,status[,3]))
(effvals[,2]=ifelse(effvals[,3]<=unbounded+obj2test,NA,effvals[,2]))
(effvals[,3]=ifelse(effvals[,3]<=unbounded+obj2test,NA,effvals[,3]))
(napick= which(is.na(effvals[,3])))
(qhat[napick]=NA)
#####
(tmp=PEERS)
(PEERS=A2SM(matrix(PEERS[,1:nrow(X)],nrow=nrow(PEERS)),SMM='RCI',NA_flag=-1E12))
(PEERS=as.data.frame(cbind(PEERS$ia,PEERS$ja,-1*PEERS$ra)))
names(PEERS)=c('dmu0','dmuz','z')
if(replaceA1==TRUE) PEERS$dmuz=ifelse(PEERS$dmuz==1,0,PEERS$dmuz)
PEERS$dmuz = ifelse(is.na(PEERS$z),NA,PEERS$dmuz)
(PEERS = unique(PEERS))
if(!is.null(dmulist0)) PEERS$dmu0=dmulist0[PEERS$dmu0]
PEERS
#######################
if(qout>=qtol|skipzprob==FALSE) {
## process peer data
(tmpq=PEERSq)
PEERSq[is.na(PEERSq)]=0
PEERSq=A2SM(matrix(PEERSq[,1:nrow(X)],nrow=nrow(PEERSq)),SMM='RCI',NA_flag=-1E12)
(PEERSq=as.data.frame(cbind(PEERSq$ia,PEERSq$ja,-1*PEERSq$ra)))
names(PEERSq)=c('dmu0','dmuz','z')
if(replaceA1==TRUE) PEERSq$dmuz=ifelse(PEERSq$dmuz==1,0,PEERSq$dmuz)
PEERSq$dmuz = ifelse(is.na(PEERSq$z),NA,PEERSq$dmuz)
(PEERSq = unique(PEERSq))
if(!is.null(dmulist0)) PEERSq$dmu0=dmulist0[PEERSq$dmu0]
PEERSq
## process D2OUT data
tmpD2=D2OUT
D2OUT[is.na(D2OUT)]=0
D2OUT=A2SM(as.matrix(D2OUT),SMM='RCI',NA_flag=-1E12)
D2OUT=as.data.frame(cbind(D2OUT$ia,D2OUT$ja,D2OUT$ra))
names(D2OUT)=c('dmu0','dmuout','dout')
if(replaceA1==TRUE) D2OUT$dmuout=ifelse(D2OUT$dmuout==1,0,D2OUT$dmuout)
D2OUT$dmuout = ifelse(is.na(D2OUT$dout),NA,D2OUT$dmuout)
D2OUT = unique(D2OUT)
D2OUT = D2OUT[,c('dmu0','dmuout')]
} # end if(qout>=qtol|skipzprob==FALSE)
#######################
###############################################################################
list(effvals=effvals,qhat=qhat,qiter=qiter,D2OUT=D2OUT,status=status,
PSOL=PSOL,PSOLq1=PSOLq1,PSOLq2=PSOLq2,PEERS=PEERS,PEERSq=PEERSq,
RCOST=RCOST,RCOSTq1=RCOSTq1,RCOSTq2=RCOSTq2,devstart=devstart,devend=devend,
LP0=LP0,LP1=LP1,LP2=LP2)
###############################################################################
} # end function qDEA_solve
###############################################################################
###############################################################################
#' qDEA_mlist: Obtain subsample qDEA results for a set (vector) of subsample sizes
#' !!!Intended to be called from qDEA function!!!!
#'
#' mcells = number of subsample sizes
#' @param XM Input levels for subsample of reference set X.
#' @param YM Output levels for subsample of reference set Y.
#' @param mpick Index of subsample dmus in original (X,Y) reference set.
#' @param qout Maximal proportion of dmu's allowed external to DEA hull.
#' @param qoutS Proportion of external points to identify using qDEA slicing.
#' @param Bdrop Index list of DMU's whose initial DDEA solutions were unbounded.
#' @param Bdropq Index list of DMU's whose initial qDEA solutions were unbounded.
#' @param X0 Inputs for set of ndmu0 dmu's to be processed.
#' @param Y0 Outputs for set of ndmu0 dmu's to be processed.
#' @param DX0 Input directions for ndmu0 dmu's in X0 and Y0.
#' @param DY0 Output directions for ndmu0 dmu's in X0 and Y0.
#' @param mlist Vector of subsample sizes
#' @param RTS Returns to scale default='CRS' options are 'CRS,'VRS,'DRS','IRS.
#' @param nqiter Maximal number of qDEA iterations.
#' @param qtol qout search tolerance with iterative qDEA.
#' @param BIGM Default Big M in RHS of qDEA stage 2 process.
#' @param eps Search tolerance in qDEA improvement tests.
#' @param skipzprob Skip qDEA if qout=0.
#' @param unbounded qDEA reported as unbounded if obj<=unbounded.
#' @param obj2test Converge tol for objective in iterative qDEA.
#' @param replaceA1 Put dmu0's data in first row of reference sets. .
#' @param baseqDEA Use basic qDEA model from EJOR article
#' @param printlog Progress of dmu's solved when (X0,Y0) >1 dmus.
#' @param prntmod Print progress every prntmod dmus.
#' @param getbootpeers Return dmu projection weights for each dmu and each bootstrap.
#' @param dmulist0 DMU0 index in originally inputs X0,Y0,DX0,and DY0
#' @param solver LP solver Default='highs'
#' @return A list containing the following components:
#' @return EFFmlist = ndmu0 by mcell matrix of subsample DDEA distances.
#' @return EFFmlistq = ndmu0 by mcell matrix of subsample qDEA distances.
#' @return XM = Input levels for subsample of reference set X.
#' @return YM = Output levels for subsample of reference set X.
#' @return qout = Maximal proportion of dmu's allowed external to DEA hull.
#' @return qoutS = Proportion of external points to identify using qDEA slicing.
#' @return X0 = Inputs for set of ndmu0 dmu's to be processed.
#' @return Y0 = Outputs for set of ndmu0 dmu's to be processed.
#' @return DX0 = Input directions for ndmu0 dmu's in X0 and Y0.
#' @return DY0 = Output directions for ndmu0 dmu's in X0 and Y0.
#' @return mlist = Vector of subsample sizes
#' @return mcells = Number of subsample sizes
#' @return peers = A data frame containing DDEA peer dmus and projection weights
#' @return peersq = A data frame containing qDEA peer dmus and projection weights
###############################################################################
# debug qDEA_mlist
###############################################################################
# XM=XM;YM=YM;X0=X0;Y0=Y0;DX0=DX0;DY0=DY0
# (qout=qout); (qoutS=qoutS)
# Bdrop=Bdrop;Bdropq=Bdropq
##############################################################################
# RTS=RTS;nqiter=nqiter;qtol=qtol;BIGM=BIGM;eps=eps;skipzprob=TRUE;unbounded=unbounded
# obj2test=obj2test;printlog=printlog;prntmod=prntmod;printtxt=printtxt;replaceA1=FALSE
# getbootpeers=getbootpeers
###############################################################################
#' @export
###############################################################################
qDEA_mlist=function(XM,YM,mpick,qout=0.10,qoutS=qout,Bdrop=NULL,Bdropq=NULL,
X0=XM,Y0=YM,DX0=XM,DY0=YM,mlist,RTS='CRS',nqiter=1,qtol=1e-6,BIGM=1E9,eps=1e-6,
skipzprob=TRUE,unbounded=(-1E3),obj2test=1e-4,replaceA1=replaceA1,baseqDEA=baseqDEA,
printlog=FALSE,prntmod=100,getbootpeers=FALSE,dmulist0=1:nrow(X0),solver='highs') {
###############################################################################
#setDTthreads(1)
###############################################################################
if(nrow(XM)!=nrow(YM)) stop('Error: nrow(X)!= nrow(Y)')
if(nrow(X0)!=nrow(Y0)) stop('Error: nrow(X0)!= nrow(Y0)')
if(ncol(XM)!=ncol(X0)) stop('Error: ncol(X)!=ncol(X0)')
if(ncol(YM)!=ncol(Y0)) stop('Error: ncol(Y)!=ncol(Y0)')
if(nrow(X0)!=nrow(DX0)) stop('Error: nrow(X0)!=nrow(DX0)')
if(nrow(Y0)!=nrow(DY0)) stop('Error: nrow(Y0)!=nrow(DY0)')
if(ncol(X0)!=ncol(DX0)) stop('Error: ncol(X0)!=ncol(DX0)')
if(ncol(Y0)!=ncol(DY0)) stop('Error: ncol(Y0)!=ncol(DY0)')
###############################################################################
(mcells=length(mlist))
(ncY=ncol(YM));(ncX=ncol(XM))
(ncYX=ncY+ncX)
(nDMU0=nrow(X0))
EFFmlist=matrix(NA,nDMU0,mcells)
EFFmlistq=matrix(NA,nDMU0,mcells)
# data frame for matching dmu index from XM with dmu index from original X
dmuset = data.frame(cbind(dmuz=mpick,dmuzm=1:length(mpick)))
#head(dmuset)
#data frame for putting unbounded dmus back in PEERSq object with NA entries
dmu0_ub = data.frame(cbind(dmu0=1:nDMU0,ub=rep(NA,nDMU0)))
#head(dmu_ub)
peers = 'NOT RUN'
peersq = 'NOT RUN'
if(getbootpeers==TRUE){
peers = list('NOT RUN')
peersq =list('NOT RUN')
}
#######################################
j=1
# debugoff
for(j in 1:mcells){
(m=mlist[j])
Ym=as.matrix(YM[1:m,]);Xm=as.matrix(XM[1:m,])
tmp=qDEA_solve(X=Xm,Y=Ym,qout=qout,qoutS=qoutS,X0=X0,Y0=Y0,DX0=DX0,DY0=DY0,
RTS=RTS,nqiter=nqiter,qtol=qtol,BIGM=BIGM,eps=eps,skipzprob=skipzprob,
unbounded=unbounded,obj2test=obj2test,printlog=printlog,
prntmod=prntmod,replaceA1=replaceA1,baseqDEA=baseqDEA,dmulist0=NULL,solver=solver)
EFFmlist[,j]=tmp$effvals[,1]
EFFmlistq[,j]=tmp$effvals[,3]
# debugoff
if(getbootpeers==TRUE){
PEERS=tmp$PEERS
names(PEERS) = c('dmu0','dmuzm','z')
PEERS=merge(PEERS,dmuset,all.x=TRUE)
PEERS$dmuz=ifelse(is.na(PEERS$dmuz),PEERS$dmu0,PEERS$dmuz)
PEERS$m = m
PEERS=PEERS[,c('m','dmu0','dmuz','z')]
PEERS=doBy::orderBy(~dmu0+dmuz,data=PEERS)
########################
PEERSq=tmp$PEERSq
names(PEERSq) = c('dmu0','dmuzm','z')
PEERSq=merge(PEERSq,dmuset,all.x=TRUE)
PEERSq$dmuz=ifelse(is.na(PEERSq$dmuz),PEERSq$dmu0,PEERSq$dmuz)
PEERSq$m = m
PEERSq=PEERSq[,c('m','dmu0','dmuz','z')]
PEERSq=doBy::orderBy(~dmu0+dmuz,data=PEERSq)
########################
peers[[j]] = PEERS
peersq[[j]] = PEERSq
#########################
} # end if(getbootpeers==TRUE)
} # end loop on j
if(getbootpeers==TRUE){
peers = dplyr::bind_rows(peers)
peersq = dplyr::bind_rows(peersq)
if(length(Bdrop)>=1){
EFFmlist[Bdrop,] = NA
pick1=match(peers$dmu0,Bdrop)
pick2=which(!is.na(pick1))
peers$dmuz[pick2] = NA
peers$z[pick2] = NA
peers=unique(peers)
} # end if(length(Bdrop)>=1)
if(length(Bdropq)>=1) {
EFFmlistq[Bdropq,] = NA
pick1=match(peersq$dmu0,Bdropq)
pick2=which(!is.na(pick1))
peersq$dmuz[pick2] = NA
peersq$z[pick2] = NA
peersq=unique(peersq)
} # end if(length(Bdropq)>=1)
} #end if(getbootpeers==TRUE)
peers=stats::na.omit(peers)
peersq=stats::na.omit(peersq)
#################################################################################
list(EFFmlist=EFFmlist,EFFmlistq=EFFmlistq,XM=XM,YM=YM,qout=qout,qoutS=qoutS,
X0=X0,Y0=Y0,DX0=DX0,DY0=DY0,mlist=mlist,mcells=mcells,peers=peers,peersq=peersq)
###############################################################################
} # end function qDEA_mlist
###############################################################################
###############################################################################
#' nCM_mpick: Select subsample size m from mlist
#' !!!!Intended to be called from qDEA function!!!!
#'
#' Pulls associated set of nboot bootstrapped values and computes bias corected
#' "s(m)" values.
#' ** Uses Simar and Wilson (2011) suggested procedure for selecting m **
#' @param stat DDEA or qDEA distance point estimate
#' @param boot nboot by mcells matrix of bootstrapped DDEA and qDEA distance
#' @param n Number of dmus in full reference set (X,Y)
#' @param mlist Vector of subsample sizes
#' @param beta Convergence rate. 0.5 => "root n" convergence
#' @param alpha Alpha level for Simar-Wilson (2011) subsample size selection procedure
#' @param CILag CILag level for Simar-Wilson (2011) subsample size selection procedure
#' @return mpick = selected subsample size,
#' @return S = nboot by mcell matrix of
#' s = stat(n) - (m/n)^beta(boot(m)-stat(n)
#' bias corrected values for each sample size m
#' @return s = vector (of length nboot) of 'bias.corected' s-values for
#' sample size 'mpick'
#' @return stat.bc = bias corrected point estimate of stat = mean(s)
###############################################################################
#' @export
###############################################################################
nCm_mpick=function(stat,boot,n,mlist,beta=0.5,alpha=0.05,CILag=1){
###############################################################################
boot=as.matrix(boot)
(mcells=length(mlist))
if(mcells!=ncol(boot)) stop('length of mlist != ncol(boot)')
###############################################################################
# Estimate confidence intervals
(m = mlist)
(nM = length(m))
S = 0 * boot
for (j in 1:nM) {
S[, j] = stat- ((m[j]/n)^beta) * (boot[, j] - stat)
}
UL = apply(S, 2, stats::quantile, probs = (1 - alpha/2),na.rm=TRUE)
LL = apply(S, 2, stats::quantile, probs = (alpha/2),na.rm=TRUE)
summary(UL)
summary(LL)
CL = stat - UL
CU = stat - LL
summary(CL)
summary(CU)
if(nM>2){
MCL = lagMat(CL, -CILag:CILag)
MCU = lagMat(CU, -CILag:CILag)
CLsd = apply(MCL, 1, stats::sd)
CUsd = apply(MCU, 1, stats::sd)
Csd = CLsd + CUsd
(mpick0 = which(Csd == suppressWarnings(min(Csd, na.rm=TRUE))))
} # end if(nM>2)
if(nM==1){
mpick0=1
} # end if(nM==1)
mpick=mpick0[1]
s=S[,mpick]
###############################
list(mpick=mpick,S=S,s=s,stat.bc=mean(s))
###############################################################################
} # end function nCm_mpick
###############################################################################
###############################################################################
#' DEAbuild: Builds DDEA LP object for use in qDEA_solve function
#'
#' @param X Reference dmu's = ndmu x number of inputs input matrix.
#' @param Y Reference dmu's = ndmu x number of outputs output matrix.
#' @param X0 Inputs for set of ndmu0 dmu's to be processed.
#' @param Y0 Outputs for set of ndmu0 dmu's to be processed.
#' @param DX0 Input directions for ndmu0 dmu's in X0 and Y0.
#' @param DY0 Output directions for ndmu0 dmu's in X0 and Y0.
#' @param dmu0 Row in (X0,Y0,DX0,DY0) to use as given dmu
#' @param RTS Returns to scale: 'CRS,'VRS,'DRS','IRS.
#' @param unbounded DEA obj restricted >= to unbounded. Default = -1E3
#' @param solver LP solver Default='highs'
#' @return Returns an LP list object containing the following elements:
#' @return LPsense = 'max' or 'min'
#' @return nnz = number of nonzero elements in 'A' matrix
#' @return nr = number or rows in 'A' matrix
#' @return nc = number of columns in 'A' matrix
#' @return obj = nc length vector of objective coefficients
#' @return ra = nonzero coefficients in 'A' matrix
#' @return ia = row indexes for non zero elements in 'A' matrix
#' @return ja = column indexes for non zero elements in 'A' matrix
#' @return SMM = Sparse matrix method: 'CMO','RMO','CRI',or 'RCI'
#' @return ZINDEX = 'zero indexing' (T) or 'one indexing' (F)
#' @return dirs = nr length vector of constraint signs ('<=','>=', or '=')
#' @return rhs = nr length vector of RHS coefficients
#' @return xlower = nc vector of lower bounds on 'x' choice variables
#' @return xupper = nc vector of upper bounds on 'x' choice variables
#' @return rlower = nr vector of Ax lower bounds i.e. rlower <= Ax
#' @return rupper = nr vector of Ax upper bounds i.e. Ax <= rupper
#' @return vartypes = nc vector of variable types: ('C','B','I') -- qDEA:rep('C',nc)
#' @return yxchngC = index of given dmu's output-input values locations in ra vector
#' (used to edit LP problem as we loop through DMU's in (X0,Y0,DX0,DY0)
#' @return dyxchngC = index of given dmu's direction values locations in ra vector
#' (used to edit LP problem as we loop through DMU's in (X0,Y0,DX0,DY0)
#' @return yxchngR = index of given dmu's obj restriction output-input values locations in ra vector
#' (used to edit LP problem as we loop through DMU's in (X0,Y0,DX0,DY0)
#' @return iyxchngC = row indexes associated with yxchngC cells
#' @return idyxchngC = row indexes associated with dyxchngC cells
#' @return iyxchngR = row indexes associated with yxchngR cells
#' @return DOBJ = matrix of obj values to change by dmu with DOBJ=cbind(-Y0,X0)
#' @return DYX = matrix of (dy,dx) values to change by dmu with DYX=cbind(DY0,DX0)
#' @return RTS = Returns to scale: 'CRS,'VRS,'DRS','IRS.
#' @return dmu0 = Row in (X0,Y0,DX0,DY0) to use as given dmu
###############################################################################
# debugging DEAbuild
###############################################################################
# debug on
# (dmu0=1); (RTS='VRS'); (SMM='CMO'); (ZINDEX=TRUE)
###############################################################################
#' @export
###############################################################################
DEAbuild=function(X,Y,X0,Y0,DX0,DY0,dmu0=1,RTS='CRS',
unbounded=-1E3,solver='highs'){
###############################################################################
# obj restriction
SMM='CMO'; ZINDEX=TRUE
(nrY=nrow(Y));(nrX=nrow(X))
(ncY=ncol(Y));(ncX=ncol(X))
if(nrY!=nrX) stop('nrY!=nrX')
DOBJ=cbind(-Y0,X0)
DYX=cbind(DY0,DX0)
y0=Y0[dmu0,];x0=X0[dmu0,];dy=DY0[dmu0,];dx=DX0[dmu0,]
(ncYX=ncY+ncX)
(nc=ncY+ncX)
(nr=nrY+2)
# Note: At this stage ja contains "columns counts"
# Include obj restriction
(M1=as.matrix(as.matrix(rbind(Y,dy,y0))))
(ra=as.vector(M1))
(ia=rep(1:(nrY+2),ncY))
# Note: At this stage ja contains "columns counts"
(ja=rep(nrY+2,ncY))
(M2=as.matrix(as.matrix(rbind(-X,dx,-x0))))
(ra=c(ra,as.vector(M2)))
(ia=c(ia,rep(1:(nrX+2),ncX)))
(ja=c(ja,rep((nrX+2),ncX)))
if(RTS!='CRS'&RTS!='crs'){
(ra=c(ra,rep(-1,nrY),-1))
(ia=c(ia,1:nrY,nrY+2))
# Note: At this stage ja contains "columns counts"
(ja=c(ja,nrY+1))
}
# compute CMO ja vector from counts
(ja=cumsum(c(1,ja)))
# change locations in constraint matrix
(yxchngC = cumsum(c(1,rep((nrY+2),(ncYX-1)))))
(dyxchngC = cumsum(c(nrY+1,rep((nrY+2),(ncYX-1)))))
(yxchngR = cumsum(c(nrY+2,rep((nrY+2),(ncYX-1)))))
###########################################################
LPsense='min'
(obj=c(-y0,x0))
(xlower=rep(0,nc))
(xupper=rep(1E8,nc))
(dirs=c(rep('<=',nrY),'=','<='))
(rhs=c(rep(0,nrY),1,-unbounded))
if(RTS!='CRS'&RTS!='crs'){
nc=nc+1
obj=c(obj,1)
if(RTS=='VRS'|RTS=='vrs'){
xlower=c(xlower,-1E8)
xupper=c(xupper,1E8)
}
if(RTS=='DRS'|RTS=='drs'){
xlower=c(xlower,0)
xupper=c(xupper,1E8)
}
if(RTS=='IRS'|RTS=='irs'){
xlower=c(xlower,0)
xupper=c(xupper,1E8)
}
} # end if(RTS!='CRS'&RTS!='crs')
rlower=rhs
rupper=rhs
rlower=ifelse(dirs=='<='|dirs=='<',-1E8,rlower)
rupper=ifelse(dirs=='<='|dirs=='<',rhs,rupper)
rlower=ifelse(dirs=='>='|dirs=='>',rhs,rlower)
rupper=ifelse(dirs=='>='|dirs=='>',1E8,rupper)
rlower=ifelse(dirs=='=='|dirs=='=',rhs,rlower)
rupper=ifelse(dirs=='=='|dirs=='=',rhs+1E-8,rupper)
rlower;rupper
###############################################################################
iyxchngC=ia[yxchngC]
idyxchngC=ia[dyxchngC]
iyxchngR=ia[yxchngR]
###############################################################################
if(ZINDEX==TRUE) {ia=ia-1;ja=ja-1}
vartypes=rep("C",nc)
###############################################################################
###############################################################################
LP=list(LPsense=LPsense,nnz=length(ra),nr=nr,nc=nc,obj=obj,ra=ra,ia=ia,ja=ja,
SMM=SMM,ZINDEX=ZINDEX,dirs=dirs,rhs=rhs,xlower=xlower,xupper=xupper,
rlower=rlower,rupper=rupper,vartypes=vartypes,
yxchngC=yxchngC,dyxchngC=dyxchngC,yxchngR=yxchngR,iyxchngC=iyxchngC,
idyxchngC=idyxchngC,iyxchngR=iyxchngR,DOBJ=DOBJ,DYX=DYX,RTS=RTS,dmu0=dmu0)
###############################################################################
LP=SM2SM(LP,SMM2='CRI',ZINDEX2=FALSE)
###############################################################################
LP
###############################################################################
} # end function DEAbuild
###############################################################################
###############################################################################
#' qDEAbuild: Builds qDEA LP object for use in qDEA_solve function
#'
#' @param X Reference dmu's = ndmu x number of inputs input matrix.
#' @param Y Reference dmu's = ndmu x number of outputs output matrix.
#' @param X0 Inputs for set of ndmu0 dmu's to be processed.
#' @param Y0 Outputs for set of ndmu0 dmu's to be processed.
#' @param DX0 Input directions for ndmu0 dmu's in X0 and Y0.
#' @param DY0 Output directions for ndmu0 dmu's in X0 and Y0.
#' @param qout Maximal proportion of dmu's allowed external to DEA hull.
#' @param dmu0 Row in (X0,Y0,DX0,DY0) to use as given dmu
#' @param RTS Returns to scale: 'CRS,'VRS,'DRS','IRS.
#' @param unbounded DEA obj restricted >= to unbounded. Default = -1E3
#' @param solver LP solver Default='highs'
#' @return Returns an LP list object containing the following elements:
#' @return LPsense = 'max' or 'min'
#' @return nnz = number of nonzero elements in 'A' matrix
#' @return nr = number or rows in 'A' matrix
#' @return nc = number of columns in 'A' matrix
#' @return obj = nc length vector of objective coefficients
#' @return ra = nonzero coefficients in 'A' matrix
#' @return ia = row indexes for non zero elements in 'A' matrix
#' @return ja = column indexes for non zero elements in 'A' matrix
#' @return SMM = Sparse matrix method: 'CMO','RMO','CRI',or 'RCI'
#' @return ZINDEX = 'zero indexing' (T) or 'one indexing' (F)
#' @return dirs = nr length vector of constaint signs ('<=','>=', or '=')
#' @return rhs = nr length vector of RHS coefficients
#' @return xlower = nc vector of lower bounds on 'x' choice variables
#' @return xupper = nc vector of upper bounds on 'x' choice variables
#' @return rlower = nr vector of Ax lower bounds i.e. rlower <= Ax
#' @return rupper = nr vector of Ax upper bounds i.e. Ax <= rupper
#' @return vartypes = nc vector of variable types: ('C','B','I') -- qDEA:rep('C',nc)
#' @return yxchngC = index of given dmu's output-input values locations in ra vector
#' (used to edit LP problem as we loop through DMU's in (X0,Y0,DX0,DY0)
#' @return dyxchngC = index of given dmu's direction values locations in ra vector
#' (used to edit LP problem as we loop through DMU's in (X0,Y0,DX0,DY0)
#' @return yxchngR = index of given dmu's output-input values locations in ra vector
#' (used to edit LP problem as we loop through DMU's in (X0,Y0,DX0,DY0)
#' @return iyxchngC = row indexes associated with yxchngC cells
#' @return idyxchngC = row indexes associated with dyxchngC cells
#' @return iyxchngR = row indexes associated with yxchngR cells
#' @return DOBJ = matrix of obj values to change by dmu with DOBJ=cbind(-Y0,X0)
#' @return DYX = matrix of (dy,dx) values to change by dmu with DYX=cbind(DY0,DX0)
#' @return qchng = index of 1/q value location in ra vector
#' (used to edit LP problem as we iterate qDEA
#' @return RTS = Returns to scale: 'CRS,'VRS,'DRS','IRS.
#' @return dmu0 = Row in (X0,Y0,DX0,DY0) to use as given dmu.
#' @return devstart = LP index for column of first qDDEA-S1 LPM "deviation" value
#' @return devend = LP index for column of last qDDEA-S1 LPM "deviation" value
#' @return tau0 = parameter used in iterative qDEA process
###############################################################################
# debug qDEAbuild
# (qout=0.10); (dmu0=1); (RTS='VRS'); (SMM='CMO'); (ZINDEX=TRUE)
###############################################################################
#' @export
###############################################################################
qDEAbuild=function(X,Y,X0,Y0,DX0,DY0,qout=0.10,dmu0=1,RTS='CRS',
unbounded=-1E3,solver='highs'){
###############################################################################
SMM='CMO'; ZINDEX=TRUE
eps = (1E-6)
DOBJ=cbind(-Y0,X0)
DYX=cbind(DY0,DX0)
(nrY=nrow(Y));(nrX=nrow(X))
(ncY=ncol(Y));(ncX=ncol(X))
if(nrY!=nrX) stop('nrY!=nrX')
(ncYX=ncY+ncX)
(nobs=nrow(Y))
(devstart=ncYX+2)
(devend=devstart+nobs-1)
y0=Y0[dmu0,];x0=X0[dmu0,];dy=DY0[dmu0,];dx=DX0[dmu0,]
q0 = qout-eps
q1 = (floor(qout*nrY)+0.95)/nrY - eps
tau0 = q1
(ncYX=ncY+ncX)
(nc=ncY+ncX+1+nrY+1)
(nr=nrY+4)
# Include obj restriction
M1=as.matrix(as.matrix(rbind(Y,dy,y0)))
ra=as.vector(M1)
ia=rep(1:(nrY+2),ncY)
# initially ja is just a count of nz elements in each column
ja=rep(nrY+2,ncY)
M2=as.matrix(rbind(-X,dx,-x0))
ra=c(ra,as.vector(M2))
ia=c(ia,rep(1:(nrX+2),ncX))
ja=c(ja,rep((nrY+2),ncX))
ra=c(ra,rep(1,nrY),-1)
ia=c(ia,1:nrY,nrY+4)
ja=c(ja,nrY+1)
ra=c(ra,rep(c(-1,1/nrY),nrY))
tmp=rep(1:nrY,each=2)
pick=seq(2,2*nrY,2)
tmp[pick]=nrY+3
ia=c(ia,tmp)
ja=c(ja,rep(2,nrY))
ra=c(ra,-1,1/tau0)
ia=c(ia,nrY+3,nrY+4)
ja=c(ja,2)
qchng=length(ra)
if(RTS!='CRS'& RTS!='crs'){
ra=c(ra,rep(-1,nrY),-1)
ia=c(ia,1:nrY,nrY+2)
ja=c(ja,nrY+1)
}
# compute CMO ja vector from counts
ja=cumsum(c(1,ja))
# change locations in constraint matrix
(yxchngC = cumsum(c(1,rep((nrY+2),(ncYX-1)))))
(dyxchngC = cumsum(c(nrY+1,rep((nrY+2),(ncYX-1)))))
(yxchngR = cumsum(c(nrY+2,rep((nrY+2),(ncYX-1)))))
###########################################################
LPsense='min'
(obj=c(-y0,x0,0,rep(0,nrY),0))
(xlower=rep(0,nc))
(xupper=rep(1E8,nc))
(dirs=c(rep('<=',nrY),'=','<=','<=','<='))
(rhs=c(rep(0,nrY),1,-unbounded,0,0))
if(RTS!='CRS'&RTS!='crs'){
nc=nc+1
obj=c(obj,1)
if(RTS=='VRS'|RTS=='vrs'){
xlower=c(xlower,-1E8)
xupper=c(xupper,1E8)
}
if(RTS=='DRS'|RTS=='drs'){
xlower=c(xlower,0)
xupper=c(xupper,1E8)
}
if(RTS=='IRS'|RTS=='irs'){
xlower=c(xlower,0)
xupper=c(xupper,1E8)
}
} # end if(RTS!='CRS'&RTS!='crs')
rlower=rhs
rupper=rhs
rlower=ifelse(dirs=='<='|dirs=='<',-1E8,rlower)
rupper=ifelse(dirs=='<='|dirs=='<',rhs,rupper)
rlower=ifelse(dirs=='>='|dirs=='>',rhs,rlower)
rupper=ifelse(dirs=='>='|dirs=='>',1E8,rupper)
rlower=ifelse(dirs=='=='|dirs=='=',rhs,rlower)
rupper=ifelse(dirs=='=='|dirs=='=',rhs+1E-8,rupper)
rlower;rupper
###############################################################################
iyxchngC=ia[yxchngC][1]
idyxchngC=ia[dyxchngC][1]
iyxchngR=ia[yxchngR][1]
###############################################################################
if(ZINDEX==TRUE) {ia=ia-1;ja=ja-1}
vartypes=rep("C",nc)
###############################################################################
###############################################################################
LP=list(LPsense=LPsense,nnz=length(ra),nr=nr,nc=nc,obj=obj,ra=ra,ia=ia,ja=ja,
SMM=SMM,ZINDEX=ZINDEX,dirs=dirs,rhs=rhs,xlower=xlower,xupper=xupper,
rlower=rlower,rupper=rupper,vartypes=vartypes,yxchngC=yxchngC,
dyxchngC=dyxchngC,yxchngR=yxchngR,iyxchngC=iyxchngC,idyxchngC=idyxchngC,
iyxchngR=iyxchngR,DOBJ=DOBJ,DYX=DYX,qchng=qchng,RTS=RTS,qout=qout,dmu0=dmu0,
devstart=devstart,devend=devend,qout=qout,tau0=tau0)
###############################################################################
LP=SM2SM(LP,SMM2='CRI',ZINDEX2=FALSE)
###############################################################################
LP
###############################################################################
} # end function qDEAbuild
###############################################################################
###############################################################################
#' my_lag: Lag a vector of values
#'
#' @param x A vector of values to be lagged
#' @param nlag Length of lag: negative value indicates an 'uplag'
#' @return = a lagged vector
###############################################################################
#' @export
###############################################################################
my_lag=function(x,nlag=1) {
nobx=length(x)
if(nlag<0) x=rev(x)
x0=rep(NA,abs(nlag))
x1=x[1:(nobx-abs(nlag))]
lagx=c(x0,x1)
if(nlag<0){
x=rev(x)
lagx=rev(lagx)
}
lagx
} # end function my_lag
###############################################################################
###############################################################################
#' lagMat: Create a matrix of lags
#'
#' @param x Vector to be 'lagged'
#' @param lags Number of lagged columns in resulting matrix
#' @param Lzero (F= do not include x[i] in row i) (T = include x[i] in row i)
#' @return = lagged matrix
###############################################################################
lagMat=function(x,lags=2,Lzero=FALSE) {
if(length(lags)==1){
if(lags[1]>0&Lzero==FALSE) lags=1:lags
if(lags[1]<0&Lzero==FALSE) lags=-1:lags
if(Lzero!=FALSE) lags=0:lags
} # end if(length(lags)==1)
nobx=length(x)
nlags=length(lags)
X=matrix(0,nobx,nlags)
for(j in 1:nlags){
X[,j]=my_lag(x,lags[j])
}
X
} # end function lagMat
###############################################################################
###############################################################################
#' my_seconds: Function to pull seconds from proc.time function
#'
#'@return = seconds from proc.time function
###############################################################################
#' @export
###############################################################################
my_seconds=function(){
time1=proc.time()
as.numeric(time1[3])
} # end function my_seconds
###############################################################################
###############################################################################
#' Write.LP
#'
#' Write a linear programming model to a csv file.
#'
#' This function converts a linear programming model to a dense
#' representation and writes it to a csv file that can be opened or solved
#' using spreadsheet based solvers.
#'
#' @param LP A linear programming object.
#' @param path Directory to write LP .csv file to. REQUIRED to use function
#' @param filename Character string specifying the csv file name.
#'
#' @return A csv file written to disk.
#'
#' @export
###############################################################################
Write.LP=function(LP,path=NULL,filename='LP.csv'){
###############################################################################
if (is.null(path)) {
stop("Please specify a directory path/folder to save the file in.")
}
###############################################################################
(fname=paste(path,'/',filename,sep=''))
###############################################################################
A = SM2A(LP)
(nc=ncol(A))
(nr=nrow(A))
if('cnames' %in% names(LP) == FALSE) LP$cnames=paste('C_',1:nc,sep='')
if('rnames' %in% names(LP) == FALSE) LP$rnames=paste('R_',1:nr,sep='')
if('rhs' %in% names(LP)==FALSE&'rupper' %in% names(LP)==FALSE){
stop('rhs and (rlower,rupper )cannot both be missing in LP object')
}
(LPsense=LP$LPsense)
obj=LP$obj
if('vartypes' %in% names(LP) == FALSE) LP$vartypes = rep('C',length(obj))
if('xlower' %in% names(LP) == FALSE) LP$xlower = rep(0,length(obj))
if('xupper' %in% names(LP) == FALSE) LP$xupper = rep(1E10,length(obj))
if('rlower' %in% names(LP) ==TRUE){if(is.na(LP$rlower[1])) LP$rlower=NULL}
if('rupper' %in% names(LP) ==TRUE){if(is.na(LP$rupper[1])) LP$rupper=NULL}
if('rlower' %in% names(LP) == FALSE) {
rlower = LP$rhs
rlower = ifelse(LP$dirs=='<='|LP$dirs=='<',-1E10,rlower)
rlower = ifelse(LP$dirs=='>='|LP$dirs=='>',LP$rhs,rlower)
LP$rlower = rlower
}
if('rupper' %in% names(LP) == FALSE) {
rupper = LP$rhs
rupper = ifelse(LP$dirs=='<='|LP$dirs=='<',LP$rhs,rupper)
rupper = ifelse(LP$dirs=='>='|LP$dirs=='>',1E10,rupper)
LP$rupper = rupper
}
solution=0*obj
if('solution' %in% names(LP) == TRUE) solution=LP$solution
cnames=LP$cnames
rnames=LP$rnames
vartypes=LP$vartypes
xlower=LP$xlower
xupper=LP$xupper
xlower[xlower<=-Inf] = -1E10
xupper[xupper>=Inf] = 1E10
dirs=LP$dirs
rhs=NA
if('rhs' %in% names(LP) == TRUE) rhs = LP$rhs
rlower=LP$rlower
rupper=LP$rupper
M=A
M=cbind(M,'|','lhsA','|')
(tmp=c(cnames,'|','LHS','|'))
if(exists('rhs')){
M=cbind(M,dirs,'|',rhs,'|')
tmp=c(tmp,'dirs','|','rhs','|')
}
if(exists('rlower')){
M=cbind(M,rlower,'|',rupper,'|')
tmp=c(tmp,'rlower','|','rupper','|')
}
M=cbind(rnames,M)
tmp=c('',tmp)
M=rbind(rep('========',ncol(M)),M)
M=rbind(tmp,M)
M=rbind(rep('========',ncol(M)),M)
M=rbind(M,rep('========',ncol(M)))
(tmp=c('OBJ',obj,'|','lhsOBJ','=',paste(LPsense)))
(tmp=c(tmp,rep('',ncol(M)-length(tmp))))
M=rbind(M,tmp)
M=rbind(M,rep('========',ncol(M)))
(tmp=c('nrows =',nr,'ncols =',nc))
M=rbind(c(tmp,rep('',ncol(M)-length(tmp))),M)
(tmp=c('xlower',xlower,'|','xlower'))
M=rbind(M,c(tmp,rep('',ncol(M)-length(tmp))))
(tmp=c('xupper',xupper,'|','xupper'))
M=rbind(M,c(tmp,rep('',ncol(M)-length(tmp))))
M=rbind(M,rep('========',ncol(M)))
(tmp=c('xvals',solution,'|','xvalues'))
M=rbind(M,c(tmp,rep('',ncol(M)-length(tmp))))
M=rbind(M,rep('========',ncol(M)))
M=rbind(rep('========',ncol(M)),M)
(tmp=c('VARTYPE',vartypes,'|','vartypes'))
M=rbind(M,c(tmp,rep('',ncol(M)-length(tmp))))
M=rbind(M,rep('========',ncol(M)))
rm(A,nr,nc)
utils::write.table(M,file=fname,row.names=FALSE,col.names=FALSE,sep=',')
###############################################################################
} # end function Write.LP
###############################################################################
#################################################################################
#' iter_delete: Function used for 'peeling' that supplements or supplants qDEA
#' 'slicing' procedure. Intended to be called from qDEA function
#'
#' @param LP0 qDEA Stage-2 LP object.
#' @param ndmus Number of dmus in reference data (X,Y).
#' @param nout0 Number of external points in current LP object's solution.
#' @param iterlim Maximal number of 'peeling' iterations.
#' @param BIGM Default Big M in RHS of qDEA stage 2 process.
#' @param unbounded LP reported as unbounded if obj<=unbounded.
#' @param eps iter_delete convergence test parameter
#' @param solver LP solver Default='highs'
#' @return obj2 = best objective level found
#' @return nout = number of external points
#' @return outlist = list of external dmus
#' @return iter = number of iterations completed
#' @return LP2 = LP object for best solution found
#' @return LP2sol = optimal LP solution
#################################################################################
# (ndmus=nrow(X))
# nout0=1
# iterlim=100
# BIGM=1E6
# unbounded=1E-3
# eps=1e-6
#################################################################################
#' @export
###############################################################################
iter_delete=function(LP0,ndmus,nout0,iterlim=max(250,nout0),BIGM=1E6,
unbounded=1E-3,eps=1e-6,solver='highs'){
###############################################################################
LP2=LP0
LP2sol=LPSOLVER(LP2,solver=solver)
(obj2=LP2sol$objval)
nobs=ndmus
lhs2=LP2sol$lhs[1:nobs]
dual2=LP2sol$dual[1:nobs]
(pick20=which(lhs2 > eps))
(nout=length(pick20))
(nleft=nout0-nout)
(pick21=which(dual2 < (-eps)))
pick22=sort(c(pick20,pick21))
(npick22=length(pick22))
go=1
iter=1
while(iter <= iterlim & npick22<=nout0){
LP2=LP0
LP2$rhs[pick22]=BIGM
LP2$rupper[pick22]=BIGM
LP2sol=LPSOLVER(LP2,solver=solver)
lhs2=LP2sol$lhs[1:nobs]
dual2=LP2sol$dual[1:nobs]
(obj2=LP2sol$objval)
(pick20=which(lhs2 > eps))
(nout=length(pick20))
(nleft=nout0-nout)
(pick21=which(dual2 < (-eps)))
(npick21=length(pick21))
(pick22=sort(c(pick20,pick21)))
(npick22=length(pick22))
(iter=iter+1)
} # end while(iter <= iterlim & npick22<=nleft){
##############################################################
iter
nout
nout0
(npick21=length(pick21))
while(iter <= iterlim & nout<nout0){
LPlist=list('NOT RUN')
objhat=0
pick20
j=1
for(j in 1:npick21){
(pick22=sort(c(pick20,pick21[j])))
LP2=LP0
LP2$rhs[pick22]=BIGM
LP2$rupper[pick22]=BIGM
LP2sol=LPSOLVER(LP2,solver=solver)
(objhat[j]=LP2sol$objval)
LPlist[[j]]=list(LP2=LP2,LP2sol=LP2sol)
} # end loop for(jon in 1:npick21)
objhat
(jpick=which(objhat==min(objhat))[1])
LP2=LPlist[[jpick]]$LP2
LP2sol=LPlist[[jpick]]$LP2sol
lhs2=LP2sol$lhs[1:nobs]
dual2=LP2sol$dual[1:nobs]
(obj2=LP2sol$objval)
(pick20=which(lhs2 > eps))
(nout=length(pick20))
(nleft=nout0-nout)
(pick21=which(dual2 < (-eps)))
(npick21=length(pick21))
(pick22=sort(c(pick20,pick21)))
(npick22=length(pick22))
(iter=iter+1)
} #end while(iter <= iterlim & nout<nout0){
list(obj2=obj2,nout=nout,outlist=pick20,iter=iter,LP2=LP2,LP2sol=LP2sol)
} # end function iter_delete
###############################################################################
###############################################################################
#require(data.table);setDTthreads(1)
###############################################################################
#' A2SM: Convert a matrix A to sparse matrix object ASM
#'
#' @param A The matrix A.
#' @param SMM Sparse matrix method with SMM='A',CRI','RCI','CMO', or 'RMO'
#' @param ZINDEX T = Use zero indexing or F = Use one indexing.
#' @param eps Value to use in non-zero test.
#' @param NA_flag Number to uses as NA flag
#'
#' @return A list object containing the matrix components:
#' @return nnz = the number of non-zero elements in ra.
#' @return nr = the number of rows in matrix A.
#' @return nc = the number of columns in matrix A.
#' @return ia = the row index.
#' @return ja = the column index.
#' @return ra = the non-zero coefficients in A
#' @return rnames = matrix row names -- may be ""
#' @return cnames = matrix column names -- may be ""
#' @return SMM = the sparse matrix type.
#' @return ZINDEX with T = Use zero indexing or F = Use one indexing.
#' @examples
#' (A = matrix(c(1,0,0,2,0,3,1,0,0,5,0,6),3,4,byrow=TRUE))
#' (ASM1 = A2SM(A,SMM='CMO',ZINDEX=TRUE))
#' (ASM2=SM2SM(ASM1,SMM2='CRI',ZINDEX2=FALSE))
#' SM2A(ASM1)
#' SM2A(ASM2)
#'
#' (A = matrix(c(1,0,0,2,0,3,1,0,0,5,0,6),3,4,byrow=TRUE)); A[2,3]=NA; A
#' (ASM1 = A2SM(A,SMM='CMO',ZINDEX=TRUE)); A; SM2A(ASM1)
#' (ASM2 = SM2SM(ASM1,SMM2='CRI',ZINDEX2=TRUE)); A ; SM2A(ASM2)
#'
#' #CMO documentation example
#' nr=5
#' nc=8
#' ra=c(3.0,5.6,1.0,2.0,1.1,1.0,-2.0,2.8,-1.0,1.0,1.0,-1.2,-1.0,1.9)
#' ia=c(0,4,0,1,1,2,0,3,0,4,2,3,0,4)
#' ja=c(0,2,4,6,8,10,11,12,14)
#' SMM='CMO'
#' ZINDEX=TRUE
#' ASM=list(nr=nr,nc=nc,ra=ra,ia=ia,ja=ja,SMM=SMM,ZINDEX=ZINDEX)
#' (A=SM2A(ASM))
###############################################################################
#' @export
###############################################################################
A2SM = function(A, SMM = 'CRI',ZINDEX=FALSE,
eps = .Machine$double.eps,NA_flag=(-1E12)){
##############################################################################
A[is.na(A)] = NA_flag
##############################################################################
(nr = nrow(A))
(nc = ncol(A))
if(is.null(rownames(A))) rownames(A)=rep('',nr)
if(is.null(colnames(A))) colnames(A)=rep('',nc)
(rnames=rownames(A))
(cnames=colnames(A))
##############################################################################
# initially uses SMM RCI
tmp=list(NULL); eps=1E-9
j=1
for(j in 1:nc){
a=as.vector(A[,j])
ipick=which(abs(a)>eps)
ra=a[ipick]
ia=ipick
ja=rep(j,length(ipick))
as.data.frame(cbind(ia,ja,ra))
tmp[[j]]=as.data.frame(cbind(ia,ja,ra))
} # end loop for(j in 1:nc)
tmp=dplyr::bind_rows(tmp)
(nnz=nrow(tmp))
ASM=list(nnz=nnz,nr=nr,nc=nc,ia=tmp$ia,ja=tmp$ja,ra=tmp$ra,
rnames=rnames,cnames=cnames,SMM='CRI',ZINDEX=FALSE)
rm(tmp)
# convert to user specified sparse matrix format
SM2SM(ASM,SMM2=SMM,ZINDEX2=ZINDEX)
} # end function A2SM
###############################################################################
###############################################################################
#' SM2A: Convert a sparse matrix object ASM into a matrix A
#'
#' ASM is a sparse matrix object containing:
#' nr = number of rows, nc = number of columns, ra = nonzero coeff,
#' ia = row indices, ja = col indices,
#' rnames = row names,cnames = column names # Note these may be missing or ""
#' SMM = space matrix method, and ZINDEX = use 0-indexing
#' @param ASM A sparse matrix object
#' @param NA_flag Number to uses as NA flag
#'
#' @return The matrix A.
#' @examples
#' (A = matrix(c(1,0,0,2,0,3,1,0,0,5,0,6),3,4,byrow=TRUE))
#' (ASM1 = A2SM(A,SMM='CMO',ZINDEX=TRUE))
#' (ASM2=SM2SM(ASM1,SMM2='CRI',ZINDEX2=FALSE))
#' SM2A(ASM1)
#' SM2A(ASM2)
#'
#' (A = matrix(c(1,0,0,2,0,3,1,0,0,5,0,6),3,4,byrow=TRUE)); A[2,3]=NA; A
#' (ASM1 = A2SM(A,SMM='CMO',ZINDEX=TRUE)); A; SM2A(ASM1)
#' (ASM2 = SM2SM(ASM1,SMM2='CRI',ZINDEX2=TRUE)); A ; SM2A(ASM2)
#'
#' #CMO documentation example
#' nr=5
#' nc=8
#' ra=c(3.0,5.6,1.0,2.0,1.1,1.0,-2.0,2.8,-1.0,1.0,1.0,-1.2,-1.0,1.9)
#' ia=c(0,4,0,1,1,2,0,3,0,4,2,3,0,4)
#' ja=c(0,2,4,6,8,10,11,12,14)
#' SMM='CMO'
#' ZINDEX=TRUE
#' ASM=list(nr=nr,nc=nc,ra=ra,ia=ia,ja=ja,SMM=SMM,ZINDEX=ZINDEX)
#' (A=SM2A(ASM))
###############################################################################
#' @export
###############################################################################
SM2A = function(ASM,NA_flag=(-1E12)){
###############################################################################
ASM$ra[is.na(ASM$ra)] = NA_flag
##
nr=ASM$nr; nc=ASM$nc
if(!exists('rnames',where=ASM)) ASM$rnames=rep('',nr)
if(!exists('cnames',where=ASM)) ASM$cnames=rep('',nc)
##
(SMM0=ASM$SMM)
if(SMM0=='A') A = ASM$A
if(SMM0!='A'){
ASM2=SM2SM(ASM,SMM2='CRI',ZINDEX2=FALSE)
nr=ASM2$nr;nc=ASM2$nc;
ia=ASM2$ia;ja=ASM2$ja;ra=ASM2$ra
(ra=as.numeric(ra))
(ia=as.integer(ia))
(ja=as.integer(ja))
(nr=as.integer(nr))
(nc=as.integer(nc))
A=matrix(0,nr,nc)
A[cbind(ia,ja)]=ra
} # end if(SMM0!='A')
##
A[A <= NA_flag] = NA
##
if(!is.null(ASM$rnames)&length(ASM$rnames)==nrow(A)) rownames(A)=ASM$rnames
if(!is.null(ASM$cnames)&length(ASM$cnames)==ncol(A)) colnames(A)=ASM$cnames
###############################################################################
A
###############################################################################
} # end function SM2A
###############################################################################
###############################################################################
#' SM2SM: Convert a sparse matrix object ASM into a different sparse matrix form.
#'
#' ASM is a sparse matrix object containing:
#' nr = number of rows, nc = number of columns, ra = nonzero coeff,
#' ia = row indices, ja = col indices, SMM = space matrix method,
#' ZINDEX=use 0-indexing
#' @param ASM A sparse matrix object
#' @param SMM2 sparse matrix method with SMM2 ='A','CRI','RCI','CMO',or'RMO'
#' @param ZINDEX2 T = Use zero indexing or F = Use one indexing.
#' @param NA_flag Number to uses as NA flag
#'
#' @return A list object containing the following sparse matrix components:
#' @return nnz = the number of non-zero elements in ra.
#' @return nr = the number of rows in matrix.
#' @return nc = the number of columns in matrix.
#' @return ia = the revised row index.
#' @return ja = the revised column index.
#' @return rnames = row names -- may be missing or ""
#' @return cnames = column names -- may be missing or ""
#' @return ra = the revised non-zero coefficients in A
#' @return SMM = the revised sparse matrix type.
#' @return ZINDEX = TRUE or F for the revised sparse form.
#' @examples
#' (A = matrix(c(1,0,0,2,0,3,1,0,0,5,0,6),3,4,byrow=TRUE))
#' (ASM1 = A2SM(A,SMM='CMO',ZINDEX=TRUE))
#' (ASM2 = SM2SM(ASM1,SMM2='CRI',ZINDEX2=FALSE))
#' SM2A(ASM1)
#' SM2A(ASM2)
#'
#' (A = matrix(c(1,0,0,2,0,3,1,0,0,5,0,6),3,4,byrow=TRUE)); A[2,3]=NA; A
#' (ASM1 = A2SM(A,SMM='CMO',ZINDEX=TRUE)); SM2A(ASM1)
#' (ASM2 = SM2SM(ASM1,SMM2='CRI',ZINDEX2=TRUE)); A ; SM2A(ASM2)
#'
#' #CMO documentation example
#' nr=5
#' nc=8
#' ra=c(3.0,5.6,1.0,2.0,1.1,1.0,-2.0,2.8,-1.0,1.0,1.0,-1.2,-1.0,1.9)
#' ia=c(0,4,0,1,1,2,0,3,0,4,2,3,0,4)
#' ja=c(0,2,4,6,8,10,11,12,14)
#' SMM='CMO'
#' ZINDEX=TRUE
#' ASM=list(nr=nr,nc=nc,ra=ra,ia=ia,ja=ja,SMM=SMM,ZINDEX=ZINDEX)
#' SM2A(ASM)
###############################################################################
#' @export
###############################################################################
SM2SM = function(ASM, SMM2='CRI', ZINDEX2=FALSE, NA_flag=(-1E12)){
###############################################################################
ASM$ra[is.na(ASM$ra)] = NA_flag
##
nnz=ASM$nnz;nr=ASM$nr; nc=ASM$nc
if(!exists('rnames',where=ASM)) ASM$rnames=rep('',nr)
if(!exists('cnames',where=ASM)) ASM$cnames=rep('',nc)
##
(rnames=ASM$rnames)
(cnames=ASM$cnames)
## original sparse form
(SMM0=ASM$SMM)
## if SMM0 indicates a matrix convert ASM to SMM='CRI' sparse form
if(SMM0=='A') {
(tmp=A2SM(ASM$A, SMM = 'CRI', ZINDEX=FALSE))
ASM=merge_lists(tmp,ASM)
} # end if(SMM0=='A')
## Pull sparse matrix elements from ASM
(nr=ASM$nr);(nc=ASM$nc);(ia=ASM$ia);(ja=ASM$ja);(ra=ASM$ra)
(SMM=ASM$SMM);(ZINDEX=ASM$ZINDEX)
## convert "zero" to "one" indexing if needed
if(ZINDEX==TRUE){
ia=ia+1
ja=ja+1
}# end if (ZINDEX == TRUE)
ia;ja
## convert CMO ja vector or RMO ia vector to complete index vectors
if(SMM=='CMO') ja=rep(1:nc,diff(ja))
if(SMM=='RMO') ia=rep(1:nr,diff(ia))
tmp=as.data.frame(cbind(ia,ja,ra))
## convert ASM to ASM2 form
if(SMM2=='CRI'|SMM2=='CMO'){
# resort to be sure in "complete and contiguous CRI form
(tmp=tmp[with(tmp,order(ja,ia)),])
# pull resorted data
ra=tmp$ra; ia=tmp$ia; ja=tmp$ja
if(SMM2=='CMO'){
(ja=c(match(1:nc,ja),nnz+1))
} # end if(SMM2=='CMO')
} # end if(SMM2=='CRI'|SMM2=='CMO')
if(SMM2=='RCI'|SMM2=='RMO'){
# resort to be sure in "complete and contiguous RCI form
(tmp=tmp[with(tmp,order(ia,ja)),])
# pull resorted data
ra=tmp$ra; ia=tmp$ia; ja=tmp$ja
if(SMM2=='RMO'){
ia=c(match(1:nr,ia),nnz+1)
} # end if(SMM2=='RMO')
} # end if(SMM2=='RCI'|SMM2=='RMO)
(ra=as.numeric(ra))
(ia=as.integer(ia))
(ja=as.integer(ja))
(nr=as.integer(nr))
(nc=as.integer(nc))
if(ZINDEX2 == TRUE) {ia = ia-1; ja = ja-1}
##
##
ASM2=ASM # inherits all objects in ASM
# put revised data into ASM2
ASM2$nnz=nnz;ASM2$nr=nr;ASM2$nc=nc;ASM2$ia=ia;ASM2$ja=ja;ASM2$ra=ra
ASM2$SMM=SMM2;ASM2$ZINDEX=ZINDEX2
if(SMM2=='A'){
ASM2=ASM # inherits all objects in original ASM
ASM2$A = SM2A(ASM)
ASM2$SMM = SMM2
}
#(ASM2 = list(nnz=nnz,nr = nr,nc = nc,ia = ia2,ja = ja2,ra = ra2,
# rnames=rnames,cnames=cnames,SMM = SMM2,ZINDEX = ZINDEX2))
###############################################################################
ASM2
###############################################################################
} # end function SM2SM
###############################################################################
# ###############################################################################
# # Examples:
# ###############################################################################
# # create matrix with missing data
# (A = matrix(c(1,0,0,2,0,3,1,0,0,5,0,6),3,4,byrow=TRUE)); A[2,3]=NA; A
# # no names for in matrix A
# c(rownames(A),colnames(A))
# # add names to A matrix
# rownames(A)=paste('r',1:nrow(A),sep='');colnames(A)=paste('c',1:ncol(A),sep='')
# A
# # Convert to sparse matrix with NA's replaced by NA_flag
# ASM1=A2SM(A,NA_flag=-1000)
# # Note below: SM2A puts " NA's" back in matrix if original matrix had them
# # or if user indicates levels in matrix that designate missing data
# SM2A(ASM1);SM2A(ASM1,NA_flag=-1E4);SM2A(ASM1,NA_flag=-1E3)
# ##############################################################################
#
# ##############################################################################
# A=t(matrix(1:15,5,3))
# A[cbind(c(3,1,2,3,1,3),c(1,2,3,3,4,5))]=0
# A
# #
# # (A=matrix(c(1,0,2,0,0,3,0,4,0,5,6,0),3,4,byrow=TRUE))
# # (A=matrix(c(1,0,0,2,0,3,4,0,0,5,0,6),3,4,byrow=TRUE))
# # (A=matrix(c(1,0,0,2,0,3,0,0,0,5,0,6),3,4,byrow=TRUE))
# # (A=matrix(c(2,3,1,4,1,2,3,4,2),3,3,byrow=TRUE))
# #
# A2SM(A,SMM='CMO',ZINDEX=TRUE)
#
# SMMlist=c('RCI','CRI','RMO','CMO')
# ZINDEXlist=c(FALSE,TRUE)
#
# SMM1 = 'RCI' # 'RCI', 'CRI', 'RMO', 'CMO'
# ZINDEX1=FALSE # T or F
# SMM2='CRI' # 'RCI', 'CRI', 'RMO', 'CMO'
# ZINDEX2=FALSE # T or F
#
# for(SMM1 in SMMlist){
# for(ZINDEX1 in ZINDEXlist) {
# for(SMM2 in SMMlist) {
# for(ZINDEX2 in ZINDEXlist) {
# print(paste(SMM1,ZINDEX1,SMM2,ZINDEX2))
# (ASM=A2SM(A,SMM=SMM1,ZINDEX=ZINDEX1))
# (ASM1=SM2SM(ASM,SMM2=SMM2,ZINDEX2=ZINDEX2))
# (ASM2=A2SM(A,SMM=SMM2,ZINDEX=ZINDEX2))
# print(summary(unlist(ASM1)==unlist(ASM2)))
# } # end for(ZINDEX2 in ZINDEXlist)
# } # end for(SMM2 in SMMlist)
# } # end for(ZINDEX1 in ZINDEXlist)
# } # end for(SMM1 in SMMlist)
#
#
# summary(SM2A(ASM1)==SM2A(ASM2))
# #############################################################
# (A=matrix(c(1,0,2,0,0,3,0,4,0,5,6,0),3,4,byrow=TRUE))
#
# (ASM=A2SM(A,SMM='A'))
# ASM2=SM2SM(ASM,SMM2=SMM2,ZINDEX2=ZINDEX2)
# (A2=SM2A(ASM2))
# summary(as.vector(A)==as.vector(A2))
#
# for(SMM2 in SMMlist) {
# for(ZINDEX2 in ZINDEXlist) {
# print(paste(SMM2,ZINDEX2))
# (ASM=A2SM(A=A,SMM='A'))
# (ASM2=SM2SM(ASM,SMM2=SMM2,ZINDEX2=ZINDEX2))
# (A2=SM2A(ASM2))
# print(summary(as.vector(A)==as.vector(A2)))
# } # end for(ZINDEX2 in ZINDEXlist)
# } # end for(SMM2 in SMMlist)
# ###########################################################
# (A=matrix(c(1,0,2,0,0,3,0,4,0,5,6,0),3,4,byrow=TRUE))
# SMM2='CRI' # 'RCI', 'CRI', 'RMO', 'CMO'
# ZINDEX2=FALSE # T or F
#
# (ASM=A2SM(A,SMM='A'))
# ASM2=SM2SM(ASM,SMM2=SMM2,ZINDEX2=ZINDEX2)
# (A2=SM2A(ASM2))
# summary(as.vector(A)==as.vector(A2))
# (A3=SM2A(ASM))
# summary(as.vector(A)==as.vector(A3))
# ##########################################################
# # CMO documentation example
# nr=5
# nc=8
# ra=c(3.0,5.6,1.0,2.0,1.1,1.0,-2.0,2.8,-1.0,1.0,1.0,-1.2,-1.0,1.9)
# ia=c(0,4,0,1,1,2,0,3,0,4,2,3,0,4)
# ja=c(0,2,4,6,8,10,11,12,14)
# SMM='CMO'
# ZINDEX=TRUE
# ASM=list(nr=nr,nc=nc,ra=ra,ia=ia,ja=ja,SMM=SMM,ZINDEX=ZINDEX)
# (A=SM2A(ASM))
# #############################################################################
###############################################################################
#' cbindSM: "column bind" two sparse matrices
#'
#' @param SM1 First sparse matrix object
#' @param SM2 Second sparse matrix object
#' @param SMM Sparse matrix method with 'CRI','RCI','CMO',or'RMO'
#' @param ZINDEX T = Use zero indexing or F = Use one indexing.
#'
#' @return A 'column bound" sparse matrix object with components"
#' @return nnz = the number of non-zero elements in ra.
#' @return nr = the number of rows in matrix A.
#' @return nc = the number of columns in matrix A.
#' @return ia = the row index.
#' @return ja = the column index.
#' @return ra = the non-zero coefficients in A
#' @return rnames = matrix row names -- may be ''
#' @return cnames = matrix column names -- may be ''
#' @return SMM = the sparse matrix type.
#' @return ZINDEX with T = Use zero indexing or F = Use one indexing.
###############################################################################
#' @export
###############################################################################
cbindSM=function(SM1,SM2,SMM='CRI',ZINDEX=FALSE){
nr1=SM1$nr; nc1=SM1$nc; nr2=SM2$nr; nc2=SM2$nc
if(nr1!=nr2) stop('STOP: row numbers not compatible')
tmp1=SM2SM(SM1,SMM2='CRI',ZINDEX2=FALSE)
tmp2=SM2SM(SM2,SMM2='CRI',ZINDEX2=FALSE)
nr=nr1; nc = nc1+nc2
ia=c(tmp1$ia,tmp2$ia)
ja=c(tmp1$ja,tmp2$ja+nc1)
ra=c(tmp1$ra,tmp2$ra)
tmp1$nnz=length(ra);tmp1$nr=nr;tmp1$nc=nc
tmp1$ia=ia;tmp1$ja=ja;tmp1$ra=ra;tmp1$SMM='CRI';tmp1$ZINDEX=FALSE
if(length(tmp1$cnames)>1) {
if(length(tmp2$cnames)>1) tmp1$cnames=c(tmp1$cnames,tmp2$cnames)
if(length(tmp2$cnames)<=1) tmp1$cnames=c(tmp1$cnames,rep('',nc2))
} # end
SM2SM(tmp1,SMM2=SMM,ZINDEX2=ZINDEX)
} # end function cbindSPM
###############################################################################
###############################################################################
#' rbindSM: "row bind" two sparse matrices
#'
#' @param SM1 First sparse matrix object
#' @param SM2 Second sparse matrix object
#' @param SMM Sparse matrix method with 'CRI','RCI','CMO',or'RMO'
#' @param ZINDEX T = Use zero indexing or F = Use one indexing.
#'
#' @return A 'row bound" sparse matrix object with components"
#' @return nnz = the number of non-zero elements in ra.
#' @return nr = the number of rows in matrix A.
#' @return nc = the number of columns in matrix A.
#' @return ia = the row index.
#' @return ja = the column index.
#' @return ra = the non-zero coefficients in A
#' @return rnames = matrix row names -- may be ''
#' @return cnames = matrix column names -- may be ''
#' @return SMM = the sparse matrix type.
#' @return ZINDEX with T = Use zero indexing or F = Use one indexing.
###############################################################################
#' @export
###############################################################################
rbindSM=function(SM1,SM2,SMM='CRI',ZINDEX=FALSE){
nr1=SM1$nr; nc1=SM1$nc; nr2=SM2$nr; nc2=SM2$nc
if(nc1!=nc2) stop('STOP: column number not compatible')
tmp1=SM2SM(SM1,SMM2='RCI',ZINDEX2=FALSE)
tmp2=SM2SM(SM2,SMM2='RCI',ZINDEX2=FALSE)
nr=nr1+nr2; nc = nc1
ia=c(tmp1$ia,tmp2$ia+nr1)
ja=c(tmp1$ja,tmp2$ja)
ra=c(tmp1$ra,tmp2$ra)
tmp1$nnz=length(ra);tmp1$nr=nr;tmp1$nc=nc
tmp1$ia=ia;tmp1$ja=ja;tmp1$ra=ra;tmp1$SMM='RCI';tmp1$ZINDEX=FALSE
if(length(tmp1$rnames)>1) {
if(length(tmp2$rnames)>1) tmp1$rnames=c(tmp1$rnames,tmp2$rnames)
if(length(tmp2$rnames)<=1) tmp1$rnames=c(tmp1$rnames,rep('',nr2))
} # end
SM2SM(tmp1,SMM2=SMM,ZINDEX2=ZINDEX)
} # end function rbindSPM
###############################################################################
#
#
#
################################################################################
# # Examples of cbindSM and rbindSM
# set.seed(1001)
# nr=7;nc=7
# (A=matrix(round(runif(nr*nc),3),nr,nc))
# A[A<0.5]=0
# colnames(A)=paste('c',1:nc,sep='')
# rownames(A)=paste('r',1:nr,sep='')
# A
#
# (A11=A[1:5,1:5]); (A12=A[1:5,6:7])
# (A21=A[6:7,1:5]); (A22=A[6:7,6:7])
#
#
# (SM11=A2SM(A11))
# (SM12=A2SM(A12))
# (SM21=A2SM(A21))
# (SM22=A2SM(A22))
#
# cbind(A11,A12)
# SM2A(cbindSM(SM11,SM12))
#
#
# rbind(A11,A21)
# SM2A(rbindSM(SM11,SM21))
#
# A
# SM2A(cbindSM(rbindSM(SM11,SM21),rbindSM(SM12,SM22)))
###############################################################################
# A2SM(A,SMM='CMO')
# cbindSM(rbindSM(SM11,SM21),rbindSM(SM12,SM22),SMM='CMO')
###############################################################################
###############################################################################
#' merge_lists: Merges two list objects.
#' This function appends any objects in list2 and not in list 1
#' to list1 with priority given to list 1 components.
#'
#' @param list1 A list object
#' @param list2 A list object
#' @return list3 A merged object
#' @examples
#' L1=list(a=1:3,b=4:6,c=7:9,e=10:12)
#' L2=list(b=13:15,d=16:18,e=19:21,f=22:24)
#' (L3=merge_lists(L1,L2))
#' (L4=merge_lists(L2,L1))
###############################################################################
#' @export
###############################################################################
merge_lists = function(list1, list2){
#####################################################################
# This function appends any objects in list2 and not in list 1
# to list1 with priority given to list 1 components.
#####################################################################
(names1 = names(list1))
(names2 = names(list2))
(list3 = list1)
(n1=length(list1))
(pick2 = which(is.na(match(names2, names1))))
if (length(pick2) > 0) {
for(j in 1:length(pick2)) list3 = c(list3,list2[pick2[j]])
}
sort_list(list3)
} # end function merge_lists
###############################################################################
###############################################################################
#' sort_list: Sorts objects in list by object name
#'
#' @param list1 A list object
#' @return list2 A resorted list object
###############################################################################
#' @export
###############################################################################
sort_list=function(list1){
(n=length(list1))
df1=data.frame(Lnames=names(list1))
df1$cell=1:n; df1$cell=as.numeric(df1$cell)
df1=df1[order(df1$Lnames),]
list2=list1[df1$cell[1]]
for(j in 2:n) list2=c(list2,list1[df1$cell[j]])
list2
} # end function sort_list
################################################################################
###############################################################################
#' LPSOLVER Function to call specified solver to sparse LP problem
#'
#' @param LP An LP object
#' @param solver LP solver Default='highs' Note: Versions for 'clp','highs','glpk' or 'rglpk' are available from author
#' @return A list containing following components:
#' @return status = status of solver
#' @return objval = objective value
#' @return solution = lp primal solution
#' @return dual = lp dual solution
#' @return rcost = lp reduced cost
#' @return lhs = lp lhs associated with primal solution
#' @return TIME = total seconds to execute LP_* function
#' @return SMtime = seconds to transform form of sparse LP (if needed)
#' @return proctime = seconds required by specified solver to solve LP
#' @return LP = LP object
###############################################################################
#' @export
###############################################################################
LPSOLVER=function(LP,solver='highs'){
LP_highs(LP)
} # end function LPSOLVER
###############################################################################
###############################################################################
#' LP_highs Function to solve sparse LP problem using highs package
#'
#' @param LP An LP object
#' @return A list containing following components:
#' @return status = status of highs solver
#' @return objval = objective value
#' @return solution = lp primal solution
#' @return dual = lp dual solution
#' @return rcost = lp reduced cost
#' @return lhs = lp lhs associated with primal solution
#' @return TIME = total seconds to execute LP_highs function
#' @return SMtime = seconds to transform form of sparse LP (if needed)
#' @return proctime = seconds required by highs to solve LP
#' @return LP = LP object (transformed to SMM='CRI', ZINDEX=FALSE if needed)
###############################################################################
#' @export
###############################################################################
LP_highs=function(LP){
###############################################################################
#require(highs);require(Matrix)
###############################################################################
if(!requireNamespace("highs", quietly=TRUE)|!requireNamespace("Matrix", quietly=TRUE)) {
stop('The packages highs and Matrix must be installed to use the option solver = highs')
}
###############################################################################
# Note: This code is set up to allow
# the user to measure the time used in
# different procedures. The following
#lines pull initial times from the system
###############################################################################
time000=my_seconds(); time00=my_seconds()
###############################################################################
#Convert LP matrix format (if necessary)
###############################################################################
# print(paste('SMM',LP$SMM))
###############################################################################
if(LP$SMM!='CRI'|LP$ZINDEX!=FALSE) LP=SM2SM(LP,SMM2='CRI',ZINDEX2=FALSE)
# Compute sparse matrix conversion time used
SMtime=my_seconds()-time00
#######################################
if('vartypes' %in% names(LP) == FALSE) LP$vartypes = rep('C',length(LP$obj))
#######################################
# Construct lower and upper bounds on x
# If xbounds are not specified in LP
# we assume 0 <= x <= Inf
#######################################
if('xlower' %in% names(LP) == FALSE) LP$xlower = rep(0,length(LP$obj))
if('xupper' %in% names(LP) == FALSE) LP$xupper = rep(Inf,length(LP$obj))
#######################################
# Construct lower and upper bounds on Ax
# If rbounds are not specified in LP
# we construct them from the LP list's
# 'rhs' and 'dirs' vectors
#######################################
if('rlower' %in% names(LP) == FALSE) {
rlower = LP$rhs
rlower = ifelse(LP$dirs=='<='|LP$dirs=='<',-Inf,rlower)
rlower = ifelse(LP$dirs=='>='|LP$dirs=='>',LP$rhs,rlower)
LP$rlower = rlower
}
if('rupper' %in% names(LP) == FALSE) {
rupper = LP$rhs
rupper = ifelse(LP$dirs=='<='|LP$dirs=='<',LP$rhs,rupper)
rupper = ifelse(LP$dirs=='>='|LP$dirs=='>',Inf,rupper)
LP$rupper = rupper
}
(LPsense = ifelse(LP$LPsense!='max'|LP$LPsense!='MAX',FALSE,TRUE))
# # highs blows up if abs of any sparse matrix ra values are <1e-9
# LP$ra[abs(LP$ra)<2e-9] = 2e-9
ASM = Matrix::sparseMatrix(LP$ia,LP$ja,x=LP$ra)
time0=my_seconds()
tmp = highs::highs_solve(L = LP$obj, lower = LP$xlower, upper = LP$xupper,
A = ASM, lhs = LP$rlower, rhs = LP$rupper,types=LP$vartypes,
maximum=LPsense,control=list(threads=1))
proctime=my_seconds()-time0
# retrieve the results
(status=tmp$status)
(objval=tmp$objective_value)
(solution=tmp$primal_solution)
rcost=tmp$solver_msg$col_dual
lhs=tmp$solver_msg$row_value
dual=tmp$solver_msg$row_dual
TIME=my_seconds()-time000
###############################################################################
list(status=status,objval=objval,solution=solution,dual=dual,rcost=rcost,
lhs=lhs,TIME=TIME,SMtime=SMtime,proctime=proctime,LP=LP)
###############################################################################
} # end function LP_highs
###############################################################################
Try the qDEA package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.