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