I want to implement a branch and bound algorithm. How to warm start dual simplex in Cylp?

[abdxyz]

Hello, I'm trying to implement my own simple branch and bound algorithm using Cylp. I want to use the basis in the parent node to warm start the child node. In the CyClpSimplex, I didn't find any API to do this. In this page, https://www.coin-or.org/Doxygen/Clp/classCoinWarmStartPrimalDual.html, I find that Clp support warm start. Do you have some idea? Thanks a lot.

[abdxyz]

I found @tkralphs has implemented a CuPPy package using the Cutting Plane algorithm to solve MIP. But the code just uses the primal simplex algorithm. It did not use the simplex information to warm start. Are there some difficulties to warm start? I'm a Ph.D. student, I want to implement my own solver to understand this subject deeply. I would appreciate it a lot if you can give me some advice.

[tkralphs]

The way simplex is called in CuPPY does in fact do the warm start.

lp.primal(startFinishOptions = 'x')

If you don't want warm-starting, you can explicitly call initialSolve().

In case you are interested, I also have Python branch-and-bound implementation here, which you may want to look at. It currently uses PuLP as the LP interface, but I was planning to change to CyLP when I get some time. I will integrate cutting plane generation to get a full branch-and-cut.

[abdxyz]

Thanks a lot. I'll first read your code and then implement my solver. I'm a little confused about the warm start. I'm taking an integer course now. It seems when adding a constraint or changing a variable bound, The primal feasibility isn't satisfied. It's better to use dual simplex. Do you know any document about CLP which explaining the details of simplex? Thanks a lot.

[tkralphs]

Yeah, you're right, it should be dual() there if you're concerned with performance. I've really only used CuPPy for solving very small MILPs, mostly for teaching, so I didn't think much about performance issues.

[abdxyz]

Ok, I'm writing the code now. I want it to be a general framework support B&B, cutting plane, and other methods like column generation and decomposition methods. Like you, I write it mainly for educational purposes. I already finished B&B. After I finished cutting plane algorithm, I will open-source it. If I have some other questions, I then ask you. Thanks a lot.

[abdxyz]

I have finished my first edition. I open-source it here https://github.com/abdxyz/mathu-solver. It is a general framework, every important function is in a python class or model. Every model is independent and can be combined together easily. I already finish the basic B&B and Cutting algorithm. And I will improve it for other advanced algorithms like Branch-and-cut and column generation.

[tkralphs]

Your framework looks nice, along the lines of the abstractions in Alps. It would be interesting to put a Python wrapper around Alps and allow the branching and bounding methods to be written in Python, while the search is done using the Alps library itself. I'm teaching integer programming next semester and since this is along the lines of what I tended to do with the Python branch-and-bound in GrUMPy, I may use and contribute to it.

I do have one request. I noticed that you put a copy of CuPPy in your repository and I think you copied and modified some stuff from it in other places. Meanwhile, your code doesn't have a license. I have been remiss in putting license and copyright information into the source files of CuPPy, as no one was really using it anyway. Since CuPPy is on Pypi and is easily installed with pip, would you mind just using it that way and not re-distributing those files without the original license? Of course, the license does allow redistribution if you include the license with the redistributed files. If you put a license on the code you have written yourself, it would be easier for me to use it with my students. Thanks for respecting that!

I should also point out DipPy, which is a Python wrapper around Dip. It already allows you to do cut and column generation in Python, while utilizing a state-of-the-art branch-cut-and-price framework underneath. That is probably not what you're going for, but thought I would point it out.

[abdxyz]

Ok, I had included an EPL license. I'm a novice to open source code. I think this license will give you freedom to do your own job. I write it mainly for educational purposes. Next, I'll use GrUMPy to make it display friendly. Then I will add other algorithms like branch-cut-and-price and column generation. I want to implement the algorithms written in the Integer Programming book. Any contribution will be welcomed.

[abdxyz]

I have tried GrUMPy. The cutting plane part works well. But I can't get the B&B tree as your image. I wrote my test result in an issue. I think it is because of the python's version.

[abdxyz]

Hello. I'm now improving my code to warm start simplex using the dual simplex algorithm. I want to use the basises information of a parent node to accelerate the child node. In my implementation, these two nodes are two simplex objects. In a child node, only one constraint is added, like an upper bound change of one variable. In my opinion, I should reuse the parent tableau and add one row. But I didn't find any API to set the tableau of the child node. I find that it may be useful to get the parent basis information and variable solution and use it to update the child node. And only do some extra work for the added constraint. I only find APIs to get the date, but it may be not supported to set the data.

    property dualConstraintSolution:
        '''
        Dual variables' solution
        :rtype: Numpy array
        '''
        def __get__(self):
            ret =  <object>self.CppSelf.getDualRowSolution()
            if self.cyLPModel:
                m = self.cyLPModel
                inds = m.inds
                d = {}
                for c in inds.constIndex.keys():
                    d[c] = ret[inds.constIndex[c]]
                ret = d
            else:
                pass
                #names = self.variableNames
                #if names:
                #    d = CyLPSolution()
                #    for i in range(len(names)):
                #        d[names[i]] = ret[i]
                #    ret = d
            return ret

