# Copyright (C) 2005, Paulo Jose da Silva e Silva

from numarray import *
from coin import *

si = OsiSolver()
si.readMps('p0033.mps')

si.setObjSense(-1.0)
si.initialSolve()
si.setObjSense(1.0)
si.enableSimplexInterface(True)
numberIterations = 0
numberColumns = si.getNumCols()
numberRows = si.getNumRows()
fakeCost = si.getObjCoefficients()
duals = zeros(numberRows, Float)
djs = zeros(numberColumns, Float)
solution = si.getColSolution()

while True:
    if numberIterations & 1 == 0:
        dj = si.getReducedCost()
        dual = si.getRowPrice()
    else:
        dj = djs
        dual = duals
        si.getReducedGradient(djs, duals, fakeCost)

    colIn = 9999
    direction = 1
    best = 1.0e-6
    bigDual = dual.argmax()
    if dual[bigDual] > best:
        best = dual[bigDual]
        direction = -1
        colIn = -bigDual - 1
    for i in range(numberColumns):
        value = dj[i]
        if value < -best and solution[i] < 1.0e-6:
            direction = 1
            best = -value
            colIn = i
        elif value > best and solution[i] > 1.0 - 1.0e-6:
            direction = -1
            best = value
            colIn = i
    if colIn == 9999: break
    out, outStatus = copy_intp(0), copy_intp(0)
    stepSize = copy_doublep(0.0)
    si.primalPivotResult(colIn, direction, out, outStatus,
                         stepSize, 'NULL')
    print 'out %d, direction %d theta %g' % (intp_value(out),
                                             intp_value(outStatus),
                                             doublep_value(stepSize))
    solution = si.getColSolution()
    numberIterations += 1
si.disableSimplexInterface()
si.resolve()
assert(not si.getIterationCount())
si.setObjSense(-1.0)
si.initialSolve()