source: sundry/stable/1.1/sundry/VolVolume.h @ 5999

Revision 5999, 21.6 KB checked in by wehart, 5 years ago (diff)

Merged revisions 5281-5998 via svnmerge from
 https://software.sandia.gov/svn/public/acro/sundry/trunk

........

r5296 | jberry | 2007-12-19 10:22:19 -0700 (Wed, 19 Dec 2007) | 2 lines


Use the multipleRandomizedSensorPlacements correctly.

........

r5376 | lafisk | 2008-05-23 13:56:37 -0600 (Fri, 23 May 2008) | 2 lines


Add link flags for static executables.

........

r5834 | wehart | 2008-11-23 20:01:45 -0700 (Sun, 23 Nov 2008) | 2 lines


Update of Copyright assertions.

........

r5853 | wehart | 2008-11-24 08:42:59 -0700 (Mon, 24 Nov 2008) | 2 lines


Update to change the license on these files to CPL.

........

r5901 | wehart | 2008-12-01 20:32:01 -0700 (Mon, 01 Dec 2008) | 2 lines


Portability fixes for gcc 4.3.1

........

Line 
1/*  _________________________________________________________________________
2 *
3 *  Acro: A Common Repository for Optimizers
4 *  Copyright (c) 2008 Sandia Corporation.
5 *  This software is distributed under the CPL License.
6 *  Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
7 *  the U.S. Government retains certain rights in this software.
8 *  For more information, see the README.txt file in the top Acro directory.
9 *  _________________________________________________________________________
10 */
11
12// Copyright (C) 2000, International Business Machines
13// Corporation and others.  All Rights Reserved.
14
15#ifndef __VOLUME_HPP__
16#define __VOLUME_HPP__
17
18#include <cfloat>
19#include <algorithm>
20#include <cstdio>
21#include <cmath>
22
23#ifndef VOL_DEBUG
24// When VOL_DEBUG is 1, we check vector indices
25#define VOL_DEBUG 0
26#endif
27
28template <class T> static inline T
29VolMax(register const T x, register const T y) {
30   return ((x) > (y)) ? (x) : (y);
31}
32
33template <class T> static inline T
34VolAbs(register const T x) {
35   return ((x) > 0) ? (x) : -(x);
36}
37
38//############################################################################
39
40#if defined(VOL_DEBUG) && (VOL_DEBUG != 0)
41#define VOL_TEST_INDEX(i, size)                 \
42{                                               \
43   if ((i) < 0 || (i) >= (size)) {              \
44      printf("bad VOL_?vector index\n");        \
45      abort();                                  \
46   }                                            \
47}
48#define VOL_TEST_SIZE(size)                     \
49{                                               \
50   if (s <= 0) {                                \
51      printf("bad VOL_?vector size\n");         \
52      abort();                                  \
53   }                                            \
54}
55#else
56#define VOL_TEST_INDEX(i, size)
57#define VOL_TEST_SIZE(size)
58#endif
59     
60//############################################################################
61
62class VOL_dvector;
63class VOL_ivector;
64class VOL_primal;
65class VOL_dual;
66class VOL_swing;
67class VOL_alpha_factor;
68class VOL_vh;
69class VOL_indc;
70class VOL_user_hooks;
71class VOL_problem;
72
73//############################################################################
74
75/**
76   This class contains the parameters controlling the Volume Algorithm
77*/
78struct VOL_parms {
79   /** initial value of lambda */
80   double lambdainit;
81   /** initial value of alpha */
82   double alphainit;
83   /** minimum value for alpha */
84   double alphamin;
85   /** when little progress is being done, we multiply alpha by alphafactor */
86   double alphafactor;
87
88   /** initial upper bound of the value of an integer solution */
89   double ubinit;
90
91   /** accept if max abs viol is less than this */
92   double primal_abs_precision;
93   /** accept if abs gap is less than this */
94   double gap_abs_precision;
95   /** accept if rel gap is less than this */
96   double gap_rel_precision; 
97   /** terminate if best_ub - lcost < granularity */
98   double granularity;
99
100   /** terminate if the relative increase in lcost through
101       <code>ascent_check_invl</code> steps is less than this */
102   double minimum_rel_ascent;
103   /** when to check for sufficient relative ascent the first time */
104   int    ascent_first_check;
105   /** through how many iterations does the relative ascent have to reach a
106       minimum */
107   int    ascent_check_invl;
108   
109   /** maximum number of iterations  */
110   int    maxsgriters;
111
112   /** controls the level of printing.
113       The flag should the the 'OR'-d value of the following options:
114       <ul>
115       <li> 0 - print nothing
116       <li> 1 - print iteration information
117       <li> 2 - add lambda information
118       <li> 4 - add number of Red, Yellow, Green iterations
119       </ul>
120       Default: 3
121   */
122   int    printflag;
123   /** controls how often do we print */
124   int    printinvl;
125   /** controls how often we run the primal heuristic */
126   int    heurinvl;
127
128   /** how many consecutive green iterations are allowed before changing
129       lambda */
130   int greentestinvl;
131   /** how many consecutive yellow iterations are allowed before changing
132       lambda */
133   int yellowtestinvl;
134   /** how many consecutive red iterations are allowed before changing
135       lambda */
136   int redtestinvl;
137
138   /** number of iterations before we check if alpha should be decreased */
139   int    alphaint;
140
141   /** name of file for saving dual solution */
142   char* temp_dualfile;
143};
144
145//############################################################################
146
147/** vector of doubles. It is used for most vector operations.
148
149    Note: If <code>VOL_DEBUG</code> is <code>#defined</code> to be 1 then each
150    time an entry is accessed in the vector the index of the entry is tested
151    for nonnegativity and for being less than the size of the vector. It's
152    good to turn this on while debugging, but in final runs it should be
153    turned off (beause of the performance hit).
154*/
155class VOL_dvector {
156public:
157   /** The array holding the vector */
158   double* v;
159   /** The size of the vector */
160   int sz;
161
162public:
163   /** Construct a vector of size s. The content of the vector is undefined. */
164   VOL_dvector(const int s) {
165      VOL_TEST_SIZE(s);
166      v = new double[sz = s];
167   }
168   /** Default constructor creates a vector of size 0. */
169   VOL_dvector() : v(0), sz(0) {}
170   /** Copy constructor makes a replica of x. */
171   VOL_dvector(const VOL_dvector& x) : v(0), sz(0) {
172      sz = x.sz;
173      if (sz > 0) {
174         v = new double[sz];
175         std::copy(x.v, x.v + sz, v);
176      }
177   }
178   /** The destructor deletes the data array. */
179   ~VOL_dvector() { delete[] v; }
180
181   /** Return the size of the vector. */
182   inline int size() const {return sz;}
183
184   /** Return a reference to the <code>i</code>-th entry. */
185   inline double& operator[](const int i) {
186      VOL_TEST_INDEX(i, sz);
187      return v[i];
188   }
189
190   /** Return the <code>i</code>-th entry. */
191   inline double operator[](const int i) const {
192      VOL_TEST_INDEX(i, sz);
193      return v[i];
194   }
195
196   /** Delete the content of the vector and replace it with a vector of length
197       0. */
198   inline void clear() {
199      delete[] v;
200      v = 0;
201      sz = 0;
202   }
203   /** Convex combination. Replace the current vector <code>v</code> with
204       <code>v = (1-gamma) v + gamma w</code>. */
205   inline void cc(const double gamma, const VOL_dvector& w) {
206      if (sz != w.sz) {
207         printf("bad VOL_dvector sizes\n");
208         abort();
209      }
210      double * p_v = v - 1;
211      const double * p_w = w.v - 1;
212      const double * const p_e = v + sz;
213      const double one_gamma = 1.0 - gamma;
214      while ( ++p_v != p_e ){
215         *p_v = one_gamma * (*p_v) + gamma * (*++p_w);
216      }
217   }
218
219   /** delete the current vector and allocate space for a vector of size
220       <code>s</code>. */
221   inline void allocate(const int s) {
222      VOL_TEST_SIZE(s);
223      delete[] v;                 
224      v = new double[sz = s];
225   }
226
227   /** swaps the vector with <code>w</code>. */
228   inline void swap(VOL_dvector& w) {
229      std::swap(v, w.v);
230      std::swap(sz, w.sz);
231   }
232
233   /** Copy <code>w</code> into the vector. */
234   VOL_dvector& operator=(const VOL_dvector& w);
235   /** Replace every entry in the vector with <code>w</code>. */
236   VOL_dvector& operator=(const double w);
237};
238
239//-----------------------------------------------------------------------------
240/** vector of ints. It's used to store indices, it has similar
241    functions as VOL_dvector.
242
243    Note: If <code>VOL_DEBUG</code> is <code>#defined</code> to be 1 then each
244    time an entry is accessed in the vector the index of the entry is tested
245    for nonnegativity and for being less than the size of the vector. It's
246    good to turn this on while debugging, but in final runs it should be
247    turned off (beause of the performance hit).
248*/
249class VOL_ivector {
250public:
251   /** The array holding the vector. */
252   int* v;
253   /** The size of the vector. */
254   int sz;
255public:
256   /** Construct a vector of size s. The content of the vector is undefined. */
257   VOL_ivector(const int s) {
258      VOL_TEST_SIZE(s);
259      v = new int[sz = s];
260   }
261   /** Default constructor creates a vector of size 0. */
262   VOL_ivector() : v(0), sz(0) {}
263   /** Copy constructor makes a replica of x. */
264   VOL_ivector(const VOL_ivector& x) {
265      sz = x.sz;
266      if (sz > 0) {
267         v = new int[sz];
268         std::copy(x.v, x.v + sz, v);
269      }
270   }
271   /** The destructor deletes the data array. */
272   ~VOL_ivector(){
273      delete [] v;
274   }
275
276   /** Return the size of the vector. */
277   inline int size() const { return sz; }
278   /** Return a reference to the <code>i</code>-th entry. */
279   inline int& operator[](const int i) {
280      VOL_TEST_INDEX(i, sz);
281      return v[i];
282   }
283
284   /** Return the <code>i</code>-th entry. */
285   inline int operator[](const int i) const {
286      VOL_TEST_INDEX(i, sz);
287      return v[i];
288   }
289
290   /** Delete the content of the vector and replace it with a vector of length
291       0. */
292   inline void clear() {
293      delete[] v;
294      v = 0;
295      sz = 0;
296   }
297
298   /** delete the current vector and allocate space for a vector of size
299       <code>s</code>. */
300   inline void allocate(const int s) {
301      VOL_TEST_SIZE(s);
302      delete[] v;
303      v = new int[sz = s];
304   }
305
306   /** swaps the vector with <code>w</code>. */
307   inline void swap(VOL_ivector& w) {
308      std::swap(v, w.v);
309      std::swap(sz, w.sz);
310   }
311
312   /** Copy <code>w</code> into the vector. */
313   VOL_ivector& operator=(const VOL_ivector& v);     
314   /** Replace every entry in the vector with <code>w</code>. */
315   VOL_ivector& operator=(const int w);
316};
317
318//############################################################################
319// A class describing a primal solution. This class is used only internally
320class VOL_primal {
321public:
322   // objective value of this primal solution
323   double value;
324   // the largest of the v[i]'s
325   double viol; 
326   // primal solution 
327   VOL_dvector x;
328   // v=b-Ax, for the relaxed constraints
329   VOL_dvector v;
330
331   VOL_primal(const int psize, const int dsize) : x(psize), v(dsize) {}
332   VOL_primal(const VOL_primal& primal) :
333      value(primal.value), viol(primal.viol), x(primal.x), v(primal.v) {}
334   ~VOL_primal() {}
335   inline VOL_primal& operator=(const VOL_primal& p) {
336      if (this == &p)
337         return *this;
338      value = p.value;
339      viol = p.viol;
340      x = p.x;
341      v = p.v;
342      return *this;
343   }
344
345   // convex combination. data members in this will be overwritten
346   // convex combination between two primal solutions
347   // x <-- alpha x + (1 - alpha) p.x
348   // v <-- alpha v + (1 - alpha) p.v
349   inline void cc(const double alpha, const VOL_primal& p) {
350      value = alpha * p.value + (1.0 - alpha) * value;
351      x.cc(alpha, p.x);
352      v.cc(alpha, p.v);
353   }
354   // find maximum of v[i]
355   void find_max_viol(const VOL_dvector& dual_lb,
356                      const VOL_dvector& dual_ub);
357};
358
359//-----------------------------------------------------------------------------
360// A class describing a dual solution. This class is used only internally
361class VOL_dual {
362public:
363   // lagrangian value
364   double lcost;
365   // reduced costs * (pstar-primal)
366   double xrc;
367   // this information is only printed
368   // dual vector
369   VOL_dvector u;
370
371   VOL_dual(const int dsize) : u(dsize) { u = 0.0;}
372   VOL_dual(const VOL_dual& dual) :
373      lcost(dual.lcost), xrc(dual.xrc), u(dual.u) {}
374   ~VOL_dual() {}
375   inline VOL_dual& operator=(const VOL_dual& p) {
376      if (this == &p)
377         return *this;
378      lcost = p.lcost;
379      xrc = p.xrc;
380      u = p.u;
381      return *this;
382   }
383   // dual step
384   void   step(const double target, const double lambda,
385               const VOL_dvector& dual_lb, const VOL_dvector& dual_ub,
386               const VOL_dvector& v);
387   double ascent(const VOL_dvector& v, const VOL_dvector& last_u) const;
388   void   compute_xrc(const VOL_dvector& pstarx, const VOL_dvector& primalx,
389                      const VOL_dvector& rc);
390
391};
392
393
394//############################################################################
395/* here we check whether an iteration is green, yellow or red. Also according
396   to this information we decide whether lambda should be changed */
397class VOL_swing {
398private:
399   VOL_swing(const VOL_swing&);
400   VOL_swing& operator=(const VOL_swing&);
401public:
402   enum condition {green, yellow, red} lastswing;
403   int lastgreeniter, lastyellowiter, lastrediter;
404   int ngs, nrs, nys;
405   int rd;
406   
407   VOL_swing() {
408      lastgreeniter = lastyellowiter = lastrediter = 0;
409      ngs = nrs = nys = 0;
410   }
411   ~VOL_swing(){}
412
413   inline void cond(const VOL_dual& dual,
414                    const double lcost, const double ascent, const int iter) {
415      double eps = 1.e-3;
416
417      if (ascent > 0.0  &&  lcost > dual.lcost + eps) {
418         lastswing = green;
419         lastgreeniter = iter;
420         ++ngs;
421         rd = 0;
422      } else {
423         if (ascent <= 0  &&  lcost > dual.lcost) {
424            lastswing = yellow;
425            lastyellowiter = iter;
426            ++nys;
427            rd = 0;
428         } else {
429            lastswing = red;
430            lastrediter = iter;
431            ++nrs;
432            rd = 1;
433         }
434      }
435   }
436
437   inline double
438   lfactor(const VOL_parms& parm, const double lambda, const int iter) {
439      double lambdafactor = 1.0;
440      double eps = 5.e-4;
441      int cons;
442
443      switch (lastswing) {
444       case green:
445         cons = iter - VolMax(lastyellowiter, lastrediter);
446         if (parm.printflag & 4)
447            printf("      G: Consecutive Gs = %3d\n\n", cons);
448         if (cons >= parm.greentestinvl && lambda < 2.0) {
449            lastgreeniter = lastyellowiter = lastrediter = iter;
450            lambdafactor = 2.0;
451            if (parm.printflag & 2)
452               printf("\n ---- increasing lamda to %g ----\n\n",
453                      lambda * lambdafactor);
454         }
455         break;
456     
457       case yellow:
458         cons = iter - VolMax(lastgreeniter, lastrediter);
459         if (parm.printflag & 4)
460            printf("      Y: Consecutive Ys = %3d\n\n", cons);
461         if (cons >= parm.yellowtestinvl) {
462            lastgreeniter = lastyellowiter = lastrediter = iter;
463            lambdafactor = 1.1;
464            if (parm.printflag & 2)
465               printf("\n **** increasing lamda to %g *****\n\n",
466                      lambda * lambdafactor);
467         }
468         break;
469     
470       case red:
471         cons = iter - VolMax(lastgreeniter, lastyellowiter);
472         if (parm.printflag & 4)
473            printf("      R: Consecutive Rs = %3d\n\n", cons);
474         if (cons >= parm.redtestinvl && lambda > eps) {
475            lastgreeniter = lastyellowiter = lastrediter = iter;
476            lambdafactor = 0.67;
477            if (parm.printflag & 2)
478               printf("\n **** decreasing lamda to %g *****\n\n",
479                      lambda * lambdafactor);
480         }
481         break;
482      }
483      return lambdafactor;
484   }
485
486   inline void
487   print() {
488      printf("**** G= %i, Y= %i, R= %i ****\n", ngs, nys, nrs);
489      ngs = nrs = nys = 0; 
490   }
491};
492
493//############################################################################
494/* alpha should be decreased if after some number of iterations the objective
495   has increased less that 1% */
496class VOL_alpha_factor {
497private:
498   VOL_alpha_factor(const VOL_alpha_factor&);
499   VOL_alpha_factor& operator=(const VOL_alpha_factor&);
500public:
501   double lastvalue;
502
503   VOL_alpha_factor() {lastvalue = -DBL_MAX;}
504   ~VOL_alpha_factor() {}
505
506   inline double factor(const VOL_parms& parm, const double lcost,
507                        const double alpha) {
508      if (alpha < parm.alphamin)
509         return 1.0;
510      register const double ll = VolAbs(lcost);
511      const double x = ll > 10 ? (lcost-lastvalue)/ll : (lcost-lastvalue);
512      lastvalue = lcost;
513      return (x <= 0.01) ? parm.alphafactor : 1.0;
514   }
515};
516
517//############################################################################
518/* here we compute the norm of the conjugate direction -hh-, the norm of the
519   subgradient -norm-, the inner product between the subgradient and the
520   last conjugate direction -vh-, and the inner product between the new
521   conjugate direction and the subgradient */
522class VOL_vh {
523private:
524   VOL_vh(const VOL_vh&);
525   VOL_vh& operator=(const VOL_vh&);
526public:
527   double hh;
528   double norm;
529   double vh;
530   double asc;
531
532   VOL_vh(const double alpha,
533          const VOL_dvector& dual_lb, const VOL_dvector& dual_ub,
534          const VOL_dvector& v, const VOL_dvector& vstar,
535          const VOL_dvector& u);
536   ~VOL_vh(){}
537};
538
539//############################################################################
540/* here we compute different parameter to be printed. v2 is the square of
541   the norm of the subgradient. vu is the inner product between the dual
542   variables and the subgradient. vabs is the maximum absolute value of
543   the violations of pstar. asc is the inner product between the conjugate
544   direction and the subgradient */
545class VOL_indc {
546private:
547   VOL_indc(const VOL_indc&);
548   VOL_indc& operator=(const VOL_indc&);
549public:
550   double v2;
551   double vu;
552   double vabs;
553   double asc;
554
555public:
556   VOL_indc(const VOL_dvector& dual_lb, const VOL_dvector& dual_ub,
557            const VOL_primal& primal, const VOL_primal& pstar,
558            const VOL_dual& dual);
559   ~VOL_indc() {}
560};
561
562//#############################################################################
563
564/** The user hooks should be overridden by the user to provide the
565    problem specific routines for the volume algorithm. The user
566    should derive a class ...
567
568    for all hooks: return value of -1 means that volume should quit
569*/
570class VOL_user_hooks {
571public:
572   virtual ~VOL_user_hooks() {}
573public:
574   // for all hooks: return value of -1 means that volume should quit
575   /** compute reduced costs   
576       @param u (IN) the dual variables
577       @param rc (OUT) the reduced cost with respect to the dual values
578   */
579   virtual int compute_rc(const VOL_dvector& u, VOL_dvector& rc) = 0;
580
581   /** Solve the subproblem for the subgradient step.
582       @param dual (IN) the dual variables
583       @param rc (IN) the reduced cost with respect to the dual values
584       @param lcost (OUT) the lagrangean cost with respect to the dual values
585       @param x (OUT) the primal result of solving the subproblem
586       @param v (OUT) b-Ax for the relaxed constraints
587       @param pcost (OUT) the primal objective value of <code>x</code>
588   */
589   virtual int solve_subproblem(const VOL_dvector& dual, const VOL_dvector& rc,
590                                double& lcost, VOL_dvector& x, VOL_dvector& v,
591                                double& pcost) = 0;
592   /** Starting from the primal vector x, run a heuristic to produce
593       an integer solution 
594       @param x (IN) the primal vector
595       @param heur_val (OUT) the value of the integer solution (return
596       <code>DBL_MAX</code> here if no feas sol was found
597   */
598   virtual int heuristics(const VOL_problem& p,
599                          const VOL_dvector& x, double& heur_val, double lb)=0;
600};
601
602//#############################################################################
603
604/** This class holds every data for the Volume Algorithm and its
605    <code>solve</code> method must be invoked to solve the problem.
606
607    The INPUT fields must be filled out completely before <code>solve</code>
608    is invoked. <code>dsol</code> have to be filled out if and only if the
609    last argument to <code>solve</code> is <code>true</code>.
610*/
611
612class VOL_problem {
613private:
614   VOL_problem(const VOL_problem&);
615   VOL_problem& operator=(const VOL_problem&);
616   void set_default_parm();
617   // ############ INPUT fields ########################
618public:
619   /**@name Constructors and destructor */
620   //@{
621   /** Default constructor. */
622   VOL_problem();
623   /** Create a a <code>VOL_problem</code> object and read in the parameters
624       from <code>filename</code>. */
625   VOL_problem(const char *filename);
626   /** Destruct the object. */
627   ~VOL_problem();
628   //@}
629
630   /**@name Method to solve the problem. */
631   //@{
632   /** Solve the problem using the <code>hooks</code>. Any information needed
633       in the hooks must be stored in the structure <code>user_data</code>
634       points to. */
635   int solve(VOL_user_hooks& hooks, const bool use_preset_dual = false);
636   //@}
637
638private:
639   /**@name Internal data (may be inquired for) */
640   //@{
641   /** value of alpha */
642   double alpha_;
643   /** value of lambda */
644   double lambda_;
645   // This union is here for padding (so that data members would be
646   // double-aligned on x86 CPU
647   union {
648      /** iteration number */
649      int iter_;
650      double __pad0;
651   };
652   //@}
653
654public:
655 
656   /**@name External data (containing the result after solve) */
657   //@{
658   /** final lagrangian value (OUTPUT) */
659   double value;
660   /** final dual solution (INPUT/OUTPUT) */
661   VOL_dvector dsol;
662   /** final primal solution (OUTPUT) */
663   VOL_dvector psol;
664   /** violations (b-Ax) for the relaxed constraints */
665   VOL_dvector viol;
666   //@}
667
668   /**@name External data (may be changed by the user before calling solve) */
669   //@{
670   /** The parameters controlling the Volume Algorithm (INPUT) */
671   VOL_parms parm;
672   /** length of primal solution (INPUT) */
673   int psize;       
674   /** length of dual solution (INPUT) */
675   int dsize;     
676   /** lower bounds for the duals (if 0 length, then filled with -inf) (INPUT)
677    */
678   VOL_dvector dual_lb;
679   /** upper bounds for the duals (if 0 length, then filled with +inf) (INPUT)
680    */
681   VOL_dvector dual_ub;
682   //@}
683
684public:
685   /**@name Methods returning final data */
686   //@{
687   /** returns the iteration number */
688   int    iter() const { return iter_; }
689   /** returns the value of alpha */
690   double alpha() const { return alpha_; }
691   /** returns the value of lambda */
692   double lambda() const { return lambda_; }
693   //@}
694
695private:
696   /**@name Private methods used internally */
697   //@{
698   /** Read in the parameters from the file <code>filename</code>. */
699   void read_params(const char* filename);
700
701   /** initializes duals, bounds for the duals, alpha, lambda */
702   int initialize(const bool use_preset_dual);
703
704   /** print volume info every parm.printinvl iterations */
705   void print_info(const int iter,
706                   const VOL_primal& primal, const VOL_primal& pstar,
707                   const VOL_dual& dual);
708
709   /** Checks if lcost is close to the target, if so it increases the target.
710       Close means that we got within 5% of the target. */
711   double readjust_target(const double oldtarget, const double lcost) const;
712
713   /** Here we decide the value of alpha1 to be used in the convex
714       combination. The new pstar will be computed as <br>
715       pstar = alpha1 * pstar + (1 - alpha1) * primal <br>
716       More details of this are in doc.ps. <br>
717       IN:  alpha, primal, pstar, dual <br>
718       @return alpha1
719   */
720   double power_heur(const VOL_primal& primal, const VOL_primal& pstar,
721                     const VOL_dual& dual) const;
722   //@}
723};
724
725#endif
Note: See TracBrowser for help on using the repository browser.