<html>
<head>
<meta content="text/html; charset=utf-8" http-equiv="Content-Type">
</head>
<body link="#0B6CDA" text="#000000" vlink="#551A8B" alink="#EE0000"
bgcolor="#ffffff">
Hello Ted,<br>
this is a revision of my previous mail with added code snippets.
They are excerpts of a copy of DIP source code version 0.92.2 to
which I added some changes, and they show the fixes I made to
address the mentioned issues. It goes without saying that I do not
claim the code to be flawless. Nevertheless, I hope it is helpful to
you.<br>
Best regards<br>
Marcus<br>
<br>
<hr size="2" width="100%"><i><br>
[</i><i><i>Implementation</i>]</i> I found Dippy to suffer from <b>memory
leaking</b>. Since it interfaces with Python, it is responsible
for maintaning the reference counters of the Python objects it deals
with. This is done in some parts of the code, but by far not all. A
crucial situation which came to my attention is the method <tt>DippyDecompApp::solveRelaxed</tt>.
The retrieved columns are converted to objects of the class <tt>DecompVar</tt>,
yet the reference counter on the original Python object is not
decreased, which prevents them from being deleted. Hence, for large
problem instances the memory floods.<br>
<br>
<u><i>DippyDecompApp.cpp:317<br>
</i></u><u><i><br>
</i></u><tt>DecompSolverStatus DippyDecompApp::solveRelaxed(const
int whichBlock,</tt><tt><br>
</tt><tt> const double* redCostX, </tt><tt><br>
</tt><tt> const double convexDual, </tt><tt><br>
</tt><tt> DecompVarList& varList)</tt><tt><br>
</tt><tt>{</tt><tt><br>
</tt><tt> if (!m_pySolveRelaxed) {</tt><tt><br>
</tt><tt> return DecompSolStatNoSolution;</tt><tt><br>
</tt><tt> }</tt><tt><br>
</tt><tt><br>
</tt><tt> PyObject* pRelaxKey = PyList_GetItem(m_relaxedKeys,
whichBlock);</tt><tt><br>
</tt><tt> PyObject* pRedCostList =
pyTupleList_FromDoubleArray(redCostX, m_colList);</tt><tt><br>
</tt><tt> PyObject* pConvexDual = PyFloat_FromDouble(convexDual);</tt><tt><br>
</tt><tt> // call solveRelaxed on DipProblem</tt><tt><br>
</tt><tt><br>
</tt><tt> PyObject* pStatandVarList = PyObject_CallMethod(m_pProb,
"solveRelaxed", "OOd", </tt><tt><br>
</tt><tt> pRelaxKey,</tt><tt><br>
</tt><tt> pRedCostList,</tt><tt><br>
</tt><tt> pConvexDual);</tt><tt><br>
<b><br>
</b><b> Py_DECREF(pRedCostList);</b><b><br>
</b><b> Py_DECREF(pConvexDual);</b><br>
</tt><tt><br>
</tt><tt> if ( (pStatandVarList == NULL) || (pStatandVarList ==
Py_None) ){</tt><tt><br>
</tt><tt> throw UtilException("Error calling method
prob.solveRelaxed()", "solveRelaxed", "DippyDecompApp");</tt><tt><br>
</tt><tt> }</tt><tt><br>
</tt><tt><br>
</tt><tt> // [status, varList] = relaxed_solver(...)</tt><tt><br>
</tt><tt> PyObject * pStatus = PyTuple_GetItem(pStatandVarList,
0);</tt><tt><br>
</tt><tt><br>
</tt><tt> int cStatus = PyInt_AsLong(pStatus);</tt><tt><br>
</tt><tt><br>
</tt><tt> DecompSolverStatus status = (DecompSolverStatus)cStatus;</tt><tt><br>
</tt><tt><br>
</tt><tt> PyObject * pVarList = PyTuple_GetItem(pStatandVarList,
1);</tt><tt><br>
</tt><tt><br>
</tt><tt> int nVars = PyObject_Length(pVarList);</tt><tt><br>
</tt><tt><br>
</tt><tt> // In the new design, we need to allow the possibility
that the user will solve</tt><tt><br>
</tt><tt> // the problem exactly, but not find any solutions with
reduced costs zero</tt><tt><br>
</tt><tt> // The below is is commented out and left in the source
for posterity</tt><tt><br>
</tt><tt> // tkr 11/11/15</tt><tt><br>
</tt><tt> //if (nVars == 0) {</tt><tt><br>
</tt><tt> // throw UtilException("Empty variable list",
"solveRelaxed", "DippyDecompApp");</tt><tt><br>
</tt><tt> //}</tt><tt><br>
</tt><tt><br>
</tt><tt> // solveRelaxed returns 3-tuples (cost, reduced cost,
dictionary of (variable, value) pairs)</tt><tt><br>
</tt><tt> // We can use these to construct a C++ DecompVar objects</tt><tt><br>
</tt><tt> double cost, rc;</tt><tt><br>
</tt><tt> PyObject* pTuple, *pDict, *pKeys, *pCol;</tt><tt><br>
</tt><tt> string name;</tt><tt><br>
</tt><tt> double value;</tt><tt><br>
</tt><tt><br>
</tt><tt> for (int j = 0; j < nVars; j++) {</tt><tt><br>
</tt><tt> pTuple = PySequence_GetItem(pVarList, j);</tt><tt><br>
</tt><tt> cost = PyFloat_AsDouble(PyTuple_GetItem(pTuple,
0));</tt><tt><br>
</tt><tt> rc = PyFloat_AsDouble(PyTuple_GetItem(pTuple,
1));</tt><tt><br>
</tt><tt><br>
</tt><tt> pDict = PyTuple_GetItem(pTuple, 2);</tt><tt><br>
</tt><tt> pKeys = PyDict_Keys(pDict);</tt><tt><br>
</tt><tt> vector<int> varInds;</tt><tt><br>
</tt><tt> vector<double> varVals;</tt><tt><br>
</tt><tt><br>
</tt><tt> for (int n = 0; n < PyObject_Length(pDict); n++) {</tt><tt><br>
</tt><tt> pCol = PyList_GetItem(pKeys, n);</tt><tt><br>
</tt><tt> value = PyFloat_AsDouble(PyDict_GetItem(pDict,
pCol));</tt><tt><br>
</tt><tt> varInds.push_back(m_colIndices[pCol]);</tt><tt><br>
</tt><tt> varVals.push_back(value);</tt><tt><br>
</tt><tt> }</tt><tt><br>
</tt><tt><b><tt><b><br>
</b><b> Py_DECREF(</b></tt></b></tt><tt><b><tt><b><tt>pKeys</tt>);</b></tt></b></tt><tt><b><tt><b><tt><b><br>
</b><b> Py_DECREF(pTuple);</b></tt></b></tt></b><b></b><br>
</tt><tt><br>
</tt><tt> DecompVar* var = new DecompVar(varInds, varVals, rc,
cost);</tt><tt><br>
</tt><tt> var->setBlockId(whichBlock);</tt><tt><br>
</tt><tt> varList.push_back(var);</tt><tt><br>
</tt><tt> }</tt><tt><br>
</tt><tt><br>
<b> Py_DECREF(pStatandVarList);</b><br>
<br>
</tt><tt> return status;</tt><tt><br>
</tt><tt>}</tt><u><i><br>
</i></u><br>
<hr size="2" width="100%"><br>
<i>[Implementation]</i> The method <tt>AlpsDecompTreeNode::getBranchedVar
</tt>seems to be meant to return the <b>variable which is branched
on</b>. However, only one of the four vectors potentially holding
branching decisions is checked. I think both, lower bounds and upper
bounds, need to be checked for at least one branch, e.g. <tt>downBranchLB</tt>_
and <tt>downBranchUB_</tt>.<br>
<br>
<i><u>AlpsDecompTreeNode.h:76</u><br>
<br>
</i><tt> int getBranchedVar() const {</tt><tt><br>
</tt><tt> if (!downBranchLB_.empty()) {</tt><tt><br>
</tt><tt> return downBranchLB_[0].first;</tt><b><tt><br>
</tt></b><b><tt> } else if (!downBranchUB_.empty()) {</tt></b><b><tt><br>
</tt></b><b><tt> return downBranchUB_[0].first;</tt></b><tt><br>
<b>// Not sure whether the following is necessary. Seems to be
hard to decide as the getBranchedVar is not "natural", i.e. it
does only consi</b><b>der one variable when there could be more.</b><b><br>
</b><b>
//</b></tt><b><tt> } else if (!upBranchLB_.empty()) {</tt></b><b><tt><br>
</tt></b><b><tt>// return upBranchLB_[0].first;</tt></b><b><tt><br>
</tt></b><b><tt>// } else if (!upBranchUB_.empty()) {</tt></b><b><tt><br>
</tt></b><b><tt>// return upBranchUB_[0].first;</tt></b><tt><br>
</tt><tt> } else {</tt><tt><br>
</tt><tt> return -1;</tt><tt><br>
</tt><tt> }</tt><tt><br>
</tt><tt> }<br>
<br>
</tt>
<hr size="2" width="100%"><br>
<i>[</i><i>Implementation</i>] In <tt>DecompAlgo::processNode</tt>,
line 1855 the processing of a node is terminated if the <b>lower
bound meets the global upper bound</b>. It takes into account
numeric inexactness. The calling method <tt>AlpsDecompTreeNode::process
</tt>repeats this check in line 309 but in an exact fashion why the
node might not be fathomed when it should be. In my oppinion this
repeated condition should match the original one. Or maybe even
better, <tt>DecompAlgo::getStopCriteria</tt> could be used?<br>
<br>
<i><u>AlpsDecompTreeNode.cpp:309</u><br>
</i><tt><br>
if (quality_ >= currentUB<b> - </b><b>DecompEpsilon</b>)
{</tt><tt><br>
</tt><tt> doFathom = true;</tt><tt><br>
</tt><tt> UTIL_DEBUG(param.msgLevel, 3,</tt><tt><br>
</tt><tt> cout << "Fathom since
thisQuality= "</tt><tt><br>
</tt><tt> << setw(10) <<
UtilDblToStr(thisQuality)</tt><tt><br>
</tt><tt> << " quality_= " <<
setw(10) << UtilDblToStr(quality_)</tt><tt><br>
</tt><tt> << " currentUB = " <<
setw(10) << UtilDblToStr(currentUB)</tt><tt><br>
</tt><tt> << " gap = " <<
setw(10) << UtilDblToStr(gap) << endl;</tt><tt><br>
</tt><tt> );</tt><tt><br>
</tt><tt> }</tt><br>
<br>
<hr size="2" width="100%"><br>
<i>[Implementation]</i> The member <tt>DecompAlgo::m_relGap</tt> is
set in <tt>DecompAlgo::updateObjBound</tt> and used in <tt>DecompAlgo::isGapTight</tt>.
<tt>DecompAlgo::m_relGap</tt>, however, is not reset when entering <tt>DecompAlgo::processNode</tt>.
Therefore, it has an <b>invalid value</b> based on a node processed
before - probably representing a tight <b>gap</b>. I think, this
might lead to stopping the processing of a node immediately as it is
mistakenly believed to have a tight gap, cf. <tt>DecompAlgo::phaseUpdate</tt>,
line 4274. I suggest to reset <tt>DecompAlgo::m_relGap </tt>appropriately
or replace it completely by calls to <tt>DecompAlgo::getNodeLPGap</tt>.<br>
<br>
<i><u>DecompAlgo.h:387</u><br>
</i><tt><br>
bool isGapTight() {</tt><tt><br>
</tt><tt> //TODO: make param</tt><tt><br>
</tt><tt> double tightGap = m_param.MasterGapLimit;</tt><tt><br>
</tt><tt><br>
</tt><tt> //printf("isGapTight m_relGap = %g\n", m_relGap);</tt><tt><br>
</tt><tt> if (m_param.LogDebugLevel >= 2) {</tt><tt><br>
</tt><tt> (*m_osLog) << "DW GAP = " <<
UtilDblToStr(m_relGap)</tt><tt><br>
</tt><tt> << " isTight = " <<
(m_relGap <= tightGap)</tt><tt><br>
</tt><tt> << "\n";</tt><tt><br>
</tt><tt> }</tt><tt><br>
</tt><tt><br>
</tt><tt> if (<b>getNodeLPGap()</b> <= tightGap) {</tt><tt><br>
</tt><tt> return true;</tt><tt><br>
</tt><tt> } else {</tt><tt><br>
</tt><tt> return false;</tt><tt><br>
</tt><tt> }</tt><tt><br>
</tt><tt> }<br>
<br>
</tt>
<hr size="2" width="100%"><br>
<i>[Functionality/Performance]</i> For the problem I consider, the
MILP formulation of the subproblems is quite huge. Therefore, I use
a dummy formulation and prevent DIP from looking at it by always
providing columns with negative reduced cost via my implementation
of <tt>DecompApp::solveRelaxed</tt> if they exist. I do so in a
two-step approach. In a first step the subproblem is solved
heuristically. If this does not lead to new columns an exact method
is applied. Since the exact method is costly, I would like to stick
to the heuristic as long as there are new columns for <i>some</i>
subproblem. This stands in contrast to the current interface of DIP
as it tries to solve subproblems without new columns using the MILP
formulation.<br>
Furthermore, the subproblems resemble each other for my problem.
This would allow to solve them in an accumulated fashion as far as
the reduced cost and branching decisions allow so. The latter is
certainly true for the root node of the branch-and-bound tree.
Hence, I would suggest to redesign the interface of DIP to enable a
solution process for <b>all the subproblems at once</b> if the user
provides it. Maybe the current treatment could act as a fallback.<br>
<br>
My approach here was to introduce the method <tt>solveRelaxedAll</tt>
for <tt>DecompApp </tt>and <tt>DecompAlgo</tt>. <tt>DecompAlgo::solveRelaxedAll
</tt>calls <tt>DecompApp::solveRelaxedAll</tt> and decides on the
return value whether it is actually implemented by the user (<tt>true</tt>)
or not (<tt>false</tt>). This result is propagated back to <tt>DecompAlgo::generateVars</tt>
which falls back on the original method if there is no
implementation. This design is not perfect at all since there are
large parts of code in <tt>DecompAlgo::solveRelaxedAll</tt> copied
from <tt>DecompAlgo::solveRelaxed.</tt><br>
<br>
<u><i>DecompApp.h</i></u><br>
<tt><br>
virtual bool solveRelaxedAll(const double* redCost, const
double* convexDuals, DecompVarList& varList,
std::vector<DecompSolverStatus>& states) {</tt><tt><br>
</tt><tt> return false;</tt><tt><br>
</tt><tt> }</tt><br>
<br>
<br>
<br>
<u><i>DecompAlgo.cpp</i></u><br>
<br>
<tt>bool DecompAlgo::solveRelaxedAll(const double* redCost,</tt><tt><br>
</tt><tt> const double* origCost,</tt><tt><br>
</tt><tt> const double* convexDuals,</tt><tt><br>
</tt><tt> const int n_origCols,</tt><tt><br>
</tt><tt> DecompSolverResult* solveResult,</tt><tt><br>
</tt><tt> DecompVarList& vars,</tt><tt><br>
</tt><tt> double timeLimit</tt><tt><br>
</tt><tt> )</tt><tt><br>
</tt><tt>{</tt><tt><br>
</tt><tt><br>
</tt><tt> UtilPrintFuncBegin(m_osLog, m_classTag,
"solveRelaxedAll()", m_param.LogDebugLevel, 2);</tt><tt><br>
</tt><tt><br>
</tt><tt> if (m_param.SubProbParallel) {</tt><tt><br>
</tt><tt> m_stats.timerOther1.reset();</tt><tt><br>
</tt><tt> } else {</tt><tt><br>
</tt><tt> m_stats.timerOther2.reset();</tt><tt><br>
</tt><tt> }</tt><tt><br>
</tt><tt><br>
</tt><b><tt> DecompVarList userVars;</tt></b><b><tt><br>
</tt></b><b><tt> vector<DecompSolverStatus>
solverStatus(m_numConvexCon, DecompSolStatNoSolution);</tt></b><b><tt><br>
</tt></b><b><tt><br>
</tt></b><b><tt> if (m_param.SolveRelaxAsIp != 1) {</tt></b><b><tt><br>
</tt></b><b><tt><br>
</tt></b><b><tt> if (!m_app->solveRelaxedAll(redCost,
convexDuals, userVars, solverStatus))</tt></b><b><tt><br>
</tt></b><b><tt> return false;</tt></b><b><tt><br>
</tt></b><b><tt><br>
</tt></b><b><tt> for (auto var : userVars)</tt></b><b><tt><br>
</tt></b><b><tt> if (var->getVarType() ==
DecompVar_Point)</tt></b><b><tt><br>
</tt></b><b><tt>
var->setReducedCost(var->getReducedCost() -
convexDuals[var->getBlockId()]);</tt></b><b><tt><br>
</tt></b><b><tt> }</tt></b><b><tt><br>
</tt></b><b><tt><br>
</tt></b><b><tt> m_isColGenExact = true;</tt></b><b><tt><br>
</tt></b><b><tt> for (int subprobIndex = 0; subprobIndex <
m_numConvexCon; subprobIndex++)</tt></b><b><tt><br>
</tt></b><b><tt> m_isColGenExact &=
solverStatus[subprobIndex] == DecompSolStatOptimal;</tt></b><b><tt><br>
</tt></b><b><tt><br>
</tt></b><b><tt> if ((m_isColGenExact || userVars.size() >
0) && (m_param.SolveRelaxAsIp != 2)) {</tt></b><b><tt><br>
</tt></b><b><tt> vars.splice(vars.end(), userVars);</tt></b><b><tt><br>
</tt></b><b><tt> } else {</tt></b><tt><br>
</tt><tt><br>
</tt><tt> for (int subprobIndex = 0; subprobIndex <
m_numConvexCon; subprobIndex++) {</tt><tt><br>
</tt><tt><br>
</tt><tt> DecompSubModel& subModel =
getModelRelax(subprobIndex);</tt><tt><br>
</tt><tt> OsiSolverInterface* subprobSI =
subModel.getOsi();</tt><tt><br>
</tt><tt> int whichBlock =
subModel.getBlockId();</tt><tt><br>
</tt><tt> bool isRoot = getNodeIndex() ?
false : true;</tt><tt><br>
</tt><tt> DecompConstraintSet* model =
subModel.getModel();</tt><tt><br>
</tt><tt><br>
</tt><tt> bool doCutoff = m_param.SubProbUseCutoff;</tt><tt><br>
</tt><tt> bool doExact = m_function ==
DecompFuncGenerateInitVars ? false : true;</tt><tt><br>
</tt><tt><br>
</tt><tt> assert(subprobSI);</tt><tt><br>
</tt><tt><br>
</tt><tt> subModel.setOsiObjCoeff(redCost);</tt><tt><br>
</tt><tt><br>
</tt><tt> if (m_param.BranchEnforceInSubProb) {</tt><tt><br>
</tt><tt> subModel.setActiveColBounds(m_colLBNode,
m_colUBNode);</tt><tt><br>
</tt><tt> }</tt><tt><br>
</tt><tt><br>
</tt><tt> if (m_param.LogDumpModel > 1) {</tt><tt><br>
</tt><tt> string baseName = "subProb_" +
subModel.getModelName();</tt><tt><br>
</tt><tt> if (m_isStrongBranch)</tt><tt><br>
</tt><tt> baseName += "_SB";</tt><tt><br>
</tt><tt> std::cout << "problem name is "
<< baseName << m_nodeStats.nodeIndex <<
m_nodeStats.cutCallsTotal << m_nodeStats.priceCallsTotal
<< whichBlock << std::endl;</tt><tt><br>
</tt><tt> printCurrentProblem(subprobSI,
baseName, m_nodeStats.nodeIndex, m_nodeStats.cutCallsTotal,
m_nodeStats.priceCallsTotal, whichBlock);</tt><tt><br>
</tt><tt> }</tt><tt><br>
</tt><tt><br>
</tt><tt> subModel.solveAsMIP(solveResult, m_param,
doExact, doCutoff, isRoot, convexDuals[subprobIndex] -
DecompEpsilon, timeLimit);</tt><tt><br>
</tt><tt><br>
</tt><tt> m_isColGenExact = solveResult->m_isOptimal;</tt><tt><br>
</tt><tt><br>
</tt><tt> if (solveResult->m_nSolutions) {</tt><tt><br>
</tt><tt> int k;</tt><tt><br>
</tt><tt> int nSol =
std::min<int>(solveResult->m_nSolutions,
m_param.SubProbNumSolLimit);</tt><tt><br>
</tt><tt> for (k = 0; k < nSol; k++) {</tt><tt><br>
</tt><tt> const double* milpSolution =
solveResult->getSolution(k);</tt><tt><br>
</tt><tt><br>
</tt><tt> vector<int> ind;</tt><tt><br>
</tt><tt> vector<double> els;</tt><tt><br>
</tt><tt> int i, c;</tt><tt><br>
</tt><tt> double varRedCost = 0.0;</tt><tt><br>
</tt><tt> double varOrigCost = 0.0;</tt><tt><br>
</tt><tt> DecompVarType varType =
!solveResult->m_isUnbounded ? DecompVar_Point : DecompVar_Ray;</tt><tt><br>
</tt><tt><br>
</tt><tt> if (model->isSparse()) {</tt><tt><br>
</tt><tt> map<int, int>::const_iterator
mcit;</tt><tt><br>
</tt><tt> const map<int, int>&
sparseToOrig = model->getMapSparseToOrig();</tt><tt><br>
</tt><tt><br>
</tt><tt> for (mcit = sparseToOrig.begin();</tt><tt><br>
</tt><tt> mcit != sparseToOrig.end();
mcit++) {</tt><tt><br>
</tt><tt> i = mcit->first;
//sparse-index</tt><tt><br>
</tt><tt> c = mcit->second;
//original-index</tt><tt><br>
</tt><tt><br>
</tt><tt> if
(!UtilIsZero(milpSolution[i], m_app->m_param.TolZero)) {</tt><tt><br>
</tt><tt> ind.push_back(c);</tt><tt><br>
</tt><tt>
els.push_back(milpSolution[i]);</tt><tt><br>
</tt><tt> varRedCost += redCost[c] *
milpSolution[i];</tt><tt><br>
</tt><tt> varOrigCost += origCost[c]
* milpSolution[i];</tt><tt><br>
</tt><tt> }</tt><tt><br>
</tt><tt> }</tt><tt><br>
</tt><tt> }</tt><tt><br>
</tt><tt> else {</tt><tt><br>
</tt><tt> for (c = 0; c < n_origCols; c++)
{</tt><tt><br>
</tt><tt> if
(!UtilIsZero(milpSolution[c], m_app->m_param.TolZero)) {</tt><tt><br>
</tt><tt> ind.push_back(c);</tt><tt><br>
</tt><tt>
els.push_back(milpSolution[c]);</tt><tt><br>
</tt><tt> varRedCost += redCost[c] *
milpSolution[c];</tt><tt><br>
</tt><tt> varOrigCost += origCost[c]
* milpSolution[c];</tt><tt><br>
</tt><tt> }</tt><tt><br>
</tt><tt> }</tt><tt><br>
</tt><tt> }</tt><tt><br>
</tt><tt><br>
</tt><tt> if (varType == DecompVar_Point) {</tt><tt><br>
</tt><tt> varRedCost -=
convexDuals[subprobIndex];</tt><tt><br>
</tt><tt> }</tt><tt><br>
</tt><tt><br>
</tt><tt> DecompVar* var = new DecompVar(ind,
els, varRedCost, varOrigCost, varType);</tt><tt><br>
</tt><tt> var->setBlockId(whichBlock);</tt><tt><br>
</tt><tt> vars.push_back(var);</tt><tt><br>
</tt><tt> }</tt><tt><br>
</tt><tt> }</tt><tt><br>
</tt><tt> }</tt><tt><br>
</tt><tt><br>
</tt><b><tt> }</tt></b><tt><br>
</tt><tt><br>
</tt><tt> if (!m_param.SubProbParallel) {</tt><tt><br>
</tt><tt>
m_stats.thisSolveRelax.push_back(m_stats.timerOther1.getRealTime());</tt><tt><br>
</tt><tt> }</tt><tt><br>
</tt><tt><br>
</tt><tt> UtilPrintFuncEnd(m_osLog, m_classTag,
"solveRelaxedAll()", m_param.LogDebugLevel, 2);</tt><tt><br>
</tt><tt><br>
</tt><tt> return true;</tt><tt><br>
</tt><tt>}<br>
</tt><u><i><br>
<br>
<br>
Decom</i></u><u><i>pAlgo.cpp:4763</i></u><tt><br>
<br>
<b> timeLimit = max(m_param.SubProbTimeLimitExact -
m_stats.timerOverall.getRealTime(), 0.0);</b><b><br>
</b><b> if (!solveRelaxedAll(redCostX, origObjective,
convexDuals, nCoreCols, &solveResult, potentialVars,
timeLimit)) {</b><br>
<br>
#ifdef _OPENMP<br>
UTIL_DEBUG(m_app->m_param.LogDebugLevel, 3,<br>
(*m_osLog) << "===== START Threaded solve of
subproblems. =====\n";<br>
);<br>
if (m_param.SubProbParallel) {<br>
omp_set_num_threads(min(m_param.NumConcurrentThreadsSubProb,
m_numConvexCon));<br>
}<br>
else {<br>
omp_set_num_threads(1);<br>
}<br>
#endif</tt><br>
<tt><br>
DecompVarList* potentialVarsT = new
DecompVarList[m_numConvexCon];<br>
CoinAssertHint(potentialVarsT, "Error: Out of
Memory");</tt><tt><br>
<br>
#pragma omp parallel for schedule(dynamic,
m_param.SubProbParallelChunksize) <br>
for (int subprobIndex = 0; subprobIndex <
m_numConvexCon; subprobIndex++) {<br>
<br>
DecompSubModel& subModel =
getModelRelax(subprobIndex);<br>
double alpha = u[nBaseCoreRows +
subprobIndex];<br>
DecompSolverResult solveResult(m_infinity);<br>
<br>
#ifdef _OPENMP<br>
UTIL_DEBUG(m_app->m_param.LogDebugLevel, 4,<br>
(*m_osLog) << "THREAD " <<
omp_get_thread_num() << " solving subproblem " <<
subprobIndex << "\n";);<br>
#else<br>
UTIL_DEBUG(m_app->m_param.LogDebugLevel, 4,<br>
(*m_osLog) << "solve relaxed model = "
<< subModel.getModelName() << endl;<br>
);<br>
#endif<br>
<br>
timeLimit = max(m_param.SubProbTimeLimitExact -
m_stats.timerOverall.getRealTime(), 0.0);<br>
solveRelaxed(redCostX, origObjective, alpha,
nCoreCols, false, subModel, &solveResult,
potentialVarsT[subprobIndex], timeLimit);<br>
<br>
if (solveResult.m_isCutoff) {<br>
mostNegRCvec[subprobIndex] =
min(mostNegRCvec[subprobIndex], 0.0);<br>
}<br>
}<br>
<br>
for (int subprobIndex = 0; subprobIndex <
m_numConvexCon; subprobIndex++) {<br>
/* printf("arg[%d].vars size=%d\n",<br>
t,
static_cast<int>(arg[t].vars->size()));<br>
*/<br>
for (it = potentialVarsT[subprobIndex].begin();<br>
it != potentialVarsT[subprobIndex].end(); it++)
{<br>
varRedCost = (*it)->getReducedCost();<br>
whichBlock = (*it)->getBlockId();<br>
<br>
if ((*it)->getVarType() == DecompVar_Point) {<br>
alpha = u[nBaseCoreRows + whichBlock];<br>
} else if ( (*it)->getVarType() ==
DecompVar_Ray) {<br>
alpha = 0;<br>
}<br>
<br>
UTIL_DEBUG(m_app->m_param.LogDebugLevel, 3,<br>
(*m_osLog)<br>
<< "alpha[block=" << whichBlock
<< "]:" << alpha<br>
<< " varRedCost: " <<
varRedCost << "\n";<br>
);<br>
}<br>
}<br>
<br>
#ifdef _OPENMP<br>
UTIL_DEBUG(m_app->m_param.LogDebugLevel, 3,<br>
(*m_osLog)<br>
<< "===== END Threaded solve of
subproblems. =====\n";);<br>
#endif<br>
<br>
for (int subprobIndex = 0; subprobIndex <
m_numConvexCon; subprobIndex++) {<br>
for (it = potentialVarsT[subprobIndex].begin(); it
!= potentialVarsT[subprobIndex].end(); it++) {<br>
potentialVars.push_back(*it);<br>
}<br>
}<br>
<br>
UTIL_DELARR(potentialVarsT);<br>
<b> }</b><b><br>
</b><b><br>
</b><b> UTIL_DELARR(convexDuals);</b><br>
<br>
</tt>
<hr size="2" width="100%"><br>
[<i>Functionality</i>] Based on the fact that my algorithm for the
subproblem is partly heuristic (see previous paragraphs), it is not
until the last iterations for solving a single node that it provides
not only feasible (<tt>DecompSolStatFeasible</tt>) but optimal
solutions (<tt>DecompSolStatOptimal</tt>). Thus, a lower bound for
the node can only be computed at the end of the solution process of
each node. This prevents me from using the <b>tailing-off mechanism</b>
provided by DIP (<tt>DecompAlgo::isTailoffLB</tt>) as it is based on
the lower bound. For that reason, I think it would be useful to
introduce an alternative tailing-off control based on the
progression of the upper bound for the relaxed problem. Would this
be a reasonable approach?<br>
<br>
<hr size="2" width="100%"><br>
<i>[Implementation/Functionality]</i> If a node is solved to
optimality in the pricing phase (<tt>PHASE_PRICE2</tt>) no more
columns are generated and the algorithm switches to <tt>PHASE_CUT</tt>
(cf. <tt>DecompAlgo::phaseUpdate</tt>, line 4215). It prevents
stopping on a tight gap as checked in <tt>DecompAlgo::phaseUpdate</tt>,
line 4274. Thus, the switch to the cutting phase it carried out.
However, there is a parameter called <tt>PCStrategy</tt>. Setting
it to<tt> favorPrice</tt> will make <tt>DecompAlgo::phaseUpdate</tt>
to immediately switch the phase from <tt>PHASE_CUT</tt> to <tt>PHASE_PRICE2</tt>
back again in line 4153 as long as the limit on the number of price
calls is not reached. This can result in alternation of the two
phases and hence an<b> infinite loop</b>. I found the remedy of
setting the <tt>RoundCutItersLimit</tt> to 0, which probably suits
my intention better. Yet, I wonder what the actual use of the <tt>PCStrategy</tt>
parameter is. At the moment it seems to be redundant as the <tt>RoundPriceItersLimit</tt>
and <tt>RoundCutItersLimit </tt>are the controlling paramerters.<br>
<br>
<hr size="2" width="100%"><br>
<i>[Performance/Implementation]</i> I found that the <b>checks for
duplicate and parallel columns</b> in DecompApp::addVarsToPool,
lines 5565 sqq. are quite expensive. First of all, I believe that
the check for parallel columns in line 5609 is redundant if the
parameter <tt>ParallelColsLimit</tt> is 1.0 (as pointed out in the
comment preceeding that line of code). Since this is the default
value of the parameter, I recommend checking for parallel columns
only if the parameter is smaller than 1.0. Apart from that in my
understanding of column generation, columns with negative reduced
cost cannot have been included before. This would render it
unnecessary to check for duplicates in the existing columns. Maybe
this is not true for some configurations like <tt>DualStab</tt>? As
mentioned in the comments in the code, the hashing of the columns is
not efficient at all. A comparison based on this hashes is even more
expensive than a direct comparison of the sparse representation. A
better hashing and conditional exact comparison would result in a
noticable speed-up, I suppose.<br>
<br>
<u><i>DecompAlgo.cpp:5609<br>
<br>
</i></u><tt> if (foundGoodCol </tt><tt><b>&&
m_param.ParallelColsLimit < 1.0</b></tt><tt> && </tt><tt><br>
</tt><tt> m_varpool.isParallel(m_vars, waitingCol,
m_param.ParallelColsLimit)) {</tt><tt><br>
</tt><tt> UTIL_DEBUG(m_app->m_param.LogDebugLevel, 3,</tt><tt><br>
</tt><tt> (*m_osLog) << "Parallel variable,
already in vars.\n";</tt><tt><br>
</tt><tt> );</tt><tt><br>
</tt><tt> waitingCol.deleteVar();</tt><tt><br>
</tt><tt> waitingCol.deleteCol();</tt><tt><br>
</tt><tt><br>
</tt><tt> if (m_algo != RELAX_AND_CUT) { //??</tt><tt><br>
</tt><tt> m_nodeStats.varsThisCall--;</tt><tt><br>
</tt><tt> m_nodeStats.varsThisRound--;</tt><tt><br>
</tt><tt> }</tt><tt><br>
</tt><tt><br>
</tt><tt> continue;</tt><tt><br>
</tt><tt> }<br>
<br>
<br>
<br>
</tt><u><i>DecompVarPool.cp</i></u><u><i>p:</i></u><u><i>150</i></u><tt><br>
<br>
bool DecompVarPool::isDuplicate(const DecompVarList& vars,<br>
const DecompWaitingCol& wcol)<br>
{<br>
const DecompVar* var = wcol.getVarPtr();<br>
const int block = var->getBlockId();<br>
const int len = var->m_s.getNumElements();<br>
const int* indices = var->m_s.getIndices();<br>
const double* values = var->m_s.getElements();<br>
// const string hash = var->getStrHash();<br>
<br>
for (DecompVarList::const_iterator vi = vars.begin(); vi !=
vars.end(); vi++) {<br>
<br>
if ((*vi)->getBlockId() != block)<br>
continue;<br>
<br>
// Could be reasonalbe if hashing is lightweight and yields an
adequate number of "buckets"<br>
// if ((*vi)->getStrHash() != hash)<br>
// continue;<br>
<br>
if ((*vi)->m_s.getNumElements() != len)<br>
continue;<br>
<br>
const int* other_indices = (*vi)->m_s.getIndices();<br>
const double* other_values = (*vi)->m_s.getElements();<br>
<br>
for (int i = 0; i < len; i++)<br>
if (other_indices[i] != indices[i] || other_values[i]
!= values[i])<br>
goto next;<br>
<br>
return true;<br>
<br>
next: {}<br>
}<br>
<br>
return false;<br>
}<br>
<br>
bool DecompVarPool::isDuplicate(const DecompWaitingCol& wcol)<br>
{<br>
</tt><tt><tt>// Repetition of code from above. Refactor?<br>
<br>
</tt> const DecompVar* var = wcol.getVarPtr();<br>
const int block = var->getBlockId();<br>
const int len = var->m_s.getNumElements();<br>
const int* indices = var->m_s.getIndices();<br>
const double* values = var->m_s.getElements();<br>
// const string hash = var->getStrHash();<br>
<br>
for (auto vi = begin(); vi != end(); vi++) {<br>
<br>
const DecompVar* other = (*vi).getVarPtr();<br>
</tt><tt><tt><br>
</tt> if (other->getBlockId() != block)<br>
continue;<br>
</tt><tt><br>
// Could be reasonalbe if hashing is lightweight and yields an
adequate number of "buckets"</tt><tt><br>
// if (other->getStrHash() != hash)<br>
// continue;<br>
<br>
if (other->m_s.getNumElements() != len)<br>
continue;<br>
<br>
const int* other_indices = other->m_s.getIndices();<br>
const double* other_values = other->m_s.getElements();<br>
<br>
for (int i = 0; i < len; i++)<br>
if (other_indices[i] != indices[i] || other_values[i]
!= values[i])<br>
goto next;<br>
<br>
return true;<br>
<br>
next: {}<br>
}<br>
<br>
return false;<br>
}<br>
</tt> <br>
<hr size="2" width="100%"><br>
<i>[Question]</i> Finally, I did not understand the usage of the
parameters <tt>BranchEnforceInMaster</tt> and <tt>BranchEnforceInSubProb</tt>
and would be grateful if you could explain the basic meaning of them
to me. They seem to control the treatment of the branching
decisions. <tt>BranchEnforceInMaster</tt> suggests to include the
decisions as new rows in the master problem and enforce them via
reduced costs, when <tt>BranchEnforceInSubProb</tt> suggests to
make it the subproblems' task to deal with them. In the latter case,
what happens to conditions which include master-only variables and
therefore cannot be treated in the subproblems? What is the best way
for getting the branching decisions for the current node when using
<tt>BranchEnforceInSubProb</tt>?<br>
<br>
<hr size="2" width="100%"><br>
<br>
<br>
<br>
<br>
<br>
<br>
<br>
<br>
<br>
<br>
<br>
<div class="moz-cite-prefix">Am 27.03.2016 um 00:53 schrieb Ted
Ralphs:<br>
</div>
<blockquote
cite="mid:CA+GYycvOXJphnpy9ZXf2ABL_OH_8-ze-H_O-1-MrrrTKQwvieA@mail.gmail.com"
type="cite">
<div dir="ltr">Hi Marcus,
<div><br>
</div>
<div>Thanks very much for the detailed feedback! Sorry for the
delay in responding. I had a quick look over your comments and
they all seem reasonable. The implementation issues look like
they can be fixed fairly quickly and easily. For the others,
we need to have a more detailed look to see what can be done
in the short and long run. In the next couple of days, I'll
try to respond with detailed comments on each of your issues.
Of course, if you want to provide patches for any of these,
that would be most welcome :).</div>
<div><br>
</div>
<div>Cheers,</div>
<div><br>
</div>
<div>Ted</div>
<div><br>
</div>
<div><br>
</div>
</div>
<div class="gmail_extra"><br>
<div class="gmail_quote">On Mon, Mar 21, 2016 at 9:59 AM, Marcus
Kaiser <span dir="ltr"><<a moz-do-not-send="true"
href="mailto:marcus.kaiser@mytum.de" target="_blank">marcus.kaiser@mytum.de</a>></span>
wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0
.8ex;border-left:1px #ccc solid;padding-left:1ex">
<div text="#000000" bgcolor="#FFFFFF"> Hello Dip/py
community,<br>
in my master's thesis im dealing with a Dantzig-Wolfe
decomposition and the resulting column generation. For the
implementation I use the DIP framework, version 0.92.2 on
Windows 7 with MSVC 2013. I appreciate that you provide
such an extensive framework. While working with it and
debbuging my code the following things came to my mind.
Some of them are performance and implementation related
others address the functionality of the algorithms in DIP.
Additionally, I have a question about the treatment of
branching decisions.<br>
<br>
<i>[</i><i><i>Implementation</i>]</i> I found Dippy to
suffer from <b>memory leaking</b>. Since it interfaces
with Python, it is responsible for maintaning the
reference counters of the Python objects it deals with.
This is done in some parts of the code, but by far not
all. A crucial situation which came to my attention is the
method <tt>DippyDecompApp::solveRelaxed</tt>. The
retrieved columns are converted to objects of the class <tt>DecompVar</tt>,
yet the reference counter on the original Python object is
not decreased, which prevents them from being deleted.
Hence, for large problem instances the memory floods.<br>
<br>
<i>[Implementation]</i> The method <tt>AlpsDecompTreeNode::getBranchedVar
</tt>seems to be meant to return the <b>variable which is
branched on</b>. However, only one of the four vectors
potentially holding branching decisions is checked. I
think both, lower bounds and upper bounds, need to be
checked for at least one branch, e.g. <tt>downBranchLB</tt>_
and <tt>downBranchUB_</tt>.<br>
<br>
<i>[</i><i>Implementation</i>] In <tt>DecompAlgo::processNode</tt>,
line 1855 the processing of a node is terminated if the <b>lower
bound meets the global upper bound</b>. It takes into
account numeric inexactness. The calling method <tt>AlpsDecompTreeNode::process
</tt>repeats this check in line 309 but in an exact
fashion why the node might not be fathomed when it should
be. In my oppinion this repeated condition should match
the original one. Or maybe even better, <tt>DecompAlgo::getStopCriteria</tt>
could be used?<br>
<br>
<i>[Implementation]</i> The member <tt>DecompAlgo::m_relGap</tt>
is set in <tt>DecompAlgo::updateObjBound</tt> and used in
<tt>DecompAlgo::isGapTight</tt>. <tt>DecompAlgo::m_relGap</tt>,
however, is not reset when entering <tt>DecompAlgo::processNode</tt>.
Therefore, it has an <b>invalid value</b> based on a node
processed before - probably representing a tight <b>gap</b>.
I think, this might lead to stopping the processing of a
node immediately as it is mistakenly believed to have a
tight gap, cf. <tt>DecompAlgo::phaseUpdate</tt>, line
4274. I suggest to reset <tt>DecompAlgo::m_relGap </tt>appropriately
or replace it completely by calls to <tt>DecompAlgo::getNodeLPGap</tt>.<br>
<br>
<i>[Functionality/Performance]</i> For the problem I
consider, the MILP formulation of the subproblems is quite
huge. Therefore, I use a dummy formulation and prevent DIP
from looking at it by always providing columns with
negative reduced cost via my implementation of <tt>DecompApp::solveRelaxed</tt>
if they exist. I do so in a two-step approach. In a first
step the subproblem is solved heuristically. If this does
not lead to new columns an exact method is applied. Since
the exact method is costly, I would like to stick to the
heuristic as long as there are new columns for <i>some</i>
subproblem. This stands in contrast to the current
interface of DIP as it tries to solve subproblems without
new columns using the MILP formulation.<br>
Furthermore, the subproblems resemble each other for my
problem. This would allow to solve them in an accumulated
fashion as far as the reduced cost and branching decisions
allow so. The latter is certainly true for the root node
of the branch-and-bound tree. Hence, I would suggest to
redesign the interface of DIP to enable a solution process
for <b>all the subproblems at once</b> if the user
provides it. Maybe the current treatment could act as a
fallback.<br>
<br>
[<i>Functionality</i>] Based on the fact that my algorithm
for the subproblem is partly heuristic (see previous
paragraphs), it is not until the last iterations for
solving a single node that it provides not only feasible (<tt>DecompSolStatFeasible</tt>)
but optimal solutions (<tt>DecompSolStatOptimal</tt>).
Thus, a lower bound for the node can only be computed at
the end of the solution process of each node. This
prevents me from using the <b>tailing-off mechanism</b>
provided by DIP (<tt>DecompAlgo::isTailoffLB</tt>) as it
is based on the lower bound. For that reason, I think it
would be useful to introduce an alternative tailing-off
control based on the progression of the upper bound for
the relaxed problem. Would this be a reasonable approach?<br>
<br>
<i>[Implementation/Functionality]</i> If a node is solved
to optimality in the pricing phase (<tt>PHASE_PRICE2</tt>)
no more columns are generated and the algorithm switches
to <tt>PHASE_CUT</tt> (cf. <tt>DecompAlgo::phaseUpdate</tt>,
line 4215). It prevents stopping on a tight gap as checked
in <tt>DecompAlgo::phaseUpdate</tt>, line 4274. Thus, the
switch to the cutting phase it carried out. However, there
is a parameter called <tt>PCStrategy</tt>. Setting it to<tt>
favorPrice</tt> will make <tt>DecompAlgo::phaseUpdate</tt>
to immediately switch the phase from <tt>PHASE_CUT</tt>
to <tt>PHASE_PRICE2</tt> back again in line 4153 as long
as the limit on the number of price calls is not reached.
This can result in alternation of the two phases and hence
an<b> infinite loop</b>. I found the remedy of setting the
<tt>RoundCutItersLimit</tt> to 0, which probably suits my
intention better. Yet, I wonder what the actual use of the
<tt>PCStrategy</tt> parameter is. At the moment it seems
to be redundant as the <tt>RoundPriceItersLimit</tt> and
<tt>RoundCutItersLimit </tt>are the controlling
paramerters.<br>
<br>
<i>[Performance/Implementation]</i> I found that the <b>checks
for duplicate and parallel columns</b> in
DecompApp::addVarsToPool, lines 5565 sqq. are quite
expensive. First of all, I believe that the check for
parallel columns in line 5609 is redundant if the
parameter <tt>ParallelColsLimit</tt> is 1.0 (as pointed
out in the comment preceeding that line of code). Since
this is the default value of the parameter, I recommend
checking for parallel columns only if the parameter is
smaller than 1.0. Apart from that in my understanding of
column generation, columns with negative reduced cost
cannot have been included before. This would render it
unnecessary to check for duplicates in the existing
columns. Maybe this is not true for some configurations
like <tt>DualStab</tt>? As mentioned in the comments in
the code, the hashing of the columns is not efficient at
all. A comparison based on this hashes is even more
expensive than a direct comparison of the sparse
representation. A better hashing and conditional exact
comparison would result in a noticable speed-up, I
suppose.<br>
<br>
<i>[Question]</i> Finally, I did not understand the usage
of the parameters <tt>BranchEnforceInMaster</tt> and <tt>BranchEnforceInSubProb</tt>
and would be grateful if you could explain the basic
meaning of them to me. They seem to control the treatment
of the branching decisions. <tt>BranchEnforceInMaster</tt>
suggests to include the decisions as new rows in the
master problem and enforce them via reduced costs, when <tt>BranchEnforceInSubProb</tt>
suggests to make it the subproblems' task to deal with
them. In the latter case, what happens to conditions which
include master-only variables and therefore cannot be
treated in the subproblems? What is the best way for
getting the branching decisions for the current node when
using <tt>BranchEnforceInSubProb</tt>?<br>
<br>
Thank you in advance,<br>
Marcus<br>
</div>
<br>
_______________________________________________<br>
Dip mailing list<br>
<a moz-do-not-send="true" href="mailto:Dip@list.coin-or.org">Dip@list.coin-or.org</a><br>
<a moz-do-not-send="true"
href="http://list.coin-or.org/mailman/listinfo/dip"
rel="noreferrer" target="_blank">http://list.coin-or.org/mailman/listinfo/dip</a><br>
</blockquote>
</div>
<br>
<br clear="all">
<div><br>
</div>
-- <br>
<div class="gmail_signature">
<div dir="ltr">Dr. Ted Ralphs<br>
Professor, Lehigh University<br>
(610) 628-1280<br>
ted 'at' lehigh 'dot' edu<br>
<a moz-do-not-send="true"
href="http://coral.ie.lehigh.edu/%7Eted" target="_blank">coral.ie.lehigh.edu/~ted</a><br>
</div>
</div>
</div>
</blockquote>
<br>
</body>
</html>