OpenGM  2.3.x
Discrete Graphical Model Library
lsatr.hxx
Go to the documentation of this file.
1 #pragma once
2 #ifndef OPENGM_LSATR_HXX
3 #define OPENGM_LSATR_HXX
4 
5 #include <algorithm>
6 #include <vector>
7 #include <string>
8 #include <iostream>
9 
10 #include "opengm/opengm.hxx"
11 #include "opengm/inference/inference.hxx"
12 #include "opengm/inference/visitors/visitors.hxx"
13 
14 #include <maxflowlib.h>
15 
16 #ifndef NOVIGRA
17 #include "vigra/multi_distance.hxx"
18 #include "vigra/multi_array.hxx"
19 #endif
20 
21 namespace opengm {
22 
36 
37 
38 
39  struct LSA_TR_WeightedEdge{
40  LSA_TR_WeightedEdge(double aw, size_t au, size_t av): w(aw), u(au), v(av){}
41  double w;
42  size_t u;
43  size_t v;
44  };
45 
46 
47  template<class LabelType>
48  class LSA_TR_HELPER{
49  public:
50  enum DISTANCE {HAMMING, EUCLIDEAN};
51 
52  LSA_TR_HELPER() { distanceType_= EUCLIDEAN;};
53  ~LSA_TR_HELPER(){ if(graph_!=NULL){delete graph_; delete changedList_;} };
54  template<class GM>
55  void init(const GM&, const std::vector<LabelType>& );
56  void set(const double);
57  void set(const std::vector<LabelType>&, const double);
58  double optimize(std::vector<LabelType>&);
59  void setDistanceType(const DISTANCE d){ distanceType_=d; };
60 
61  double eval(const std::vector<LabelType>&) const;
62  double evalAprox(const std::vector<LabelType>&,const std::vector<LabelType>&, const double) const;
63  void evalBoth(const std::vector<LabelType>& label, const std::vector<LabelType>& workingPoint, const double lambda, double& value, double& valueAprox) const;
64 
65  private:
66  typedef maxflowLib::Graph<double,double,double> graph_type;
67  typedef maxflowLib::Block<typename graph_type::node_id> block_type;
68 
69  void updateDistance();
70 
71  size_t numVar_;
72  double lambda_;
73  double constTerm_;
74  double constTermApproximation_;
75  double constTermTrustRegion_;
76  std::vector<LabelType> workingPoint_;
77  std::vector<double> distance_;
78  std::vector<double> unaries_;
79  std::vector<double> approxUnaries_;
80  std::vector< LSA_TR_WeightedEdge> supEdges_;
81  std::vector< LSA_TR_WeightedEdge> subEdges_;
82  graph_type* graph_;
83  block_type* changedList_;
84  bool solved_;
85  DISTANCE distanceType_;
86  std::vector<size_t> shape_;
87  };
88 
89 
90  template<class GM, class ACC>
91  class LSA_TR : public Inference<GM, ACC>
92  {
93  public:
94  typedef ACC AccumulationType;
95  typedef GM GraphicalModelType;
96  OPENGM_GM_TYPE_TYPEDEFS;
97  typedef opengm::visitors::VerboseVisitor<LSA_TR<GM,ACC> > VerboseVisitorType;
98  typedef opengm::visitors::EmptyVisitor<LSA_TR<GM,ACC> > EmptyVisitorType;
99  typedef opengm::visitors::TimingVisitor<LSA_TR<GM,ACC> > TimingVisitorType;
100 
101  class Parameter {
102  public:
103  enum DISTANCE {HAMMING, EUCLIDEAN};
104  size_t randSeed_;
105  double maxLambda_;
106  double initialLambda_;
107  double precisionLambda_;
108  double lambdaMultiplier_;
109  double reductionRatio_;
110  DISTANCE distance_;
111 
112  Parameter(){
113  randSeed_ = 42;
114  maxLambda_ = 1e5;
115  initialLambda_ = 0.1;
116  precisionLambda_ = 1e-9; // used to compare GEO lambda in parametric maxflow
117  lambdaMultiplier_ = 2; // used for jumps in backtracking;
118  reductionRatio_ = 0.25; // used to decide whether to increase or decrease lambda using the multiplier
119  distance_ = EUCLIDEAN;
120  }
121  };
122 
123  LSA_TR(const GraphicalModelType&);
124  LSA_TR(const GraphicalModelType&, const Parameter&);
125  ~LSA_TR();
126  std::string name() const;
127  const GraphicalModelType& graphicalModel() const;
128  InferenceTermination infer();
129  void reset();
130  template<class VisitorType>
131  InferenceTermination infer(VisitorType&);
132  void setStartingPoint(typename std::vector<LabelType>::const_iterator);
133  virtual InferenceTermination arg(std::vector<LabelType>&, const size_t = 1) const ;
134  virtual ValueType value()const{ return gm_.evaluate(curState_);}
135 
136  private:
137  void init();
138  double findMinimalChangeBrakPoint(const double lambda, const std::vector<LabelType>& workingPoint);
139 
140  LSA_TR_HELPER<LabelType> helper_;
141  const GraphicalModelType& gm_;
142  Parameter param_;
143  std::vector<LabelType> curState_;
144  size_t numVar_;
145 
146  ValueType constTerm_;
147  std::vector<ValueType> unaries_;
148  std::vector<LSA_TR_WeightedEdge> subEdges_;
149  std::vector<LSA_TR_WeightedEdge> supEdges_;
150  std::vector<ValueType> approxUnaries_;
151 
152  };
153 
155  template<class LabelType>
156  template<class GM>
157  void LSA_TR_HELPER<LabelType>::init(const GM& gm, const std::vector<LabelType>& workingPoint){
158  typedef size_t IndexType;
159  solved_ = false;
160  numVar_ = gm.numberOfVariables();
161  workingPoint_ = workingPoint;
162  lambda_ = 0.2;
163  constTerm_ = 0;
164  unaries_.resize(numVar_,0);
165  distance_.resize(numVar_,0);
166 
167  const LabelType label00[] = {0,0};
168  const LabelType label01[] = {0,1};
169  const LabelType label10[] = {1,0};
170  const LabelType label11[] = {1,1};
171  for(IndexType f=0; f<gm.numberOfFactors();++f){
172  OPENGM_ASSERT(gm[f].numberOfVariables() <= 2);
173  if(gm[f].numberOfVariables() == 0){
174  constTerm_ += gm[f](label00);
175  }
176  else if(gm[f].numberOfVariables() == 1){
177  const double v0 = gm[f](label00);
178  const double v1 = gm[f](label11);
179  const IndexType var0 = gm[f].variableIndex(0);
180  constTerm_ += v0;
181  unaries_[var0] += v1-v0;
182  }
183  else if(gm[f].numberOfVariables() == 2){
184  const double v00 = gm[f](label00);
185  const double v01 = gm[f](label01);
186  const double v10 = gm[f](label10);
187  const double v11 = gm[f](label11);
188  const IndexType var0 = gm[f].variableIndex(0);
189  const IndexType var1 = gm[f].variableIndex(1);
190  constTerm_ += v00;
191  const double D = 0.5*(v11-v00);
192  const double M = 0.5*(v10-v01);
193  unaries_[var0] += D+M;
194  unaries_[var1] += D-M;
195  const double V = v10-v00-D-M;
196  if(V>0){//submodular
197  subEdges_.push_back( LSA_TR_WeightedEdge(V,var0,var1));
198  }
199  else if(V<0){//supermodular
200  unaries_[var0] += V;
201  unaries_[var1] += V;
202  supEdges_.push_back( LSA_TR_WeightedEdge(-2*V,var0,var1));
203  }
204  }
205  }
206  std::cout << std::endl;
207  std::cout << subEdges_.size() <<" submodular edges."<<std::endl;
208  std::cout << supEdges_.size() <<" supermodular edges."<<std::endl;
209 
210  graph_ = new graph_type(gm.numberOfVariables(),subEdges_.size()+1);
211  changedList_ = new block_type(gm.numberOfVariables());
212 
213  graph_->add_node(numVar_);
214  for(size_t i=0; i<subEdges_.size(); ++i){
215  graph_->add_edge( subEdges_[i].u, subEdges_[i].v, subEdges_[i].w, subEdges_[i].w);
216  }
217  approxUnaries_.assign(unaries_.begin(),unaries_.end());
218  for(size_t i=0; i<supEdges_.size(); ++i){
219  const size_t var0 = supEdges_[i].u;
220  const size_t var1 = supEdges_[i].v;
221  const double w = supEdges_[i].w;
222  if(workingPoint[var0]==1)
223  approxUnaries_[var1] += w;
224  if(workingPoint[var1]==1)
225  approxUnaries_[var0] += w;
226  if(workingPoint[var0]==1 && workingPoint[var1]==1)
227  constTermApproximation_ -= w;
228  }
229 
230 
231  shape_.resize(1,numVar_);
232  std::vector<size_t> neigbor_count(numVar_,0);
233  for(size_t i=0; i<supEdges_.size(); ++i){
234  ++neigbor_count[supEdges_[i].u];
235  ++neigbor_count[supEdges_[i].v];
236  }
237  for(size_t i=0; i<subEdges_.size(); ++i){
238  ++neigbor_count[subEdges_[i].u];
239  ++neigbor_count[subEdges_[i].v];
240  }
241  size_t min_deg = *std::min_element(neigbor_count.begin(),neigbor_count.end());
242  std::vector<size_t> corners;
243  for(size_t i=0; i<neigbor_count.size(); ++i)
244  if (neigbor_count[i] == min_deg)
245  corners.push_back(i);
246  if(corners.size()==4){
247  if( !(corners[1]-corners[0] != corners[3]-corners[2])&&
248  !(corners[0] != 0 || corners[3] != numVar_-1) ){
249  shape_.resize(2);
250  shape_[0] = corners[1]-corners[0]+1;
251  shape_[1] = numVar_ / shape_[0];
252  }
253  }
254  if(shape_.size() ==1 && distanceType_ == EUCLIDEAN)
255  std::cout << "Warning : Shape of labeling is 1 and Euclidean distance does not make sense! Maybe autodetection of shape fails ..." <<std::endl;
256 
257 
258 
259  updateDistance();
260  constTermTrustRegion_ = 0;
261  for(int i=0; i<approxUnaries_.size(); ++i){
262  approxUnaries_[i] += lambda_*distance_[i];
263  graph_->add_tweights( i, 0, approxUnaries_[i]);
264  if(distance_[i]<0)
265  constTermTrustRegion_-=lambda_*distance_[i];
266  }
267  };
268 
269  template<class LabelType>
270  void LSA_TR_HELPER<LabelType>::updateDistance() {
271  if (distanceType_==HAMMING){
272  for(int i=0; i<numVar_; ++i){
273  if(workingPoint_[i]==0){
274  distance_[i] = 1;
275  }
276  else{
277  distance_[i] = -1;
278  }
279  }
280  }
281 #ifdef NOVIGRA
282  else if(distanceType_==EUCLIDEAN){
283  std::cout << "Warning : The useage of euclidean distance requires VIGRA!" <<std::endl;
284  std::cout << " Vigra is disabled -> Switch to Hamming distance!" <<std::endl;
285  distanceType_=HAMMING;
286  for(int i=0; i<numVar_; ++i){
287  if(workingPoint_[i]==0){
288  distance_[i] = 1;
289  }
290  else{
291  distance_[i] = -1;
292  }
293  }
294  }
295 #else
296  else if(distanceType_==EUCLIDEAN){
297  std::vector<size_t> s = shape_;
298  std::vector<double> dist0(numVar_,0);
299  std::vector<double> dist1(numVar_,0);
300  if(s.size()==1){
301  typedef vigra::MultiArrayView<1, LabelType> ArrayType;
302  typedef vigra::MultiArrayView<1, double> DArrayType;
303  typedef typename ArrayType::difference_type ShapeType;
304  ShapeType shape(s[0]);
305  ShapeType stride(1);
306 
307 
308  ArrayType source( shape, stride, &workingPoint_[0] );
309  DArrayType dest0( shape, stride, &dist0[0] );
310  DArrayType dest1( shape, stride, &dist1[0] );
311 
312  vigra::separableMultiDistance(source, dest0, false);
313  vigra::separableMultiDistance(source, dest1, true);
314  for(int i=0; i<numVar_; ++i){
315  if(workingPoint_[i]==0){
316  distance_[i] = (dist1[i]-0.5);
317  }
318  else{
319  distance_[i] = -(dist0[i]-0.5);
320  }
321  }
322  }
323  else if(s.size()==2){
324  typedef vigra::MultiArrayView<2, LabelType> ArrayType;
325  typedef vigra::MultiArrayView<2, double> DArrayType;
326  typedef typename ArrayType::difference_type ShapeType;
327  ShapeType shape(s[0],s[1]);
328  ShapeType stride(1,s[0]);
329 
330 
331  ArrayType source( shape, stride, &workingPoint_[0] );
332  DArrayType dest0( shape, stride, &dist0[0] );
333  DArrayType dest1( shape, stride, &dist1[0] );
334 
335  vigra::separableMultiDistance(source, dest0, false);
336  vigra::separableMultiDistance(source, dest1, true);
337  for(int i=0; i<numVar_; ++i){
338  if(workingPoint_[i]==0){
339  distance_[i] = (dist1[i]-0.5);
340  }
341  else{
342  distance_[i] = -(dist0[i]-0.5);
343  }
344  }
345  }
346  else if(s.size()==3){
347  typedef vigra::MultiArrayView<3, LabelType> ArrayType;
348  typedef vigra::MultiArrayView<3, double> DArrayType;
349  typedef typename ArrayType::difference_type ShapeType;
350  ShapeType shape(s[0],s[1],s[2]);
351  ShapeType stride(1,s[0],s[0]*s[1]);
352 
353 
354  ArrayType source( shape, stride, &workingPoint_[0] );
355  DArrayType dest0( shape, stride, &dist0[0] );
356  DArrayType dest1( shape, stride, &dist1[0] );
357 
358  vigra::separableMultiDistance(source, dest0, false);
359  vigra::separableMultiDistance(source, dest1, true);
360  for(int i=0; i<numVar_; ++i){
361  if(workingPoint_[i]==0){
362  distance_[i] = (dist1[i]-0.5);
363  }
364  else{
365  distance_[i] = -(dist0[i]-0.5);
366  }
367  }
368  }
369  }//end EUCLIDEAN
370 #endif
371  else{
372  std::cout <<"Unknown distance"<<std::endl;
373  }
374  return;
375  }
376 
377  template<class LabelType>
378  double LSA_TR_HELPER<LabelType>::optimize(std::vector<LabelType>& label){
379  double value;
380  //std::cout << lambda_ <<std::endl;
381  if(solved_){ //use warmstart
382  value = graph_->maxflow(true,changedList_);
383  for (typename graph_type::node_id* ptr = changedList_->ScanFirst(); ptr; ptr = changedList_->ScanNext()) {
384  typename graph_type::node_id var = *ptr;
385  OPENGM_ASSERT(var>=0 && var<numVar_);
386  graph_->remove_from_changed_list(var);
387  }
388 
389  for(size_t var=0; var<numVar_; ++var) {
390  if (graph_->what_segment(var) == graph_type::SOURCE) { label[var]=1;}
391  else { label[var]=0;}
392  }
393  changedList_->Reset();
394  }
395  else{ //first round without warmstart
396  value = graph_->maxflow();
397  for(size_t var=0; var<numVar_; ++var) {
398  if (graph_->what_segment(var) == graph_type::SOURCE) { label[var]=1;}
399  else { label[var]=0;}
400  }
401  solved_=true;
402  }
403  return value + constTerm_ + constTermApproximation_ + constTermTrustRegion_;
404  }
405 
406  template<class LabelType>
407  void LSA_TR_HELPER<LabelType>::set(const double newLambda){
408  if( newLambda == lambda_ ) return;
409  double difLambda = newLambda - lambda_;
410  lambda_ = newLambda;
411  constTermTrustRegion_ = 0;
412  if(solved_){
413  for(int i=0; i<approxUnaries_.size(); ++i){
414  double oldcap = graph_->get_trcap(i);
415  approxUnaries_[i] += difLambda*distance_[i];
416  graph_->add_tweights( i, 0, difLambda*distance_[i] );
417  if(distance_[i]<0)
418  constTermTrustRegion_ -= difLambda*distance_[i];
419  double newcap = graph_->get_trcap(i);
420  if (!((newcap > 0 && oldcap > 0)||(newcap < 0 && oldcap < 0))){
421  graph_->mark_node(i);
422  }
423  }
424  }else{
425  for(int i=0; i<approxUnaries_.size(); ++i){
426  approxUnaries_[i] += difLambda*distance_[i];
427  graph_->add_tweights( i, 0, difLambda*distance_[i] );
428  if(distance_[i]<0)
429  constTermTrustRegion_ -= difLambda*distance_[i];
430  }
431  }
432  }
433 
434  template<class LabelType>
435  void LSA_TR_HELPER<LabelType>::set(const std::vector<LabelType>& newWorkingPoint, const double newLambda){
436  workingPoint_ = newWorkingPoint;
437  lambda_ = newLambda;
438  constTermTrustRegion_ = 0;
439  constTermApproximation_ = 0;
440 
441  updateDistance();
442 
443  std::vector<double> newApproxUnaries = unaries_;
444  for(size_t i=0; i<supEdges_.size(); ++i){
445  const size_t var0 = supEdges_[i].u;
446  const size_t var1 = supEdges_[i].v;
447  const double w = supEdges_[i].w;
448  if(workingPoint_[var0]==1)
449  newApproxUnaries[var1] += w;
450  if(workingPoint_[var1]==1)
451  newApproxUnaries[var0] += w;
452  if(workingPoint_[var0]==1 && workingPoint_[var1]==1)
453  constTermApproximation_ -= w;
454  }
455  if(solved_){
456  for(int i=0; i<numVar_; ++i){
457  double oldcap = graph_->get_trcap(i);
458  newApproxUnaries[i] += lambda_*distance_[i];
459  graph_->add_tweights( i, 0, newApproxUnaries[i]-approxUnaries_[i] );
460  if(distance_[i]<0)
461  constTermTrustRegion_ -= lambda_*distance_[i];
462  double newcap = graph_->get_trcap(i);
463  if (!((newcap > 0 && oldcap > 0)||(newcap < 0 && oldcap < 0))){
464  graph_->mark_node(i);
465  }
466  }
467  }else{
468  for(int i=0; i<numVar_; ++i){
469  newApproxUnaries[i] += lambda_*distance_[i];
470  graph_->add_tweights( i, 0, newApproxUnaries[i]-approxUnaries_[i]);
471  if(distance_[i]<0)
472  constTermTrustRegion_ -= lambda_*distance_[i];
473  }
474  }
475  approxUnaries_.assign(newApproxUnaries.begin(),newApproxUnaries.end());
476 
477  }
478 
479 
480  template<class LabelType>
481  double LSA_TR_HELPER<LabelType>::eval(const std::vector<LabelType>& label) const
482  {
483  typedef double ValueType;
484  ValueType v = constTerm_;
485  for(size_t var=0; var<numVar_;++var)
486  if(label[var]==1)
487  v += unaries_[var];
488  for(size_t i=0; i<subEdges_.size(); ++i)
489  if(label[subEdges_[i].u] != label[subEdges_[i].v])
490  v += subEdges_[i].w;
491  for(size_t i=0; i<supEdges_.size(); ++i)
492  if(label[supEdges_[i].u] == 1 && label[supEdges_[i].v] == 1)
493  v += supEdges_[i].w;
494  return v;
495  }
496 
497  template<class LabelType>
498  double LSA_TR_HELPER<LabelType>::evalAprox(const std::vector<LabelType>& label, const std::vector<LabelType>& workingPoint, const double lambda) const
499  {
500  typedef double ValueType;
501  ValueType v = constTerm_;
502  for(size_t var=0; var<numVar_;++var)
503  if(label[var]==1)
504  v += unaries_[var];
505  for(size_t i=0; i<subEdges_.size(); ++i)
506  if(label[subEdges_[i].u] != label[subEdges_[i].v])
507  v += subEdges_[i].w;
508  for(size_t i=0; i<supEdges_.size(); ++i){
509  if(label[supEdges_[i].u] == 1 && workingPoint[supEdges_[i].v] == 1 )
510  v += supEdges_[i].w;
511  if(workingPoint[supEdges_[i].u] == 1 && label[supEdges_[i].v] == 1 )
512  v += supEdges_[i].w;
513  if(workingPoint[supEdges_[i].u] == 1 && workingPoint[supEdges_[i].v] == 1 )
514  v -= supEdges_[i].w;
515  }
516  for(size_t i=0; i<numVar_; ++i){
517  if(label[i] != workingPoint[i])
518  v += lambda * std::fabs(distance_[i]);
519  }
520  return v;
521  }
522 
523  template<class LabelType>
524  void LSA_TR_HELPER<LabelType>::evalBoth(const std::vector<LabelType>& label, const std::vector<LabelType>& workingPoint, const double lambda, double& value, double& valueAprox) const
525  {
526  value = constTerm_;
527  for(size_t var=0; var<numVar_;++var)
528  if(label[var]==1)
529  value += unaries_[var];
530  for(size_t i=0; i<subEdges_.size(); ++i)
531  if(label[subEdges_[i].u] != label[subEdges_[i].v])
532  value += subEdges_[i].w;
533  valueAprox = value;
534  for(size_t i=0; i<supEdges_.size(); ++i){
535  if(label[supEdges_[i].u] == 1 && label[supEdges_[i].v] == 1 )
536  value += supEdges_[i].w;
537  if(label[supEdges_[i].u] == 1 && workingPoint[supEdges_[i].v] == 1 )
538  valueAprox += supEdges_[i].w;
539  if(workingPoint[supEdges_[i].u] == 1 && label[supEdges_[i].v] == 1 )
540  valueAprox += supEdges_[i].w;
541  if(workingPoint[supEdges_[i].u] == 1 && workingPoint[supEdges_[i].v] == 1 )
542  valueAprox -= supEdges_[i].w;
543  }
544  for(size_t i=0; i<numVar_; ++i){
545  if(label[i] != workingPoint[i])
546  valueAprox += lambda * std::fabs(distance_[i]);
547  }
548  }
549 
550 
551 
553 
554  template<class GM, class ACC>
555  LSA_TR<GM, ACC>::~LSA_TR(){}
556 
557  template<class GM, class ACC>
558  inline
559  LSA_TR<GM, ACC>::LSA_TR
560  (
561  const GraphicalModelType& gm
562  )
563  : gm_(gm),
564  param_(Parameter())
565  {
566  init();
567  }
568 
569  template<class GM, class ACC>
570  LSA_TR<GM, ACC>::LSA_TR
571  (
572  const GraphicalModelType& gm,
573  const Parameter& parameter
574  )
575  : gm_(gm),
576  param_(parameter)
577  {
578  init();
579  }
580 
581  template<class GM, class ACC>
582  void LSA_TR<GM, ACC>::init()
583  {
584  srand(param_.randSeed_);
585  numVar_ = gm_.numberOfVariables();
586  curState_.resize(numVar_,1);
587  for (size_t i=0; i<numVar_; ++i) curState_[i]= rand()%2;
588  helper_.init(gm_, curState_);
589  if(param_.distance_ == Parameter::HAMMING)
590  helper_.setDistanceType(LSA_TR_HELPER<LabelType>::HAMMING);
591  else if(param_.distance_ == Parameter::EUCLIDEAN)
592  helper_.setDistanceType(LSA_TR_HELPER<LabelType>::EUCLIDEAN);
593  else
594  std::cout << "Warning: Unknown distance type !"<<std::endl;
595  }
596 
597 
598  template<class GM, class ACC>
599  inline void
600  LSA_TR<GM, ACC>::reset()
601  {
602  curState_.resize(numVar_,1);
603  }
604 
605  template<class GM, class ACC>
606  inline void
607  LSA_TR<GM,ACC>::setStartingPoint(typename std::vector<typename LSA_TR<GM,ACC>::LabelType>::const_iterator begin) {
608  curState_.assign(begin, begin+numVar_);
609  }
610 
611  template<class GM, class ACC>
612  inline std::string
613  LSA_TR<GM, ACC>::name() const
614  {
615  return "LSA_TR";
616  }
617 
618  template<class GM, class ACC>
619  inline const typename LSA_TR<GM, ACC>::GraphicalModelType&
620  LSA_TR<GM, ACC>::graphicalModel() const
621  {
622  return gm_;
623  }
624 
625  template<class GM, class ACC>
626  inline InferenceTermination
627  LSA_TR<GM,ACC>::infer()
628  {
629  EmptyVisitorType v;
630  return infer(v);
631  }
632 
633 
634  template<class GM, class ACC>
635  template<class VisitorType>
636  InferenceTermination LSA_TR<GM,ACC>::infer
637  (
638  VisitorType& visitor
639  )
640  {
641  const ValueType tau1 = 0;
642  const ValueType tau2 = param_.reductionRatio_;
643  bool exitInf=false;
644  std::vector<LabelType> label(numVar_);
645  std::vector<ValueType> energies;
646  std::vector<ValueType> energiesAprox;
647  double lambda = param_.initialLambda_;
648  helper_.set(curState_,lambda);
649  visitor.begin(*this);
650 
651  ValueType curr_value_aprox = helper_.evalAprox(curState_, curState_, lambda);
652  ValueType curr_value = helper_.eval(curState_);
653  bool changedWorkingpoint = false;
654  bool changedLambda = false;
655  ValueType value_after;
656  ValueType value_after_aprox;
657 
658  OPENGM_ASSERT(std::fabs(curr_value-gm_.evaluate(curState_))<0.0001);
659  for (size_t i=0; i<10000 ; ++i){
660  //std::cout << "round "<<i<<" (lambda = "<<lambda<<"): " <<std::endl;
661  if(lambda>param_.maxLambda_) break;
662 
663  if(changedWorkingpoint)
664  helper_.set(curState_,lambda);
665  else if(changedLambda)
666  helper_.set(lambda);
667  changedWorkingpoint = false;
668  changedLambda = false;
669  helper_.optimize(label);
670  helper_.evalBoth(label, curState_, lambda, value_after, value_after_aprox);
671 
672  //if(std::fabs(curr_value_aprox-curr_value)>0.0001)
673  // std::cout << "WARNING : "<< helper_.evalAprox(curState_, curState_, lambda) << " != " << helper_.eval(curState_) << " == " <<gm_.evaluate(curState_)<<std::endl;
674  OPENGM_ASSERT(std::fabs(helper_.eval(curState_)-gm_.evaluate(curState_))<0.0001);
675  OPENGM_ASSERT(helper_.eval(curState_) == curr_value);
676  OPENGM_ASSERT(std::fabs(helper_.eval(label)-gm_.evaluate(label))<0.0001);
677 
678  const ValueType P = curr_value_aprox - value_after_aprox;
679  const ValueType R = curr_value - value_after;
680 
681 
682 
683  //std::cout <<P << " " <<curr_value_aprox << " " << value_after_aprox <<std::endl;
684  if(P==0){
685  // ** Search for smallest possible step (largest penalty that give progress)
686  //std::cout << "Approximation does not improve energy ... searching for better lambda ... "<< std::flush;
687  lambda = findMinimalChangeBrakPoint(lambda, curState_);
688  helper_.set(lambda);
689  helper_.optimize(label);
690  //std::cout<<"set lambda to "<< lambda <<std::endl;
691 
692  helper_.evalBoth(label, curState_, lambda, value_after, value_after_aprox);
693  const ValueType P = curr_value_aprox - value_after_aprox;
694  const ValueType R = curr_value - value_after;
695  if(R<=0){
696  visitor(*this);
697  break;
698  }else if(R>0){
699  // ** Update Working Point
700  //std::cout<<"Update Working Point"<<std::endl;
701  curState_.assign(label.begin(),label.end());
702  changedWorkingpoint = true;
703  curr_value = value_after;
704  curr_value_aprox = value_after;
705  //OPENGM_ASSERT(std::fabs( curr_value_aprox-helper_.evalAprox(curState_, curState_, lambda) )<0.0001);
706  }
707  }
708  else{
709  if(P<0) std::cout << "WARNING : "<< curr_value_aprox << " < " << value_after_aprox << std::endl;
710  if(R>tau1){
711  // ** Update Working Point
712  //std::cout<<"Update Working Point"<<std::endl;
713  curState_.assign(label.begin(),label.end());
714  changedWorkingpoint = true;
715  //helper_.set(curState_,lambda);
716  curr_value = value_after;
717  curr_value_aprox = value_after;
718  //OPENGM_ASSERT(std::fabs( curr_value_aprox-helper_.evalAprox(curState_, curState_, lambda) )<0.0001);
719  }
720  }
721 
722  // ** Update trust region term
723  if(R/P>tau2){ ;
724  lambda = lambda / param_.lambdaMultiplier_;
725  changedLambda = true;
726  //helper_.set(lambda);
727  //std::cout<<"Decrease TR to "<< lambda <<std::endl;
728  }
729  else{
730  lambda = lambda * param_.lambdaMultiplier_;
731  changedLambda = true;
732  //helper_.set(lambda);
733  //std::cout<<"Increase TR to "<< lambda<<std::endl;
734  }
735 
736  // ** Store values
737  energies.push_back (curr_value);
738  energiesAprox.push_back(value_after_aprox);
739 
740  // ** Call Visitor
741  if( visitor(*this) != visitors::VisitorReturnFlag::ContinueInf ) break;
742  }
743 
744  visitor.end(*this);
745  return NORMAL;
746  }
747 
748  template<class GM, class ACC>
749  inline InferenceTermination
750  LSA_TR<GM,ACC>::arg
751  (
752  std::vector<LabelType>& x,
753  const size_t N
754  ) const
755  {
756  if(N==1) {
757  //x.resize(gm_.numberOfVariables());
758  //for (size_t i=0; i<gm_.numberOfVariables(); ++i)
759  // x[i] = curState_[i];
760  x.assign(curState_.begin(), curState_.end());
761  return NORMAL;
762  }
763  else {
764  return UNKNOWN;
765  }
766  }
767 
768 
769  template<class GM, class ACC>
770  double LSA_TR<GM,ACC>::findMinimalChangeBrakPoint(const double lambda, const std::vector<LabelType>& workingPoint){
771 
772  ValueType topLambda = lambda;
773  ValueType bottomLambda = param_.precisionLambda_;
774  std::vector<LabelType> topLabel(numVar_);
775  std::vector<LabelType> bottomLabel(numVar_);
776  std::vector<LabelType> label(numVar_);
777  // upper bound for best lambda
778  while(true){
779  helper_.set(topLambda);
780  helper_.optimize(topLabel);
781  if(!std::equal(topLabel.begin(),topLabel.end(),workingPoint.begin()))
782  topLambda = topLambda * 2;
783  else
784  break;
785  }
786 
787  // lower bound for lambda
788  helper_.set(bottomLambda);
789  helper_.optimize(bottomLabel);
790 
791  // binary search for minimal change point
792  while(true){
793  double middleLambda = (topLambda + bottomLambda)/2.0;
794  //std::cout <<"test "<< bottomLambda<<" < "<<middleLambda<<" < "<<topLambda<<std::endl;
795  helper_.set(middleLambda);
796  helper_.optimize(label);
797 
798  if(!std::equal(label.begin(),label.end(),topLabel.begin())){
799  bottomLambda = middleLambda;
800  bottomLabel = label;
801  }
802  else if(!std::equal(label.begin(),label.end(),bottomLabel.begin())){
803  topLambda = middleLambda;
804  topLabel = label;
805  }
806  else{
807  return bottomLambda;
808  }
809  if((topLambda-bottomLambda) < param_.precisionLambda_){
810  return bottomLambda;
811  }
812  }
813  }
814 
815 } // namespace opengm
816 
817 #endif // #ifndef OPENGM_LSATR_HXX
opengm::LSA_TR::Parameter::randSeed_
size_t randSeed_
Definition: lsatr.hxx:104
opengm::LSA_TR_HELPER::DISTANCE
DISTANCE
Definition: lsatr.hxx:50
opengm::LSA_TR::Parameter::initialLambda_
double initialLambda_
Definition: lsatr.hxx:106
opengm
The OpenGM namespace.
Definition: config.hxx:43
opengm::LSA_TR::GraphicalModelType
GM GraphicalModelType
Definition: lsatr.hxx:95
opengm::LSA_TR::EmptyVisitorType
opengm::visitors::EmptyVisitor< LSA_TR< GM, ACC > > EmptyVisitorType
Definition: lsatr.hxx:98
opengm::LSA_TR_HELPER::LSA_TR_HELPER
LSA_TR_HELPER()
Definition: lsatr.hxx:52
opengm::LSA_TR::Parameter::Parameter
Parameter()
Definition: lsatr.hxx:112
opengm::LSA_TR_HELPER::set
void set(const double)
Definition: lsatr.hxx:407
opengm::LSA_TR::Parameter::maxLambda_
double maxLambda_
Definition: lsatr.hxx:105
opengm::LSA_TR::Parameter::EUCLIDEAN
Definition: lsatr.hxx:103
opengm::LSA_TR::reset
void reset()
Definition: lsatr.hxx:600
opengm::LSA_TR::~LSA_TR
~LSA_TR()
Definition: lsatr.hxx:555
opengm::LSA_TR::name
std::string name() const
Definition: lsatr.hxx:613
opengm::LSA_TR::AccumulationType
ACC AccumulationType
Definition: lsatr.hxx:94
OPENGM_ASSERT
#define OPENGM_ASSERT(expression)
Definition: opengm.hxx:77
opengm::LSA_TR::Parameter
Definition: lsatr.hxx:101
opengm::LSA_TR_HELPER::evalAprox
double evalAprox(const std::vector< LabelType > &, const std::vector< LabelType > &, const double) const
Definition: lsatr.hxx:498
opengm::LSA_TR_HELPER::init
void init(const GM &, const std::vector< LabelType > &)
Definition: lsatr.hxx:157
opengm::LSA_TR
Definition: lsatr.hxx:91
opengm::LSA_TR_HELPER::setDistanceType
void setDistanceType(const DISTANCE d)
Definition: lsatr.hxx:59
opengm::LSA_TR_WeightedEdge::u
size_t u
Definition: lsatr.hxx:42
opengm::LSA_TR::TimingVisitorType
opengm::visitors::TimingVisitor< LSA_TR< GM, ACC > > TimingVisitorType
Definition: lsatr.hxx:99
opengm::LSA_TR::setStartingPoint
void setStartingPoint(typename std::vector< LabelType >::const_iterator)
set initial labeling
Definition: lsatr.hxx:607
opengm::LSA_TR_HELPER::optimize
double optimize(std::vector< LabelType > &)
Definition: lsatr.hxx:378
opengm::LSA_TR_WeightedEdge::w
double w
Definition: lsatr.hxx:41
opengm::Inference::ValueType
GraphicalModelType::ValueType ValueType
Definition: inference.hxx:41
opengm::LSA_TR_HELPER::EUCLIDEAN
Definition: lsatr.hxx:50
opengm::LSA_TR_WeightedEdge
Local Submodular Approximation with Trust Region regularization .
Definition: lsatr.hxx:39
opengm::Inference
Inference algorithm interface.
Definition: inference.hxx:34
opengm::LSA_TR::Parameter::lambdaMultiplier_
double lambdaMultiplier_
Definition: lsatr.hxx:108
opengm::LSA_TR::Parameter::reductionRatio_
double reductionRatio_
Definition: lsatr.hxx:109
opengm::LSA_TR_HELPER::~LSA_TR_HELPER
~LSA_TR_HELPER()
Definition: lsatr.hxx:53
opengm::LSA_TR::Parameter::precisionLambda_
double precisionLambda_
Definition: lsatr.hxx:107
opengm::LSA_TR::OPENGM_GM_TYPE_TYPEDEFS
OPENGM_GM_TYPE_TYPEDEFS
Definition: lsatr.hxx:96
opengm::LSA_TR_WeightedEdge::v
size_t v
Definition: lsatr.hxx:43
opengm::LSA_TR_HELPER
Definition: lsatr.hxx:48
opengm::LSA_TR_WeightedEdge::LSA_TR_WeightedEdge
LSA_TR_WeightedEdge(double aw, size_t au, size_t av)
Definition: lsatr.hxx:40
opengm::LSA_TR_HELPER::evalBoth
void evalBoth(const std::vector< LabelType > &label, const std::vector< LabelType > &workingPoint, const double lambda, double &value, double &valueAprox) const
Definition: lsatr.hxx:524
opengm::Inference::LabelType
GraphicalModelType::LabelType LabelType
Definition: inference.hxx:39
opengm::LSA_TR_HELPER::eval
double eval(const std::vector< LabelType > &) const
Definition: lsatr.hxx:481
opengm::LSA_TR::graphicalModel
const GraphicalModelType & graphicalModel() const
Definition: lsatr.hxx:620
opengm::LSA_TR::Parameter::HAMMING
Definition: lsatr.hxx:103
opengm::LSA_TR::Parameter::distance_
DISTANCE distance_
Definition: lsatr.hxx:110
opengm::LSA_TR::VerboseVisitorType
opengm::visitors::VerboseVisitor< LSA_TR< GM, ACC > > VerboseVisitorType
Definition: lsatr.hxx:97
opengm::LSA_TR::value
virtual ValueType value() const
return the solution (value)
Definition: lsatr.hxx:134
opengm::LSA_TR::arg
virtual InferenceTermination arg(std::vector< LabelType > &, const size_t=1) const
output a solution
Definition: lsatr.hxx:751
opengm::LSA_TR::Parameter::DISTANCE
DISTANCE
Definition: lsatr.hxx:103
opengm::LSA_TR_HELPER::HAMMING
Definition: lsatr.hxx:50
opengm::LSA_TR::LSA_TR
LSA_TR(const GraphicalModelType &)
Definition: lsatr.hxx:560
opengm::LSA_TR::infer
InferenceTermination infer()
Definition: lsatr.hxx:627
opengm::InferenceTermination
InferenceTermination
Definition: inference.hxx:24
Generated by   doxygen 1.8.9.1