OpenGM  2.3.x
Discrete Graphical Model Library
mqpbo.hxx
Go to the documentation of this file.
1 #pragma once
2 #ifndef OPENGM_MQPBO_HXX
3 #define OPENGM_MQPBO_HXX
4 
5 #include <vector>
6 #include <string>
7 #include <iostream>
8 
9 #include "opengm/opengm.hxx"
10 #include "opengm/inference/visitors/visitors.hxx"
11 #include "opengm/inference/inference.hxx"
12 #include <opengm/utilities/metaprogramming.hxx>
13 #include "opengm/utilities/tribool.hxx"
14 #include <opengm/inference/messagepassing/messagepassing.hxx>
15 #include <opengm/functions/view_fix_variables_function.hxx>
16 
17 //#define MQPBOHotFixOutPutPartialOPtimalityMap
18 #ifdef MQPBOHotFixOutPutPartialOPtimalityMap
19 #include <opengm/datastructures/marray/marray_hdf5.hxx>
20 #endif
21 
22 #include "opengm/inference/external/qpbo.hxx"
23 
24 namespace opengm {
25 
36  template<class GM, class ACC>
37  class MQPBO : public Inference<GM, ACC>
38  {
39  public:
40  typedef ACC AccumulationType;
41  typedef GM GmType;
42  typedef GM GraphicalModelType;
43  OPENGM_GM_TYPE_TYPEDEFS;
44  typedef visitors::VerboseVisitor<MQPBO<GM, ACC> > VerboseVisitorType;
45  typedef visitors::EmptyVisitor<MQPBO<GM, ACC> > EmptyVisitorType;
46  typedef visitors::TimingVisitor<MQPBO<GM, ACC> > TimingVisitorType;
47  typedef ValueType GraphValueType;
48 
49  enum PermutationType {NONE, RANDOM, MINMARG};
50 
51  class Parameter{
52  public:
53  Parameter(): useKovtunsMethod_(true), probing_(false), strongPersistency_(false), rounds_(0), permutationType_(NONE) {};
54  std::vector<LabelType> label_;
55  bool useKovtunsMethod_;
56  const bool probing_; //do not use this!
57  bool strongPersistency_;
58  size_t rounds_;
59  PermutationType permutationType_;
60  };
61 
62  MQPBO(const GmType&, const Parameter& = Parameter());
63  ~MQPBO();
64  std::string name() const;
65  const GmType& graphicalModel() const;
66  InferenceTermination infer();
67  void reset();
68  typename GM::ValueType bound() const;
69  typename GM::ValueType value() const;
70  template<class VisitorType>
71  InferenceTermination infer(VisitorType&);
72  InferenceTermination testQuess(std::vector<LabelType> &guess);
73  InferenceTermination testPermutation(PermutationType permutationType);
74  void setStartingPoint(typename std::vector<typename GM::LabelType>::const_iterator);
75  virtual InferenceTermination arg(std::vector<LabelType>&, const size_t = 1) const ;
76 
77  const std::vector<opengm::Tribool>& partialOptimality(IndexType var) const {return partialOptimality_[var];}
78  bool partialOptimality(IndexType var, LabelType& l) const {l=label_[var]; return optimal_[var];}
79 
80  double optimalityV() const;
81  double optimality() const;
82  private:
83  InferenceTermination testQuess(LabelType guess);
84  void AddUnaryTerm(int var, ValueType v0, ValueType v1);
85  void AddPairwiseTerm(int var0, int var1,ValueType v00,ValueType v01,ValueType v10,ValueType v11);
86 
87  const GmType& gm_;
88  Parameter param_;
89 
90  kolmogorov::qpbo::QPBO<GraphValueType>* qpbo_;
91  ValueType constTerm_;
92  ValueType bound_;
93  //int* label_;
94  //int* defaultLabel_;
95 
96  std::vector<std::vector<LabelType> > permutation_; // org -> new
97  std::vector<std::vector<LabelType> > inversePermutation_; // new -> org
98  std::vector<std::vector<opengm::Tribool> > partialOptimality_;
99  std::vector<bool> optimal_;
100  std::vector<LabelType> label_;
101  std::vector<size_t> variableOffset_;
102 
103  size_t numNodes_;
104  size_t numEdges_;
105 
106  GraphValueType scale;
107 
108  };
110 
111 
112  template<class GM, class ACC>
113  MQPBO<GM,ACC>::MQPBO
114  (
115  const GmType& gm,
116  const Parameter& parameter
117  )
118  : gm_(gm),
119  param_(parameter),
120  scale(1)
121  {
122  for(size_t j = 0; j < gm_.numberOfFactors(); ++j) {
123  if(gm_[j].numberOfVariables() > 2) {
124  throw RuntimeError("This implementation of MQPBO supports only factors of order <= 2.");
125  }
126  }
127 
128  //Allocate Memory for Permutations
129  permutation_.resize(gm_.numberOfVariables());
130  inversePermutation_.resize(gm_.numberOfVariables());
131  for(IndexType var=0; var<gm_.numberOfVariables(); ++var){
132  permutation_[var].resize(gm_.numberOfLabels(var));
133  inversePermutation_[var].resize(gm_.numberOfLabels(var));
134  }
135 
136  //Set Default Optimality
137  partialOptimality_.resize(gm_.numberOfVariables());
138  optimal_.resize(gm_.numberOfVariables(),false);
139  label_.resize(gm_.numberOfVariables());
140  for(IndexType var=0; var<gm_.numberOfVariables(); ++var){
141  partialOptimality_[var].resize(gm_.numberOfLabels(var),opengm::Tribool::Maybe);
142  }
143 
144  //Calculated number of nodes and edges
145  numNodes_=0;
146  numEdges_=0;
147  size_t numSOF=0;
148  variableOffset_.resize(gm_.numberOfVariables(),0);
149  for(IndexType var=0; var<gm_.numberOfVariables(); ++var){
150  variableOffset_[var] = numNodes_;
151  numNodes_ += gm_.numberOfLabels(var)-1;
152  }
153  for(IndexType var=1; var<gm_.numberOfVariables(); ++var){
154  OPENGM_ASSERT( variableOffset_[var-1]< variableOffset_[var]);
155  }
156 
157  for(IndexType f=0; f<gm_.numberOfFactors(); ++f){
158  if(gm_[f].numberOfVariables()==1)
159  numEdges_ += gm_[f].numberOfLabels(0)-2;
160  if(gm_[f].numberOfVariables()==2){
161  numEdges_ += (gm_[f].numberOfLabels(0)-1);//*(gm_[f].numberOfLabels(1)-1);
162  ++numSOF;
163  }
164  }
165 
166  if(param_.rounds_>0){
167  std::cout << "Large" <<std::endl;
168  qpbo_ = new kolmogorov::qpbo::QPBO<GraphValueType > (numNodes_, numEdges_); // max number of nodes & edges
169  qpbo_->AddNode(numNodes_);
170  }
171  else{
172  std::cout << "Small" <<std::endl;
173  qpbo_ = new kolmogorov::qpbo::QPBO<GraphValueType > (gm_.numberOfVariables(), numSOF); // max number of nodes & edges
174  qpbo_->AddNode(gm_.numberOfVariables());
175  }
176  }
177 
178  template<class GM, class ACC>
179  MQPBO<GM,ACC>::~MQPBO
180  (
181  )
182  {
183  delete qpbo_;
184  }
185 
188  template<class GM, class ACC>
189  inline void
190  MQPBO<GM,ACC>::reset()
191  {
193  }
194 
196  template<class GM, class ACC>
197  inline void
198  MQPBO<GM,ACC>::setStartingPoint
199  (
200  typename std::vector<typename GM::LabelType>::const_iterator begin
201  )
202  {
204  }
205 
206  template<class GM, class ACC>
207  inline std::string
208  MQPBO<GM,ACC>::name() const
209  {
210  return "MQPBO";
211  }
212 
213  template<class GM, class ACC>
214  inline const typename MQPBO<GM,ACC>::GmType&
215  MQPBO<GM,ACC>::graphicalModel() const
216  {
217  return gm_;
218  }
219 
220 
221  template<class GM, class ACC>
222  inline void
223  MQPBO<GM,ACC>::AddUnaryTerm(int var, ValueType v0, ValueType v1){
224  qpbo_->AddUnaryTerm(var, scale*v0, scale*v1);
225  }
226 
227  template<class GM, class ACC>
228  inline void
229  MQPBO<GM,ACC>::AddPairwiseTerm(int var0, int var1,ValueType v00,ValueType v01,ValueType v10,ValueType v11){
230  qpbo_->AddPairwiseTerm(var0, var1,scale*v00,scale*v01,scale*v10,scale*v11);
231  }
232 
233  template<class GM, class ACC>
234  inline InferenceTermination
235  MQPBO<GM,ACC>::testQuess(LabelType guess)
236  {
237  qpbo_->Reset();
238  qpbo_->AddNode(gm_.numberOfVariables());
239  for(size_t f = 0; f < gm_.numberOfFactors(); ++f) {
240  if(gm_[f].numberOfVariables() == 0) {
241  ;
242  }
243  else if(gm_[f].numberOfVariables() == 1) {
244  const LabelType numLabels = gm_[f].numberOfLabels(0);
245  const IndexType var = gm_[f].variableIndex(0);
246 
247  ValueType v0 = gm_[f](&guess);
248  ValueType v1; ACC::neutral(v1);
249  for(LabelType i=0; i<guess; ++i)
250  ACC::op(gm_[f](&i),v1);
251  for(LabelType i=guess+1; i<numLabels; ++i)
252  ACC::op(gm_[f](&i),v1);
253  AddUnaryTerm(var, v0, v1);
254  }
255  else if(gm_[f].numberOfVariables() == 2) {
256  const IndexType var0 = gm_[f].variableIndex(0);
257  const IndexType var1 = gm_[f].variableIndex(1);
258 
259  LabelType c[2] = {guess,guess};
260  LabelType c2[2] = {0,1};
261 
262  ValueType v00 = gm_[f](c);
263  ValueType v01 = gm_[f](c2);
264  ValueType v10 = v01;
265  ValueType v11 = std::min(v00,v01);
266  AddPairwiseTerm(var0, var1,v00,v01,v10,v11);
267  }
268  }
269  qpbo_->MergeParallelEdges();
270  qpbo_->Solve();
271  for(IndexType var=0; var<gm_.numberOfVariables();++var){
272  if(qpbo_->GetLabel(var)==0){
273  for(LabelType l=0; l<gm_.numberOfLabels(var); ++l){
274  partialOptimality_[var][l] =opengm::Tribool::False;
275  }
276  partialOptimality_[var][guess] =opengm::Tribool::True;
277  optimal_[var]=true;
278  label_[var]=guess;
279  }
280  }
281  return NORMAL;
282  }
283 
284 
285  template<class GM, class ACC>
286  inline InferenceTermination
287  MQPBO<GM,ACC>::testQuess(std::vector<LabelType> &guess)
288  {
289  qpbo_->Reset();
290  qpbo_->AddNode(gm_.numberOfVariables());
291 
292  for(size_t var=0; var<gm_.numberOfVariables(); ++var){
293  std::vector<ValueType> v(gm_.numberOfLabels(var),0);
294  for(size_t i=0; i<gm_.numberOfFactors(var); ++i){
295  size_t f = gm_.factorOfVariable(var, i);
296  if(gm_[f].numberOfVariables()==1){
297  for(size_t j=0; j<v.size(); ++j){
298  v[j] += gm_[f](&j);
299  }
300 
301  }
302  else if(gm_[f].numberOfVariables() == 2) {
303  LabelType c[] = {guess[gm_[f].variableIndex(0)],guess[gm_[f].variableIndex(1)]};
304  if(gm_[f].variableIndex(0)==var){
305  for(c[0]=0; c[0]<guess[var]; ++c[0]){
306  v[c[0]] += gm_[f](c);
307  }
308  for(c[0]=guess[var]+1; c[0]<v.size(); ++c[0]){
309  v[c[0]] += gm_[f](c);
310  }
311  }
312  else if(gm_[f].variableIndex(1)==var){
313  for(c[1]=0; c[1]<guess[var]; ++c[1]){
314  v[c[1]] += gm_[f](c);
315  }
316  for(c[1]=guess[var]+1; c[1]<v.size(); ++c[1]){
317  v[c[1]] += gm_[f](c);
318  }
319  }
320  else{
321  OPENGM_ASSERT(false);
322  }
323  }
324  }
325  ValueType v0 = v[guess[var]];
326  ValueType v1; ACC::neutral(v1);
327  for(size_t j=0; j<guess[var]; ++j){
328  ACC::op(v[j],v1);
329  }
330  for(size_t j=guess[var]+1; j<v.size(); ++j){
331  ACC::op(v[j],v1);
332  }
333  AddUnaryTerm(var, v0, v1);
334  }
335 
336 
337  for(size_t f = 0; f < gm_.numberOfFactors(); ++f) {
338  if(gm_[f].numberOfVariables() < 2) {
339  continue;
340  }
341  else if(gm_[f].numberOfVariables() == 2) {
342  const IndexType var0 = gm_[f].variableIndex(0);
343  const IndexType var1 = gm_[f].variableIndex(1);
344 
345  LabelType c[2] = {guess[var0],guess[var1]};
346  LabelType c0[2] = {guess[var0],guess[var1]};
347  LabelType c1[2] = {guess[var0],guess[var1]};
348  ValueType v00 = gm_[f](c);
349  ValueType v01 = 0;
350  ValueType v10 = 0;
351  ValueType v11; ACC::neutral(v11);
352 
353  for(c[0]=0; c[0]<gm_[f].numberOfLabels(0); ++c[0]){
354  for(c[1]=0; c[1]<gm_[f].numberOfLabels(1); ++c[1]){
355  if(c[0]==guess[var0] || c[1]==guess[var1]){
356  continue;
357  }
358  else{
359  c0[0]=c[0];
360  c1[1]=c[1];
361  ValueType v = gm_[f](c) - gm_[f](c0) - gm_[f](c1);
362  ACC::op(v,v11);
363  }
364  }
365  }
366  AddPairwiseTerm(var0, var1,v00,v01,v10,v11);
367  }
368  }
369  qpbo_->MergeParallelEdges();
370  qpbo_->Solve();
371  for(IndexType var=0; var<gm_.numberOfVariables();++var){
372  if(qpbo_->GetLabel(var)==0){
373  for(LabelType l=0; l<gm_.numberOfLabels(var); ++l){
374  partialOptimality_[var][l] =opengm::Tribool::False;
375  }
376  partialOptimality_[var][guess[var]] =opengm::Tribool::True;
377  optimal_[var]=true;
378  label_[var]=guess[var];
379  }
380  }
381  return NORMAL;
382  }
383 
384  template<class GM, class ACC>
385  inline InferenceTermination
386  MQPBO<GM,ACC>::testPermutation(PermutationType permutationType)
387  {
388  //Set up MQPBO for current partial optimality
389  std::vector<IndexType> var2VarR(gm_.numberOfVariables());
390  std::vector<IndexType> varR2Var;
391  std::vector<size_t> varROffset;
392  size_t numBVar=0;
393  for(size_t var = 0; var < gm_.numberOfVariables(); ++var) {
394  if(optimal_[var]){
395  ;//do nothing
396  }
397  else{
398  varROffset.push_back(numBVar);
399  numBVar = numBVar + gm_.numberOfLabels(var)-1;
400  var2VarR[var]=varR2Var.size();
401  varR2Var.push_back(var);
402  }
403  }
404  std::cout << varR2Var.size() <<" / "<<gm_.numberOfVariables()<<std::endl;
405 
406  //Find Permutation
407  if(permutationType==NONE){
408  for(IndexType var=0; var<gm_.numberOfVariables(); ++var){
409  for(LabelType l=0; l<gm_.numberOfLabels(var); ++l){
410  permutation_[var][l]=l;
411  }
412  }
413  }
414  else if(permutationType==RANDOM){
415  srand ( unsigned ( time (NULL) ) );
416  for(IndexType var=0; var<gm_.numberOfVariables(); ++var){
417  LabelType numStates = gm_.numberOfLabels(var);
418  //IDENTYTY PERMUTATION
419  for(LabelType i=0; i<numStates;++i){
420  permutation_[var][i]=i;
421  }
422  //SHUFFEL PERMUTATION
423  std::random_shuffle(permutation_[var].begin(),permutation_[var].end());
424  }
425  }
426  else if(permutationType==MINMARG){
427 
428  typedef typename opengm::GraphicalModel<ValueType, OperatorType, opengm::ViewFixVariablesFunction<GM>, DiscreteSpace<IndexType, LabelType> > SUBGM;
429 
430 
431 
432  std::vector<LabelType> numberOfLabels(varR2Var.size());
433  for(size_t i=0; i<varR2Var.size(); ++i)
434  numberOfLabels[i] = gm_.numberOfLabels(varR2Var[i]);
435  typename SUBGM::SpaceType subspace(numberOfLabels.begin(),numberOfLabels.end());
436  SUBGM gm(subspace);
437  for(IndexType f=0; f<gm_.numberOfFactors();++f){
438  std::vector<PositionAndLabel<IndexType, LabelType> > fixed;
439  std::vector<IndexType> vars;
440  for(IndexType i=0; i<gm_[f].numberOfVariables();++i){
441  const IndexType var = gm_[f].variableIndex(i);
442  if(optimal_[var]){
443  fixed.push_back(PositionAndLabel<IndexType, LabelType>(i,label_[var]));
444  }
445  else{
446  vars.push_back(var2VarR[var]);
447  }
448  }
449  opengm::ViewFixVariablesFunction<GM> func(gm_[f], fixed);
450  gm.addFactor(gm.addFunction(func),vars.begin(),vars.end());
451  }
452 
453  typedef typename opengm::MessagePassing<SUBGM, ACC,opengm::BeliefPropagationUpdateRules<SUBGM,ACC>, opengm::MaxDistance> LBP;
454  typename LBP::Parameter para;
455  para.maximumNumberOfSteps_ = 100;
456  para.damping_ = 0.5;
457  LBP bp(gm,para);
458  bp.infer();
459 
460  //for(IndexType var=0; var<gm_.numberOfVariables(); ++var){
461  for(IndexType varR=0; varR<gm.numberOfVariables(); ++varR){
462  IndexType var = varR2Var[varR];
463  LabelType numStates = gm_.numberOfLabels(var);
464  typename GM::IndependentFactorType marg;
465  bp.marginal(varR, marg);
466 
467  //SHUFFEL PERMUTATION
468  std::vector<LabelType> list(numStates);
469  for(LabelType i=0; i<numStates;++i){
470  list[i]=i;
471  }
472  LabelType t;
473  for(LabelType i=0; i<numStates;++i){
474  for(LabelType j=i+1; i<numStates;++i){
475  if(marg(&list[j])<marg(&list[i])){
476  t = list[i];
477  list[i]=list[j];
478  list[j]=t;
479  }
480  }
481  }
482  for(LabelType i=0; i<numStates;++i){
483  permutation_[var][i] = list[i];
484  }
485  }
486  }
487  else{
488  throw RuntimeError("Error: Unknown Permutation!");
489  }
490  //CALCULATE INVERSE PERMUTATION
491  for(IndexType var=0; var<gm_.numberOfVariables(); ++var){
492  for(LabelType l=0; l<gm_.numberOfLabels(var); ++l){
493  inversePermutation_[var][permutation_[var][l]]=l;
494  }
495  }
496 
497 
498  //Build Graph
499  ValueType constValue = 0;
500  qpbo_->Reset();
501  qpbo_->AddNode(numBVar);
502  //qpbo_->AddNode(numNodes_);
503 
504  for(IndexType varR = 0; varR < varR2Var.size(); ++varR) {
505  IndexType var = varR2Var[varR];
506  for(size_t l = 0; l+1<gm_.numberOfLabels(var); ++l){
507  AddUnaryTerm((int) (varROffset[varR]+l), 0.0, 0.0);
508  }
509  for(LabelType l=1; l+1<gm_.numberOfLabels(var); ++l){
510  AddPairwiseTerm((int) (varROffset[varR]+l-1), (int) (varROffset[varR]+l), 0.0, 1e30, 0.0, 0.0);
511  }
512  }
513  /*
514  for(size_t var = 0; var < gm_.numberOfVariables(); ++var) {
515  for(size_t l = 0; l+1<gm_.numberOfLabels(var); ++l){
516  AddUnaryTerm((int) (variableOffset_[var]+l), 0.0, 0.0);
517  }
518  for(LabelType l=1; l+1<gm_.numberOfLabels(var); ++l){
519  AddPairwiseTerm((int) (variableOffset_[var]+l-1), (int) (variableOffset_[var]+l), 0.0, 1000000.0, 0.0, 0.0);
520  }
521  }
522  */
523 
524  for(size_t f = 0; f < gm_.numberOfFactors(); ++f) {
525  if(gm_[f].numberOfVariables() == 0) {
526  const LabelType l = 0;
527  constValue += gm_[f](&l);
528  }
529  else if(gm_[f].numberOfVariables() == 1) {
530  const LabelType numLabels = gm_[f].numberOfLabels(0);
531  const IndexType var = gm_[f].variableIndex(0);
532  if(optimal_[var]){
533  constValue += gm_[f](&(label_[var]));
534  }
535  else{
536  LabelType l0 = inversePermutation_[var][0];
537  LabelType l1;
538  constValue += gm_[f](&l0);
539  const IndexType varR = var2VarR[var];
540  for(LabelType i=1 ; i<numLabels; ++i){
541  l0 = inversePermutation_[var][i-1];
542  l1 = inversePermutation_[var][i];
543  AddUnaryTerm((int) (varROffset[varR]+i-1), 0.0, gm_[f](&l1)-gm_[f](&l0));
544  //AddUnaryTerm((int) (variableOffset_[var]+i-1), 0.0, gm_[f](&l1)-gm_[f](&l0));
545  }
546  }
547  }
548  else if(gm_[f].numberOfVariables() == 2) {
549  const IndexType var0 = gm_[f].variableIndex(0);
550  const IndexType var1 = gm_[f].variableIndex(1);
551  const IndexType varR0 = var2VarR[var0];
552  const IndexType varR1 = var2VarR[var1];
553 
554  if(optimal_[var0]&&optimal_[var1]){
555  LabelType l[2] = { label_[var0], label_[var1]};
556  constValue += gm_[f](l);
557  }
558  else if(optimal_[var0]){
559  const LabelType numLabels = gm_[f].numberOfLabels(1);
560  LabelType l0[2] = { label_[var0], inversePermutation_[var1][0]};
561  LabelType l1[2] = { label_[var0], 0};
562  constValue += gm_[f](l0);
563  for(LabelType i=1 ; i<numLabels; ++i){
564  l0[1] = inversePermutation_[var1][i-1];
565  l1[1] = inversePermutation_[var1][i];
566  //AddUnaryTerm((int) (variableOffset_[var1]+i-1), 0.0, gm_[f](l1)-gm_[f](l0));
567  AddUnaryTerm((int) (varROffset[varR1]+i-1), 0.0, gm_[f](l1)-gm_[f](l0));
568  }
569  }
570  else if(optimal_[var1]){
571  const LabelType numLabels = gm_[f].numberOfLabels(0);
572  LabelType l0[2] = { inversePermutation_[var0][0], label_[var1]};
573  LabelType l1[2] = { 0, label_[var1]};
574  constValue += gm_[f](l0);
575  for(LabelType i=1 ; i<numLabels; ++i){
576  l0[0] = inversePermutation_[var0][i-1];
577  l1[0] = inversePermutation_[var0][i];
578  AddUnaryTerm((int) (varROffset[varR0]+i-1), 0.0, gm_[f](l1)-gm_[f](l0));
579  //AddUnaryTerm((int) (variableOffset_[var0]+i-1), 0.0, gm_[f](l1)-gm_[f](l0));
580  }
581  }
582  else{
583  {
584  const LabelType l[2]={inversePermutation_[var0][0],inversePermutation_[var1][0]};
585  constValue += gm_[f](l);
586  }
587  for(size_t i=1; i<gm_[f].numberOfLabels(0);++i){
588  const LabelType l1[2]={inversePermutation_[var0][i] ,inversePermutation_[var1][0]};
589  const LabelType l2[2]={inversePermutation_[var0][i-1],inversePermutation_[var1][0]};
590  AddUnaryTerm((int) (varROffset[varR0]+i-1), 0.0, gm_[f](l1)-gm_[f](l2));
591  //AddUnaryTerm((int) (variableOffset_[var0]+i-1), 0.0, gm_[f](l1)-gm_[f](l2));
592  }
593  for(size_t i=1; i<gm_[f].numberOfLabels(1);++i){
594  const LabelType l1[2]={inversePermutation_[var0][0],inversePermutation_[var1][i]};
595  const LabelType l2[2]={inversePermutation_[var0][0],inversePermutation_[var1][i-1]};
596  AddUnaryTerm((int) (varROffset[varR1]+i-1), 0.0, gm_[f](l1)-gm_[f](l2));
597  //AddUnaryTerm((int) (variableOffset_[var1]+i-1), 0.0, gm_[f](l1)-gm_[f](l2));
598  }
599  for(size_t i=1; i<gm_[f].numberOfLabels(0);++i){
600  for(size_t j=1; j<gm_[f].numberOfLabels(1);++j){
601  const int node0 = varROffset[varR0]+i-1;
602  const int node1 = varROffset[varR1]+j-1;
603  //const int node0 = variableOffset_[var0]+i-1;
604  //const int node1 = variableOffset_[var1]+j-1;
605  ValueType v = 0;
606  int l[2] = {(int)inversePermutation_[var0][i],(int)inversePermutation_[var1][j]}; v += gm_[f](l);
607  l[0]=inversePermutation_[var0][i-1]; v -= gm_[f](l);
608  l[1]=inversePermutation_[var1][j-1]; v += gm_[f](l);
609  l[0]=inversePermutation_[var0][i]; v -= gm_[f](l);
610  if(v!=0.0)
611  AddPairwiseTerm(node0, node1,0.0,0.0,0.0,v);
612  }
613  }
614  }
615  }
616  }
617  qpbo_->MergeParallelEdges();
618 
619  //Optimize
620 
621  qpbo_->Solve();
622  if(!param_.strongPersistency_)
623  qpbo_->ComputeWeakPersistencies();
624  // if(!parameter_.strongPersistency_) {
625  // qpbo_->ComputeWeakPersistencies();
626  // }
627 
628  bound_ = constValue + 0.5 * qpbo_->ComputeTwiceLowerBound();
629 
630  /*PROBEING*/
631  if(param_.probing_) {
632  std::cout << "Start Probing ..."<<std::endl;
633  // Initialize mapping for probe
634  int *mapping = new int[numBVar];
635  //int *mapping = new int[numNodes_];
636  for(int i = 0; i < static_cast<int>(numBVar); ++i) {
637  //for(int i = 0; i < static_cast<int>(numNodes_); ++i) {
638  qpbo_->SetLabel(i, qpbo_->GetLabel(i));
639  mapping[i] = i * 2;
640  }
641  typename kolmogorov::qpbo::QPBO<GraphValueType>::ProbeOptions options;
642  options.C = 1000000000;
643  if(!param_.strongPersistency_)
644  options.weak_persistencies = 1;
645  else
646  options.weak_persistencies = 0;
647  qpbo_->Probe(mapping, options);
648  if(!param_.strongPersistency_)
649  qpbo_->ComputeWeakPersistencies();
650 
651  for(IndexType var=0; var<gm_.numberOfVariables();++var){
652  if(optimal_[var]) continue;
653  IndexType varR = var2VarR[var];
654  //Lable==0
655  {
656  int l = qpbo_->GetLabel(mapping[varROffset[varR]]/2);
657  if(l>=0) l = (l + mapping[varROffset[varR]]) % 2;
658  //int l = qpbo_->GetLabel(mapping[variableOffset_[var]]/2);
659  //if(l>=0) l = (l + mapping[variableOffset_[var]]) % 2;
660  if(l==0) {partialOptimality_[var][inversePermutation_[var][0]]&=opengm::Tribool::True;}
661  else if(l==1){partialOptimality_[var][inversePermutation_[var][0]]&=opengm::Tribool::False;}
662  else {partialOptimality_[var][inversePermutation_[var][0]]&=opengm::Tribool::Maybe;}
663  }
664  //Label==max
665  {
666  int l = qpbo_->GetLabel(mapping[varROffset[varR]+gm_.numberOfLabels(var)-2]/2);
667  if(l>=0) l = (l + mapping[varROffset[varR]+gm_.numberOfLabels(var)-2]) % 2;
668  //int l = qpbo_->GetLabel(mapping[variableOffset_[var]+gm_.numberOfLabels(var)-2]/2);
669  //if(l>=0) l = (l + mapping[variableOffset_[var]+gm_.numberOfLabels(var)-2]) % 2;
670  if(l==0) {partialOptimality_[var][inversePermutation_[var][gm_.numberOfLabels(var)-1]]&=opengm::Tribool::False;}
671  else if(l==1){partialOptimality_[var][inversePermutation_[var][gm_.numberOfLabels(var)-1]]&=opengm::Tribool::True;}
672  else {partialOptimality_[var][inversePermutation_[var][gm_.numberOfLabels(var)-1]]&=opengm::Tribool::Maybe;}
673  }
674  //ELSE
675 
676  for(LabelType l=1; l+1<gm_.numberOfLabels(var);++l)
677  {
678  int l1 = qpbo_->GetLabel(mapping[varROffset[varR]+l-1]/2);
679  int l2 = qpbo_->GetLabel(mapping[varROffset[varR]+l]/2);
680  if(l1>=0) l1 = (l1 + mapping[varROffset[varR]+l-1]) % 2;
681  if(l2>=0) l2 = (l2 + mapping[varROffset[varR]+l]) % 2;
682  //int l1 = qpbo_->GetLabel(mapping[variableOffset_[var]+l-1]/2);
683  //int l2 = qpbo_->GetLabel(mapping[variableOffset_[var]+l]/2);
684  //if(l1>=0) l1 = (l1 + mapping[variableOffset_[var]+l-1]) % 2;
685  //if(l2>=0) l2 = (l2 + mapping[variableOffset_[var]+l]) % 2;
686 
687  OPENGM_ASSERT(!(l1==0 && l2==1));
688  if(l1==1 && l2==0) {partialOptimality_[var][inversePermutation_[var][l]]&=opengm::Tribool::True;}
689  else if(l2==1) {partialOptimality_[var][inversePermutation_[var][l]]&=opengm::Tribool::False;}
690  else if(l1==0) {partialOptimality_[var][inversePermutation_[var][l]]&=opengm::Tribool::False;}
691  //else {partialOptimality_[var][inversePermutation_[var][l]]&=opengm::Tribool::Maybe;}
692  }
693  }
694  delete mapping;
695  }
696  else{
697  for(IndexType var=0; var<gm_.numberOfVariables();++var){
698  if(optimal_[var]) continue;
699  IndexType varR = var2VarR[var];
700  //Lable==0
701  {
702  int l = qpbo_->GetLabel(varROffset[varR]);
703  //int l = qpbo_->GetLabel(variableOffset_[var]);
704  if(l==0){
705  OPENGM_ASSERT(!(partialOptimality_[var][inversePermutation_[var][0]]==opengm::Tribool::False));
706  partialOptimality_[var][inversePermutation_[var][0]]&=opengm::Tribool::True;
707  }
708  else if(l==1){
709  OPENGM_ASSERT(!(partialOptimality_[var][inversePermutation_[var][0]]==opengm::Tribool::True));
710  partialOptimality_[var][inversePermutation_[var][0]]&=opengm::Tribool::False;
711  }
712  // else {partialOptimality_[var][permutation_[var][0]]&=opengm::Tribool::Maybe;}
713  }
714  //Label==max
715  {
716  int l = qpbo_->GetLabel(varROffset[varR]+gm_.numberOfLabels(var)-2);
717  //int l = qpbo_->GetLabel(variableOffset_[var]+gm_.numberOfLabels(var)-2);
718  if(l==0){
719  OPENGM_ASSERT(!(partialOptimality_[var][inversePermutation_[var][gm_.numberOfLabels(var)-1]]==opengm::Tribool::True));
720  partialOptimality_[var][inversePermutation_[var][gm_.numberOfLabels(var)-1]]&=opengm::Tribool::False;
721  }
722  else if(l==1){
723  OPENGM_ASSERT(!(partialOptimality_[var][inversePermutation_[var][gm_.numberOfLabels(var)-1]]==opengm::Tribool::False));
724  partialOptimality_[var][inversePermutation_[var][gm_.numberOfLabels(var)-1]]&=opengm::Tribool::True;
725  }
726  //else {partialOptimality_[var][permutation_[var][gm_.numberOfLabels(var)-1]]&=opengm::Tribool::Maybe;}
727  }
728  //ELSE
729 
730  for(LabelType l=1; l+1<gm_.numberOfLabels(var);++l)
731  {
732  int l1 = qpbo_->GetLabel(varROffset[varR]+l-1);
733  int l2 = qpbo_->GetLabel(varROffset[varR]+l);
734  //int l1 = qpbo_->GetLabel(variableOffset_[var]+l-1);
735  //int l2 = qpbo_->GetLabel(variableOffset_[var]+l);
736  OPENGM_ASSERT(!(l1==0 && l2==1));
737  if(l1==1 && l2==0) {
738  OPENGM_ASSERT(!(partialOptimality_[var][inversePermutation_[var][l]]==opengm::Tribool::False));
739  partialOptimality_[var][inversePermutation_[var][l]]&=opengm::Tribool::True;
740  }
741  else if(l2==1){
742  OPENGM_ASSERT(!(partialOptimality_[var][inversePermutation_[var][l]]==opengm::Tribool::True));
743  partialOptimality_[var][inversePermutation_[var][l]]&=opengm::Tribool::False;
744  }
745  else if(l1==0){
746  OPENGM_ASSERT(!(partialOptimality_[var][inversePermutation_[var][l]]==opengm::Tribool::True));
747  partialOptimality_[var][inversePermutation_[var][l]]&=opengm::Tribool::False;
748  }
749  //else{
750  // partialOptimality_[var][permutation_[var][l]]&=opengm::Tribool::Maybe;
751  //}
752  }
753  }
754  }
755  for(IndexType var=0; var<gm_.numberOfVariables();++var){
756  if(optimal_[var]) continue;
757  LabelType countTRUE = 0;
758  LabelType countFALSE = 0;
759  for(LabelType l=1; l+1<gm_.numberOfLabels(var);++l){
760  if(partialOptimality_[var][l]==opengm::Tribool::True)
761  ++countTRUE;
762  if(partialOptimality_[var][l]==opengm::Tribool::False)
763  ++countFALSE;
764  }
765  if(countTRUE==1){
766  optimal_[var]=true;
767  for(LabelType l=1; l+1<gm_.numberOfLabels(var);++l){
768  if(partialOptimality_[var][l]==opengm::Tribool::True)
769  label_[var]=l;
770  else
771  partialOptimality_[var][l]=opengm::Tribool::False;
772  }
773  }
774  if(countFALSE+1==gm_.numberOfLabels(var)){
775  optimal_[var]=true;
776  for(LabelType l=1; l+1<gm_.numberOfLabels(var);++l){
777  if(partialOptimality_[var][l]==opengm::Tribool::Maybe){
778  label_[var]=l;
779  partialOptimality_[var][l]=opengm::Tribool::True;
780  }
781  }
782  }
783  }
784  return NORMAL;
785  }
786 
787  template<class GM, class ACC>
788  InferenceTermination MQPBO<GM,ACC>::infer
789  ()
790  {
791  EmptyVisitorType visitor;
792  return infer(visitor);
793  }
794 
795  template<class GM, class ACC>
796  template<class VisitorType>
797  InferenceTermination MQPBO<GM,ACC>::infer
798  (
799  VisitorType& visitor
800  )
801  {
802  visitor.addLog("optimality");
803  visitor.addLog("optimalityV");
804  if(param_.rounds_>1 && param_.strongPersistency_==false)
805  std::cout << "WARNING: Using weak persistency and several rounds may lead to wrong results if solution is not unique!"<<std::endl;
806 
807  LabelType maxNumberOfLabels = 0;
808  for(IndexType var=0; var<gm_.numberOfVariables();++var){
809  maxNumberOfLabels = std::max(maxNumberOfLabels, gm_.numberOfLabels(var));
810  }
811  bool isPotts = true;
812 
813  for(IndexType f=0; f< gm_.numberOfFactors(); ++f){
814  if(gm_[f].numberOfVariables()<2) continue;
815  isPotts &= gm_[f].isPotts();
816  if(!isPotts) break;
817  }
818 
819  visitor.begin(*this);
820 
821  if(param_.useKovtunsMethod_){
822  if(isPotts){
823  std::cout << "Use Kovtuns method for potts"<<std::endl;
824  for(LabelType l=0; l<maxNumberOfLabels; ++l) {
825  testQuess(l);
826  double xoptimality = optimality();
827  double xoptimalityV = optimalityV();
828  visitor(*this);
829  visitor.log("optimality",xoptimality);
830  visitor.log("optimalityV",xoptimalityV);
831 
832  //std::cout << "partialOptimality : " << optimality() << std::endl;
833  }
834  }
835  else{
836  std::cout << "Use Kovtuns method for non-potts is not supported yet"<<std::endl;
837  /*
838  for(LabelType l=0; l<maxNumberOfLabels; ++l){
839  std::vector<LabelType> guess(gm_.numberOfVariables(),l);
840  for(IndexType var=0; var<gm_.numberOfVariables();++var){
841  if(l>=gm_.numberOfLabels(var)){
842  guess[var]=l-1;
843  }
844  }
845  testQuess(guess);
846  double xoptimality = optimality();
847  visitor(*this,this->value(),bound(),"partialOptimality",xoptimality);
848  //std::cout << "partialOptimality : " << optimality() << std::endl;
849  }
850  */
851  }
852  }
853 
854  if(param_.rounds_>0){
855  std::cout << "Start "<<param_.rounds_ << " of multilabel QPBO for different permutations" <<std::endl;
856  for(size_t rr=0; rr<param_.rounds_;++rr){
857  testPermutation(param_.permutationType_);
858  double xoptimality = optimality();
859  double xoptimalityV = optimalityV();
860  visitor(*this);
861  visitor.log("optimality",xoptimality);
862  visitor.log("optimalityV",xoptimalityV);
863 
864  //std::cout << "partialOptimality : " << optimality() << std::endl;
865  }
866  }
867 
868 #ifdef MQPBOHotFixOutPutPartialOPtimalityMap
869  hid_t fid = marray::hdf5::createFile("mqpbotmp.h5");
870  std::vector<double> optimal;
871  for(size_t i=0; i<optimal_.size();++i)
872  optimal.push_back((double)(optimal_[i]));
873  marray::hdf5::save(fid, "popt", optimal);
874  marray::hdf5::closeFile(fid);
875 #endif
876 
877  visitor.end(*this);
878 
879  return NORMAL;
880  }
881 
882  template<class GM, class ACC>
883  double
884  MQPBO<GM,ACC>::optimality
885  () const
886  {
887  size_t labeled = 0;
888  size_t unlabeled = 0;
889  for(IndexType var=0; var<gm_.numberOfVariables();++var){
890  for(LabelType l=0; l<gm_.numberOfLabels(var);++l){
891  if(partialOptimality_[var][l]==opengm::Tribool::Maybe)
892  ++unlabeled;
893  else
894  ++labeled;
895  }
896  }
897  return labeled*1.0/(labeled+unlabeled);
898  }
899 
900  template<class GM, class ACC>
901  double
902  MQPBO<GM,ACC>::optimalityV
903  () const
904  {
905  size_t labeled = 0;
906  for(IndexType var=0; var<gm_.numberOfVariables();++var){
907  for(LabelType l=0; l<gm_.numberOfLabels(var);++l){
908  if(partialOptimality_[var][l]==opengm::Tribool::True){
909  ++labeled;
910  continue;
911  }
912  }
913  }
914  return labeled*1.0/gm_.numberOfVariables();
915  }
916 
917  template<class GM, class ACC>
918  typename GM::ValueType
919  MQPBO<GM,ACC>::bound
920  () const
921  {
922  return bound_;
923  }
924 
925  template<class GM, class ACC>
926  typename GM::ValueType MQPBO<GM,ACC>::value() const {
927  std::vector<LabelType> states;
928  arg(states);
929  return gm_.evaluate(states);
930  }
931 
932  template<class GM, class ACC>
933  inline InferenceTermination
934  MQPBO<GM,ACC>::arg
935  (
936  std::vector<LabelType>& x,
937  const size_t N
938  ) const
939  {
940  if(N==1){
941  x.resize(gm_.numberOfVariables(),0);
942 
943  for(IndexType var=0; var<gm_.numberOfVariables(); ++var){
944  size_t countTrue = 0;
945  size_t countFalse = 0;
946  size_t countMaybe = 0;
947  x[var]=0;
948  for(LabelType l=0; l<gm_.numberOfLabels(var);++l){
949  if(partialOptimality_[var][l]==opengm::Tribool::Maybe){
950  x[var] = l;
951  ++countMaybe;
952  }
953  if(partialOptimality_[var][l]==opengm::Tribool::True){
954  x[var] = l;
955  ++countTrue;
956  }
957  if(partialOptimality_[var][l]==opengm::Tribool::False){
958  ++countFalse;
959  }
960  }
961  OPENGM_ASSERT(countTrue+countFalse+countMaybe == gm_.numberOfLabels(var));
962  OPENGM_ASSERT(countTrue<2);
963  OPENGM_ASSERT(countFalse<gm_.numberOfLabels(var));
964  }
965  return NORMAL;
966  }
967  else {
968  return UNKNOWN;
969  }
970  }
971 } // namespace opengm
972 
973 #endif // #ifndef OPENGM_MQPBO_HXX
opengm::MQPBO::AccumulationType
ACC AccumulationType
Definition: mqpbo.hxx:40
opengm::MQPBO::setStartingPoint
void setStartingPoint(typename std::vector< typename GM::LabelType >::const_iterator)
set starting point
Definition: mqpbo.hxx:199
opengm
The OpenGM namespace.
Definition: config.hxx:43
marray::hdf5::createFile
hid_t createFile(const std::string &, HDF5Version=DEFAULT_HDF5_VERSION)
Create an HDF5 file.
Definition: marray_hdf5.hxx:747
opengm::DiscreteSpace
Discrete space in which variables can have differently many labels.
Definition: discretespace.hxx:17
opengm::MQPBO::value
GM::ValueType value() const
return the solution (value)
Definition: mqpbo.hxx:926
opengm::MQPBO::name
std::string name() const
Definition: mqpbo.hxx:208
opengm::MQPBO::reset
void reset()
reset assumes that the structure of the graphical model has not changed
Definition: mqpbo.hxx:190
opengm::MessagePassing
A framework for message passing algorithms Cf. F. R. Kschischang, B. J. Frey and H...
Definition: messagepassing.hxx:50
opengm::MQPBO::MINMARG
Definition: mqpbo.hxx:49
opengm::MQPBO::NONE
Definition: mqpbo.hxx:49
opengm::MQPBO::testPermutation
InferenceTermination testPermutation(PermutationType permutationType)
Definition: mqpbo.hxx:386
opengm::MQPBO
[class mqpbo] Multilabel QPBO (MQPBO) Implements the algorithms described in i) Ivan Kovtun: Partial ...
Definition: mqpbo.hxx:37
OPENGM_ASSERT
#define OPENGM_ASSERT(expression)
Definition: opengm.hxx:77
opengm::MQPBO::GraphValueType
ValueType GraphValueType
Definition: mqpbo.hxx:47
opengm::MQPBO::OPENGM_GM_TYPE_TYPEDEFS
OPENGM_GM_TYPE_TYPEDEFS
Definition: mqpbo.hxx:43
view_fix_variables_function.hxx
opengm::Inference::IndexType
GraphicalModelType::IndexType IndexType
Definition: inference.hxx:40
opengm::MQPBO::MQPBO
MQPBO(const GmType &, const Parameter &=Parameter())
[class mqpbo]
Definition: mqpbo.hxx:114
opengm::MQPBO::arg
virtual InferenceTermination arg(std::vector< LabelType > &, const size_t=1) const
output a solution
Definition: mqpbo.hxx:935
opengm::MQPBO::optimality
double optimality() const
Definition: mqpbo.hxx:885
opengm::Inference::ValueType
GraphicalModelType::ValueType ValueType
Definition: inference.hxx:41
opengm::Inference
Inference algorithm interface.
Definition: inference.hxx:34
opengm::MQPBO::VerboseVisitorType
visitors::VerboseVisitor< MQPBO< GM, ACC > > VerboseVisitorType
Definition: mqpbo.hxx:44
marray::hdf5::closeFile
void closeFile(const hid_t &)
Close an HDF5 file.
Definition: marray_hdf5.hxx:810
opengm::MQPBO::optimalityV
double optimalityV() const
Definition: mqpbo.hxx:903
marray::hdf5::save
void save(const hid_t &, const std::string &, const Marray< T > &)
Save an Marray as an HDF5 dataset.
Definition: marray_hdf5.hxx:265
opengm::MQPBO::testQuess
InferenceTermination testQuess(std::vector< LabelType > &guess)
Definition: mqpbo.hxx:287
opengm::MQPBO::partialOptimality
bool partialOptimality(IndexType var, LabelType &l) const
Definition: mqpbo.hxx:78
opengm::MQPBO::Parameter::probing_
const bool probing_
Definition: mqpbo.hxx:56
opengm::MQPBO::bound
GM::ValueType bound() const
return a bound on the solution
Definition: mqpbo.hxx:920
opengm::MQPBO::PermutationType
PermutationType
Definition: mqpbo.hxx:49
opengm::ViewFixVariablesFunction
Funcion that refers to a factor of another GraphicalModel in which some variables are fixed...
Definition: view_fix_variables_function.hxx:22
marray_hdf5.hxx
opengm::MQPBO::Parameter::label_
std::vector< LabelType > label_
Definition: mqpbo.hxx:53
opengm::DiscreteSpace::numberOfLabels
LabelType numberOfLabels(const IndexType) const
Definition: discretespace.hxx:124
opengm::Inference::LabelType
GraphicalModelType::LabelType LabelType
Definition: inference.hxx:39
opengm::MQPBO::TimingVisitorType
visitors::TimingVisitor< MQPBO< GM, ACC > > TimingVisitorType
Definition: mqpbo.hxx:46
opengm::MQPBO::RANDOM
Definition: mqpbo.hxx:49
opengm::MQPBO::EmptyVisitorType
visitors::EmptyVisitor< MQPBO< GM, ACC > > EmptyVisitorType
Definition: mqpbo.hxx:45
opengm::MQPBO::GraphicalModelType
GM GraphicalModelType
Definition: mqpbo.hxx:42
opengm::MQPBO::Parameter::permutationType_
PermutationType permutationType_
Definition: mqpbo.hxx:59
opengm::RuntimeError
OpenGM runtime error.
Definition: opengm.hxx:100
opengm::MQPBO::infer
InferenceTermination infer()
Definition: mqpbo.hxx:789
opengm::MQPBO::graphicalModel
const GmType & graphicalModel() const
Definition: mqpbo.hxx:215
opengm::InferenceTermination
InferenceTermination
Definition: inference.hxx:24
opengm::MQPBO::partialOptimality
const std::vector< opengm::Tribool > & partialOptimality(IndexType var) const
Definition: mqpbo.hxx:77
Generated by   doxygen 1.8.9.1