R version 2.9.0 (2009-04-17) Copyright (C) 2009 The R Foundation for Statistical Computing ISBN 3-900051-07-0 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. Natural language support but running in an English locale R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > # smplx.r simplex algorithm > # J F Monahan, July 2009 > # > # first the code > # > # phase 1 / phase 2 algorithm > p1p2 <- function(A,b,c,m,n) { + p <- n+(1:m) # indices of starting basic feasible solution + print(dim(A)) + Anew <- cbind(A*sign(b),diag(rep(1,m))) # append identity + print(dim(Anew)) + cnew <- c(rep(0,n),rep(1,m)) # new cost vector + feas <- simplex(p,Anew,abs(b),cnew,m,n+m) # solve lp problem + print("******donewithphase1****************") + print(feas) + print("phase1resultsabove") + if( (feas[[3]] == 0)*(feas[[4]] == 0 ) ) { # found feasible soln + p <- feas[[1]] + p1p2 <- simplex(p,A,b,c,m,n) } + else { # no feasible soln + print("nofeasiblesolution") + p1p2 <- list(feas[[1]],feas[[2]],3,0) } + } # end of function p1p2 > # > # simplex algorithm > simplex <- function(p,A,b,c,m,n) { # dumb simplex algorithm + # p is vector of indices for basic, feasible soln + # A matrix, b rhs vector, c costs + # m constraints, n variables + #p# print(dim(A)) + #p# print(p) + B <- A[,p] # matrix of starting solution + BiAb <- solve(B,cbind(A,b)) + # start loop + for( i in (1:100)) { # begin loop + #p# print(p) + pn <- (1:n)[-p] + #p# print(pn) + cp <- c[p] + cn <- c[pn] + x <- rep(0,n) + x[p] <- BiAb[,n+1] + #p# print(BiAb[,n+1]) + #p# print(sum(cp*BiAb[,n+1])) # current obj fn value + # find feasible directions + zjmcj <- cp %*% BiAb[,pn] - cn + #p# print(zjmcj) + flist <- (1:(n-m))*(zjmcj > 1.e-10) + print(flist) + #p print((BiAb[,flist]>1.e-10)) + #p print(apply((BiAb[,flist]> 1.e-10)+0,2,sum)) + if( sum(flist) == 0 ) { print("done") + code <- 0 + break } + ibig <- which.max(zjmcj) # ibig is in, what is out? + #p# print(c(ibig,max(zjmcj))) + ibig <- pn[ibig] + #p# print(ibig) + if( max(BiAb[,ibig]) < 1.e-10 ) { + print("unbound") + code <- 2 + break } + v <- BiAb[,n+1]/BiAb[,ibig] + v[ (BiAb[,ibig]<1.e-10) ] <- 1.e30 + lmin <- which.min( v ) # lmin is out + ### print(c(lmin,v[lmin])) + v <- BiAb[,ibig]/BiAb[lmin,ibig] + v[lmin] <- (BiAb[lmin,ibig]-1)/BiAb[lmin,ibig] + BiAb <- BiAb - outer(v,BiAb[lmin,]) + ### print(BiAb) + p[lmin] <- ibig + code <- 1 # too many iterations if leave with this + } # end of loop + x <- rep(0,n) + x[p] <- BiAb[,n+1] + # return indices, x, return code, obj fn value + simplex <- list(p,x,code,sum(c*x)) } # done with simplex > # call phase 1/ phase2 algorithm > # > # first problem from Wagner, p103, 119-120 > m <- 3 > n <- 7 > A <- matrix(c(1,7,3,1,5,5,1,3,10,1,2,15,1,0,0,0,1,0,0,0,1),m,n) > b <- c(15,120,100) > c <- c(-4,-5,-9,-11,0,0,0) > A [,1] [,2] [,3] [,4] [,5] [,6] [,7] [1,] 1 1 1 1 1 0 0 [2,] 7 5 3 2 0 1 0 [3,] 3 5 10 15 0 0 1 > b [1] 15 120 100 > c [1] -4 -5 -9 -11 0 0 0 > prob1 <- p1p2(A,b,c,m,n) [1] 3 7 [1] 3 10 [1,] 1 2 3 4 5 6 7 [1,] 1 2 3 4 5 0 0 [1,] 0 0 0 4 5 0 0 [1,] 0 0 0 0 0 0 0 [1] "done" [1] "******donewithphase1****************" [[1]] [1] 1 6 4 [[2]] [1] 10.416667 0.000000 0.000000 4.583333 0.000000 37.916667 0.000000 [8] 0.000000 0.000000 0.000000 [[3]] [1] 0 [[4]] [1] 0 [1] "phase1resultsabove" [1,] 0 2 0 0 [1,] 0 0 0 0 [1] "done" > # return indices, x, return code, obj fn value > prob1 [[1]] [1] 1 6 3 [[2]] [1] 7.142857 0.000000 7.857143 0.000000 0.000000 46.428571 0.000000 [[3]] [1] 0 [[4]] [1] -99.28571 > # > # second problem from VJ Bowman > m <- 6 > n <- 10 > A <- cbind( + matrix(c(1,2,5,0,0,0,4,3,1,0,0,0,5,0,0,1,0,-3,2,0,0,0,1,-4),6,4), + diag(6) ) # block identity in cols 5-10 > b <- c(7,6,5,4,3,-12) > c <- c(-1,-8,-5,-6,rep(0,6)) > A [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [1,] 1 4 5 2 1 0 0 0 0 0 [2,] 2 3 0 0 0 1 0 0 0 0 [3,] 5 1 0 0 0 0 1 0 0 0 [4,] 0 0 1 0 0 0 0 1 0 0 [5,] 0 0 0 1 0 0 0 0 1 0 [6,] 0 0 -3 -4 0 0 0 0 0 1 > b [1] 7 6 5 4 3 -12 > c [1] -1 -8 -5 -6 0 0 0 0 0 0 > prob2 <- p1p2(A,b,c,m,n) [1] 6 10 [1] 6 16 [1,] 1 2 3 4 5 6 7 8 9 0 [1,] 1 2 3 0 5 6 7 8 0 0 [1,] 0 2 0 4 0 6 7 0 0 0 [1,] 0 0 0 4 5 6 7 0 0 0 [1,] 0 0 0 4 5 6 0 0 0 0 [1,] 0 0 0 4 5 0 0 0 0 0 [1,] 0 0 0 4 0 0 0 0 0 0 [1,] 0 2 0 0 0 0 0 0 0 0 [1,] 0 0 0 0 0 0 0 0 0 0 [1] "done" [1] "******donewithphase1****************" [[1]] [1] 4 6 1 8 7 3 [[2]] [1] 1.000000e+00 0.000000e+00 1.314768e-31 3.000000e+00 0.000000e+00 [6] 4.000000e+00 8.881784e-15 4.000000e+00 0.000000e+00 0.000000e+00 [11] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 [16] 0.000000e+00 [[3]] [1] 0 [[4]] [1] 0 [1] "phase1resultsabove" [1,] 1 0 0 0 [1,] 0 0 0 0 [1] "done" > # return indices, x, return code, obj fn value > prob2 [[1]] [1] 4 6 2 8 7 3 [[2]] [1] 0.00 0.25 0.00 3.00 0.00 5.25 4.75 4.00 0.00 0.00 [[3]] [1] 0 [[4]] [1] -20 > # > # third problem from Wagner -- unbounded > m <- 2 > n <- 4 > A <- matrix(c(1,-1,-1,1,1,0,0,1),m,n) > b <- c(1,1) > c <- c(-1,-1,0,0) > A [,1] [,2] [,3] [,4] [1,] 1 -1 1 0 [2,] -1 1 0 1 > b [1] 1 1 > c [1] -1 -1 0 0 > prob3 <- p1p2(A,b,c,m,n) [1] 2 4 [1] 2 6 [1,] 0 0 3 4 [1,] 0 2 3 0 [1,] 0 0 0 0 [1] "done" [1] "******donewithphase1****************" [[1]] [1] 3 2 [[2]] [1] 0 1 2 0 0 0 [[3]] [1] 0 [[4]] [1] 0 [1] "phase1resultsabove" [1,] 1 0 [1] "unbound" > # return indices, x, return code, obj fn value > prob3 [[1]] [1] 3 2 [[2]] [1] 0 1 2 0 [[3]] [1] 2 [[4]] [1] -1 > # > # fourth problem from Wagner -- no solution > # same A, c, but negate b > b <- -b > A [,1] [,2] [,3] [,4] [1,] 1 -1 1 0 [2,] -1 1 0 1 > b [1] -1 -1 > c [1] -1 -1 0 0 > prob4 <- p1p2(A,b,c,m,n) [1] 2 4 [1] 2 6 [1,] 0 0 0 0 [1] "done" [1] "******donewithphase1****************" [[1]] [1] 5 6 [[2]] [1] 0 0 0 0 1 1 [[3]] [1] 0 [[4]] [1] 2 [1] "phase1resultsabove" [1] "nofeasiblesolution" > # return indices, x, return code, obj fn value > prob4 [[1]] [1] 5 6 [[2]] [1] 0 0 0 0 1 1 [[3]] [1] 3 [[4]] [1] 0 > # > # fifth problem is L1 regression > m <- 6 # number of obs > n <- 16 # 2*nobs + 2*(cols of X) > X <- cbind(rep(1,m),(1:m)) # design matrix > A <- cbind(diag(rep(1,m)),diag(rep(-1,m)),X,-X) > b <- c(4.5,5.5,6.5,8,10,12) # y's > c <- c(rep(1,2*m),0,0,0,0) > A [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [1,] 1 0 0 0 0 0 -1 0 0 0 0 0 1 1 [2,] 0 1 0 0 0 0 0 -1 0 0 0 0 1 2 [3,] 0 0 1 0 0 0 0 0 -1 0 0 0 1 3 [4,] 0 0 0 1 0 0 0 0 0 -1 0 0 1 4 [5,] 0 0 0 0 1 0 0 0 0 0 -1 0 1 5 [6,] 0 0 0 0 0 1 0 0 0 0 0 -1 1 6 [,15] [,16] [1,] -1 -1 [2,] -1 -2 [3,] -1 -3 [4,] -1 -4 [5,] -1 -5 [6,] -1 -6 > b [1] 4.5 5.5 6.5 8.0 10.0 12.0 > c [1] 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 > p <- (1:m) # we have feasible soln, b positive > prob5 <- simplex(p,A,b,c,m,n) [1,] 0 0 0 0 0 0 7 8 0 0 [1,] 0 0 0 0 5 0 0 8 0 0 [1,] 0 0 0 0 0 6 0 8 0 0 [1,] 0 0 0 0 0 6 0 0 0 0 [1,] 0 0 0 0 0 0 0 0 0 0 [1] "done" > # return indices, x, return code, obj fn value > prob5 [[1]] [1] 1 9 13 14 10 6 [[2]] [1] 0.5 0.0 0.0 0.0 0.0 0.5 0.0 0.0 0.5 0.5 0.0 0.0 2.5 1.5 0.0 0.0 [[3]] [1] 0 [[4]] [1] 2 > # done > rm(list=ls()) > q()