from math import sqrt, exp
from numpy import matrix, bmat, ones, zeros, identity, diagflat, asscalar
from numpy.linalg import solve, inv, det

""" Returns an n-dimension column vector of zeros """
def z(n):
  return zeros([n,1])

""" Returns an n-dimension column vector of ones """
def e(n):
  return ones([n,1])

""" Checks that A, b, c have consistent dimensions, and prints a message."""  
def check_input(A, b, c, function_name):
  (m, n) = A.shape
  assert b.shape == (m, 1)
  assert c.shape == (n, 1)
  print ">", function_name, "called for model with", n, "variables and", m, "constraints"
  return (m, n)

""" Solves the model min c'x subject to Ax <= b, x >= 0.
    This function returns a tuple (obj, x, y), where obj is the optimal 
    objective value, x is the optimal primal solution, and y is the 
    optimal dual solution."""
def intpm_leq(A, b, c, alpha):
  (m, n) = check_input(A, b, c, "intpm_leq")

  # Introduce slack variables
  A = bmat([[A, identity(m)]])
  c = bmat([[c], [z(m)]])
  
  # Run the interior path algorithm for models with equality constraints
  return intpm_eq(A, b, c, alpha)


""" Solves the model min c'x subject to Ax = b, x >= 0. 
    This function returns a tuple (obj, x, y), where obj is the optimal 
    objective value, x is the optimal primal solution, and y is the 
    optimal dual solution."""
def intpm_eq(A, b, c, alpha):
  (m, n) = check_input(A, b, c, "intpm_eq")
  
  # Construct an auxiliary problem for which it is easy to find an intiial 
  # solution. See Section 6.4.2.
  MP = alpha * alpha
  MD = alpha * alpha * (n+1) - alpha * sum(c)

  A = bmat([
    [A, b-alpha*A*e(n), z(m)], 
    [(alpha*e(n)-c).T, z(1), alpha*e(1)]
  ])
  b = bmat([[b], [MD*e(1)]])
  c = bmat([[c], [MP*e(1)], [z(1)]])
  x0 = bmat([[alpha*e(n)], [e(1)], [alpha*e(1)]])
  mu = alpha * alpha
  
  # Run the interior path algorithm
  return intpm(A, b, c, x0, mu)

""" Solves the model min c'x subject to Ax = b, x >= 0, by using the given
    initial strictly feasible point x0, interior path parameter mu, and 
    precision t.     
    This function returns a tuple (obj, x, y), where obj is the optimal 
    objective value, x is the optimal primal solution, and y is the 
    optimal dual solution."""
def intpm(A, b, c, x0, mu0, t = 10):
  (m, n) = check_input(A, b, c, "intpm")
  
  # Calculate theta and the target value of mu
  theta = 1. / (6. * sqrt(n))
  mu_final = exp(-t)/n
  x = x0
  mu = mu0
  iterations = 0

  # Perform interior path iterations
  print "Starting interior path algorithm"
  while mu > mu_final:
    iterations += 1
    
    # Update mu
    mu = (1 - theta) * mu
    
    # Calculate next point
    # NOTE: solve(AX * AX.T, AX) calculates inv(AX * AX.T) * AX in a more 
    # efficient and numerically stable manner.
    X = diagflat(x) 
    AX = A * X
    PAX = identity(n) - AX.T * solve(AX * AX.T, AX)
    x = x + X * PAX * (e(n) - X * c / mu)

  # Calculate the corresponding dual solution and the optimal objective value
  w = mu * solve(X, e(n))
  y = solve(A * A.T, A * (c - w))
  obj = asscalar(c.T * x)
  
  # Print out some information
  print "Interior path algorithm terminated after", iterations, "iterations"
  print "Objective value =", obj
  print "Duality gap =", asscalar(b.T * y - c.T * x)
  print "Primal solution:"
  for (i, value) in enumerate(x):
    print "x%-2d = %16.8f" % (i, value)
  print "Dual solution:"
  for (i, value) in enumerate(y):
    print "y%-2d = %16.8f" % (i, value)
  return (obj, x, w)


@@ -1,12 +1,1143 @@
+# Interior path method implementation%0A# Copyright (C) 2013-2015 Y. Zwols and G. Sierksma%0A# %0A# Permission is hereby granted, free of charge, to any person obtaining%0A# a copy of this software and associated documentation files (the%0A# %22Software%22), to deal in the Software without restriction, including%0A# without limitation the rights to use, copy, modify, merge, publish,%0A# distribute, sublicense, and/or sell copies of the Software, and to%0A# permit persons to whom the Software is furnished to do so, subject to%0A# the following conditions:%0A# %0A# The above copyright notice and this permission notice shall be%0A# included in all copies or substantial portions of the Software.%0A# %0A# THE SOFTWARE IS PROVIDED %22AS IS%22, WITHOUT WARRANTY OF ANY KIND,%0A# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF%0A# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.%0A# IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY CLAIM, DAMAGES OR%0A# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,%0A# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR%0A# OTHER DEALINGS IN THE SOFTWARE.%0A%0A
 from math im

Update intpm.py