The feasiblity pump

#
# Feasibility pump (prototype)
# using the Xpress Python interface
#

from __future__ import print_function

import xpress as xp

def getRoundedSol (sol, I):
    rsol = sol[:]
    for i in I:
        rsol[i] = round (sol[i])
    return rsol

def computeViol (p, viol, rtype, m):
    for i in range (m):
        if rtype[i] == 'L':
            viol[i] = -viol[i]
        elif rtype[i] == 'E':
            viol[i] = abs (viol[i])
    return max (viol[:m])

p = xp.problem ()
p.read ('test.lp')

n = p.attributes.cols # number of columns
m = p.attributes.rows # number of rows
N = range (n)

vtype = []                   # create empty vector
p.getcoltype (vtype, 0, n-1) # obtain variable type ('C', for continuous)

I = [i for i in N if vtype[i] != 'C'] # discrete variables

V = p.getVariable ()

p.lpoptimize()
sol = []
p.getlpsol (x = sol)
roundsol = getRoundedSol (sol, I)
slack = []
p.calcslacks (roundsol, slack)
rtype = []
p.getrowtype (rtype, 0, m - 1)

#
# If x_I is the vector of integer variables and x_I* is its LP value,
# the auxiliary problem is
#
# min |x_I - x_I*|_1
# s.t. x in X
#
# where X is the original feasible set. Add new variables y and set
# their sum as the objective, then define the l1 norm with the
# constraints
#
# y_i >=    x_i - x_i*
# y_i >= - (x_i - x_i*)
#

y = [xp.var () for i in I]

p.addVariable (y) # add variables for new objective
p.setObjective (xp.Sum (y)) # objective 

defPenaltyPos = [y[i] >=   V[I[i]] for i in range (len (I))] # rhs to be configured later
defPenaltyNeg = [y[i] >= - V[I[i]] for i in range (len (I))]

p.addConstraint (defPenaltyPos, defPenaltyNeg)

slackTol = 1e-4

while computeViol (p, slack, rtype, m) > slackTol:

    # modify definition of penalty variable
    p.chgrhs (defPenaltyPos, [- roundsol[i] for i in I])
    p.chgrhs (defPenaltyNeg, [  roundsol[i] for i in I])

    # reoptimize
    p.lpoptimize ()
    p.getlpsol (x = sol)

    roundsol = getRoundedSol (sol, I)
    p.calcslacks (roundsol, slack)

# Found feasible solution

print ('feasible solution:', roundsol[:n])