# 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
     print(dim(A))
     print(p)
     B <- A[,p]     # matrix of starting solution
     BiAb <- solve(B,cbind(A,b))
   # start loop
    for( i in (1:100)) {      # begin loop
       print(p)
     pn <- (1:n)[-p]
       print(pn)
     cp <- c[p]
     cn <- c[pn]
     x <- rep(0,n)
     x[p] <- BiAb[,n+1]
       print(BiAb[,n+1])
       print(sum(cp*BiAb[,n+1]))        # current obj fn value
  # find feasible directions     
     zjmcj <- cp %*% BiAb[,pn] - cn
       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?
       print(c(ibig,max(zjmcj)))
     ibig <- pn[ibig]
       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]
     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 
  alldat <- scan("smplx11.dat")
  Ab <- matrix(alldat,m,n+1,byrow=TRUE)
  A <- Ab[,1:n]
  b <- Ab[,n+1]
  c <- alldat[m*(n+1)+(1:n)]
  A
  b
  c
   prob1 <- p1p2(A,b,c,m,n)
   prob1
 #
 # second problem from VJ Bowman
  m <- 6
  n <- 10 
  alldat <- scan("smplxy.dat")
  Ab <- matrix(alldat,m,n+1,byrow=TRUE)
  A <- Ab[,1:n]
  b <- Ab[,n+1]
  c <- alldat[m*(n+1)+(1:n)]
  A
  b
  c
   prob2 <- p1p2(A,b,c,m,n)
   prob2
 #
 # 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
  b
  c
   prob3 <- p1p2(A,b,c,m,n)
   prob3
 #
 # fourth problem from Wagner -- no solution
 #  same A, c, but negate b
  b <- -b
  A
  b
  c
   prob4 <- p1p2(A,b,c,m,n)
   prob4
 #
 # 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
  b
  c
  p <- (1:m)                  # we have feasible soln, b positive
   prob5 <- simplex(p,A,b,c,m,n)
   prob5
 # done
  rm(list=ls())
  q()