OpenGM  2.3.x
Discrete Graphical Model Library
multicut.hxx
Go to the documentation of this file.
1 #pragma once
2 #ifndef OPENGM_MULTICUT_HXX
3 #define OPENGM_MULTICUT_HXX
4 
5 #include <algorithm>
6 #include <vector>
7 #include <queue>
8 #include <utility>
9 #include <string>
10 #include <iostream>
11 #include <fstream>
12 #include <typeinfo>
13 #include <limits>
14 #ifdef WITH_BOOST
15 #include <boost/unordered_map.hpp>
16 #include <boost/unordered_set.hpp>
17 #else
18 #include <ext/hash_map>
19 #include <ext/hash_set>
20 #endif
21 
22 #include "opengm/datastructures/marray/marray.hxx"
23 #include "opengm/opengm.hxx"
24 #include "opengm/inference/inference.hxx"
25 #include "opengm/inference/visitors/visitors.hxx"
26 #include "opengm/utilities/timer.hxx"
27 #include "opengm/utilities/queues.hxx"
28 #include "opengm/utilities/partitions.hxx"
29 
30 #include <ilcplex/ilocplex.h>
31 //ILOSTLBEGIN
32 
33 namespace opengm {
34 
36 class HigherOrderTerm
37 {
38 public:
39  size_t factorID_;
40  bool potts_;
41  size_t valueIndex_;
42  std::vector<size_t> lpIndices_;
43  HigherOrderTerm(size_t factorID, bool potts, size_t valueIndex)
44  : factorID_(factorID), potts_(potts),valueIndex_(valueIndex) {}
45  HigherOrderTerm()
46  : factorID_(0), potts_(false),valueIndex_(0) {}
47 };
49 
72 struct ParamHeper{
73 enum MWCRounding {NEAREST,DERANDOMIZED,PSEUDODERANDOMIZED};
74 };
75 
76 template<class GM, class ACC>
77 class Multicut : public Inference<GM, ACC>
78 {
79 public:
80  typedef ACC AccumulationType;
81  typedef GM GraphicalModelType;
82  OPENGM_GM_TYPE_TYPEDEFS;
83  typedef size_t LPIndexType;
84  typedef visitors::VerboseVisitor<Multicut<GM,ACC> > VerboseVisitorType;
85  typedef visitors::EmptyVisitor<Multicut<GM,ACC> > EmptyVisitorType;
86  typedef visitors::TimingVisitor<Multicut<GM,ACC> > TimingVisitorType;
87 
88 
89 #ifdef WITH_BOOST
90  typedef boost::unordered_map<IndexType, LPIndexType> EdgeMapType;
91  typedef boost::unordered_set<IndexType> MYSET;
92 #else
93  typedef __gnu_cxx::hash_map<IndexType, LPIndexType> EdgeMapType;
94  typedef __gnu_cxx::hash_set<IndexType> MYSET;
95 #endif
96 
97 
98  template<class GM_, class ACC_>
99  struct rebind{
100  typedef Multicut<GM_, ACC_> type;
101  };
102 
103  struct Parameter : public ParamHeper{
104  public:
105 
106 
107  int numThreads_;
108  bool verbose_;
109  bool verboseCPLEX_;
110  double cutUp_;
111  double timeOut_;
112  std::string workFlow_;
113  size_t maximalNumberOfConstraintsPerRound_;
114  double edgeRoundingValue_;
115  ParamHeper::MWCRounding MWCRounding_;
116  size_t reductionMode_;
117  std::vector<bool> allowCutsWithin_;
118  bool useOldPriorityQueue_;
119  bool useChordalSearch_;
120  bool useBufferedStates_;
121  bool initializeWith3Cycles_;
122 
125  Parameter
126  (
127  int numThreads=0,
128  double cutUp=1.0e+75
129  )
130  : numThreads_(numThreads), verbose_(false),verboseCPLEX_(false), cutUp_(cutUp),
131  timeOut_(36000000), maximalNumberOfConstraintsPerRound_(1000000),
132  edgeRoundingValue_(0.00000001),MWCRounding_(NEAREST), reductionMode_(3),useOldPriorityQueue_(false), useChordalSearch_(false), useBufferedStates_(false),
133  initializeWith3Cycles_(false)
134  {};
135 
136  template<class OTHER_PARAM>
137  Parameter
138  (
139  const OTHER_PARAM & p
140  )
141  : numThreads_(p.numThreads_), verbose_(p.verbose_),verboseCPLEX_(p.verboseCPLEX_), cutUp_(p.cutUp_),
142  timeOut_(p.timeOut_), maximalNumberOfConstraintsPerRound_(p.maximalNumberOfConstraintsPerRound_),
143  edgeRoundingValue_(p.edgeRoundingValue_),MWCRounding_(p.MWCRounding_), reductionMode_(p.reductionMode_),
144  useOldPriorityQueue_(p.useOldPriorityQueue_), useChordalSearch_(p.useChordalSearch_),
145  initializeWith3Cycles_(false)
146  {};
147  };
148 
149  virtual ~Multicut();
150  Multicut(const GraphicalModelType&, Parameter para=Parameter());
151  Multicut(const size_t, const std::map<UInt64Type, ValueType> & accWeights, const Parameter & para=Parameter());
152 
153 
154  virtual std::string name() const {return "Multicut";}
155  const GraphicalModelType& graphicalModel() const;
156  virtual InferenceTermination infer();
157  template<class VisitorType> InferenceTermination infer(VisitorType&);
158  virtual InferenceTermination arg(std::vector<LabelType>&, const size_t = 1) const;
159  ValueType bound() const;
160  ValueType value() const;
161  ValueType calcBound(){ return 0; }
162  ValueType evaluate(std::vector<LabelType>&) const;
163 
164  template<class LPVariableIndexIterator, class CoefficientIterator>
165  void addConstraint(LPVariableIndexIterator, LPVariableIndexIterator,
166  CoefficientIterator, const ValueType&, const ValueType&);
167  std::vector<double> getEdgeLabeling() const;
168  std::vector<size_t> getSegmentation() const;
169 
170  template<class IT>
171  size_t getLPIndex(IT a, IT b) { return neighbours[a][b]; };
172 
173  size_t inferenceState_;
174  size_t constraintCounter_;
175 private:
176  enum ProblemType {INVALID, MC, MWC};
177 
178  const GraphicalModelType& gm_;
179  ProblemType problemType_;
180  Parameter parameter_;
181  double constant_;
182  double bound_;
183  double bufferedValue_;
184  std::vector<LabelType> bufferedStates_;
185  const double infinity_;
186  LabelType numberOfTerminals_;
187  IndexType numberOfNodes_;
188  LPIndexType numberOfTerminalEdges_;
189  LPIndexType numberOfInternalEdges_;
190  LPIndexType terminalOffset_;
191  LPIndexType numberOfHigherOrderValues_;
192  LPIndexType numberOfInterTerminalEdges_;
193 
194  std::vector<std::vector<size_t> > workFlow_;
195  std::vector<std::pair<IndexType,IndexType> > edgeNodes_;
196 
199  std::vector<EdgeMapType > neighbours;
200 
201  IloEnv env_;
202  IloModel model_;
203  IloNumVarArray x_;
204  IloRangeArray c_;
205  IloObjective obj_;
206  IloNumArray sol_;
207  IloCplex cplex_;
208 
209  bool integerMode_;
210  const double EPS_; //small number: for numerical issues constraints are still valid if the not up to EPS_
211  Partitions<size_t,LabelType> P_;
212 
213  void initCplex();
214 
215  size_t findCycleConstraints(IloRangeArray&, bool = true, bool = true);
216  size_t findIntegerCycleConstraints(IloRangeArray&, bool = true);
217  size_t findTerminalTriangleConstraints(IloRangeArray&);
218  size_t findIntegerTerminalTriangleConstraints(IloRangeArray&, std::vector<LabelType>& conf);
219  size_t findMultiTerminalConstraints(IloRangeArray&);
220  size_t findOddWheelConstraints(IloRangeArray&);
221  size_t removeUnusedConstraints(); //TODO: implement
222  size_t enforceIntegerConstraints();
223  size_t add3CycleConstraints();
224 
225  bool readWorkFlow(std::string);
226 
227  InferenceTermination partition(std::vector<LabelType>&, std::vector<std::list<size_t> >&, double = 0.5) const;
228  ProblemType setProblemType();
229  LPIndexType getNeighborhood(const LPIndexType, std::vector<EdgeMapType >&,std::vector<std::pair<IndexType,IndexType> >&, std::vector<HigherOrderTerm>&);
230 
231  template <class DOUBLEVECTOR>
232  double shortestPath(const IndexType, const IndexType, const std::vector<EdgeMapType >&, const DOUBLEVECTOR&, std::vector<IndexType>&, const double = std::numeric_limits<double>::infinity(), bool = true) const;
233  template <class DOUBLEVECTOR>
234  double shortestPath2(const IndexType, const IndexType, const std::vector<EdgeMapType >&, const DOUBLEVECTOR&, std::vector<IndexType>&,
235  std::vector<IndexType>&, opengm::ChangeablePriorityQueue<double>&,
236  const double = std::numeric_limits<double>::infinity(), bool = true) const;
237 
238  InferenceTermination derandomizedRounding(std::vector<LabelType>&) const;
239  InferenceTermination pseudoDerandomizedRounding(std::vector<LabelType>&, size_t = 1000) const;
240  double derandomizedRoundingSubProcedure(std::vector<LabelType>&,const std::vector<LabelType>&, const double) const;
241 
242  //PROTOCOLATION
243 
244  enum{
245  Protocol_ID_Solve = 0,
246  Protocol_ID_AddConstraints = 1,
247  Protocol_ID_RemoveConstraints = 2,
248  Protocol_ID_IntegerConstraints = 3,
249  Protocol_ID_CC = 4,
250  Protocol_ID_TTC = 5,
251  Protocol_ID_MTC = 6,
252  Protocol_ID_OWC = 7,
253  Protocol_ID_Unknown = 8
254  };
255 
256  enum{
257  Action_ID_RemoveConstraints = 0,
258  Action_ID_IntegerConstraints = 1,
259  Action_ID_CC = 10,
260  Action_ID_CC_I = 11,
261  Action_ID_CC_IFD = 12,
262  Action_ID_CC_FD = 13,
263  Action_ID_CC_B = 14,
264  Action_ID_CC_FDB = 15,
265  Action_ID_TTC = 20,
266  Action_ID_TTC_I = 21,
267  Action_ID_MTC = 30,
268  Action_ID_OWC = 40
269  };
270 
271  std::vector<std::vector<double> > protocolateTiming_;
272  std::vector<std::vector<size_t> > protocolateConstraints_;
273 
274 };
275 
276 
277 
278 template<class GM, class ACC>
279 Multicut<GM, ACC>::Multicut
280 (
281  const size_t numNodes,
282  const std::map<UInt64Type, ValueType> & accWeights,
283  const Parameter & para
284  ) : gm_(GM()), parameter_(para) , bound_(-std::numeric_limits<double>::infinity()), infinity_(1e8), integerMode_(false),
285  EPS_(1e-8)
286 {
287  if(typeid(ACC) != typeid(opengm::Minimizer) || typeid(OperatorType) != typeid(opengm::Adder)) {
288  throw RuntimeError("This implementation does only supports Min-Plus-Semiring.");
289  }
290  if(parameter_.reductionMode_<0 ||parameter_.reductionMode_>3) {
291  throw RuntimeError("Reduction Mode has to be 1, 2 or 3!");
292  }
293 
294  //Set Problem Type
295  problemType_ = MC;
296  numberOfTerminalEdges_ = 0;
297  numberOfTerminals_ = 0;
298  numberOfInterTerminalEdges_ = 0;
299  numberOfHigherOrderValues_ = 0;
300  numberOfNodes_ = numNodes;
301  size_t numEdges = accWeights.size();
302  //Calculate Neighbourhood
303  neighbours.resize(numberOfNodes_);
304  numberOfInternalEdges_=0;
305  LPIndexType numberOfAdditionalInternalEdges=0;
306 
307  if(para.useBufferedStates_){
308  bufferedValue_ = std::numeric_limits<double>::infinity();
309  bufferedStates_.resize(numNodes,0);
310  }
311 
312 
313  typedef std::map<IndexType, ValueType> MapType;
314  typedef typename MapType::const_iterator MapIter;
315 
316  // cplex stuff
317  IloInt N = numEdges;
318  model_ = IloModel(env_);
319  x_ = IloNumVarArray(env_);
320  c_ = IloRangeArray(env_);
321  obj_ = IloMinimize(env_);
322  sol_ = IloNumArray(env_,N);
323  // set variables and objective
324  x_.add(IloNumVarArray(env_, N, 0, 1, ILOFLOAT));
325 
326  IloNumArray obj(env_,N);
327  // add edges
328 
329 
330  for(MapIter i = accWeights.begin(); i!=accWeights.end(); ++i){
331  const UInt64Type key = i->first;
332  const ValueType weight = i->second;
333  const UInt64Type u = key/numberOfNodes_;
334  const UInt64Type v = key - u*numberOfNodes_;
335  if(neighbours[u].find(v)==neighbours[u].end()) {
336  neighbours[u][v] = numberOfInternalEdges_;
337  neighbours[v][u] = numberOfInternalEdges_;
338  edgeNodes_.push_back(std::pair<IndexType,IndexType>(v,u));
339  obj[numberOfInternalEdges_] = weight;
340  ++numberOfInternalEdges_;
341 
342  }
343  else{
344  OPENGM_CHECK_OP(true,==,false,"");
345  }
346  }
347 
348  obj_.setLinearCoefs(x_,obj);
349  model_.add(obj_);
350  if(para.initializeWith3Cycles_){
351  add3CycleConstraints();
352  }
353  // initialize solver
354  cplex_ = IloCplex(model_);
355 }
356 
357 
358 
359 template<class GM, class ACC>
360 typename Multicut<GM, ACC>::LPIndexType Multicut<GM, ACC>::getNeighborhood
361 (
362  const LPIndexType numberOfTerminalEdges,
363  std::vector<EdgeMapType >& neighbours,
364  std::vector<std::pair<IndexType,IndexType> >& edgeNodes,
365  std::vector<HigherOrderTerm>& higherOrderTerms
366  )
367 {
368  //Calculate Neighbourhood
369  neighbours.resize(gm_.numberOfVariables());
370  LPIndexType numberOfInternalEdges=0;
371  LPIndexType numberOfAdditionalInternalEdges=0;
372  // Add edges that have to be included
373  for(size_t f=0; f<gm_.numberOfFactors(); ++f) {
374  if(gm_[f].numberOfVariables()==2) { // Second Order Potts
375  IndexType u = gm_[f].variableIndex(1);
376  IndexType v = gm_[f].variableIndex(0);
377  if(neighbours[u].find(v)==neighbours[u].end()) {
378  neighbours[u][v] = numberOfTerminalEdges+numberOfInternalEdges;
379  neighbours[v][u] = numberOfTerminalEdges+numberOfInternalEdges;
380  edgeNodes.push_back(std::pair<IndexType,IndexType>(v,u));
381  ++numberOfInternalEdges;
382  }
383  }
384  }
385  for(size_t f=0; f<gm_.numberOfFactors(); ++f) {
386  if(gm_[f].numberOfVariables()>2 && !gm_[f].isPotts()){ // Generalized Potts
387  higherOrderTerms.push_back(HigherOrderTerm(f, false, 0));
388  for(size_t i=0; i<gm_[f].numberOfVariables();++i) {
389  for(size_t j=0; j<i;++j) {
390  IndexType u = gm_[f].variableIndex(i);
391  IndexType v = gm_[f].variableIndex(j);
392  if(neighbours[u].find(v)==neighbours[u].end()) {
393  neighbours[u][v] = numberOfTerminalEdges+numberOfInternalEdges;
394  neighbours[v][u] = numberOfTerminalEdges+numberOfInternalEdges;
395  edgeNodes.push_back(std::pair<IndexType,IndexType>(v,u));
396  ++numberOfInternalEdges;
397  ++numberOfAdditionalInternalEdges;
398  }
399  }
400  }
401  }
402  }
403  //Add for higher order potts term only neccesary edges
404  for(size_t f=0; f<gm_.numberOfFactors(); ++f) {
405  if(gm_[f].numberOfVariables()>2 && gm_[f].isPotts()) { //Higher order Potts
406  higherOrderTerms.push_back(HigherOrderTerm(f, true, 0));
407  std::vector<LPIndexType> lpIndexVector;
408  //Find spanning tree vor the variables nb(f) using edges that already exist.
409  std::vector<bool> variableInSpanningTree(gm_.numberOfVariables(),true);
410  for(size_t i=0; i<gm_[f].numberOfVariables();++i) {
411  variableInSpanningTree[gm_[f].variableIndex(i)]=false;
412  }
413  size_t connection = 2;
414  // 1 = find a spanning tree and connect higher order auxilary variable to this
415  // 2 = find a spanning subgraph including at least all eges in the subset and connect higher order auxilary variable to this
416  if(connection==2){
417  // ADD ALL
418  for(size_t i=0; i<gm_[f].numberOfVariables();++i) {
419  const IndexType u = gm_[f].variableIndex(i);
420  for(typename EdgeMapType::const_iterator it=neighbours[u].begin() ; it != neighbours[u].end(); ++it){
421  const IndexType v = (*it).first;
422  if(variableInSpanningTree[v] == false && u<v){
423  lpIndexVector.push_back((*it).second);
424  }
425  }
426  }
427  }
428  else if(connection==1){
429  // ADD TREE
430  for(size_t i=0; i<gm_[f].numberOfVariables();++i) {
431  const IndexType u = gm_[f].variableIndex(i);
432  for(typename EdgeMapType::const_iterator it=neighbours[u].begin() ; it != neighbours[u].end(); ++it){
433  const IndexType v = (*it).first;
434  if(variableInSpanningTree[v] == false){
435  variableInSpanningTree[v] = true;
436  lpIndexVector.push_back((*it).second);
437  }
438  }
439  }
440  }
441  else{
442  OPENGM_ASSERT(false);
443  }
444  higherOrderTerms.back().lpIndices_=lpIndexVector;
445 
446  // Check if edges need to be added to have a spanning subgraph
447  //TODO
448  }
449  }
450  //std::cout << "Additional Edges: "<<numberOfAdditionalInternalEdges<<std::endl;
451  return numberOfInternalEdges;
452 }
453 
454 template<class GM, class ACC>
455 Multicut<GM, ACC>::~Multicut() {
456  env_.end();
457 }
458 
459 template<class GM, class ACC>
460 Multicut<GM, ACC>::Multicut
461 (
462  const GraphicalModelType& gm,
463  Parameter para
464  ) : gm_(gm), parameter_(para) , bound_(-std::numeric_limits<double>::infinity()), infinity_(1e8), integerMode_(false),
465  EPS_(1e-7)
466 {
467  if(typeid(ACC) != typeid(opengm::Minimizer) || typeid(OperatorType) != typeid(opengm::Adder)) {
468  throw RuntimeError("This implementation does only supports Min-Plus-Semiring.");
469  }
470  if(parameter_.reductionMode_<0 ||parameter_.reductionMode_>3) {
471  throw RuntimeError("Reduction Mode has to be 1, 2 or 3!");
472  }
473 
474  //Set Problem Type
475  setProblemType();
476  if(problemType_ == INVALID)
477  throw RuntimeError("Invalid Model for Multicut-Solver! Solver requires a generalized potts model!");
478 
479  //Calculate Neighbourhood
480  std::vector<double> valuesHigherOrder;
481  std::vector<HigherOrderTerm> higherOrderTerms;
482  numberOfInternalEdges_ = getNeighborhood(numberOfTerminalEdges_, neighbours, edgeNodes_ ,higherOrderTerms);
483  numberOfNodes_ = gm_.numberOfVariables();
484 
485  if(parameter_.useBufferedStates_){
486  bufferedValue_ = std::numeric_limits<double>::infinity();
487  bufferedStates_.resize(numberOfNodes_,0);
488  }
489 
490  // Display some info
491  if(parameter_.verbose_ == true) {
492  std::cout << "** Multicut Info" << std::endl;
493  if(problemType_==MC)
494  std::cout << " problemType_: Multicut" << std::endl;
495  if(problemType_==MWC)
496  std::cout << " problemType_: Multiway Cut" << std::endl;
497  std::cout << " numberOfInternalEdges_: " << numberOfInternalEdges_ << std::endl;
498  std::cout << " numberOfNodes_: " << numberOfNodes_ << std::endl;
499  std::cout << " allowCutsWithin_: ";
500  if(problemType_==MWC && parameter_.allowCutsWithin_.size() == numberOfTerminals_){
501  for(size_t i=0; i<parameter_.allowCutsWithin_.size(); ++i)
502  if(parameter_.allowCutsWithin_[i]) std::cout<<i<<" ";
503  }
504  else{
505  std::cout<<"none";
506  }
507  std::cout << std::endl;
508  std::cout << " higherOrderTerms.size(): " << higherOrderTerms.size() << std::endl;
509  std::cout << " numberOfTerminals_: " << numberOfTerminals_ << std::endl;
510  }
511 
512  //Build Objective Value
513  constant_=0;
514  size_t valueSize;
515  if(numberOfTerminals_==0) valueSize = numberOfInternalEdges_;
516  else valueSize = numberOfTerminalEdges_+numberOfInternalEdges_+numberOfInterTerminalEdges_;
517  std::vector<double> values (valueSize,0);
518 
519  for(size_t f=0; f<gm_.numberOfFactors(); ++f) {
520  if(gm_[f].numberOfVariables() == 0) {
521  LabelType l = 0;
522  constant_ += gm_[f](&l);
523  }
524  else if(gm_[f].numberOfVariables() == 1) {
525  IndexType node = gm_[f].variableIndex(0);
526  for(LabelType i=0; i<gm_.numberOfLabels(node); ++i) {
527  for(LabelType j=0; j<gm_.numberOfLabels(node); ++j) {
528  if(i==j) values[node*numberOfTerminals_+i] += (1.0/(numberOfTerminals_-1)-1) * gm_[f](&j);
529  else values[node*numberOfTerminals_+i] += (1.0/(numberOfTerminals_-1)) * gm_[f](&j);
530  }
531  }
532  }
533  else if(gm_[f].numberOfVariables() == 2) {
534  if(gm_[f].numberOfLabels(0)==2 && gm_[f].numberOfLabels(1)==2){
535  IndexType node0 = gm_[f].variableIndex(0);
536  IndexType node1 = gm_[f].variableIndex(1);
537  LabelType cc[] = {0,0}; ValueType a = gm_[f](cc);
538  cc[0]=1;cc[1]=1; ValueType b = gm_[f](cc);
539  cc[0]=0;cc[1]=1; ValueType c = gm_[f](cc);
540  cc[0]=1;cc[1]=0; ValueType d = gm_[f](cc);
541 
542  values[neighbours[gm_[f].variableIndex(0)][gm_[f].variableIndex(1)]] += ((c+d-a-a) - (b-a))/2.0;
543  values[node0*numberOfTerminals_+0] += ((b-a)-(-d+c))/2.0;
544  values[node1*numberOfTerminals_+0] += ((b-a)-( d-c))/2.0;
545  constant_ += a;
546  }else{
547  LabelType cc0[] = {0,0};
548  LabelType cc1[] = {0,1};
549  values[neighbours[gm_[f].variableIndex(0)][gm_[f].variableIndex(1)]] += gm_[f](cc1) - gm_[f](cc0);
550  constant_ += gm_[f](cc0);
551  }
552  }
553  }
554  for(size_t h=0; h<higherOrderTerms.size();++h){
555  if(higherOrderTerms[h].potts_) {
556  const IndexType f = higherOrderTerms[h].factorID_;
557  higherOrderTerms[h].valueIndex_= valuesHigherOrder.size();
558  OPENGM_ASSERT(gm_[f].numberOfVariables() > 2);
559  std::vector<LabelType> cc0(gm_[f].numberOfVariables(),0);
560  std::vector<LabelType> cc1(gm_[f].numberOfVariables(),0);
561  cc1[0] = 1;
562  valuesHigherOrder.push_back(gm_[f](cc1.begin()) - gm_[f](cc0.begin()) );
563  constant_ += gm_[f](cc0.begin());
564  }
565  else{
566  const IndexType f = higherOrderTerms[h].factorID_;
567  higherOrderTerms[h].valueIndex_= valuesHigherOrder.size();
568  if(gm_[f].numberOfVariables() == 3) {
569  size_t i[] = {0, 1, 2 };
570  valuesHigherOrder.push_back(gm_[f](i));
571  i[0]=0; i[1]=0; i[2]=1;
572  valuesHigherOrder.push_back(gm_[f](i));
573  i[0]=0; i[1]=1; i[2]=0;
574  valuesHigherOrder.push_back(gm_[f](i));
575  i[0]=1; i[1]=0; i[2]=0;
576  valuesHigherOrder.push_back(gm_[f](i));
577  i[0]=0; i[1]=0; i[2]=0;
578  valuesHigherOrder.push_back(gm_[f](i));
579  }
580  else if(gm_[f].numberOfVariables() == 4) {
581  size_t i[] = {0, 1, 2, 3 };//0
582  if(numberOfTerminals_>=4){
583  valuesHigherOrder.push_back(gm_[f](i));
584  }else{
585  valuesHigherOrder.push_back(0.0);
586  }
587  if(numberOfTerminals_>=3){
588  i[0]=0; i[1]=0; i[2]=1; i[3] = 2;//1
589  valuesHigherOrder.push_back(gm_[f](i));
590  i[0]=0; i[1]=1; i[2]=0; i[3] = 2;//2
591  valuesHigherOrder.push_back(gm_[f](i));
592  i[0]=0; i[1]=1; i[2]=1; i[3] = 2;//4
593  valuesHigherOrder.push_back(gm_[f](i));
594  }else{
595  valuesHigherOrder.push_back(0.0);
596  valuesHigherOrder.push_back(0.0);
597  valuesHigherOrder.push_back(0.0);
598  }
599  i[0]=0; i[1]=0; i[2]=0; i[3] = 1;//7
600  valuesHigherOrder.push_back(gm_[f](i));
601  i[0]=0; i[1]=1; i[2]=2; i[3] = 0;//8
602  valuesHigherOrder.push_back(gm_[f](i));
603  i[0]=0; i[1]=1; i[2]=1; i[3] = 0;//12
604  valuesHigherOrder.push_back(gm_[f](i));
605  i[0]=1; i[1]=0; i[2]=2; i[3] = 0;//16
606  valuesHigherOrder.push_back(gm_[f](i));
607  i[0]=1; i[1]=0; i[2]=1; i[3] = 0;//18
608  valuesHigherOrder.push_back(gm_[f](i));
609  i[0]=0; i[1]=0; i[2]=1; i[3] = 0;//25
610  valuesHigherOrder.push_back(gm_[f](i));
611  i[0]=1; i[1]=2; i[2]=0; i[3] = 0;//32
612  valuesHigherOrder.push_back(gm_[f](i));
613  i[0]=1; i[1]=1; i[2]=0; i[3] = 0;//33
614  valuesHigherOrder.push_back(gm_[f](i));
615  i[0]=0; i[1]=1; i[2]=0; i[3] = 0;//42
616  valuesHigherOrder.push_back(gm_[f](i));
617  i[0]=1; i[1]=0; i[2]=0; i[3] = 0;//52
618  valuesHigherOrder.push_back(gm_[f](i));
619  i[0]=0; i[1]=0; i[2]=0; i[3] = 0;//63
620  valuesHigherOrder.push_back(gm_[f](i));
621  }
622  else{
623  const IndexType f = higherOrderTerms[h].factorID_;
624  higherOrderTerms[h].valueIndex_= valuesHigherOrder.size();
625  P_.resize(gm_[f].numberOfVariables());
626  std::vector<LabelType> l(gm_[f].numberOfVariables());
627  for(size_t i=0; i<P_.BellNumber(gm_[f].numberOfVariables()); ++i){
628  P_.getPartition(i,l);
629  valuesHigherOrder.push_back(gm_[f](l.begin()));
630  }
631  //throw RuntimeError("Generalized Potts Terms of an order larger than 4 a currently not supported. If U really need them let us know!");
632  }
633  }
634  }
635 
636  //count auxilary variables
637  numberOfHigherOrderValues_ = valuesHigherOrder.size();
638 
639  // build LP
640  //std::cout << "Higher order auxilary variables " << numberOfHigherOrderValues_ << std::endl;
641  //std::cout << "TerminalEdges " << numberOfTerminalEdges_ << std::endl;
642  OPENGM_ASSERT( numberOfTerminalEdges_ == gm_.numberOfVariables()*numberOfTerminals_ );
643  //std::cout << "InternalEdges " << numberOfInternalEdges_ << std::endl;
644 
645  OPENGM_ASSERT(values.size() == numberOfTerminalEdges_+numberOfInternalEdges_+numberOfInterTerminalEdges_);
646  IloInt N = values.size() + numberOfHigherOrderValues_;
647  model_ = IloModel(env_);
648  x_ = IloNumVarArray(env_);
649  c_ = IloRangeArray(env_);
650  obj_ = IloMinimize(env_);
651  sol_ = IloNumArray(env_,N);
652 
653  // set variables and objective
654  x_.add(IloNumVarArray(env_, N, 0, 1, ILOFLOAT));
655 
656  IloNumArray obj(env_,N);
657  for (size_t i=0; i< values.size();++i) {
658  if(values[i]==0)
659  obj[i] = 0.0;//1e-50; //for numerical reasons
660  else
661  obj[i] = values[i];
662  }
663  {
664  size_t count =0;
665  for (size_t i=0; i<valuesHigherOrder.size();++i) {
666  obj[values.size()+count++] = valuesHigherOrder[i];
667  }
668  OPENGM_ASSERT(count == numberOfHigherOrderValues_);
669  }
670  obj_.setLinearCoefs(x_,obj);
671 
672  // set constraints
673  size_t constraintCounter = 0;
674  // multiway cut constraints
675  if(problemType_ == MWC) {
676  // From each internal-node only one terminal-edge should be 0
677  for(IndexType var=0; var<gm_.numberOfVariables(); ++var) {
678  c_.add(IloRange(env_, numberOfTerminals_-1, numberOfTerminals_-1));
679  for(LabelType i=0; i<gm_.numberOfLabels(var); ++i) {
680  c_[constraintCounter].setLinearCoef(x_[var*numberOfTerminals_+i],1);
681  }
682  ++constraintCounter;
683  }
684  // Inter-terminal-edges have to be 1
685  for(size_t i=0; i<(size_t)(numberOfTerminals_*(numberOfTerminals_-1)/2); ++i) {
686  c_.add(IloRange(env_, 1, 1));
687  c_[constraintCounter].setLinearCoef(x_[numberOfTerminalEdges_+numberOfInternalEdges_+i],1);
688  ++constraintCounter;
689  }
690  }
691 
692 
693  // higher order constraints
694  size_t count = 0;
695  for(size_t i=0; i<higherOrderTerms.size(); ++i) {
696  size_t factorID = higherOrderTerms[i].factorID_;
697  size_t numVar = gm_[factorID].numberOfVariables();
698  OPENGM_ASSERT(numVar>2);
699 
700  if(higherOrderTerms[i].potts_) {
701  double b = higherOrderTerms[i].lpIndices_.size();
702 
703 
704  // Add only one constraint is sufficient with {0,1} constraints
705  // ------------------------------------------------------------
706  // ** -|E|+1 <= -|E|*y_H+\sum_{e\in H} y_e <= 0
707  if(parameter_.reductionMode_ % 2 == 1){
708  c_.add(IloRange(env_, -b+1 , 0));
709  for(size_t i1=0; i1<higherOrderTerms[i].lpIndices_.size();++i1) {
710  const LPIndexType edgeID = higherOrderTerms[i].lpIndices_[i1];
711  c_[constraintCounter].setLinearCoef(x_[edgeID],1);
712  }
713  c_[constraintCounter].setLinearCoef(x_[values.size()+count],-b);
714  constraintCounter += 1;
715  }
716  // In general this additional contraints and more local constraints leeds to tighter relaxations
717  // ---------------------------------------------------------------------------------------------
718  if(parameter_.reductionMode_ % 4 >=2){
719  // ** y_H <= sum_{e \in H} y_e
720  c_.add(IloRange(env_, -2.0*b, 0));
721  for(size_t i1=0; i1<higherOrderTerms[i].lpIndices_.size();++i1) {
722  const LPIndexType edgeID = higherOrderTerms[i].lpIndices_[i1];
723  c_[constraintCounter].setLinearCoef(x_[edgeID],-1);
724  }
725  c_[constraintCounter].setLinearCoef(x_[values.size()+count],1);
726  constraintCounter += 1;
727 
728  // ** forall e \in H : y_H>=y_e
729  for(size_t i1=0; i1<higherOrderTerms[i].lpIndices_.size();++i1) {
730  const LPIndexType edgeID = higherOrderTerms[i].lpIndices_[i1];
731  c_.add(IloRange(env_, 0 , 1));
732  c_[constraintCounter].setLinearCoef(x_[edgeID],-1);
733  c_[constraintCounter].setLinearCoef(x_[values.size()+count],1);
734  constraintCounter += 1;
735  }
736  }
737  count++;
738  }else{
739  if(numVar==3) {
740  OPENGM_ASSERT(higherOrderTerms[i].valueIndex_<=valuesHigherOrder.size());
741  LPIndexType edgeIDs[3];
742  edgeIDs[0] = neighbours[gm_[factorID].variableIndex(0)][gm_[factorID].variableIndex(1)];
743  edgeIDs[1] = neighbours[gm_[factorID].variableIndex(0)][gm_[factorID].variableIndex(2)];
744  edgeIDs[2] = neighbours[gm_[factorID].variableIndex(1)][gm_[factorID].variableIndex(2)];
745 
746  const unsigned int P[] = {0,1,2,4,7,8,12,16,18,25,32,33,42,52,63};
747  double c[3];
748 
749  c_.add(IloRange(env_, 1, 1));
750  size_t lvc=0;
751  for(size_t p=0; p<5; p++){
752  if(true || valuesHigherOrder[higherOrderTerms[i].valueIndex_+p]!=0){
753  c_[constraintCounter].setLinearCoef(x_[values.size()+count+lvc],1);
754  ++lvc;
755  }
756  }
757  ++constraintCounter;
758 
759  for(size_t p=0; p<5; p++){
760  if(true || valuesHigherOrder[higherOrderTerms[i].valueIndex_+p]!=0){
761  double ub = 2.0;
762  double lb = 0.0;
763  unsigned int mask = 1;
764  for(size_t n=0; n<3; n++){
765  if(P[p] & mask){
766  c[n] = -1.0;
767  ub--;
768  lb--;
769  }
770  else{
771  c[n] = 1.0;
772  }
773  mask = mask << 1;
774  }
775  c_.add(IloRange(env_, lb, ub));
776  for(size_t n=0; n<3; n++){
777  c_[constraintCounter].setLinearCoef(x_[edgeIDs[n]],c[n]);
778  }
779  c_[constraintCounter].setLinearCoef(x_[values.size()+count],-1);
780  ++constraintCounter;
781 
782  for(size_t n=0; n<3; n++){
783  if(c[n]>0){
784  c_.add(IloRange(env_, 0, 1));
785  c_[constraintCounter].setLinearCoef(x_[edgeIDs[n]],1);
786  c_[constraintCounter].setLinearCoef(x_[values.size()+count],-1);
787  ++constraintCounter;
788  }else{
789  c_.add(IloRange(env_, -1, 0));
790  c_[constraintCounter].setLinearCoef(x_[edgeIDs[n]],-1);
791  c_[constraintCounter].setLinearCoef(x_[values.size()+count],-1);
792  ++constraintCounter;
793  }
794  }
795  ++count;
796  }
797  }
798  }
799  else if(numVar==4) {
800  OPENGM_ASSERT(higherOrderTerms[i].valueIndex_<=valuesHigherOrder.size());
801  LPIndexType edgeIDs[6];
802  edgeIDs[0] = neighbours[gm_[factorID].variableIndex(0)][gm_[factorID].variableIndex(1)];
803  edgeIDs[1] = neighbours[gm_[factorID].variableIndex(0)][gm_[factorID].variableIndex(2)];
804  edgeIDs[2] = neighbours[gm_[factorID].variableIndex(1)][gm_[factorID].variableIndex(2)];
805  edgeIDs[3] = neighbours[gm_[factorID].variableIndex(0)][gm_[factorID].variableIndex(3)];
806  edgeIDs[4] = neighbours[gm_[factorID].variableIndex(1)][gm_[factorID].variableIndex(3)];
807  edgeIDs[5] = neighbours[gm_[factorID].variableIndex(2)][gm_[factorID].variableIndex(3)];
808 
809 
810  const unsigned int P[] = {0,1,2,4,7,8,12,16,18,25,32,33,42,52,63};
811  double c[6];
812 
813  c_.add(IloRange(env_, 1, 1));
814  size_t lvc=0;
815  for(size_t p=0; p<15; p++){
816  if(true ||valuesHigherOrder[higherOrderTerms[i].valueIndex_+p]!=0){
817  c_[constraintCounter].setLinearCoef(x_[values.size()+count+lvc],1);
818  ++lvc;
819  }
820  }
821  ++constraintCounter;
822 
823 
824  for(size_t p=0; p<15; p++){
825  double ub = 5.0;
826  double lb = 0.0;
827  unsigned int mask = 1;
828  for(size_t n=0; n<6; n++){
829  if(P[p] & mask){
830  c[n] = -1.0;
831  ub--;
832  lb--;
833  }
834  else{
835  c[n] = 1.0;
836  }
837  mask = mask << 1;
838  }
839  c_.add(IloRange(env_, lb, ub));
840  for(size_t n=0; n<6; n++){
841  c_[constraintCounter].setLinearCoef(x_[edgeIDs[n]],c[n]);
842  }
843  c_[constraintCounter].setLinearCoef(x_[values.size()+count],-1);
844  ++constraintCounter;
845 
846  for(size_t n=0; n<6; n++){
847  if(c[n]>0){
848  c_.add(IloRange(env_, 0, 1));
849  c_[constraintCounter].setLinearCoef(x_[edgeIDs[n]],1);
850  c_[constraintCounter].setLinearCoef(x_[values.size()+count],-1);
851  ++constraintCounter;
852  }else{
853  c_.add(IloRange(env_, -1, 0));
854  c_[constraintCounter].setLinearCoef(x_[edgeIDs[n]],-1);
855  c_[constraintCounter].setLinearCoef(x_[values.size()+count],-1);
856  ++constraintCounter;
857  }
858  }
859  ++count;
860  }
861  }
862  else{
863  std::vector<LPIndexType> edgeIDs(P_.BellNumber(numVar));
864  {
865  size_t cc=0;
866  for(size_t v1=1; v1<numVar; ++v1){
867  for(size_t v2=0; v2<v1; ++v2){
868  edgeIDs[cc] = neighbours[gm_[factorID].variableIndex(v2)][gm_[factorID].variableIndex(v1)];
869  ++cc;
870  }
871  }
872  }
873  c_.add(IloRange(env_, 1, 1));
874  size_t lvc=0;
875  for(size_t p=0; p<P_.BellNumber(numVar); p++){
876  if(true || valuesHigherOrder[higherOrderTerms[i].valueIndex_+p]!=0){
877  c_[constraintCounter].setLinearCoef(x_[values.size()+count+lvc],1);
878  ++lvc;
879  }
880  }
881  ++constraintCounter;
882 
883  std::vector<double> c(numVar*(numVar-1)/2,0);
884  for(size_t p=0; p<P_.BellNumber(numVar); p++){
885  double ub = numVar*(numVar-1)/2 -1;
886  double lb = 0.0;
887  unsigned int mask = 1;
888  size_t el = P_.getPartition(p);
889  for(size_t n=0; n<numVar*(numVar-1)/2; n++){
890  if(el & mask){
891  c[n] = -1.0;
892  ub--;
893  lb--;
894  }
895  else{
896  c[n] = 1.0;
897  }
898  mask = mask << 1;
899  }
900  c_.add(IloRange(env_, lb, ub));
901  for(size_t n=0; n<numVar*(numVar-1)/2; n++){
902  c_[constraintCounter].setLinearCoef(x_[edgeIDs[n]],c[n]);
903  }
904  c_[constraintCounter].setLinearCoef(x_[values.size()+count],-1);
905  ++constraintCounter;
906 
907  for(size_t n=0; n<numVar*(numVar-1)/2; n++){
908  if(c[n]>0){
909  c_.add(IloRange(env_, 0, 1));
910  c_[constraintCounter].setLinearCoef(x_[edgeIDs[n]],1);
911  c_[constraintCounter].setLinearCoef(x_[values.size()+count],-1);
912  ++constraintCounter;
913  }else{
914  c_.add(IloRange(env_, -1, 0));
915  c_[constraintCounter].setLinearCoef(x_[edgeIDs[n]],-1);
916  c_[constraintCounter].setLinearCoef(x_[values.size()+count],-1);
917  ++constraintCounter;
918  }
919  }
920  ++count;
921  }
922  }
923  }
924  }
925 
926 
927  model_.add(obj_);
928  if(constraintCounter>0) {
929  model_.add(c_);
930  }
931  if(para.initializeWith3Cycles_){
932  add3CycleConstraints();
933  }
934 
935  // initialize solver
936  cplex_ = IloCplex(model_);
937 
938 }
939 
940 template<class GM, class ACC>
941 typename Multicut<GM, ACC>::ProblemType Multicut<GM, ACC>::setProblemType() {
942  problemType_ = MC;
943  for(size_t f=0; f<gm_.numberOfFactors();++f) {
944  if(gm_[f].numberOfVariables()==1) {
945  problemType_ = MWC;
946  }
947  if(gm_[f].numberOfVariables()>1) {
948  for(size_t i=0; i<gm_[f].numberOfVariables();++i) {
949  if(gm_[f].numberOfLabels(i)<gm_.numberOfVariables()) {
950  problemType_ = MWC;
951  }
952  }
953  }
954  if(gm_[f].numberOfVariables()==2 && gm_[f].numberOfLabels(0)==2 && gm_[f].numberOfLabels(1)==2){
955  LabelType l00[] = {0,0};
956  LabelType l01[] = {0,1};
957  LabelType l10[] = {1,0};
958  LabelType l11[] = {1,1};
959  if(gm_[f](l00)!=gm_[f](l11) || gm_[f](l01)!=gm_[f](l10))
960  problemType_ = MWC; //OK - can be reparmetrized
961  }
962  else if(gm_[f].numberOfVariables()>1 && !gm_[f].isGeneralizedPotts()) {
963  problemType_ = INVALID;
964  break;
965  }
966  }
967 
968  // set member variables
969  if(problemType_ == MWC) {
970  numberOfTerminals_ = gm_.numberOfLabels(0);
971  numberOfInterTerminalEdges_ = (numberOfTerminals_*(numberOfTerminals_-1))/2;
972  numberOfTerminalEdges_ = 0;
973  for(IndexType i=0; i<gm_.numberOfVariables(); ++i) {
974  for(LabelType j=0; j<gm_.numberOfLabels(i); ++j) {
975  ++numberOfTerminalEdges_;
976  }
977  }
978  }
979  else{
980  numberOfTerminalEdges_ = 0;
981  numberOfTerminals_ = 0;
982  numberOfInterTerminalEdges_ = 0;
983  }
984 
985  return problemType_;
986 }
987 
988 //**********************************************
989 //**
990 //** Functions that find violated Constraints
991 //**
992 //**********************************************
993 
994 template<class GM, class ACC>
995 size_t Multicut<GM, ACC>::removeUnusedConstraints()
996 {
997  std::cout << "Not Implemented " <<std::endl ;
998  return 0;
999 }
1000 
1001 
1002 
1003 template<class GM, class ACC>
1004 size_t Multicut<GM, ACC>::enforceIntegerConstraints()
1005 {
1006  bool enforceIntegerConstraintsOnTerminalEdges = true;
1007  bool enforceIntegerConstraintsOnInternalEdges = false;
1008 
1009  if(numberOfTerminalEdges_ == 0 || parameter_.allowCutsWithin_.size() == numberOfTerminals_) {
1010  enforceIntegerConstraintsOnInternalEdges = true;
1011  }
1012 
1013  size_t N = 0;
1014  if (enforceIntegerConstraintsOnTerminalEdges)
1015  N += numberOfTerminalEdges_;
1016  if (enforceIntegerConstraintsOnInternalEdges)
1017  N += numberOfInternalEdges_;
1018 
1019  for(size_t i=0; i<N; ++i)
1020  model_.add(IloConversion(env_, x_[i], ILOBOOL));
1021 
1022  for(size_t i=0; i<numberOfHigherOrderValues_; ++i)
1023  model_.add(IloConversion(env_, x_[numberOfTerminalEdges_+numberOfInternalEdges_+numberOfInterTerminalEdges_+i], ILOBOOL));
1024 
1025  integerMode_ = true;
1026 
1027  return N+numberOfHigherOrderValues_;
1028 }
1029 
1036 template<class GM, class ACC>
1037 size_t Multicut<GM, ACC>::findTerminalTriangleConstraints(IloRangeArray& constraint)
1038 {
1039  OPENGM_ASSERT(problemType_ == MWC);
1040  if(!(problemType_ == MWC)) return 0;
1041  size_t tempConstrainCounter = constraintCounter_;
1042 
1043  size_t u,v;
1044  if(parameter_.allowCutsWithin_.size()!=numberOfTerminals_){
1045  for(size_t i=0; i<numberOfInternalEdges_;++i) {
1046  u = edgeNodes_[i].first;//[0];
1047  v = edgeNodes_[i].second;//[1];
1048  for(size_t l=0; l<numberOfTerminals_;++l) {
1049  if(-sol_[numberOfTerminalEdges_+i]+sol_[u*numberOfTerminals_+l]+sol_[v*numberOfTerminals_+l]<-EPS_) {
1050  constraint.add(IloRange(env_, 0 , 2));
1051  constraint[constraintCounter_].setLinearCoef(x_[numberOfTerminalEdges_+i],-1);
1052  constraint[constraintCounter_].setLinearCoef(x_[u*numberOfTerminals_+l],1);
1053  constraint[constraintCounter_].setLinearCoef(x_[v*numberOfTerminals_+l],1);
1054  ++constraintCounter_;
1055  }
1056  if(sol_[numberOfTerminalEdges_+i]-sol_[u*numberOfTerminals_+l]+sol_[v*numberOfTerminals_+l]<-EPS_) {
1057  constraint.add(IloRange(env_, 0 , 2));
1058  constraint[constraintCounter_].setLinearCoef(x_[numberOfTerminalEdges_+i],1);
1059  constraint[constraintCounter_].setLinearCoef(x_[u*numberOfTerminals_+l],-1);
1060  constraint[constraintCounter_].setLinearCoef(x_[v*numberOfTerminals_+l],1);
1061  ++constraintCounter_;
1062  }
1063  if(sol_[numberOfTerminalEdges_+i]+sol_[u*numberOfTerminals_+l]-sol_[v*numberOfTerminals_+l]<-EPS_) {
1064  constraint.add(IloRange(env_, 0 , 2));
1065  constraint[constraintCounter_].setLinearCoef(x_[numberOfTerminalEdges_+i],1);
1066  constraint[constraintCounter_].setLinearCoef(x_[u*numberOfTerminals_+l],1);
1067  constraint[constraintCounter_].setLinearCoef(x_[v*numberOfTerminals_+l],-1);
1068  ++constraintCounter_;
1069  }
1070  }
1071  if(constraintCounter_-tempConstrainCounter >= parameter_.maximalNumberOfConstraintsPerRound_)
1072  break;
1073  }
1074  }
1075  else{
1076  for(size_t i=0; i<numberOfInternalEdges_;++i) {
1077  u = edgeNodes_[i].first;//[0];
1078  v = edgeNodes_[i].second;//[1];
1079  for(size_t l=0; l<numberOfTerminals_;++l) {
1080  if(parameter_.allowCutsWithin_[l])
1081  continue;
1082  if(-sol_[numberOfTerminalEdges_+i]+sol_[u*numberOfTerminals_+l]+sol_[v*numberOfTerminals_+l]<-EPS_) {
1083  constraint.add(IloRange(env_, 0 , 2));
1084  constraint[constraintCounter_].setLinearCoef(x_[numberOfTerminalEdges_+i],-1);
1085  constraint[constraintCounter_].setLinearCoef(x_[u*numberOfTerminals_+l],1);
1086  constraint[constraintCounter_].setLinearCoef(x_[v*numberOfTerminals_+l],1);
1087  ++constraintCounter_;
1088  }
1089  if(sol_[numberOfTerminalEdges_+i]-sol_[u*numberOfTerminals_+l]+sol_[v*numberOfTerminals_+l]<-EPS_) {
1090  constraint.add(IloRange(env_, 0 , 2));
1091  constraint[constraintCounter_].setLinearCoef(x_[numberOfTerminalEdges_+i],1);
1092  constraint[constraintCounter_].setLinearCoef(x_[u*numberOfTerminals_+l],-1);
1093  constraint[constraintCounter_].setLinearCoef(x_[v*numberOfTerminals_+l],1);
1094  ++constraintCounter_;
1095  }
1096  if(sol_[numberOfTerminalEdges_+i]+sol_[u*numberOfTerminals_+l]-sol_[v*numberOfTerminals_+l]<-EPS_) {
1097  constraint.add(IloRange(env_, 0 , 2));
1098  constraint[constraintCounter_].setLinearCoef(x_[numberOfTerminalEdges_+i],1);
1099  constraint[constraintCounter_].setLinearCoef(x_[u*numberOfTerminals_+l],1);
1100  constraint[constraintCounter_].setLinearCoef(x_[v*numberOfTerminals_+l],-1);
1101  ++constraintCounter_;
1102  }
1103  }
1104  if(constraintCounter_-tempConstrainCounter >= parameter_.maximalNumberOfConstraintsPerRound_)
1105  break;
1106  }
1107  }
1108  return constraintCounter_-tempConstrainCounter;
1109 }
1110 
1117 template<class GM, class ACC>
1118 size_t Multicut<GM, ACC>::findMultiTerminalConstraints(IloRangeArray& constraint)
1119 {
1120  OPENGM_ASSERT(problemType_ == MWC);
1121  if(!(problemType_ == MWC)) return 0;
1122  size_t tempConstrainCounter = constraintCounter_;
1123  if(parameter_.allowCutsWithin_.size()==numberOfTerminals_){
1124  for(size_t i=0; i<parameter_.allowCutsWithin_.size();++i)
1125  if(parameter_.allowCutsWithin_[i])
1126  return 0; //Can not gurantee that Multi Terminal Constraints are valid cuts
1127  }
1128 
1129  size_t u,v;
1130  for(size_t i=0; i<numberOfInternalEdges_;++i) {
1131  u = edgeNodes_[i].first;//[0];
1132  v = edgeNodes_[i].second;//[1];
1133  std::vector<size_t> terminals1;
1134  std::vector<size_t> terminals2;
1135  double sum1 = 0;
1136  double sum2 = 0;
1137  for(size_t l=0; l<numberOfTerminals_;++l) {
1138  if(sol_[u*numberOfTerminals_+l]-sol_[v*numberOfTerminals_+l] > EPS_) {
1139  terminals1.push_back(l);
1140  sum1 += sol_[u*numberOfTerminals_+l]-sol_[v*numberOfTerminals_+l];
1141  }
1142  if(sol_[v*numberOfTerminals_+l]-sol_[u*numberOfTerminals_+l] > EPS_) {
1143  terminals2.push_back(l);
1144  sum2 +=sol_[v*numberOfTerminals_+l]-sol_[u*numberOfTerminals_+l];
1145  }
1146  }
1147  if(sol_[numberOfTerminalEdges_+i]-sum1<-EPS_) {
1148  constraint.add(IloRange(env_, 0 , 200000));
1149  constraint[constraintCounter_].setLinearCoef(x_[numberOfTerminalEdges_+i],1);
1150  for(size_t k=0; k<terminals1.size(); ++k) {
1151  constraint[constraintCounter_].setLinearCoef(x_[u*numberOfTerminals_+terminals1[k]],-1);
1152  constraint[constraintCounter_].setLinearCoef(x_[v*numberOfTerminals_+terminals1[k]],1);
1153  }
1154  ++constraintCounter_;
1155  }
1156  if(sol_[numberOfTerminalEdges_+i]-sum2<-EPS_) {
1157  constraint.add(IloRange(env_, 0 , 200000));
1158  constraint[constraintCounter_].setLinearCoef(x_[numberOfTerminalEdges_+i],1);
1159  for(size_t k=0; k<terminals2.size(); ++k) {
1160  constraint[constraintCounter_].setLinearCoef(x_[u*numberOfTerminals_+terminals2[k]],1);
1161  constraint[constraintCounter_].setLinearCoef(x_[v*numberOfTerminals_+terminals2[k]],-1);
1162  }
1163  ++constraintCounter_;
1164  }
1165  if(constraintCounter_-tempConstrainCounter >= parameter_.maximalNumberOfConstraintsPerRound_)
1166  break;
1167  }
1168  return constraintCounter_-tempConstrainCounter;
1169 }
1170 
1177 template<class GM, class ACC>
1178 size_t Multicut<GM, ACC>::findIntegerTerminalTriangleConstraints(IloRangeArray& constraint, std::vector<LabelType>& conf)
1179 {
1180  OPENGM_ASSERT(integerMode_);
1181  OPENGM_ASSERT(problemType_ == MWC);
1182  if(!(problemType_ == MWC)) return 0;
1183  size_t tempConstrainCounter = constraintCounter_;
1184 
1185  size_t u,v;
1186  if(parameter_.allowCutsWithin_.size()!=numberOfTerminals_){
1187  for(size_t i=0; i<numberOfInternalEdges_;++i) {
1188  u = edgeNodes_[i].first;//[0];
1189  v = edgeNodes_[i].second;//[1];
1190  if(sol_[numberOfTerminalEdges_+i]<EPS_ && (conf[u]!=conf[v]) ) {
1191  if(sol_[numberOfTerminalEdges_+i]-sol_[u*numberOfTerminals_+conf[u]]+sol_[v*numberOfTerminals_+conf[u]]<=0) {
1192  constraint.add(IloRange(env_, 0 , 10));
1193  constraint[constraintCounter_].setLinearCoef(x_[numberOfTerminalEdges_+i],1);
1194  constraint[constraintCounter_].setLinearCoef(x_[u*numberOfTerminals_+conf[u]],-1);
1195  constraint[constraintCounter_].setLinearCoef(x_[v*numberOfTerminals_+conf[u]],1);
1196  ++constraintCounter_;
1197  }
1198  if(sol_[numberOfTerminalEdges_+i]-sol_[u*numberOfTerminals_+conf[u]]+sol_[v*numberOfTerminals_+conf[u]]<=0) {
1199  constraint.add(IloRange(env_, 0 , 10));
1200  constraint[constraintCounter_].setLinearCoef(x_[numberOfTerminalEdges_+i],1);
1201  constraint[constraintCounter_].setLinearCoef(x_[u*numberOfTerminals_+conf[u]],1);
1202  constraint[constraintCounter_].setLinearCoef(x_[v*numberOfTerminals_+conf[u]],-1);
1203  ++constraintCounter_;
1204  }
1205  if(sol_[numberOfTerminalEdges_+i]-sol_[u*numberOfTerminals_+conf[v]]+sol_[v*numberOfTerminals_+conf[v]]<=0) {
1206  constraint.add(IloRange(env_, 0 , 10));
1207  constraint[constraintCounter_].setLinearCoef(x_[numberOfTerminalEdges_+i],1);
1208  constraint[constraintCounter_].setLinearCoef(x_[u*numberOfTerminals_+conf[v]],-1);
1209  constraint[constraintCounter_].setLinearCoef(x_[v*numberOfTerminals_+conf[v]],1);
1210  ++constraintCounter_;
1211  }
1212  if(sol_[numberOfTerminalEdges_+i]+sol_[u*numberOfTerminals_+conf[v]]-sol_[v*numberOfTerminals_+conf[v]]<=0) {
1213  constraint.add(IloRange(env_, 0 , 10));
1214  constraint[constraintCounter_].setLinearCoef(x_[numberOfTerminalEdges_+i],1);
1215  constraint[constraintCounter_].setLinearCoef(x_[u*numberOfTerminals_+conf[v]],1);
1216  constraint[constraintCounter_].setLinearCoef(x_[v*numberOfTerminals_+conf[v]],-1);
1217  ++constraintCounter_;
1218  }
1219  }
1220  if(sol_[numberOfTerminalEdges_+i]>1-EPS_ && (conf[u]==conf[v]) ) {
1221  constraint.add(IloRange(env_, 0 , 10));
1222  constraint[constraintCounter_].setLinearCoef(x_[numberOfTerminalEdges_+i],-1);
1223  constraint[constraintCounter_].setLinearCoef(x_[u*numberOfTerminals_+conf[u]],1);
1224  constraint[constraintCounter_].setLinearCoef(x_[v*numberOfTerminals_+conf[v]],1);
1225  ++constraintCounter_;
1226  }
1227  if(constraintCounter_-tempConstrainCounter >= parameter_.maximalNumberOfConstraintsPerRound_)
1228  break;
1229  }
1230  }
1231  else{ // Allow Cuts within classes
1232  for(size_t i=0; i<numberOfInternalEdges_;++i) {
1233  u = edgeNodes_[i].first;//[0];
1234  v = edgeNodes_[i].second;//[1];
1235  if(sol_[numberOfTerminalEdges_+i]<EPS_ && (conf[u]!=conf[v]) ) {
1236  if(sol_[numberOfTerminalEdges_+i]-sol_[u*numberOfTerminals_+conf[u]]+sol_[v*numberOfTerminals_+conf[u]]<=0) {
1237  constraint.add(IloRange(env_, 0 , 10));
1238  constraint[constraintCounter_].setLinearCoef(x_[numberOfTerminalEdges_+i],1);
1239  constraint[constraintCounter_].setLinearCoef(x_[u*numberOfTerminals_+conf[u]],-1);
1240  constraint[constraintCounter_].setLinearCoef(x_[v*numberOfTerminals_+conf[u]],1);
1241  ++constraintCounter_;
1242  }
1243  if(sol_[numberOfTerminalEdges_+i]-sol_[u*numberOfTerminals_+conf[u]]+sol_[v*numberOfTerminals_+conf[u]]<=0) {
1244  constraint.add(IloRange(env_, 0 , 10));
1245  constraint[constraintCounter_].setLinearCoef(x_[numberOfTerminalEdges_+i],1);
1246  constraint[constraintCounter_].setLinearCoef(x_[u*numberOfTerminals_+conf[u]],1);
1247  constraint[constraintCounter_].setLinearCoef(x_[v*numberOfTerminals_+conf[u]],-1);
1248  ++constraintCounter_;
1249  }
1250  if(sol_[numberOfTerminalEdges_+i]-sol_[u*numberOfTerminals_+conf[v]]+sol_[v*numberOfTerminals_+conf[v]]<=0) {
1251  constraint.add(IloRange(env_, 0 , 10));
1252  constraint[constraintCounter_].setLinearCoef(x_[numberOfTerminalEdges_+i],1);
1253  constraint[constraintCounter_].setLinearCoef(x_[u*numberOfTerminals_+conf[v]],-1);
1254  constraint[constraintCounter_].setLinearCoef(x_[v*numberOfTerminals_+conf[v]],1);
1255  ++constraintCounter_;
1256  }
1257  if(sol_[numberOfTerminalEdges_+i]+sol_[u*numberOfTerminals_+conf[v]]-sol_[v*numberOfTerminals_+conf[v]]<=0) {
1258  constraint.add(IloRange(env_, 0 , 10));
1259  constraint[constraintCounter_].setLinearCoef(x_[numberOfTerminalEdges_+i],1);
1260  constraint[constraintCounter_].setLinearCoef(x_[u*numberOfTerminals_+conf[v]],1);
1261  constraint[constraintCounter_].setLinearCoef(x_[v*numberOfTerminals_+conf[v]],-1);
1262  ++constraintCounter_;
1263  }
1264  }
1265  if(sol_[numberOfTerminalEdges_+i]>1-EPS_ && (conf[u]==conf[v]) && !parameter_.allowCutsWithin_[conf[u]] ) {
1266  constraint.add(IloRange(env_, 0 , 10));
1267  constraint[constraintCounter_].setLinearCoef(x_[numberOfTerminalEdges_+i],-1);
1268  constraint[constraintCounter_].setLinearCoef(x_[u*numberOfTerminals_+conf[u]],1);
1269  constraint[constraintCounter_].setLinearCoef(x_[v*numberOfTerminals_+conf[v]],1);
1270  ++constraintCounter_;
1271  }
1272  if(constraintCounter_-tempConstrainCounter >= parameter_.maximalNumberOfConstraintsPerRound_)
1273  break;
1274  }
1275  }
1276 
1277  return constraintCounter_-tempConstrainCounter;
1278 }
1283 template<class GM, class ACC>
1284 size_t Multicut<GM, ACC>::add3CycleConstraints()
1285 {
1286  //TODO: The search can be made faster, on should consider this later.
1287  IloRangeArray constraint = IloRangeArray(env_);
1288  size_t constraintCounter =0;
1289  LPIndexType edge1,edge2,edge3;
1290  typename EdgeMapType::const_iterator it2;
1291  typename EdgeMapType::const_iterator it3;
1292  typename EdgeMapType::const_iterator it4;
1293  for(IndexType node1=0; node1<numberOfNodes_; ++node1){
1294  for(it2=neighbours[node1].begin() ; it2 != neighbours[node1].end(); ++it2) {
1295  const IndexType node2=(*it2).first;
1296  edge1 = (*it2).second;
1297  if(node2<=node1) continue;
1298  for(it3=neighbours[node1].begin() ; it3 != neighbours[node1].end(); ++it3) {
1299  const IndexType node3=(*it3).first;
1300  edge2 = (*it3).second;
1301  if(node3<=node1) continue;
1302  if(node3<=node2) continue;
1303  it4 = neighbours[node2].find(node3);
1304  if(it4 != neighbours[node2].end()) {
1305  edge3 = (*it4).second;
1306  //found 3cycle -> add it.
1307  constraint.add(IloRange(env_, 0 , 1000000000));
1308  constraint[constraintCounter].setLinearCoef(x_[edge1],-1);
1309  constraint[constraintCounter].setLinearCoef(x_[edge2],1);
1310  constraint[constraintCounter].setLinearCoef(x_[edge3],1);
1311  ++constraintCounter;
1312  //
1313  constraint.add(IloRange(env_, 0 , 1000000000));
1314  constraint[constraintCounter].setLinearCoef(x_[edge1],1);
1315  constraint[constraintCounter].setLinearCoef(x_[edge2],-1);
1316  constraint[constraintCounter].setLinearCoef(x_[edge3],1);
1317  ++constraintCounter;
1318  //
1319  constraint.add(IloRange(env_, 0 , 1000000000));
1320  constraint[constraintCounter].setLinearCoef(x_[edge1],1);
1321  constraint[constraintCounter].setLinearCoef(x_[edge2],1);
1322  constraint[constraintCounter].setLinearCoef(x_[edge3],-1);
1323  ++constraintCounter;
1324  }
1325  }
1326  }
1327  }
1328  if(constraintCounter>0){
1329  std::cout << "Add "<<constraintCounter<<" constraints for the initial relaxation"<<std::endl;
1330  model_.add(constraint);
1331  }
1332  return constraintCounter;
1333 }
1334 
1339 template<class GM, class ACC>
1340 size_t Multicut<GM, ACC>::findCycleConstraints(
1341  IloRangeArray& constraint,
1342  bool addOnlyFacetDefiningConstraints,
1343  bool usePreBounding
1344 )
1345 {
1346  std::vector<LabelType> partit;
1347  std::vector<std::list<size_t> > neighbours0;
1348 
1349  size_t tempConstrainCounter = constraintCounter_;
1350 
1351  if(usePreBounding){
1352  partition(partit,neighbours0,1-EPS_);
1353  }
1354 
1355  std::map<std::pair<IndexType,IndexType>,size_t> counter;
1356 
1357  if(!parameter_.useOldPriorityQueue_){
1358  std::vector<IndexType> prev(neighbours.size());
1359  opengm::ChangeablePriorityQueue<double> openNodes(neighbours.size());
1360  std::vector<IndexType> path;
1361  path.reserve(neighbours.size());
1362  for(size_t i=0; i<numberOfInternalEdges_;++i) {
1363 
1364  IndexType u = edgeNodes_[i].first;//[0];
1365  IndexType v = edgeNodes_[i].second;//[1];
1366 
1367  if(usePreBounding && partit[u] != partit[v])
1368  continue;
1369 
1370  OPENGM_ASSERT(i+numberOfTerminalEdges_ == neighbours[u][v]);
1371  OPENGM_ASSERT(i+numberOfTerminalEdges_ == neighbours[v][u]);
1372 
1373  if(sol_[numberOfTerminalEdges_+i]>EPS_){
1374  //search for cycle
1375  double pathLength;
1376  pathLength = shortestPath2(u,v,neighbours,sol_,path,prev,openNodes,sol_[numberOfTerminalEdges_+i],addOnlyFacetDefiningConstraints && parameter_.useChordalSearch_);
1377  // pathLength = shortestPath(u,v,neighbours,sol_,path,sol_[numberOfTerminalEdges_+i],addOnlyFacetDefiningConstraints);
1378 
1379  if(sol_[numberOfTerminalEdges_+i]-EPS_>pathLength){
1380  bool postChordlessCheck = addOnlyFacetDefiningConstraints && !parameter_.useChordalSearch_;
1381  bool chordless = true;
1382  if(postChordlessCheck){
1383  for(size_t n1=0;n1<path.size();++n1){
1384  for(size_t n2=n1+2;n2<path.size();++n2){
1385  if(path[n1]==v && path[n2]==u) continue;
1386  if(neighbours[path[n2]].find(path[n1])!=neighbours[path[n2]].end()) {
1387  chordless = false; // do not update node if path is chordal
1388  break;
1389  }
1390  }
1391  if(!chordless)
1392  break;
1393  }
1394  }
1395  if(chordless){
1396  OPENGM_ASSERT(path.size()>2);
1397  constraint.add(IloRange(env_, 0 , 1000000000));
1398  //negative zero seemed to be required for numerical reasons, even CPlex handel this by its own, too.
1399  constraint[constraintCounter_].setLinearCoef(x_[numberOfTerminalEdges_+i],-1);
1400  for(size_t n=0;n<path.size()-1;++n){
1401  constraint[constraintCounter_].setLinearCoef(x_[neighbours[path[n]][path[n+1]]],1);
1402  }
1403  ++constraintCounter_;
1404  }
1405  }
1406  }
1407  if(constraintCounter_-tempConstrainCounter >= parameter_.maximalNumberOfConstraintsPerRound_)
1408  break;
1409  }
1410  }
1411  else{
1412  for(size_t i=0; i<numberOfInternalEdges_;++i) {
1413 
1414  IndexType u = edgeNodes_[i].first;//[0];
1415  IndexType v = edgeNodes_[i].second;//[1];
1416 
1417  if(usePreBounding && partit[u] != partit[v])
1418  continue;
1419 
1420  OPENGM_ASSERT(i+numberOfTerminalEdges_ == neighbours[u][v]);
1421  OPENGM_ASSERT(i+numberOfTerminalEdges_ == neighbours[v][u]);
1422 
1423  if(sol_[numberOfTerminalEdges_+i]>EPS_){
1424  //search for cycle
1425  std::vector<IndexType> path;
1426  double pathLength;
1427  pathLength = shortestPath(u,v,neighbours,sol_,path,sol_[numberOfTerminalEdges_+i],addOnlyFacetDefiningConstraints);
1428  if(sol_[numberOfTerminalEdges_+i]-EPS_>pathLength){
1429  OPENGM_ASSERT(path.size()>2);
1430  constraint.add(IloRange(env_, 0 , 1000000000));
1431  //negative zero seemed to be required for numerical reasons, even CPlex handel this by its own, too.
1432  constraint[constraintCounter_].setLinearCoef(x_[numberOfTerminalEdges_+i],-1);
1433  for(size_t n=0;n<path.size()-1;++n){
1434  constraint[constraintCounter_].setLinearCoef(x_[neighbours[path[n]][path[n+1]]],1);
1435  }
1436  ++constraintCounter_;
1437  }
1438  }
1439  if(constraintCounter_-tempConstrainCounter >= parameter_.maximalNumberOfConstraintsPerRound_)
1440  break;
1441  }
1442  }
1443  return constraintCounter_-tempConstrainCounter;
1444 }
1445 
1446 template<class GM, class ACC>
1447 size_t Multicut<GM, ACC>::findOddWheelConstraints(IloRangeArray& constraints){
1448  size_t tempConstrainCounter = constraintCounter_;
1449  std::vector<IndexType> var2node(gm_.numberOfVariables(),std::numeric_limits<IndexType>::max());
1450  for(size_t center=0; center<gm_.numberOfVariables();++center){
1451  var2node.assign(gm_.numberOfVariables(),std::numeric_limits<IndexType>::max());
1452  size_t N = neighbours[center].size();
1453  std::vector<IndexType> node2var(N);
1454  std::vector<EdgeMapType> E(2*N);
1455  std::vector<double> w;
1456  typename EdgeMapType::const_iterator it;
1457  size_t id=0;
1458  for(it=neighbours[center].begin() ; it != neighbours[center].end(); ++it) {
1459  IndexType var = (*it).first;
1460  node2var[id] = var;
1461  var2node[var] = id++;
1462  }
1463 
1464  for(it=neighbours[center].begin() ; it != neighbours[center].end(); ++it) {
1465  const IndexType var1 = (*it).first;
1466  const LPIndexType u = var2node[var1];
1467  typename EdgeMapType::const_iterator it2;
1468  for(it2=neighbours[var1].begin() ; it2 != neighbours[var1].end(); ++it2) {
1469  const IndexType var2 = (*it2).first;
1470  const LPIndexType v = var2node[var2];
1471  if( v != std::numeric_limits<IndexType>::max()){
1472  if(u<v){
1473  E[2*u][2*v+1]=w.size();
1474  E[2*v+1][2*u]=w.size();
1475  E[2*u+1][2*v]=w.size();
1476  E[2*v][2*u+1]=w.size();
1477  double weight = 0.5-sol_[neighbours[var1][var2]]+0.5*(sol_[neighbours[center][var1]]+sol_[neighbours[center][var2]]);
1478  //OPENGM_ASSERT(weight>-1e-7);
1479  if(weight<0) weight=0;
1480  w.push_back(weight);
1481 
1482  }
1483  }
1484  }
1485  }
1486 
1487  //Search for odd wheels
1488  for(size_t n=0; n<N;++n) {
1489  std::vector<IndexType> path;
1490  double pathLength;
1491  if(!parameter_.useOldPriorityQueue_){
1492  std::vector<IndexType> prev(2*N);
1493  opengm::ChangeablePriorityQueue<double> openNodes(2*N);
1494  pathLength = shortestPath2(2*n,2*n+1,E,w,path,prev,openNodes,1e22,false);
1495  }else{
1496  pathLength = shortestPath(2*n,2*n+1,E,w,path,1e22,false);
1497  }
1498  if(pathLength<0.5-EPS_*path.size()){// && (path.size())>3){
1499  OPENGM_ASSERT((path.size())>3);
1500  OPENGM_ASSERT(((path.size())/2)*2 == path.size() );
1501 
1502  constraints.add(IloRange(env_, -100000 , (path.size()-2)/2 ));
1503  for(size_t k=0;k<path.size()-1;++k){
1504  constraints[constraintCounter_].setLinearCoef(x_[neighbours[center][node2var[path[k]/2]]],-1);
1505  }
1506  for(size_t k=0;k<path.size()-1;++k){
1507  const IndexType u= node2var[path[k]/2];
1508  const IndexType v= node2var[path[k+1]/2];
1509  constraints[constraintCounter_].setLinearCoef(x_[neighbours[u][v]],1);
1510  }
1511  ++constraintCounter_;
1512  }
1513  if(constraintCounter_-tempConstrainCounter >= parameter_.maximalNumberOfConstraintsPerRound_)
1514  break;
1515  }
1516 
1517  //Reset var2node
1518  for(it=neighbours[center].begin() ; it != neighbours[center].end(); ++it) {
1519  var2node[(*it).first] = std::numeric_limits<IndexType>::max();
1520  }
1521 
1522  }
1523 
1524  return constraintCounter_-tempConstrainCounter;
1525 }
1526 
1532 template<class GM, class ACC>
1533 size_t Multicut<GM, ACC>::findIntegerCycleConstraints(
1534  IloRangeArray& constraint,
1535  bool addFacetDefiningConstraintsOnly
1536 )
1537 {
1538  OPENGM_ASSERT(integerMode_);
1539  std::vector<LabelType> partit(numberOfNodes_,0);
1540  std::vector<std::list<size_t> > neighbours0;
1541  partition(partit,neighbours0);
1542  size_t tempConstrainCounter = constraintCounter_;
1543 
1544  //Find Violated Cycles inside a Partition
1545  size_t u,v;
1546  opengm::FiFoQueue<size_t> nodeList(numberOfNodes_);
1547  std::vector<size_t> path(numberOfNodes_,std::numeric_limits<size_t>::max());
1548  for(size_t i=0; i<numberOfInternalEdges_;++i) {
1549  u = edgeNodes_[i].first;//[0];
1550  v = edgeNodes_[i].second;//[1];
1551  OPENGM_ASSERT(partit[u] >= 0);
1552  if(sol_[numberOfTerminalEdges_+i]>=EPS_ && partit[u] == partit[v]) {
1553  //find shortest path from u to v by BFS
1554  //std::queue<size_t> nodeList;
1555  //opengm::FiFoQueue<size_t> nodeList(numberOfNodes_);
1556  //std::vector<size_t> path(numberOfNodes_,std::numeric_limits<size_t>::max());
1557  nodeList.clear();
1558  std::fill(path.begin(),path.end(),std::numeric_limits<size_t>::max());
1559  size_t n = u;
1560  const bool preChordalcheck =addFacetDefiningConstraintsOnly && parameter_.useChordalSearch_;
1561  while(n!=v) {
1562  std::list<size_t>::iterator it;
1563  for(it=neighbours0[n].begin() ; it != neighbours0[n].end(); ++it) {
1564  if(path[*it]==std::numeric_limits<size_t>::max()) {
1565  //Check if this path is chordless
1566  if(preChordalcheck) {
1567  bool isChordless = true;
1568  size_t s = n;
1569  const size_t c = *it;
1570  while(s!=u){
1571  s = path[s];
1572  if(s==u && c==v) continue;
1573  if(neighbours[c].find(s)!=neighbours[c].end()) {
1574  isChordless = false;
1575  break;
1576  }
1577  }
1578  if(isChordless){
1579  path[*it]=n;
1580  nodeList.push(*it);
1581  }
1582  }
1583  else{
1584  path[*it]=n;
1585  nodeList.push(*it);
1586  }
1587  }
1588  }
1589  //if(nodeList.size()==0)
1590  if(nodeList.empty())
1591  break;
1592  n = nodeList.front(); nodeList.pop();
1593  }
1594  if(path[v] != std::numeric_limits<size_t>::max()){
1595  if(!integerMode_){//check if it is realy violated
1596  double w=0;
1597  while(n!=u) {
1598  w += sol_[neighbours[n][path[n]]];
1599  n=path[n];
1600  }
1601  if(sol_[neighbours[u][v]]-EPS_<w)//constraint is not violated
1602  continue;
1603  }
1604  const bool postChordlessCheck = addFacetDefiningConstraintsOnly && !parameter_.useChordalSearch_;
1605  bool chordless = true;
1606  if(postChordlessCheck){
1607  size_t s = v;
1608  while(s!=u){
1609  if(path[s]==u)
1610  break;
1611  size_t t = path[path[s]];
1612  while(true){
1613  if(s==v && t==u) break;
1614  if(neighbours[t].find(s)!=neighbours[t].end()) {
1615  chordless = false;
1616  break;
1617  }
1618  if(t==u) break;
1619  t=path[t];
1620  }
1621  if(!chordless)
1622  break;
1623  s=path[s];
1624  }
1625  }
1626  if(chordless){
1627  constraint.add(IloRange(env_, 0 , 1000000000));
1628  constraint[constraintCounter_].setLinearCoef(x_[neighbours[u][v]],-1);
1629  while(n!=u) {
1630  constraint[constraintCounter_].setLinearCoef(x_[neighbours[n][path[n]]],1);
1631  n=path[n];
1632  }
1633  ++constraintCounter_;
1634  }
1635  }
1636  if(constraintCounter_-tempConstrainCounter >= parameter_.maximalNumberOfConstraintsPerRound_)
1637  break;
1638  }
1639  if(constraintCounter_-tempConstrainCounter >= parameter_.maximalNumberOfConstraintsPerRound_)
1640  break;
1641  }
1642  return constraintCounter_-tempConstrainCounter;
1643 }
1644 //************************************************
1645 
1646 template <class GM, class ACC>
1647 InferenceTermination
1648 Multicut<GM,ACC>::infer()
1649 {
1650  EmptyVisitorType mcv;
1651  return infer(mcv);
1652 }
1653 
1654 template <class GM, class ACC>
1655 void
1656 Multicut<GM,ACC>::initCplex()
1657 {
1658 
1659  cplex_.setParam(IloCplex::Threads, parameter_.numThreads_);
1660  cplex_.setParam(IloCplex::CutUp, parameter_.cutUp_);
1661  cplex_.setParam(IloCplex::MIPDisplay,0);
1662  cplex_.setParam(IloCplex::BarDisplay,0);
1663  cplex_.setParam(IloCplex::NetDisplay,0);
1664  cplex_.setParam(IloCplex::SiftDisplay,0);
1665  cplex_.setParam(IloCplex::SimDisplay,0);
1666 
1667  cplex_.setParam(IloCplex::EpOpt,1e-9);
1668  cplex_.setParam(IloCplex::EpRHS,1e-8); //setting this to 1e-9 seemed to be to agressive!
1669  cplex_.setParam(IloCplex::EpInt,0);
1670  cplex_.setParam(IloCplex::EpAGap,0);
1671  cplex_.setParam(IloCplex::EpGap,0);
1672 
1673  if(parameter_.verbose_ == true && parameter_.verboseCPLEX_) {
1674  cplex_.setParam(IloCplex::MIPDisplay,2);
1675  cplex_.setParam(IloCplex::BarDisplay,1);
1676  cplex_.setParam(IloCplex::NetDisplay,1);
1677  cplex_.setParam(IloCplex::SiftDisplay,1);
1678  cplex_.setParam(IloCplex::SimDisplay,1);
1679  }
1680 }
1681 
1682 
1683 template <class GM, class ACC>
1684 template<class VisitorType>
1685 InferenceTermination
1686 Multicut<GM,ACC>::infer(VisitorType& mcv)
1687 {
1688  std::vector<LabelType> conf(gm_.numberOfVariables());
1689  initCplex();
1690  //cplex_.setParam(IloCplex::RootAlg, IloCplex::Primal);
1691 
1692  if(problemType_ == INVALID){
1693  throw RuntimeError("Error: Model can not be solved!");
1694  }
1695  else if(!readWorkFlow(parameter_.workFlow_)){//Use given workflow if posible
1696  if(parameter_.workFlow_.size()>0){
1697  std::cout << "Warning: can not parse workflow : " << parameter_.workFlow_ <<std::endl;
1698  std::cout << "Using default workflow ";
1699  }
1700  if(problemType_ == MWC){
1701  //std::cout << "(TTC)(MTC)(IC)(CC-IFD,TTC-I)" <<std::endl;
1702  readWorkFlow("(TTC)(MTC)(IC)(CC-IFD,TTC-I)");
1703  }
1704  else if(problemType_ == MC){
1705  //std::cout << "(CC-FDB)(IC)(CC-I)" <<std::endl;
1706  readWorkFlow("(CC-FDB)(IC)(CC-I)");
1707  }
1708  else{
1709  throw RuntimeError("Error: Model can not be solved!");
1710  }
1711  }
1712 
1713  Timer timer,timer2;
1714  timer.tic();
1715  mcv.begin(*this);
1716  size_t workingState = 0;
1717  while(workingState<workFlow_.size()){
1718  protocolateTiming_.push_back(std::vector<double>(20,0));
1719  protocolateConstraints_.push_back(std::vector<size_t>(20,0));
1720  std::vector<double>& T = protocolateTiming_.back();
1721  std::vector<size_t>& C = protocolateConstraints_.back();
1722 
1723  // Check for timeout
1724  timer.toc();
1725  if(timer.elapsedTime()>parameter_.timeOut_) {
1726  break;
1727  }
1728  //check for integer constraints
1729  for (size_t it=1; it<10000000000; ++it) {
1730  cplex_.setParam(IloCplex::Threads, parameter_.numThreads_);
1731  cplex_.setParam(IloCplex::TiLim, parameter_.timeOut_-timer.elapsedTime());
1732  timer2.tic();
1733  if(!cplex_.solve()) {
1734  if(cplex_.getStatus() != IloAlgorithm::Unbounded){
1735  std::cout << "failed to optimize. " <<cplex_.getStatus()<< std::endl;
1736  //Serious problem -> exit
1737  mcv(*this);
1738  return NORMAL;
1739  }
1740  else{
1741  //undbounded ray - most likely numerical problems
1742  }
1743  }
1744  if(cplex_.getStatus()!= IloAlgorithm::Unbounded){
1745  if(!integerMode_)
1746  bound_ = cplex_.getObjValue()+constant_;
1747  else{
1748  bound_ = cplex_.getBestObjValue()+constant_;
1749  if(!cplex_.solveFixed()) {
1750  std::cout << "failed to fixed optimize." << std::endl;
1751  mcv(*this);
1752  return NORMAL;
1753  }
1754  }
1755  }
1756  else{
1757  //bound is not set - todo
1758  }
1759  try{ cplex_.getValues(sol_, x_);}
1760  catch(IloAlgorithm::NotExtractedException e) {
1761  //std::cout << "UPS: solution not extractable due to unbounded dual ... solving this"<<std::endl;
1762  // The following code is very ugly -> todo: using rays instead
1763  sol_.clear();
1764  for(IndexType v=0; v<numberOfTerminalEdges_+numberOfInternalEdges_+numberOfInterTerminalEdges_ + numberOfHigherOrderValues_; ++v){
1765  try{
1766  sol_.add(cplex_.getValue(x_[v]));
1767  } catch(IloAlgorithm::NotExtractedException e) {
1768  sol_.add(0);
1769  }
1770  }
1771  }
1772  if(parameter_.useBufferedStates_){
1773  std::vector<LabelType> s(gm_.numberOfVariables());
1774  parameter_.useBufferedStates_ = false;
1775  arg(s);
1776  parameter_.useBufferedStates_ = true;
1777  ValueType v = gm_.evaluate(s);
1778  if(bufferedValue_ > v){
1779  bufferedValue_ = v;
1780  bufferedStates_.assign(s.begin(), s.end());
1781  }
1782  }
1783 
1784  timer2.toc();
1785  T[Protocol_ID_Solve] += timer2.elapsedTime();
1786  if(mcv(*this)!=0){
1787  workingState = workFlow_.size(); // go to the end of the workflow
1788  break;
1789  }
1790 
1791  //std::cout << "... done."<<std::endl;
1792 
1793  //Find Violated Constraints
1794  IloRangeArray constraint = IloRangeArray(env_);
1795  constraintCounter_ = 0;
1796 
1797  //std::cout << "Search violated constraints ..." <<std::endl;
1798 
1799 
1800  size_t cycleConstraints = std::numeric_limits<size_t>::max();
1801  bool constraintAdded = false;
1802  for(std::vector<size_t>::iterator it=workFlow_[workingState].begin() ; it != workFlow_[workingState].end(); it++ ){
1803  size_t n = 0;
1804  size_t protocol_ID = Protocol_ID_Unknown;
1805  timer2.tic();
1806  if(*it == Action_ID_TTC){
1807  if(parameter_.verbose_) std::cout << "* Add terminal triangle constraints: " << std::flush;
1808  n = findTerminalTriangleConstraints(constraint);
1809  if(parameter_.verbose_) std::cout << n << std::endl;
1810  protocol_ID = Protocol_ID_TTC;
1811  }
1812  else if(*it == Action_ID_TTC_I){
1813  if(!integerMode_){
1814  throw RuntimeError("Error: Calling integer terminal triangle constraint (TTC-I) seperation provedure before switching in integer mode (IC)");
1815  }
1816  if(parameter_.verbose_) std::cout << "* Add integer terminal triangle constraints: " << std::flush;
1817  arg(conf);
1818  n = findIntegerTerminalTriangleConstraints(constraint, conf);
1819  if(parameter_.verbose_) std::cout << n << std::endl;
1820  protocol_ID = Protocol_ID_TTC;
1821 
1822  }
1823  else if(*it == Action_ID_MTC){
1824  if(parameter_.verbose_) std::cout << "* Add multi terminal constraints: " << std::flush;
1825  n = findMultiTerminalConstraints(constraint);
1826  if(parameter_.verbose_) std::cout << n << std::endl;
1827  protocol_ID = Protocol_ID_MTC;
1828 
1829  }
1830  else if(*it == Action_ID_CC_I){
1831  if(!integerMode_){
1832  throw RuntimeError("Error: Calling integer cycle constraints (CC-I) seperation provedure before switching in integer mode (IC)");
1833  }
1834  if(parameter_.verbose_) std::cout << "Add integer cycle constraints: " << std::flush;
1835  n = findIntegerCycleConstraints(constraint, false);
1836  if(parameter_.verbose_) std::cout << n << std::endl;
1837  protocol_ID = Protocol_ID_CC;
1838  }
1839  else if(*it == Action_ID_CC_IFD){
1840  if(!integerMode_){
1841  throw RuntimeError("Error: Calling facet defining integer cycle constraints (CC-IFD) seperation provedure before switching in integer mode (IC)");
1842  }
1843  if(parameter_.verbose_) std::cout << "Add facet defining integer cycle constraints: " << std::flush;
1844  n = findIntegerCycleConstraints(constraint, true);
1845  if(parameter_.verbose_) std::cout << n << std::endl;
1846  protocol_ID = Protocol_ID_CC;
1847  }
1848  else if(*it == Action_ID_CC){
1849  if(parameter_.verbose_) std::cout << "Add cycle constraints: " << std::flush;
1850  n = findCycleConstraints(constraint, false, false);
1851  cycleConstraints=n;
1852  if(parameter_.verbose_) std::cout << n << std::endl;
1853  protocol_ID = Protocol_ID_CC;
1854  }
1855  else if(*it == Action_ID_CC_FD){
1856  if(parameter_.verbose_) std::cout << "Add facet defining cycle constraints: " << std::flush;
1857  n = findCycleConstraints(constraint, true, false);
1858  cycleConstraints=n;
1859  if(parameter_.verbose_) std::cout << n << std::endl;
1860  protocol_ID = Protocol_ID_CC;
1861  }
1862  else if(*it == Action_ID_CC_FDB){
1863  if(parameter_.verbose_) std::cout << "Add facet defining cycle constraints (with bounding): " << std::flush;
1864  n = findCycleConstraints(constraint, true, true);
1865  cycleConstraints=n;
1866  if(parameter_.verbose_) std::cout << n << std::endl;
1867  protocol_ID = Protocol_ID_CC;
1868  }
1869  else if(*it == Action_ID_CC_B){
1870  if(parameter_.verbose_) std::cout << "Add cycle constraints (with bounding): " << std::flush;
1871  n = findCycleConstraints(constraint, false, true);
1872  cycleConstraints=n;
1873  if(parameter_.verbose_) std::cout << n << std::endl;
1874  protocol_ID = Protocol_ID_CC;
1875  }
1876  else if(*it == Action_ID_RemoveConstraints){
1877  if(parameter_.verbose_) std::cout << "Remove unused constraints: " << std::flush;
1878  n = removeUnusedConstraints();
1879  if(parameter_.verbose_) std::cout << n << std::endl;
1880  protocol_ID = Protocol_ID_RemoveConstraints;
1881  }
1882  else if(*it == Action_ID_IntegerConstraints){
1883  if(integerMode_) continue;
1884  if(parameter_.verbose_) std::cout << "Add integer constraints: " << std::flush;
1885  n = enforceIntegerConstraints();
1886  if(parameter_.verbose_) std::cout << n << std::endl;
1887  protocol_ID = Protocol_ID_IntegerConstraints;
1888  }
1889  else if(*it == Action_ID_OWC){
1890  if(cycleConstraints== std::numeric_limits<size_t>::max()){
1891  std::cout << "WARNING: using Odd Wheel Contraints without Cycle Constrains! -> Use CC,OWC!"<<std::endl;
1892  n=0;
1893  }
1894  else if(cycleConstraints==0){
1895  if(parameter_.verbose_) std::cout << "Add odd wheel constraints: " << std::flush;
1896  n = findOddWheelConstraints(constraint);
1897  if(parameter_.verbose_) std::cout << n << std::endl;
1898  }
1899  else{
1900  //since cycle constraints are found we have to search for more violated cycle constraints first
1901  }
1902  protocol_ID = Protocol_ID_OWC;
1903  }
1904  else{
1905  std::cout <<"Unknown Inference State "<<*it<<std::endl;
1906  }
1907  timer2.toc();
1908  T[protocol_ID] += timer2.elapsedTime();
1909  C[protocol_ID] += n;
1910  if(n>0) constraintAdded = true;
1911  }
1912  //std::cout <<"... done!"<<std::endl;
1913 
1914 
1915 
1916  if(!constraintAdded){
1917  //Switch to next working state
1918  ++workingState;
1919  if(workingState<workFlow_.size())
1920  if(parameter_.verbose_) std::cout <<std::endl<< "** Switching into next state ( "<< workingState << " )**" << std::endl;
1921  break;
1922  }
1923  else{
1924  timer2.tic();
1925  //Add Constraints
1926  model_.add(constraint);
1927  //cplex_.addCuts(constraint);
1928  timer2.toc();
1929  T[Protocol_ID_AddConstraints] += timer2.elapsedTime();
1930  }
1931 
1932  // Check for timeout
1933  timer.toc();
1934  if(timer.elapsedTime()>parameter_.timeOut_) {
1935  break;
1936  }
1937 
1938  } //end inner loop over one working state
1939  } //end loop over all working states
1940 
1941  mcv.end(*this);
1942  if(parameter_.verbose_){
1943  std::cout << "end normal"<<std::endl;
1944  std::cout <<"Protokoll:"<<std::endl;
1945  std::cout <<"=========="<<std::endl;
1946  std::cout << " i | SOLVE | ADD | CC | OWC | TTC | MTV | IC " <<std::endl;
1947  std::cout << "-----+----------+----------+----------+----------+----------+----------+-----------" <<std::endl;
1948  }
1949  std::vector<size_t> IDS;
1950  IDS.push_back(Protocol_ID_Solve);
1951  IDS.push_back(Protocol_ID_AddConstraints);
1952  IDS.push_back(Protocol_ID_CC);
1953  IDS.push_back(Protocol_ID_OWC);
1954  IDS.push_back(Protocol_ID_TTC);
1955  IDS.push_back(Protocol_ID_MTC);
1956  IDS.push_back(Protocol_ID_IntegerConstraints);
1957 
1958  if(parameter_.verbose_){
1959  for(size_t i=0; i<protocolateTiming_.size(); ++i){
1960  std::cout << std::setw(5)<< i ;
1961  for(size_t n=0; n<IDS.size(); ++n){
1962  std::cout << "|"<< std::setw(10) << setiosflags(std::ios::fixed)<< std::setprecision(4) << protocolateConstraints_[i][IDS[n]];
1963  }
1964  std::cout << std::endl;
1965  std::cout << " " ;
1966  for(size_t n=0; n<IDS.size(); ++n){
1967  std::cout << "|"<< std::setw(10) << setiosflags(std::ios::fixed)<< std::setprecision(4) << protocolateTiming_[i][IDS[n]];
1968  }
1969  std::cout << std::endl;
1970  std::cout << "-----+----------+----------+----------+----------+----------+----------+-----------" <<std::endl;
1971  }
1972  std::cout << "sum ";
1973  double t_all=0;
1974  for(size_t n=0; n<IDS.size(); ++n){
1975  double t_one=0;
1976  for(size_t i=0; i<protocolateTiming_.size(); ++i){
1977  t_one += protocolateTiming_[i][IDS[n]];
1978  }
1979  std::cout << "|"<< std::setw(10) << setiosflags(std::ios::fixed)<< std::setprecision(4) << t_one;
1980  t_all += t_one;
1981  }
1982  std::cout << " | " <<t_all <<std::endl;
1983  std::cout << "-----+----------+----------+----------+----------+----------+----------+-----------" <<std::endl;
1984  }
1985  return NORMAL;
1986 }
1987 
1988 template <class GM, class ACC>
1989 InferenceTermination
1990 Multicut<GM,ACC>::arg
1991 (
1992  std::vector<typename Multicut<GM,ACC>::LabelType>& x,
1993  const size_t N
1994  ) const
1995 {
1996  if(N!=1) {
1997  return UNKNOWN;
1998  }
1999  else{
2000  if(parameter_.useBufferedStates_){
2001  x.assign(bufferedStates_.begin(),bufferedStates_.end());
2002  return NORMAL;
2003  }
2004 
2005  if(problemType_ == MWC) {
2006  if(parameter_.MWCRounding_== parameter_.NEAREST){
2007  x.resize(gm_.numberOfVariables());
2008  for(IndexType node = 0; node<gm_.numberOfVariables(); ++node) {
2009  double v = sol_[numberOfTerminals_*node+0];
2010  x[node] = 0;
2011  for(LabelType i=0; i<gm_.numberOfLabels(node); ++i) {
2012  if(sol_[numberOfTerminals_*node+i]<v) {
2013  x[node]=i;
2014  }
2015  }
2016  }
2017  return NORMAL;
2018  }
2019  else if(parameter_.MWCRounding_==Parameter::DERANDOMIZED){
2020  return derandomizedRounding(x);
2021  }
2022  else if(parameter_.MWCRounding_==Parameter::PSEUDODERANDOMIZED){
2023  return pseudoDerandomizedRounding(x,1000);
2024  }
2025  else{
2026  return UNKNOWN;
2027  }
2028  }
2029  else{
2030  std::vector<std::list<size_t> > neighbours0;
2031  InferenceTermination r = partition(x, neighbours0,parameter_.edgeRoundingValue_);
2032  return r;
2033  }
2034  }
2035 }
2036 template <class GM, class ACC>
2037 std::vector<size_t>
2038 Multicut<GM,ACC>::getSegmentation() const {
2039  std::vector<size_t> seg;
2040  std::vector<std::list<size_t> > neighbours0;
2041  partition(seg, neighbours0, 0.3);
2042  return seg;
2043 }
2044 
2045 
2046 template <class GM, class ACC>
2047 InferenceTermination
2048 Multicut<GM,ACC>::pseudoDerandomizedRounding
2049 (
2050  std::vector<typename Multicut<GM,ACC>::LabelType>& x,
2051  size_t bins
2052  ) const
2053 {
2054  std::vector<bool> processedBins(bins+1,false);
2055  std::vector<LabelType> temp;
2056  double value = 1000000000000.0;
2057  std::vector<LabelType> labelorder1(numberOfTerminals_,numberOfTerminals_-1);
2058  std::vector<LabelType> labelorder2(numberOfTerminals_,numberOfTerminals_-1);
2059  for(LabelType i=0; i<numberOfTerminals_-1;++i){
2060  labelorder1[i]=i;
2061  labelorder2[i]=numberOfTerminals_-2-i;
2062  }
2063  for(size_t i=0; i<numberOfTerminals_*gm_.numberOfVariables();++i){
2064  size_t bin;
2065  double t,d;
2066  if(sol_[i]<=0){
2067  bin=0;
2068  t=0;
2069  }
2070  else if(sol_[i]>=1){
2071  bin=bins;
2072  t=1;
2073  }
2074  else{
2075  bin = (size_t)(sol_[i]*bins)%bins;
2076  t = sol_[i];
2077  }
2078  if(!processedBins[bin]){
2079  processedBins[bin]=true;
2080  if(value>(d=derandomizedRoundingSubProcedure(temp,labelorder1,sol_[i]))){
2081  value=d;
2082  x=temp;
2083  }
2084  if(value>(d=derandomizedRoundingSubProcedure(temp,labelorder2,sol_[i]))){
2085  value=d;
2086  x=temp;
2087  }
2088  }
2089  }
2090  return NORMAL;
2091 }
2092 
2093 template <class GM, class ACC>
2094 InferenceTermination
2095 Multicut<GM,ACC>::derandomizedRounding
2096 (
2097  std::vector<typename Multicut<GM,ACC>::LabelType>& x
2098  ) const
2099 {
2100  std::vector<LabelType> temp;
2101  double value = 1000000000000.0;
2102  std::vector<LabelType> labelorder1(numberOfTerminals_,numberOfTerminals_-1);
2103  std::vector<LabelType> labelorder2(numberOfTerminals_,numberOfTerminals_-1);
2104  for(LabelType i=0; i<numberOfTerminals_-1;++i){
2105  labelorder1[i]=i;
2106  labelorder2[i]=numberOfTerminals_-2-i;
2107  }
2108  // Test primitives
2109  double d;
2110  if(value>(d=derandomizedRoundingSubProcedure(temp,labelorder1,1e-8))){
2111  value=d;
2112  x=temp;
2113  }
2114  if(value>(d=derandomizedRoundingSubProcedure(temp,labelorder2,1e-8))){
2115  value=d;
2116  x=temp;
2117  }
2118  if(value>(d=derandomizedRoundingSubProcedure(temp,labelorder1,1-1e-8))){
2119  value=d;
2120  x=temp;
2121  }
2122  if(value>(d=derandomizedRoundingSubProcedure(temp,labelorder2,1-1e-8))){
2123  value=d;
2124  x=temp;
2125  }
2126  for(size_t i=0; i<numberOfTerminals_*gm_.numberOfVariables();++i){
2127  if( sol_[i]>1e-8 && sol_[i]<1-1e-8){
2128  if(value>(d=derandomizedRoundingSubProcedure(temp,labelorder1,sol_[i]))){
2129  value=d;
2130  x=temp;
2131  }
2132  if(value>(d=derandomizedRoundingSubProcedure(temp,labelorder2,sol_[i]))){
2133  value=d;
2134  x=temp;
2135  }
2136  }
2137  }
2138  return NORMAL;
2139 }
2140 
2141 template <class GM, class ACC>
2142 double
2143 Multicut<GM,ACC>::derandomizedRoundingSubProcedure
2144 (
2145  std::vector<typename Multicut<GM,ACC>::LabelType>& x,
2146  const std::vector<typename Multicut<GM,ACC>::LabelType>& labelorder,
2147  const double threshold
2148  ) const
2149 {
2150  const LabelType lastLabel = labelorder.back();
2151  const IndexType numVar = gm_.numberOfVariables();
2152 
2153  x.assign(numVar,lastLabel);
2154 
2155  for(size_t i=0; i<labelorder.size()-1; ++i){
2156  for(IndexType v=0; v<numVar; ++v){
2157  if(x[v]==lastLabel && sol_[numberOfTerminals_*v+i]<=threshold){
2158  x[v]=labelorder[i];
2159  }
2160  }
2161  }
2162  return gm_.evaluate(x);
2163 }
2164 
2165 
2166 
2167 
2168 template <class GM, class ACC>
2169 InferenceTermination
2170 Multicut<GM,ACC>::partition
2171 (
2172  std::vector<LabelType>& partit,
2173  std::vector<std::list<size_t> >& neighbours0,
2174  double threshold
2175  ) const
2176 {
2177 
2178  partit.resize(numberOfNodes_,0);
2179  neighbours0.resize(numberOfNodes_, std::list<size_t>());
2180 
2181  size_t counter=0;
2182  for(size_t i=0; i<numberOfInternalEdges_; ++i) {
2183  IndexType u = edgeNodes_[i].first;//variableIndex(0);
2184  IndexType v = edgeNodes_[i].second;//variableIndex(1);
2185  if(sol_[numberOfTerminalEdges_+counter] <= threshold) {
2186  neighbours0[u].push_back(v);
2187  neighbours0[v].push_back(u);
2188  }
2189  ++counter;
2190  }
2191 
2192  LabelType p=0;
2193  std::vector<bool> assigned(numberOfNodes_,false);
2194  for(size_t node=0; node<numberOfNodes_; ++node) {
2195  if(assigned[node])
2196  continue;
2197  else{
2198  std::list<size_t> nodeList;
2199  partit[node] = p;
2200  assigned[node] = true;
2201  nodeList.push_back(node);
2202  while(!nodeList.empty()) {
2203  size_t n=nodeList.front(); nodeList.pop_front();
2204  std::list<size_t>::const_iterator it;
2205  for(it=neighbours0[n].begin() ; it != neighbours0[n].end(); ++it) {
2206  if(!assigned[*it]) {
2207  partit[*it] = p;
2208  assigned[*it] = true;
2209  nodeList.push_back(*it);
2210  }
2211  }
2212  }
2213  ++p;
2214  }
2215  }
2216  return NORMAL;
2217 }
2218 
2219 
2220 template <class GM, class ACC>
2221 typename Multicut<GM,ACC>::ValueType
2222 Multicut<GM,ACC>::evaluate
2223 (
2224  std::vector<LabelType>& conf
2225  ) const
2226 {
2227 
2228  return gm_.evaluate(conf);
2229 }
2230 
2231 template<class GM, class ACC>
2232 inline const typename Multicut<GM, ACC>::GraphicalModelType&
2233 Multicut<GM, ACC>::graphicalModel() const
2234 {
2235  return gm_;
2236 }
2237 
2238 template<class GM, class ACC>
2239 typename GM::ValueType Multicut<GM, ACC>::bound() const
2240 {
2241  return bound_;
2242 }
2243 
2244 template<class GM, class ACC>
2245 typename GM::ValueType Multicut<GM, ACC>::value() const
2246 {
2247  std::vector<LabelType> c;
2248  arg(c);
2249  ValueType value = gm_.evaluate(c);
2250  return value;
2251 }
2252 
2253 template<class GM, class ACC>
2254 bool Multicut<GM, ACC>::readWorkFlow(std::string s)
2255 {
2256  if(s[0]!='(' || s[s.size()-1] !=')')
2257  return false;
2258  workFlow_.clear();
2259  size_t n=0;
2260  std::string sepString;
2261  if(s.size()<2)
2262  return false;
2263  while(n<s.size()){
2264  if(s[n]==',' || s[n]==')'){//End of sepString
2265  if(sepString.compare("CC")==0)
2266  workFlow_.back().push_back(Action_ID_CC);
2267  else if(sepString.compare("CC-I")==0)
2268  workFlow_.back().push_back(Action_ID_CC_I);
2269  else if(sepString.compare("CC-IFD")==0)
2270  workFlow_.back().push_back(Action_ID_CC_IFD);
2271  else if(sepString.compare("CC-B")==0)
2272  workFlow_.back().push_back(Action_ID_CC_B);
2273  else if(sepString.compare("CC-FDB")==0)
2274  workFlow_.back().push_back(Action_ID_CC_FDB);
2275  else if(sepString.compare("CC-FD")==0)
2276  workFlow_.back().push_back(Action_ID_CC_FD);
2277  else if(sepString.compare("TTC")==0)
2278  workFlow_.back().push_back(Action_ID_TTC);
2279  else if(sepString.compare("TTC-I")==0)
2280  workFlow_.back().push_back(Action_ID_TTC_I);
2281  else if(sepString.compare("MTC")==0)
2282  workFlow_.back().push_back(Action_ID_MTC);
2283  else if(sepString.compare("OWC")==0)
2284  workFlow_.back().push_back(Action_ID_OWC);
2285  else if(sepString.compare("IC")==0)
2286  workFlow_.back().push_back(Action_ID_IntegerConstraints);
2287  else
2288  std::cout << "WARNING: Unknown Seperation Procedure ' "<<sepString<<"' is skipped!"<<std::endl;
2289  sepString.clear();
2290  }
2291  else if(s[n]=='('){//New Round
2292  workFlow_.push_back(std::vector<size_t>());
2293  }
2294  else{
2295  sepString += s[n];
2296  }
2297  ++n;
2298  }
2299  return true;
2300 }
2301 
2302 
2308 template<class GM, class ACC>
2309 template<class DOUBLEVECTOR>
2310 inline double Multicut<GM, ACC>::shortestPath(
2311  const IndexType startNode,
2312  const IndexType endNode,
2313  const std::vector<EdgeMapType >& E, //E[n][i].first/.second are the i-th neighbored node and weight-index (for w), respectively.
2314  const DOUBLEVECTOR& w,
2315  std::vector<IndexType>& shortestPath,
2316  const double maxLength,
2317  bool chordless
2318 ) const
2319 {
2320 
2321  const IndexType numberOfNodes = E.size();
2322  const double inf = std::numeric_limits<double>::infinity();
2323  const IndexType nonePrev = endNode;
2324 
2325  std::vector<IndexType> prev(numberOfNodes,nonePrev);
2326  std::vector<double> dist(numberOfNodes,inf);
2327  MYSET openNodes;
2328 
2329  openNodes.insert(startNode);
2330  dist[startNode]=0.0;
2331 
2332  while(!openNodes.empty()){
2333  IndexType node;
2334  //Find smallest open node
2335  {
2336  typename MYSET::iterator it, itMin;
2337  double v = std::numeric_limits<double>::infinity();
2338  for(it = openNodes.begin(); it!= openNodes.end(); ++it){
2339  if( dist[*it]<v ){
2340  v = dist[*it];
2341  itMin = it;
2342  }
2343  }
2344  node = *itMin;
2345  openNodes.erase(itMin);
2346  }
2347  // Check if target is reached
2348  if(node == endNode)
2349  break;
2350  // Update all neigbors of node
2351  {
2352  typename EdgeMapType::const_iterator it;
2353  for(it=E[node].begin() ; it != E[node].end(); ++it) {
2354  const IndexType node2 = (*it).first; //second edge-node
2355  const LPIndexType weighId = (*it).second; //index in weigh-vector w
2356  double cuttedWeight = std::max(0.0,w[weighId]); //cut up negative edge-weights
2357  const double weight2 = dist[node]+cuttedWeight;
2358 
2359 
2360  if(dist[node2] > weight2 && weight2 < maxLength){
2361  //check chordality
2362  bool updateNode = true;
2363  if(chordless) {
2364  IndexType s = node;
2365  while(s!=startNode){
2366  s= prev[s];
2367  if(s==startNode && node2==endNode) continue;
2368  if(neighbours[node2].find(s)!=neighbours[node2].end()) {
2369  updateNode = false; // do not update node if path is chordal
2370  break;
2371  }
2372  }
2373  }
2374  if(updateNode){
2375  prev[node2] = node;
2376  dist[node2] = weight2;
2377  openNodes.insert(node2);
2378  }
2379  }
2380  }
2381  }
2382  }
2383 
2384  if(prev[endNode] != nonePrev){//find path?
2385  shortestPath.clear();
2386  shortestPath.push_back(endNode);
2387  IndexType n = endNode;
2388  do{
2389  n=prev[n];
2390  shortestPath.push_back(n);
2391  }while(n!=startNode);
2392  OPENGM_ASSERT(shortestPath.back() == startNode);
2393  }
2394  else{
2395  // std::cout <<"ERROR" <<std::flush;
2396  }
2397 // std::cout <<"*"<< shortestPath.size()<<"-"<<dist[endNode]<<std::flush;
2398  return dist[endNode];
2399 }
2400 
2401 
2402 template<class GM, class ACC>
2403 template<class DOUBLEVECTOR>
2404 inline double Multicut<GM, ACC>::shortestPath2(
2405  const IndexType startNode,
2406  const IndexType endNode,
2407  const std::vector<EdgeMapType >& E, //E[n][i].first/.second are the i-th neighbored node and weight-index (for w), respectively.
2408  const DOUBLEVECTOR& w,
2409  std::vector<IndexType>& shortestPath,
2410  std::vector<IndexType>& prev,
2411  opengm::ChangeablePriorityQueue<double>& openNodes,
2412  const double maxLength,
2413  bool chordless
2414 ) const
2415 {
2416 
2417  const IndexType numberOfNodes = E.size();
2418  const IndexType nonePrev = endNode;
2419 
2420 // std::vector<IndexType> prev(numberOfNodes,nonePrev);
2421 // opengm::ChangeablePriorityQueue<double> openNodes(numberOfNodes);
2422  openNodes.reset();
2423  openNodes.setPriorities(10000);
2424  openNodes.push(startNode,0.0);
2425  prev[endNode] = nonePrev;
2426 
2427  IndexType node;
2428  double priority;
2429  // std::cout <<"1"<<std::flush;
2430  while(!openNodes.empty()){
2431  //Find smallest open node
2432  node = openNodes.top();
2433  priority = openNodes.topPriority();
2434  openNodes.pop();
2435  // std::cout << node << "("<< priority << ") "<<std::flush;
2436 
2437  // Check if target is reached
2438  if(node == endNode)
2439  break;
2440  // Update all neigbors of node
2441  {
2442  typename EdgeMapType::const_iterator it;
2443  for(it=E[node].begin() ; it != E[node].end(); ++it) {
2444  const IndexType node2 = (*it).first; //second edge-node
2445  const LPIndexType weighId = (*it).second; //index in weigh-vector w
2446  double cuttedWeight = std::max(0.0,w[weighId]); //cut up negative edge-weights
2447  const double weight2 = priority+cuttedWeight;
2448 
2449 
2450  //if((!openNodes.contains(node2) || openNodes.priority(node2) > weight2) && weight2 < maxLength ){
2451  if(openNodes.priority(node2) > weight2 && weight2 < maxLength ){
2452  //check chordality
2453  bool updateNode = true;
2454  if(chordless) {
2455  IndexType s = node;
2456  while(s!=startNode){
2457  s= prev[s];
2458  if(s==startNode && node2==endNode) continue;
2459  if(neighbours[node2].find(s)!=neighbours[node2].end()) {
2460  updateNode = false; // do not update node if path is chordal
2461  break;
2462  }
2463  }
2464  }
2465  if(updateNode){
2466  openNodes.push(node2,weight2);
2467  prev[node2] = node;
2468  }
2469  }
2470  }
2471  }
2472  }
2473 // std::cout <<"2"<<std::flush;
2474  if(prev[endNode] != nonePrev){//find path?
2475  shortestPath.resize(0);
2476  shortestPath.push_back(endNode);
2477  IndexType n = endNode;
2478  do{
2479  n=prev[n];
2480  shortestPath.push_back(n);
2481  }while(n!=startNode);
2482  OPENGM_ASSERT(shortestPath.back() == startNode);
2483  }
2484  else{
2485  }
2486 // std::cout <<"*"<< shortestPath.size()<<"-"<<priority<<std::flush;
2487  return openNodes.priority(endNode);
2488 }
2489 
2490 
2491 template<class GM, class ACC>
2492 std::vector<double>
2493 Multicut<GM, ACC>::getEdgeLabeling
2494 () const
2495 {
2496  std::vector<double> l(numberOfInternalEdges_,0);
2497  for(size_t i=0; i<numberOfInternalEdges_; ++i) {
2498  l[i] = sol_[numberOfTerminalEdges_+i];
2499  }
2500  return l;
2501 }
2502 
2503 // WARNING: this function is considered experimental.
2504 // variable indices refer to variables of the LP set up
2505 // in the constructor of the class template LPCplex,
2506 // NOT to the variables of the graphical model.
2507 template<class GM, class ACC>
2508 template<class LPVariableIndexIterator, class CoefficientIterator>
2509 void Multicut<GM, ACC>::addConstraint
2510 (
2511  LPVariableIndexIterator viBegin,
2512  LPVariableIndexIterator viEnd,
2513  CoefficientIterator coefficient,
2514  const ValueType& lowerBound,
2515  const ValueType& upperBound)
2516 {
2517  IloRange constraint(env_, lowerBound, upperBound);
2518  while(viBegin != viEnd) {
2519  constraint.setLinearCoef(x_[*viBegin], *coefficient);
2520  ++viBegin;
2521  ++coefficient;
2522  }
2523  model_.add(constraint);
2524  // this update of model_ does not require a
2525  // re-initialization of the object cplex_.
2526  // cplex_ is initialized in the constructor.
2527 }
2528 
2529 } // end namespace opengm
2530 
2531 #endif
opengm::ChangeablePriorityQueue::reset
void reset()
reset heap - priorities are not changed
Definition: queues.hxx:98
opengm::Multicut::infer
virtual InferenceTermination infer()
Definition: multicut.hxx:1648
opengm::Multicut::bound
ValueType bound() const
return a bound on the solution
Definition: multicut.hxx:2239
opengm
The OpenGM namespace.
Definition: config.hxx:43
opengm::Partitions::getPartition
EdgeLabelType getPartition(size_t n)
Definition: partitions.hxx:19
opengm::Adder
Addition as a binary operation.
Definition: adder.hxx:10
opengm::ChangeablePriorityQueue::pop
void pop()
Remove the current top element.
Definition: queues.hxx:153
opengm::Multicut::constraintCounter_
size_t constraintCounter_
Definition: multicut.hxx:174
opengm::Multicut::arg
virtual InferenceTermination arg(std::vector< LabelType > &, const size_t=1) const
output a solution
Definition: multicut.hxx:1991
opengm::ParamHeper
Multicut Algorithm [1] J. Kappes, M. Speth, B. Andres, G. Reinelt and C. Schnoerr, "Globally Optimal Image Partitioning by Multicuts", EMMCVPR 2011 [2] J. Kappes, M. Speth, G. Reinelt and C. Schnoerr, "Higher-order Segmentation via Multicuts", Technical Report (http://ipa.iwr.uni-heidelberg.de/ipabib/Papers/kappes-2013-multicut.pdf) .
Definition: multicut.hxx:72
opengm::Partitions::resize
void resize(size_t N)
Definition: partitions.hxx:17
opengm::Multicut::graphicalModel
const GraphicalModelType & graphicalModel() const
Definition: multicut.hxx:2233
opengm::ChangeablePriorityQueue::priority
priority_type priority(const value_type i) const
returns the value associated with index i
Definition: queues.hxx:162
opengm::ParamHeper::PSEUDODERANDOMIZED
Definition: multicut.hxx:73
opengm::Multicut::MYSET
__gnu_cxx::hash_set< IndexType > MYSET
Definition: multicut.hxx:94
opengm::Multicut::getEdgeLabeling
std::vector< double > getEdgeLabeling() const
Definition: multicut.hxx:2494
opengm::Partitions::BellNumber
size_t BellNumber(size_t N)
Definition: partitions.hxx:18
opengm::Timer::toc
void toc()
Definition: timer.hxx:109
size_t
opengm::ChangeablePriorityQueue
Heap-based changable priority queue with a maximum number of elemements.
Definition: queues.hxx:66
opengm::Timer
Platform-independent runtime measurements.
Definition: timer.hxx:29
opengm::Multicut::calcBound
ValueType calcBound()
Definition: multicut.hxx:161
OPENGM_ASSERT
#define OPENGM_ASSERT(expression)
Definition: opengm.hxx:77
opengm::Multicut::Parameter::MWCRounding_
ParamHeper::MWCRounding MWCRounding_
Definition: multicut.hxx:115
opengm::ChangeablePriorityQueue::empty
bool empty() const
check if the PQ is empty
Definition: queues.hxx:105
opengm::Multicut::getLPIndex
size_t getLPIndex(IT a, IT b)
Definition: multicut.hxx:171
opengm::UInt64Type
detail_types::UInt64Type UInt64Type
uint64
Definition: config.hxx:300
opengm::Inference::OperatorType
GraphicalModelType::OperatorType OperatorType
Definition: inference.hxx:42
opengm::ChangeablePriorityQueue::push
void push(const value_type i, const priority_type p)
Insert a index with a given priority.
Definition: queues.hxx:126
opengm::ParamHeper::DERANDOMIZED
Definition: multicut.hxx:73
opengm::Inference::IndexType
GraphicalModelType::IndexType IndexType
Definition: inference.hxx:40
opengm::Multicut::name
virtual std::string name() const
Definition: multicut.hxx:154
opengm::Multicut::Parameter::allowCutsWithin_
std::vector< bool > allowCutsWithin_
Definition: multicut.hxx:117
opengm::Timer::elapsedTime
double elapsedTime() const
Definition: timer.hxx:130
opengm::Timer::tic
void tic()
Definition: timer.hxx:96
opengm::ParamHeper::NEAREST
Definition: multicut.hxx:73
opengm::Multicut::EdgeMapType
__gnu_cxx::hash_map< IndexType, LPIndexType > EdgeMapType
Definition: multicut.hxx:93
opengm::ChangeablePriorityQueue::topPriority
priority_type topPriority() const
get top priority
Definition: queues.hxx:147
opengm::Multicut::TimingVisitorType
visitors::TimingVisitor< Multicut< GM, ACC > > TimingVisitorType
Definition: multicut.hxx:86
opengm::Inference::ValueType
GraphicalModelType::ValueType ValueType
Definition: inference.hxx:41
opengm::Inference
Inference algorithm interface.
Definition: inference.hxx:34
opengm::Multicut::~Multicut
virtual ~Multicut()
Definition: multicut.hxx:455
opengm::Multicut::addConstraint
void addConstraint(LPVariableIndexIterator, LPVariableIndexIterator, CoefficientIterator, const ValueType &, const ValueType &)
Definition: multicut.hxx:2510
opengm::Multicut::getSegmentation
std::vector< size_t > getSegmentation() const
Definition: multicut.hxx:2038
opengm::Multicut::Parameter::maximalNumberOfConstraintsPerRound_
size_t maximalNumberOfConstraintsPerRound_
Definition: multicut.hxx:113
opengm::Multicut::LPIndexType
size_t LPIndexType
Definition: multicut.hxx:83
OPENGM_CHECK_OP
#define OPENGM_CHECK_OP(A, OP, B, TXT)
Definition: submodel2.hxx:24
opengm::ParamHeper::MWCRounding
MWCRounding
Definition: multicut.hxx:73
opengm::ChangeablePriorityQueue::setPriorities
void setPriorities(T newPriority)
set all priorities to the given value
Definition: queues.hxx:91
opengm::Multicut::rebind::type
Multicut< GM_, ACC_ > type
Definition: multicut.hxx:100
opengm::Multicut::EmptyVisitorType
visitors::EmptyVisitor< Multicut< GM, ACC > > EmptyVisitorType
Definition: multicut.hxx:85
opengm::Multicut::evaluate
ValueType evaluate(std::vector< LabelType > &) const
Definition: multicut.hxx:2223
opengm::Multicut::Parameter::Parameter
Parameter(int numThreads=0, double cutUp=1.0e+75)
Definition: multicut.hxx:126
opengm::Inference::LabelType
GraphicalModelType::LabelType LabelType
Definition: inference.hxx:39
opengm::Minimizer
Minimization as a unary accumulation.
Definition: minimizer.hxx:12
opengm::Multicut::VerboseVisitorType
visitors::VerboseVisitor< Multicut< GM, ACC > > VerboseVisitorType
Definition: multicut.hxx:84
opengm::RuntimeError
OpenGM runtime error.
Definition: opengm.hxx:100
opengm::Multicut::value
ValueType value() const
return the solution (value)
Definition: multicut.hxx:2245
opengm::Multicut::Multicut
Multicut(const GraphicalModelType &, Parameter para=Parameter())
Definition: multicut.hxx:461
opengm::Multicut::inferenceState_
size_t inferenceState_
Definition: multicut.hxx:171
opengm::ChangeablePriorityQueue::top
const_reference top() const
get index with top priority
Definition: queues.hxx:141
opengm::InferenceTermination
InferenceTermination
Definition: inference.hxx:24
Generated by   doxygen 1.8.9.1