The following code is from the example of CLP.

     model2.transposeTimes(1.0, randomArray, model2.objective());
     delete [] randomArray;
     // solve
     model2.primal();
     // first check okay if solution values back
     memcpy(model.primalColumnSolution(), model2.primalColumnSolution(),
            originalNumberColumns * sizeof(double));
     memcpy(model.primalRowSolution(), model2.primalRowSolution(),
            numberRows * sizeof(double));
     int iColumn;
     for (iColumn = 0; iColumn < originalNumberColumns; iColumn++)
          model.setColumnStatus(iColumn, model2.getColumnStatus(iColumn));

     for (iRow = 0; iRow < numberRows; iRow++) {
          if (model2.getRowStatus(iRow) == ClpSimplex::basic) {
               model.setRowStatus(iRow, ClpSimplex::basic);
          } else {
               model.setRowStatus(iRow, model2.getColumnStatus(iRow + originalNumberColumns));
          }
     }
     model.primal(0);
     // and now without solution values

     for (iColumn = 0; iColumn < originalNumberColumns; iColumn++)
          model3.setColumnStatus(iColumn, model2.getColumnStatus(iColumn));

     for (iRow = 0; iRow < numberRows; iRow++)
          model3.setRowStatus(iRow, model2.getColumnStatus(iRow + originalNumberColumns));
     model3.primal(0);

I think it may be useful. Thanks a lot. Hope you are interested to improve CyLP so that we can implement a really mip solver.

[abdxyz]

I add the following code in the CyClpSimplex.pyx

from libc.stdlib cimport malloc, free, memcpy

def set_dualvariable(self,np.ndarray[np.double_t,ndim=1] value):
        print(123)
        memcpy(<double*> self.CppSelf.getDualRowSolution(),<double*> value.data,self.nConstraints*8)

But it doesn't work. It output

AttributeError: 'cylp.cy.CyClpSimplex.CyClpSimplex' object has no attribute 'set_dualvariable'

It is very strange if I change print(123) to print(123 I can still install CyLP without any error. So I think it may be some error from Cython and pip building process. Do you have some time to try this. Thanks a lot.

[abdxyz]

I adopted the following steps. It works fine.

from libc.string cimport memcpy

property dualConstraintSolution:
        '''
        Dual variables' solution

        :rtype: Numpy array
        '''
        def __get__(self):
            ret =  <object>self.CppSelf.getDualRowSolution()
            if self.cyLPModel:
                m = self.cyLPModel
                inds = m.inds
                d = {}
                for c in inds.constIndex.keys():
                    d[c] = ret[inds.constIndex[c]]
                ret = d
            else:
                pass
                #names = self.variableNames
                #if names:
                #    d = CyLPSolution()
                #    for i in range(len(names)):
                #        d[names[i]] = ret[i]
                #    ret = d
            return ret
        def __set__(self,np.ndarray[np.double_t,ndim=1] value):
                print(123)
                memcpy(<void *>self.CppSelf.getDualRowSolution(),<void*>value.data,self.nConstraints*8)

Then I use

./regenerate_cpp_files.sh 

Then I call setBasisStatus and directly assign the dual variable values with the optimal value. It works normally and says optimal arrived. The LP iterations equal 0. Following is my test script.

import numpy as np
from cylp.cy import CyClpSimplex
from cylp.py.modeling.CyLPModel import CyLPArray

A = [[50, 31], [-3, 2, ]]
b = [250, 4]
c = [-1, -0.64]
l = [1, 0]
u = [10000, 10000]

model = CyClpSimplex()
x = model.addVariable('x', len(c))
A = np.matrix(A)
b = np.matrix(b)
l = CyLPArray(l)
u = CyLPArray(u)
c = CyLPArray(c)
model += (A * x <= b)
model += l <= x <= u
model.objective = c * x

model.setBasisStatus(np.array([1, 1]).astype(np.int32), np.array([3, 3]).astype(np.int32))
model.dualConstraintSolution = np.array([-0.02031088, -0.00518135])

model.dual()
print(model.getBasisStatus())
print(model.primalConstraintSolution)  # row value
print(model.primalVariableSolutionAll)  # row slack
print(model.iteration)
print(model.dualConstraintSolution)  # dual variable
print(model.dualVariableSolution)  # dual slack

Following is the output

123
Clp0000I Optimal - objective value -5.0984456
(array([1, 1], dtype=int32), array([3, 3], dtype=int32))
{'R_1': array([250.,   4.])}
[1.94818653e+000 4.92227979e+000 1.10638848e+200 1.63041663e-322]
0
{'R_1': array([-0.02031088, -0.00518135])}
{'x': array([0., 0.])}