*! version 1.2.0 18JUN2011 capture program drop dea program define dea, rclass version 10.0 /** Terms Description: * ----------------------------------------------------------------------------- * RTS: Return To Scale * CRS: Constant Return to Scale * VRS: Variant Return to Scale * IRS: Increasing Return to Scale * DRS: Decreasing Retrn to Scale * ----------------------------------------------------------------------------- */ // syntax checking and validation----------------------------------------------- // rts - return to scale, ort - orientation // input varlist = output varlist // example: // dea Employee Area = Sales Profits, rts(CRS) ort(IN) tol1(1e-14) // dea Employee Area = Sales Profits, sav // ----------------------------------------------------------------------------- // returns 1 if the first nonblank character of local macro `0' is a comma, // or if `0' is empty. if replay() { dis as err "ivars and ovars required." exit 198 } // get and check invarnames gettoken word 0 : 0, parse(" =:,") while `"`word'"' != ":" & `"`word'"' != "=" { if `"`word'"' == "," | `"`word'"'=="" { error 198 } local invarnames `invarnames' `word' gettoken word 0 : 0, parse(" :,") } unab invarnames : `invarnames' #del ; syntax varlist(min=1) [if] [in] [using/] [, RTS(string) // ignore case sensitive,[{CRS|CCR}|{BCC|VRS}|DRS|IRS] ORT(string) // ignore case sensitive,[{IN|INPUT}|{OUT|OUTPUT}] STAGE(integer 2) // stage 1 or 2 TOL1(real 1e-14) // entering or leaving value tolerance TOL2(real 1e-8) // B inverse tolerance: 2.22e-12 MINSUBScript // minimal subscript pivot TRACE // Whether or not to do the log SAVing(string) // log file name REPLACE // Whether or not to replace the log file ]; #del cr local num_invar : word count `invarnames' local i 1 while (`i'<=`num_invar') { local invarn : word `i' of `invarnames' local junk : subinstr local invarnames "`invarn'" "", /// word all count(local j) if `j' > 1 { di as error /// "cannot specify the same input variable more than once" error 498 } local i = `i' + 1 } // default model - CRS(Constant Return to Scale) if ("`rts'" == "") local rts = "CRS" else { local rts = upper("`rts'") if ("`rts'" == "CCR") local rts = "CRS" else if ("`rts'" == "BCC") local rts = "VRS" else if (~("`rts'" == "CRS" | "`rts'" == "VRS" | /// "`rts'" == "DRS" | "`rts'" == "IRS")) { di as err "option rts allow for case-insensitive " _c di as err "CRS (eq CCR) or VRS (eq BCC) or DRS or IRS or nothing." exit 198 } } // default orientation - Input Oriented if ("`ort'" == "") local ort = "IN" else { local ort = upper("`ort'") if ("`ort'" == "I" | "`ort'" == "IN" | "`ort'" == "INPUT") { local ort = "IN" } else if ("`ort'" == "O" | "`ort'" == "OUT" | "`ort'" == "OUTPUT") { local ort = "OUT" } else { di as err "option ort allowed for case-insensitive " _c di as err "(i|in|input|o|out|output) or nothing." exit 198 } } // default stage - 1 if (~("`stage'" == "1" | "`stage'" == "2")) { dis as err "option stage allows 1 or 2 or nothing exclusively." exit 198 } if ("`using'" != "") use "`using'", clear if (~(`c(N)' > 0 & `c(k)' > 0)) { dis as err "dataset required!" exit 198 } // end of syntax checking and validation --------------------------------------- set more off capture log close dealog log using "dea.log", replace text name(dealog) preserve if ("`if'" != "" | "`in'" != "") { qui keep `in' `if' // filtering : keep in range [if exp] } deanormal, ivars(`invarnames') ovars(`varlist') /// rts(`rts') ort(`ort') stage(`stage') /// tol1(`tol1') tol2(`tol2') `minsubscript' /// `trace' saving(`saving') `replace' return add restore, preserve log close dealog end ******************************************************************************** * DEA Normal - Data Envelopment Analysis Normal ******************************************************************************** program define deanormal, rclass #del ; syntax , IVARS(string) OVARS(string) RTS(string) ORT(string) [ STAGE(integer 2) TOL1(real 1e-3) TOL2(real 1e-12) MINSUBScript TRACE SAVing(string) REPLACE ]; #del cr preserve // di _n as input /// // "options: RTS(`rts') ORT(`ort') STAGE(`stage') SAVing(`saving')" // di as input "Input Data:" // list // ------------------------------------------------------------------------- tempname dmuIn dmuOut frameMat deamainrslt dearslt vrsfrontier crslambda mkDmuMat `ivars', dmumat(`dmuIn') sprefix("i") mkDmuMat `ovars', dmumat(`dmuOut') sprefix("o") local dmuCount = colsof(`dmuIn') if ("`minsubscript'" == "") local minsubscript = "nominsubscript" if ("`ort'" == "OUT") local minrank = 1 else local minrank = 0 mata: _mkframemat("`frameMat'", "`dmuIn'", "`dmuOut'", /// "`rts'", "`ort'", `stage') deamain `dmuIn' `dmuOut' `frameMat' `rts' `ort' `stage' /// `tol1' `tol2' `minsubscript' `trace' matrix `deamainrslt' = r(deamainrslt) mata: _dmurank("`deamainrslt'", /// `=rowsof(`dmuIn')', `=rowsof(`dmuOut')', `minrank', `tol2') matrix `dearslt' = r(rank), `deamainrslt' // use mata function '_setup_dearslt_names' because the maximum string // variable length needs to be kept under the 244 for all the time mata: _setup_dearslt_names("`dearslt'", "`dmuIn'", "`dmuOut'") if ("`rts'" == "VRS") { // caution : join order! (CRS -> VRS -> DRS) // 1. VRS TE(VRS Technical Efficiency) matrix `deamainrslt' = r(deamainrslt) matrix `vrsfrontier' = `deamainrslt'[1...,1] // 2. CRS TE(CRS Technical Efficiency) local stageForVRS = 2 mata: _mkframemat("`frameMat'", "`dmuIn'", "`dmuOut'", /// "CRS", "`ort'", `stageForVRS') deamain `dmuIn' `dmuOut' `frameMat' "CRS" `ort' `stageForVRS' /// `tol1' `tol2' `minsubscript' `trace' matrix `deamainrslt' = r(deamainrslt) matrix `vrsfrontier' = `deamainrslt'[1...,1], `vrsfrontier' matrix `crslambda' = `deamainrslt'[1...,2..(`dmuCount' + 1)] // 3. DRS TE(DRS Technical Efficiency) local stageForVRS = 1 mata: _mkframemat("`frameMat'", "`dmuIn'", "`dmuOut'", /// "DRS", "`ort'", `stageForVRS') deamain `dmuIn' `dmuOut' `frameMat' "DRS" `ort' `stageForVRS' /// `tol1' `tol2' `minsubscript' `trace' matrix `deamainrslt' = r(deamainrslt) matrix `vrsfrontier' = `vrsfrontier', `deamainrslt'[1...,1] // VRS mata: _roundmat("`vrsfrontier'", 1e-14) // round off // matrix `vrsfrontier' = `vrsfrontier', J(rowsof(`vrsfrontier'), 2, 0) matrix rownames `vrsfrontier' = `: colfullnames `dmuIn'' matrix colnames `vrsfrontier' = CRS_TE VRS_TE DRS_TE SCALE RTS if ("`trace'" == "trace") matrix list `crslambda' forvalues i = 1/`=rowsof(`vrsfrontier')' { matrix `vrsfrontier'[`i',4] = /* */ float(`vrsfrontier'[`i',1]/`vrsfrontier'[`i',2]) /******************************************************************* * if CRS(CCR) == VRS(BCC) then CRS * else * if the sum of ¥ë equals to 1, then mark CRS * if sum of ¥ë is greater than 1, then mark DRS * if sum of ¥ë is less than 1, then mark IRS ******************************************************************/ // (-1:drs, 0:crs, 1:irs) if (`vrsfrontier'[`i',1] == `vrsfrontier'[`i',2]) { matrix `vrsfrontier'[`i',5] = 0 } else { local sumlambda = 0 forvalues j = 1/`dmuCount' { if (`crslambda'[`i', `j'] < .) { local sumlambda = `sumlambda' + `crslambda'[`i', `j'] } } if (`sumlambda' < 1) { matrix `vrsfrontier'[`i',5] = 1 // irs mark } else { // if (`sumlambda' >= 1) { matrix `vrsfrontier'[`i',5] = -1 // drs mark } } } } // ------------------------------------------------------------------------- // REPORT // ------------------------------------------------------------------------- di as result "" di as input "options: RTS(`rts') ORT(`ort') STAGE(`stage')" di as result "`rts'-`ort'PUT Oriented DEA Efficiency Results:" matrix list `dearslt', noblank nohalf noheader f(%9.6g) if ("`saving'" != "") { // if the saved file exists and doesn't replace, make the backup file. if ("`replace'" == "") { local dotpos = strpos("`saving'",".") if (`dotpos' > 0) { mata: _file_exists("`saving'") } else { mata: _file_exists("`saving'.dta") } if r(fileexists) { local curdt = subinstr("`c(current_date)'", " ", "", .) + /* */ subinstr("`c(current_time)'", ":", "", .) if (`dotpos' > 0) { #del ; local savefn = substr("`saving'", 1, `dotpos' - 1) + "_bak_`curdt'" + substr("`saving'",`dotpos', .); #del cr qui copy "`saving'" "`savefn'", replace } else { local savefn = "`saving'_bak_`curdt'" + ".dta" qui copy "`saving'.dta" "`savefn'", replace } } } if ("`rts'" != "VRS") { restore, preserve svmat `dearslt', names(eqcol) capture { renpfix _ renpfix ref ref_ renpfix islack is_ renpfix oslack os_ } capture save `saving', replace // di as result "" // di as result "DEA Result file:" // list } } return matrix dearslt = `dearslt' if ("`rts'" == "VRS") { if ("`trace'" == "trace") { di _n(2) as result "CRS lambda(¥ë)" matrix list `crslambda' , noblank nohalf noheader f(%9.6f) } // to verify: if you don't have to verify, you can comment this sentence. di _n(2) as result "VRS Frontier(-1:drs, 0:crs, 1:irs)" matrix list `vrsfrontier', noblank nohalf noheader f(%9.6f) di _n(2) as result "VRS Frontier:" restore, preserve svmat `vrsfrontier', names(col) rename RTS RTSNUM qui generate RTS = "drs" if RTSNUM == -1 qui replace RTS = "-" if RTSNUM == 0 qui replace RTS = "irs" if RTSNUM == 1 drop DRS_TE RTSNUM format CRS_TE VRS_TE SCALE %9.6f list if ("`saving'"!="") { capture save `saving', replace } } restore, preserve end ******************************************************************************** * DEA Main - Data Envelopment Analysis Main ******************************************************************************** program define deamain, rclass args dmuIn dmuOut frameMat rts ort stage tol1 tol2 minsubscript trace tempname efficientVec deamainrslt // stage step 1. if ("`trace'" == "trace") { di _n(2) as txt "RTS(`rts') ORT(`ort') 1st stage." } mata: _dealp("`frameMat'", "`dmuIn'", "`dmuOut'", "`rts'", "`ort'", /// 1, `tol1', `tol2', "`minsubscript'", "`efficientVec'", "`trace'") matrix `deamainrslt' = r(dealprslt) // stage step 2. if ("`stage'" == "2") { if ("`trace'" == "trace") { di _n(2) as txt "RTS(`rts') ORT(`ort') 2nd stage." } matrix `efficientVec' = `deamainrslt'[1...,1] mata: _dealp("`frameMat'", "`dmuIn'", "`dmuOut'", "`rts'", "`ort'", /// 2, `tol1', `tol2', "`minsubscript'", "`efficientVec'", "`trace'") matrix `deamainrslt' = r(dealprslt) } return matrix deamainrslt = `deamainrslt' end ******************************************************************************** * Data Import and Conversion ******************************************************************************** // Make DMU Matrix ------------------------------------------------------------- program define mkDmuMat #del ; syntax varlist(numeric) [if] [in], DMUmat(name) [ SPREfix(string) ]; #del cr qui ds // variable not found error if ("`varlist'" == "") { di as err "variable not found" exit 111 } // make matrix mkmat `varlist' `if' `in', matrix(`dmumat') rownames(dmu) matrix roweq `dmumat' = "dmu" matrix coleq `dmumat' = `=lower("`sprefix'") + "slack"' matrix `dmumat' = `dmumat'' end // Start of the MATA Definition Area ------------------------------------------- version 10 mata: mata set matastrict on /** * Declare the variable's boundary condition structure. */ struct BoundCond { real scalar val, lower, upper, free } /** * Declare the LP's parameter structure. */ struct LpParamStruct { string scalar rts // return to scale(CRS|VRS|IRS|DRS) real scalar isin // if 1 then 'in', other then 'out' real scalar stagestep // stage step. 1 or 2 real scalar dmus // number of dmus real scalar slacks // number of slacks real scalar artificials // number of artificials real scalar tol1 // tolerance 1 real scalar tol2 // tolerance 2 real scalar isminsubscript // whether min subscript or not. string scalar trace // whether trace or not. } /** * Declare the LP's result structure. */ struct LpResultStruct { real scalar xVal // value of objective. real matrix XB // basic feasible solution. real scalar rc // return code(zero means success) string scalar rmsg // return message } /** * make frame matrix and set matrix value at the param frameMat * rts - return to scale, ort - orientation */ function _mkframemat( string scalar frameMat, string scalar dmuIn, string scalar dmuOut, string scalar rts, string scalar ort, real scalar stage ) { real matrix F, DI, DO real scalar row, col, sig, isin real scalar dmus, slackins, slackouts, slacks real scalar frows, fcols DI = st_matrix(dmuIn) DO = st_matrix(dmuOut) if (cols(DI) != cols(DO)) { _error(3200, "in and out count of dmu is not match!") } // basic value setting for artificial variabels isin = (ort == "IN"); sig = (isin ? -1 : 1) dmus = cols(DI) // or cols(DO), because cols(DI) == cols(DO) slackins = rows(DI); slackouts = rows(DO) if (rts == "CRS") { slacks = slackins + slackouts // target coefficient\slackins\slackouts frows = 1 + slackins + slackouts // target coefficient,theta,dmus,slackins,slackouts,rhs fcols = 1 + 1 + dmus + slackins + slackouts + 1 } else if (rts == "VRS") { slacks = slackins + slackouts // target coefficient\slackins\slackouts\sum of lamda frows = 1 + slackins + slackouts + 1 // target coefficient,theta,dmus,slackins,slackouts,rhs fcols = 1 + 1 + dmus + slackins + slackouts + 1 } else if (rts == "IRS") { slacks = slackins + slackouts + 1 // target coefficient\slackins\slackouts\sum of lamda frows = 1 + slackins + slackouts + 1 // target coefficient,theta,dmus,slackins,slackouts,sum of lamda,rhs fcols = 1 + 1 + dmus + slackins + slackouts + 1 + 1 } else if (rts == "DRS") { slacks = slackins + slackouts + 1 // target coefficient\slackins\slackouts\sum of lamda frows = 1 + slackins + slackouts + 1 // target coefficient,theta,dmus,slackins,slackouts,sum of lamda,rhs fcols = 1 + 1 + dmus + slackins + slackouts + 1 + 1 } else { _error(3498, "invalid rts optoin.") } // make frame matrix for CRS(CCR) F = J(frows, fcols, 0) F[1, 1] = 1 replacesubmat(F, 2, 3, sig * DI) replacesubmat(F, 2 + slackins, 3, -sig * DO) replacesubmat(F, 2, 3 + dmus, sig * I(slacks)) // adjustment if (rts == "VRS") { replacesubmat(F, frows, 3, J(1, dmus, 1)) F[frows,fcols] = 1 } else if (rts == "IRS") { replacesubmat(F, frows, 3, J(1, dmus, 1)) F[frows,2 + dmus + slacks] = -1 F[frows,fcols] = 1 } else if (rts == "DRS") { replacesubmat(F, frows, 3, J(1, dmus, 1)) F[frows,2 + dmus + slacks] = 1 F[frows,fcols] = 1 } // return result st_matrix(frameMat, F) } /** * DEA Loop - Data Envelopment Analysis Loop for DMUs */ function _dealp ( string scalar frameMat, string scalar dmuIn, string scalar dmuOut, string scalar rts, string scalar ort, real scalar stagestep, real scalar tol1, real scalar tol2, string scalar minsubscript, string scalar efficientVec, string scalar trace ) { real matrix F, DI, DO, M, VARS, LPRSLT, DEALPRSLT, ARTIF real scalar dmus, slackins, slackouts, slacks, artificials, artificialrow real scalar frows, fcols, isin, isminsubscript real colvector effvec, slackidx, skipdmu string scalar tracename struct BoundCond matrix boundF, boundM struct LpParamStruct scalar param F = st_matrix(frameMat) DI = st_matrix(dmuIn) DO = st_matrix(dmuOut) if (cols(DI) != cols(DO)) { _error(3200, "in and out count of dmu is not match!") } if (!(rts == "CRS" || rts == "VRS" || rts == "IRS" || rts == "DRS")) { _error(3498, "rts must be one of CRS, VRS, IRS, DRS") } // basic value setting for artificial variabels isin = (ort == "IN") isminsubscript = (minsubscript != "nominsubscript") frows = rows(F); fcols = cols(F) dmus = cols(DI) // or cols(DO), because cols(DI) == cols(DO) slackins = rows(DI); slackouts = rows(DO) tracename = rts + "-" + (isin ? "IN" : "OUT") tracename = tracename + "-" + (stagestep == 1 ? "SI" : "SII") // ------------------------------------------------------------------------- // define number of slacks by rts if (rts == "CRS" || rts == "VRS") slacks = slackins + slackouts else if (rts == "IRS" || rts == "DRS") slacks = slackins + slackouts + 1 // define number of artificials by rts, ort, stage if (rts == "CRS" || rts == "DRS") { if (stagestep == 1) { if (isin) { artificials = slackins+slackouts; artificialrow = 2; } else artificials = 0 } else { artificials = slackouts; artificialrow = 2+slackins; } } else if (rts == "VRS" || rts == "IRS") { if (stagestep == 1) { if (isin) { artificials = slackins+slackouts+1; artificialrow = 2; } else { artificials = 1; artificialrow = frows //== 2+slackins+slackouts } } else { artificials = slackouts+1; artificialrow = 2+slackins } } if (artificials > 0) { ARTIF = J(1, artificials, 1) \ J(frows-1, artificials, 0) replacesubmat(ARTIF, artificialrow, 1, I(artificials)) F = F[,1..fcols-1], ARTIF, F[,fcols] frows = rows(F); fcols = cols(F) // revise frows, fcols } // ------------------------------------------------------------------------- // put artificial for slackins on inactive status if (stagestep == 2) { effvec = st_matrix(efficientVec) skipdmu = (effvec :== .) if (isin) { replacesubmat(F, 2, 3, -F[2..1+slackins,3::2+dmus+slackins]) } else { replacesubmat(F, 2+slackins, 3, -F[2+slackins..1+slackins+slackouts,3::2+dmus+slacks]) } } else skipdmu = J(1, dmus, 0) // ------------------------------------------------------------------------- boundF = J(1, fcols, BoundCond()); // set the boundary for the efficiency variable(theta, eta): // -INFINITE <= efficiency <= INFINITE boundF[1,2].val = 0; boundF[1,2].lower = 0; boundF[1,2].upper = . // set boundary value for the weight variable(lamda, mu): // 0 <= weight <= INFINITE for (i=3; i 0) LPRSLT[1, VARS[1,j+1]] =lpresult.XB[j, 1] } if (param.stagestep == 1 && LPRSLT[1, 1] <= 0) { LPRSLT[1, 1] = lpresult.xVal } } if (param.trace == "trace") { msg = sprintf("[DMUi=%g]%s-FINAL", dmui, tracename); printf("\n%s: original VARS.\n", msg); orgVARS printf("\n%s: VARS.\n", msg); VARS printf("\n%s: XB.\n", msg); lpresult.XB printf("\n%s: LPRSLT.\n", msg); LPRSLT } return(LPRSLT) } /** * Declare the LP's tableau structure. */ struct LpTableauStruct { pointer(real matrix) scalar CB, CNj pointer(real matrix) scalar B, Nj, b pointer(real matrix) scalar Bi,CBBi, rawXB, XB } /** * return 0: sucess * return 1: B inverse error * return 2: XB has negative value. */ real scalar function decompsition(real matrix M, struct BoundCond matrix boundM, real scalar mrows, real scalar mcols, real scalar slacks, struct LpTableauStruct scalar tbl, struct LpParamStruct scalar param ) { real matrix CB, CNj real matrix B, Nj, b real matrix Bi,CBBi, rawXB, XB, BiNjXj real scalar Njcols, result // set the tableau. tbl.CB = &CB; tbl.CNj = &CNj; tbl.B = &B; tbl.Nj = &Nj; tbl.b = &b tbl.Bi = &Bi; tbl.CBBi = &CBBi tbl.rawXB = &rawXB; tbl.XB = &XB; CB = M[1,2::1+slacks]; CNj = M[1,2+slacks::mcols-1] B = M[2..mrows,2::1+slacks]; Nj = M[2..mrows,2+slacks::mcols-1] b = M[2..mrows,mcols] Bi = lusolve(B, I(rows(B)), 1e-14) if (any(Bi :== .)) { // B is singular matrix. return (result = 1); // Bi = svsolve(B, I(rows(B)), 1e-14) // if (any(Bi :== .)) return (result = 1); } CBBi = CB*Bi // BFS(basic feasible solution) Njcols = cols(Nj); BiNjXj = J(rows(Nj), 1, 0) for (j=1; j<=Njcols; j++) { BiNjXj = BiNjXj :+ (Bi*Nj[.,1]*(boundM[1,1+slacks+j].val)) } rawXB = Bi*b - BiNjXj XB = edittozerotol(rawXB, param.tol2) // BFS(basic feasible solution) // BFS(basic feasible solution) must be nonnegative. if (any(XB :< 0)) return (result = 2); return (result = 0); // sucess } struct LpResultStruct function lp ( real matrix M, struct BoundCond matrix boundM, real matrix VARS, real scalar dmui, real scalar phase, string scalar tracename, struct LpParamStruct scalar param ) { real matrix B, CB, Bi, b, XB, rawXB, BiNjXj, CBBi, Nj, CNj, Aj, alpha real matrix T, TH1, TH2, valT, lowerT, upperT, LVi, V real scalar mrows, mcols, enteringVar, leavingVar, calci, maxiter real scalar existArtificial, xVal, alphaVal real scalar minYn, tcols, tempVal, minVal, maxVal struct BoundCond matrix boundT struct LpResultStruct scalar lpresult struct LpTableauStruct scalar tbl real colvector enterings, leavings real scalar enteringi, leavingi real scalar slacks, isin, tol1, tol2, isminsubscript string scalar trace // ------------------------------------------------------------------------- slacks = rows(M) - 1 // number of basic feasible solution isin = param.isin tol1 = param.tol1; tol2 = param.tol2 isminsubscript = param.isminsubscript trace = param.trace minYn = 0 if (param.stagestep == 1) minYn = (phase == 1) ? 1 : isin else minYn = (phase == 1) // ------------------------------------------------------------------------- mrows = rows(M); mcols = cols(M) LVi = J(slacks, 1, .) // leaving variable index matrix. // ------------------------------------------------------------------------- if (trace == "trace") { displayas("txt"); msg = "initial tableau in the LP." printf("\n[DMUi=%g]%s: %s\n", dmui, tracename, msg); M } lpresult.rc = 0; lpresult.rmsg = "" existArtificial = (phase == 2 && any(VARS[,2::1+slacks] :< 0)); maxiter = st_numscalar("c(maxiter)") for (calci=1 ; calci<=maxiter ; calci++) { // prevent infinite loop if (trace == "trace") { printf("\n[DMUi=%g]%s-LOOP[%g] Start...\n", dmui, tracename, calci) } B = M[2..mrows,2::1+slacks]; CB = M[1,2::1+slacks] Nj = M[2..mrows,2+slacks::mcols-1]; CNj = M[1,2+slacks::mcols-1] b = M[2..mrows,mcols] Bi = lusolve(B, I(rows(B)), 1e-14) if (any(Bi :== .)) { // B is singular matrix. Bi = svsolve(B, I(rows(B)), 1e-14) if (any(Bi :== .)) { lpresult.rc = 3498; lpresult.rmsg = sprintf("%s[DMUi=%g][LOOP=%g]%s", "No Solution(BFS's inverse is not exist):", dmui, calci, tracename) break; /* display("B:");B display("rank(B) : det(B)"); rank(B), det(B) _error(3498, "No Solution(BFS's inverse is not exist):" + "[DMUi=" + strofreal(dmui) + "]" + "[LOOP=" + strofreal(calci) + "]" + tracename) */ } } CBBi = CB*Bi // BFS(basic feasible solution) Njcols = cols(Nj); BiNjXj = J(rows(Nj), 1, 0) for (j=1; j<=Njcols; j++) { BiNjXj = BiNjXj :+ (Bi*Nj[.,1]*(boundM[1,1+slacks+j].val)) } rawXB = Bi*b - BiNjXj XB = edittozerotol(rawXB, tol2) // BFS(basic feasible solution) if (any(XB :== .)) { lpresult.rc = 3498; lpresult.rmsg = sprintf("%s[DMUi=%g][LOOP=%g]%s", "No Solution(XB contains missing value):", dmui, calci, tracename) break; } // BFS(basic feasible solution) must be nonnegative. boundT = boundM[1,2+slacks::mcols-1] tcols = cols(boundT); valT = J(1,tcols,.) for(j=1; j<=tcols; j++) { valT[1,j] = boundT[1,j].val } xVal = edittozerotol((CB*rawXB + CNj*valT'), tol2) // objective funtion value T = VARS[1,2::1+rows(XB)] // rows(XB) equal to number of slacks. if (trace == "trace") { printf("\n[DMUi=%g]%s-LOOP[%g]: CBBi * b = %g\n", dmui, tracename, calci, xVal); display("Entered index(if value is negative, that's artificial):"); T; display("XB = Bi*b - BiNjXj:"); rawXB; } // --------------------------------------------------------------------- // loop terminated condition. // --------------------------------------------------------------------- if (phase == 1) { // objective function value is zero if (xVal == 0 ) { // if all artificals are out or remaining artificals are at zero, // stop and go phaseII T = (T :< 0) // artificial remain or not? if (allof(T, 0) || allof(select(XB, T'), 0)) break; // If remaining artificails are not zero, No Solution. lpresult.rc = 3498 lpresult.rmsg = sprintf("%s[DMUi=%g][LOOP=%g]%s", "No Solution(Remaining artificails are not zero):", dmui, calci, tracename) break; } } // --------------------------------------------------------------------- // Select entering variable. // --------------------------------------------------------------------- enteringVar = 0; tempVal = 0; minVal = 0; maxVal = 0; boundi = 0 Njcols = cols(Nj); T = J(1, Njcols, .) for (j=1; j<=Njcols; j++) { tempVal = CBBi * Nj[,j] - CNj[1,j] // if (abs(tempVal) < tol1) continue; boundi = 1 + slacks + j if (boundM[1,boundi].val == boundM[1,boundi].lower) { // lower bound T[1,j] = tempVal } else { // upper bound T[1,j] = -tempVal } } // end of for if (!minYn) { // maximization. T = T :/ (T :< 0) if (!allof(T, .)) { minindex(T, 1, enterings, w); enteringi = 1 enteringVar = enterings[enteringi] evj = 1+slacks+enteringVar } } else { // minimization. T = T :/ (T :> 0) if (!allof(T, .)) { maxindex(T, 1, enterings, w); enteringi = 1 enteringVar = enterings[enteringi] evj = 1+slacks+enteringVar } } // No more candidate for entering variable. if (enteringVar == 0) { if (trace == "trace") { printf("\n[DMUi=%g]%s-LOOP[%g]:", dmui, tracename, calci) printf("No more candidate for entering variable.\n:(CB*Bi*Nj)-Cj\n");T } if (phase == 1) { lpresult.rc = 3498; lpresult.rmsg = sprintf("%s[DMUi=%g][LOOP=%g]%s", "No Solution(No more candidate for entering variable):", dmui, calci, tracename) } break } if (trace == "trace") { displayas("txt"); msg = "Select entering variable." printf("\n[DMUi=%g]%s-LOOP[%g]: %s(%g:%g)\n:(CB*Bi*Nj)-Cj\n", dmui, tracename, calci, msg, enteringVar, T[enteringVar]); T } // --------------------------------------------------------------------- // Select leaving variable. // --------------------------------------------------------------------- leavingVar = 0 Aj = Nj[,enteringVar] alpha = edittozerotol(Bi*Aj, tol1) if (existArtificial) { T = VARS[1,2::1+slacks]; tcols = cols(T) for (j=1; j<=tcols; j++) { if (T[1,j] < 0 && alpha[j,1] != 0) { leavingVar = j; lvj = 1+leavingVar break; } } } if (leavingVar == 0) { boundT = boundM[1,2::1+slacks] tcols = cols(boundT) lowerT = upperT = J(1,tcols,.) for(j=1; j<=tcols; j++) { lowerT[1,j] = boundT[1,j].lower upperT[1,j] = boundT[1,j].upper } // XB=Bi*b leavingCase = 0 if (boundM[1,evj].val == boundM[1,evj].lower) { minVal = .; // 1. alpha's positive min value TH1 = boundM[1,evj].lower :+ ((rawXB :- lowerT') :/ (alpha :* (alpha :> 0))) // TH1 = edittozerotol(TH1, tol1) if (any(TH1 :< minVal)) { minindex(TH1, 1, mi, w) leavingVar = mi[1]; lvj = 1+leavingVar if (phase == 1 && w[1,2] >= 2 && VARS[1,lvj] > 0) { // if phase 1 and same min ratio test result, // artificial variable must leave first. tcols = w[1,2] for (j=2; j<=tcols; j++) { if (VARS[1,mi[j]+1] < 0) { leavingVar = mi[j]; lvj = 1+leavingVar break; } } } minVal = TH1[leavingVar,1]; leavingCase = 1 } // 2. alpha's negative min value TH2 = boundM[1,evj].lower :+ ((rawXB :- upperT') :/ (alpha :* (alpha :< 0))) // TH2 = edittozerotol(TH2, tol1) if (any(TH2 :< minVal)) { minindex(TH2, 1, mi, w) leavingVar = mi[1]; lvj = 1+leavingVar minVal = TH1[leavingVar,1]; leavingCase = 2 } // 3. get the enteringVar's upper value if (boundM[1,evj].upper < minVal) { minVal = boundM[1,evj].upper; leavingCase = 3 } if (trace == "trace") { displayas("txt"); msg = "Select leaving variable.[MinVal]" printf("\n[DMUi=%g]%s-LOOP[%g]: %s(%g:%g)\n", dmui, tracename, calci, msg, leavingVar, minVal) display("XB | alpha:(XB=Bi*b, alpha=Bi*Aj):");rawXB,alpha display("[MinVal]enteringVar's upper | theta1 | theta2") printf("\n[boundM[1,%g].upper:%g][leavingCase:%g]\n", evj, boundM[1,evj].upper, leavingCase); TH1,TH2 } if (leavingCase == 1) { boundM[1,lvj].val = boundM[1,lvj].lower } else if (leavingCase == 2) { boundM[1,lvj].val = boundM[1,lvj].upper } else { // if (leavingCase == 3) boundM[1,evj].val = boundM[1,evj].upper; continue; } } else { // if (boundM[1,evj].val == boundM[1,evj].upper) maxVal = 0; // 1. alpha's positive min value TH1 = boundM[1,evj].upper :+ ((rawXB :- upperT') :/ (alpha :* (alpha :> 0))) // TH1 = edittozerotol(TH1, tol1) if (any(TH1 :> maxVal)) { maxindex(TH1, 1, mi, w) leavingVar = mi[1]; lvj = 1+leavingVar maxVal = TH1[leavingVar,1]; leavingCase = 1 } // 2. alpha's negative min value TH2 = boundM[1,evj].upper :+ ((rawXB :- lowerT') :/ (alpha :* (alpha :< 0))) // TH2 = edittozerotol(TH2, tol1) if (any(TH2 :> maxVal)) { maxindex(TH2, 1, mi, w) leavingVar = mi[1]; lvj = 1+leavingVar maxVal = TH1[leavingVar,1]; leavingCase = 2 } // 3. get the enteringVar's lower value if (boundM[1,evj].lower > maxVal) { maxVal = boundM[1,evj].lower; leavingCase = 3 } if (trace == "trace") { displayas("txt"); msg = "Select leaving variable.[MaxVal]" printf("\n[DMUi=%g]%s-LOOP[%g]: %s(%g:%g)\n", dmui, tracename, calci, msg, leavingVar, maxVal) display("XB | alpha:(XB=Bi*b, alpha=Bi*Aj):");XB,alpha display("[MaxVal]enteringVar's lower | theta1 | theta2") printf("\n[boundM[1,%g].lower:%g][leavingCase:%g]\n", evj, boundM[1,evj].lower, leavingCase); TH1,TH2 } if (leavingCase == 1) { boundM[1,lvj].val = boundM[1,lvj].upper } else if (leavingCase == 2) { boundM[1,lvj].val = boundM[1,lvj].lower } else { // if (leavingCase == 3) boundM[1,evj].val = boundM[1,evj].lower; continue; } } } // If no leaving variable exits if (leavingVar == 0) { if (trace == "trace") display("Break: No more candidate for leaving variable.") lpresult.rc = 3498; lpresult.rmsg = sprintf("%s[DMUi=%g][LOOP=%g]%s", "No Solution(No more candidate for leaving variable):", dmui, calci, tracename) break } // theta is not leaving at phaseII. if (phase == 2 && VARS[,lvj] == 1) { if (trace == "trace") display("Break: theta(¥è) is not leaving.") break } // --------------------------------------------------------------------- // reply calculatation result. // --------------------------------------------------------------------- LVi[leavingVar,1] = VARS[,evj] _swapcols(M, lvj, evj) _swapcols(boundM, lvj, evj) _swapcols(VARS, lvj, evj) // Clear artificial variable if (VARS[,evj] < 0) { T = J(1, cols(VARS), 1); T[1,evj] = 0 VARS = select(VARS, T) M = select(M, T) boundM = select(boundM, T) mcols = cols(M) // if exist artificial(phase II) if (existArtificial) { existArtificial = any(VARS[,2::1+slacks] :< 0) } } if (trace == "trace") { printf("\n[DMUi=%g]%s-LOOP[%g]: updated tableau.[%g(%g) <--> %g(%g)]\n", dmui, tracename, calci, lvj, leavingVar, evj, enteringVar); M display("LVi: Entered variable's VARS index value."); LVi display("VARS: Variable's index."); VARS } } //end of main for // return lpresult if (calci > maxiter) { lpresult.rc = 3498; lpresult.rmsg = sprintf("%s[DMUi=%g][LOOP=%g]%s", "No Solution(LOOP greater than maxiter):", dmui, calci, tracename) } if(lpresult.rc) display(lpresult.rmsg) lpresult.xVal = xVal lpresult.XB = XB return(lpresult) } /* A[.,lvj] <--> A[.,evj] */ function _swapcols( transmorphic matrix A, real scalar lvj, real scalar evj ) { transmorphic colvector v v = A[., lvj] A[., lvj] = A[., evj] A[., evj] = v } function replacesubmat ( transmorphic matrix M, real scalar row, real scalar col, transmorphic matrix T ) { M[|row,col\row + rows(T) - 1, col + cols(T) - 1|] = T } function _setup_dearslt_names(string scalar dearsltmat, string scalar dmuinmat, string scalar dmuoutmat ) { string matrix DMU_CS // dmu in matrix column stripes string matrix DEARSLT_CS // dea result matrix column stripes string matrix DEARSLT_RS // dea result matrix row stripes real matrix M real scalar mcols, cnt M = st_matrix(dearsltmat) mcols = cols(M) DMU_CS = st_matrixcolstripe(dmuinmat) for (i = 1; i <= rows(DMU_CS); i++) { DMU_CS[i, 1] = "ref" } DEARSLT_CS = ("","rank"\"","theta")\DMU_CS\ // column join st_matrixrowstripe(dmuinmat)\st_matrixrowstripe(dmuoutmat) if (mcols - rows(DEARSLT_CS) > 0) { cnt = 0 for (i = rows(DEARSLT_CS)+1 ; i <= mcols ; i++) { DEARSLT_CS = DEARSLT_CS \ ("slack", "slack_" + strofreal(++cnt)) } } DEARSLT_RS = st_matrixcolstripe(dmuinmat) // name the row and column of dea result matrix st_matrixrowstripe(dearsltmat, DEARSLT_RS) st_matrixcolstripe(dearsltmat, DEARSLT_CS) } /** * deamat - dmucount x ( 1(theta) + dmu count + slcak(in, out) count) */ function _dmurank( string scalar deamat, real scalar dmuincount, real scalar dmuoutcount, real scalar minrank, real scalar tol ) { real matrix M real rowvector v, vv, retvec, slcaksum real scalar m, mm, row M = st_matrix(deamat) v = round(M[,1], tol) if (minrank) minindex(v, rows(v), i, w) else maxindex(v, rows(v), i, w) retvec = J(rows(v), 1, .) if (allof(w[,2], 1)) { retvec[i[1::rows(v)]] = (1::rows(v)) } else { // rank correction for ties slcaksum = rowsum(M[|1,cols(M) - (dmuincount + dmuoutcount - 1)\.,.|]) for (m = 1; m <= rows(w); m++) { if (w[m,2] >= 2) { vv = i[w[m,1]::(w[m,1] + w[m,2] - 1)] if (minrank) maxindex(slcaksum[vv], w[m,2], ii, ww) else minindex(slcaksum[vv], w[m,2], ii, ww) for (mm = 1; mm <= rows(ww); mm++) { for (row = ww[mm,1]; row < ww[mm,1] + ww[mm,2]; row++) { retvec[vv[ii[row]]] = w[m,1] + ww[mm,1] - 1 } } } else { retvec[i[w[m,1]]] = w[m,1] // row = w[m,1] } } } st_matrix("r(rank)", retvec) } function maxvecindex( string scalar vecname ) { real matrix A A = st_matrix(vecname) maxindex(A, 1, i, w) st_numscalar("r(maxval)", A[i[1]]) st_numscalar("r(maxindex)", i[1]) st_matrix("r(maxindexes)", i) } function minvecindex( string scalar vecname ) { real matrix A A = st_matrix(vecname) if (sum(A :< .) > 0) { minindex(A, 1, i, w) st_numscalar("r(minval)", A[i[1]]) st_numscalar("r(minindex)", i[1]) st_matrix("r(minindexes)", i) } // if overall missing value. else { st_numscalar("r(minval)", .) st_numscalar("r(minindex)", 0) st_matrix("r(minindexes)", 0) } } function _roundmat( string scalar matname, real scalar tol ) { real matrix A A = round(st_matrix(matname), tol) st_matrix(matname, A) } function _uniqrowmat( string scalar matname, string scalar varname ) { st_matrix(matname, sort(uniqrows(st_data(., varname)), 1)) } function _file_exists( string scalar fn ) { st_numscalar("r(fileexists)", fileexists(fn)) } end // End of the MATA Definition Area ---------------------------------------------