Main Page
Classes
Files
File List
File Members
opt
build
coinor-vol
coinor-vol-1.1.7
Vol
src
VolVolume.hpp
Go to the documentation of this file.
1
// Copyright (C) 2000, International Business Machines
2
// Corporation and others. All Rights Reserved.
3
#ifndef __VOLUME_HPP__
4
#define __VOLUME_HPP__
5
6
#include <cfloat>
7
#include <algorithm>
8
#include <cstdio>
9
#include <cmath>
10
11
#ifndef VOL_DEBUG
12
// When VOL_DEBUG is 1, we check vector indices
13
#define VOL_DEBUG 0
14
#endif
15
16
template
<
class
T>
static
inline
T
17
VolMax
(
register
const
T x,
register
const
T y) {
18
return
((x) > (y)) ? (x) : (y);
19
}
20
21
template
<
class
T>
static
inline
T
22
VolAbs
(
register
const
T x) {
23
return
((x) > 0) ? (x) : -(x);
24
}
25
26
//############################################################################
27
28
#if defined(VOL_DEBUG) && (VOL_DEBUG != 0)
29
#define VOL_TEST_INDEX(i, size) \
30
{ \
31
if ((i) < 0 || (i) >= (size)) { \
32
printf("bad VOL_?vector index\n"); \
33
abort(); \
34
} \
35
}
36
#define VOL_TEST_SIZE(size) \
37
{ \
38
if (s <= 0) { \
39
printf("bad VOL_?vector size\n"); \
40
abort(); \
41
} \
42
}
43
#else
44
#define VOL_TEST_INDEX(i, size)
45
#define VOL_TEST_SIZE(size)
46
#endif
47
48
//############################################################################
49
50
class
VOL_dvector
;
51
class
VOL_ivector
;
52
class
VOL_primal
;
53
class
VOL_dual
;
54
class
VOL_swing
;
55
class
VOL_alpha_factor
;
56
class
VOL_vh
;
57
class
VOL_indc
;
58
class
VOL_user_hooks
;
59
class
VOL_problem
;
60
61
//############################################################################
62
66
struct
VOL_parms
{
68
double
lambdainit
;
70
double
alphainit
;
72
double
alphamin
;
74
double
alphafactor
;
75
77
double
ubinit
;
78
80
double
primal_abs_precision
;
82
double
gap_abs_precision
;
84
double
gap_rel_precision
;
86
double
granularity
;
87
90
double
minimum_rel_ascent
;
92
int
ascent_first_check
;
95
int
ascent_check_invl
;
96
98
int
maxsgriters
;
99
110
int
printflag
;
112
int
printinvl
;
114
int
heurinvl
;
115
118
int
greentestinvl
;
121
int
yellowtestinvl
;
124
int
redtestinvl
;
125
127
int
alphaint
;
128
130
char
*
temp_dualfile
;
131
};
132
133
//############################################################################
134
143
class
VOL_dvector
{
144
public
:
146
double
*
v
;
148
int
sz
;
149
150
public
:
152
VOL_dvector
(
const
int
s) {
153
VOL_TEST_SIZE
(s);
154
v
=
new
double
[
sz
= s];
155
}
157
VOL_dvector
() :
v
(0),
sz
(0) {}
159
VOL_dvector
(
const
VOL_dvector
& x) :
v
(0),
sz
(0) {
160
sz
= x.
sz
;
161
if
(
sz
> 0) {
162
v
=
new
double
[
sz
];
163
std::copy(x.
v
, x.
v
+
sz
,
v
);
164
}
165
}
167
~VOL_dvector
() {
delete
[]
v
; }
168
170
inline
int
size
()
const
{
return
sz
;}
171
173
inline
double
&
operator[]
(
const
int
i) {
174
VOL_TEST_INDEX
(i,
sz
);
175
return
v
[i];
176
}
177
179
inline
double
operator[]
(
const
int
i)
const
{
180
VOL_TEST_INDEX
(i,
sz
);
181
return
v
[i];
182
}
183
186
inline
void
clear
() {
187
delete
[]
v
;
188
v
= 0;
189
sz
= 0;
190
}
193
inline
void
cc
(
const
double
gamma,
const
VOL_dvector
& w) {
194
if
(
sz
!= w.
sz
) {
195
printf(
"bad VOL_dvector sizes\n"
);
196
abort();
197
}
198
double
* p_v =
v
- 1;
199
const
double
* p_w = w.
v
- 1;
200
const
double
*
const
p_e =
v
+
sz
;
201
const
double
one_gamma = 1.0 - gamma;
202
while
( ++p_v != p_e ){
203
*p_v = one_gamma * (*p_v) + gamma * (*++p_w);
204
}
205
}
206
209
inline
void
allocate
(
const
int
s) {
210
VOL_TEST_SIZE
(s);
211
delete
[]
v
;
212
v
=
new
double
[
sz
= s];
213
}
214
216
inline
void
swap
(
VOL_dvector
& w) {
217
std::swap(
v
, w.
v
);
218
std::swap(
sz
, w.
sz
);
219
}
220
222
VOL_dvector
&
operator=
(
const
VOL_dvector
& w);
224
VOL_dvector
&
operator=
(
const
double
w);
225
};
226
227
//-----------------------------------------------------------------------------
237
class
VOL_ivector
{
238
public
:
240
int
*
v
;
242
int
sz
;
243
public
:
245
VOL_ivector
(
const
int
s) {
246
VOL_TEST_SIZE
(s);
247
v
=
new
int
[
sz
= s];
248
}
250
VOL_ivector
() :
v
(0),
sz
(0) {}
252
VOL_ivector
(
const
VOL_ivector
& x) {
253
sz
= x.
sz
;
254
if
(
sz
> 0) {
255
v
=
new
int
[
sz
];
256
std::copy(x.
v
, x.
v
+
sz
,
v
);
257
}
258
}
260
~VOL_ivector
(){
261
delete
[]
v
;
262
}
263
265
inline
int
size
()
const
{
return
sz
; }
267
inline
int
&
operator[]
(
const
int
i) {
268
VOL_TEST_INDEX
(i,
sz
);
269
return
v
[i];
270
}
271
273
inline
int
operator[]
(
const
int
i)
const
{
274
VOL_TEST_INDEX
(i,
sz
);
275
return
v
[i];
276
}
277
280
inline
void
clear
() {
281
delete
[]
v
;
282
v
= 0;
283
sz
= 0;
284
}
285
288
inline
void
allocate
(
const
int
s) {
289
VOL_TEST_SIZE
(s);
290
delete
[]
v
;
291
v
=
new
int
[
sz
= s];
292
}
293
295
inline
void
swap
(
VOL_ivector
& w) {
296
std::swap(
v
, w.
v
);
297
std::swap(
sz
, w.
sz
);
298
}
299
301
VOL_ivector
&
operator=
(
const
VOL_ivector
&
v
);
303
VOL_ivector
&
operator=
(
const
int
w);
304
};
305
306
//############################################################################
307
// A class describing a primal solution. This class is used only internally
308
class
VOL_primal
{
309
public
:
310
// objective value of this primal solution
311
double
value
;
312
// the largest of the v[i]'s
313
double
viol
;
314
// primal solution
315
VOL_dvector
x
;
316
// v=b-Ax, for the relaxed constraints
317
VOL_dvector
v
;
318
319
VOL_primal
(
const
int
psize,
const
int
dsize) :
x
(psize),
v
(dsize) {}
320
VOL_primal
(
const
VOL_primal
& primal) :
321
value
(primal.
value
),
viol
(primal.
viol
),
x
(primal.
x
),
v
(primal.
v
) {}
322
~VOL_primal
() {}
323
inline
VOL_primal
&
operator=
(
const
VOL_primal
& p) {
324
if
(
this
== &p)
325
return
*
this
;
326
value
= p.
value
;
327
viol
= p.
viol
;
328
x
= p.
x
;
329
v
= p.
v
;
330
return
*
this
;
331
}
332
333
// convex combination. data members in this will be overwritten
334
// convex combination between two primal solutions
335
// x <-- alpha x + (1 - alpha) p.x
336
// v <-- alpha v + (1 - alpha) p.v
337
inline
void
cc
(
const
double
alpha,
const
VOL_primal
& p) {
338
value
= alpha * p.
value
+ (1.0 - alpha) *
value
;
339
x
.
cc
(alpha, p.
x
);
340
v
.
cc
(alpha, p.
v
);
341
}
342
// find maximum of v[i]
343
void
find_max_viol
(
const
VOL_dvector
& dual_lb,
344
const
VOL_dvector
& dual_ub);
345
};
346
347
//-----------------------------------------------------------------------------
348
// A class describing a dual solution. This class is used only internally
349
class
VOL_dual
{
350
public
:
351
// lagrangian value
352
double
lcost
;
353
// reduced costs * (pstar-primal)
354
double
xrc
;
355
// this information is only printed
356
// dual vector
357
VOL_dvector
u
;
358
359
VOL_dual
(
const
int
dsize) :
u
(dsize) {
u
= 0.0;}
360
VOL_dual
(
const
VOL_dual
& dual) :
361
lcost
(dual.
lcost
),
xrc
(dual.
xrc
),
u
(dual.
u
) {}
362
~VOL_dual
() {}
363
inline
VOL_dual
&
operator=
(
const
VOL_dual
& p) {
364
if
(
this
== &p)
365
return
*
this
;
366
lcost
= p.
lcost
;
367
xrc
= p.
xrc
;
368
u
= p.
u
;
369
return
*
this
;
370
}
371
// dual step
372
void
step
(
const
double
target,
const
double
lambda,
373
const
VOL_dvector
& dual_lb,
const
VOL_dvector
& dual_ub,
374
const
VOL_dvector
& v);
375
double
ascent
(
const
VOL_dvector
& v,
const
VOL_dvector
& last_u)
const
;
376
void
compute_xrc
(
const
VOL_dvector
& pstarx,
const
VOL_dvector
& primalx,
377
const
VOL_dvector
& rc);
378
379
};
380
381
382
//############################################################################
383
/* here we check whether an iteration is green, yellow or red. Also according
384
to this information we decide whether lambda should be changed */
385
class
VOL_swing
{
386
private
:
387
VOL_swing
(
const
VOL_swing
&);
388
VOL_swing
&
operator=
(
const
VOL_swing
&);
389
public
:
390
enum
condition
{
green
,
yellow
,
red
}
lastswing
;
391
int
lastgreeniter
,
lastyellowiter
,
lastrediter
;
392
int
ngs
,
nrs
,
nys
;
393
int
rd
;
394
395
VOL_swing
() {
396
lastgreeniter
=
lastyellowiter
=
lastrediter
= 0;
397
ngs
=
nrs
=
nys
= 0;
398
}
399
~VOL_swing
(){}
400
401
inline
void
cond
(
const
VOL_dual
& dual,
402
const
double
lcost,
const
double
ascent,
const
int
iter) {
403
double
eps = 1.e-3;
404
405
if
(ascent > 0.0 && lcost > dual.
lcost
+ eps) {
406
lastswing
=
green
;
407
lastgreeniter
= iter;
408
++
ngs
;
409
rd
= 0;
410
}
else
{
411
if
(ascent <= 0 && lcost > dual.
lcost
) {
412
lastswing
=
yellow
;
413
lastyellowiter
= iter;
414
++
nys
;
415
rd
= 0;
416
}
else
{
417
lastswing
=
red
;
418
lastrediter
= iter;
419
++
nrs
;
420
rd
= 1;
421
}
422
}
423
}
424
425
inline
double
426
lfactor
(
const
VOL_parms
& parm,
const
double
lambda,
const
int
iter) {
427
double
lambdafactor = 1.0;
428
double
eps = 5.e-4;
429
int
cons;
430
431
switch
(
lastswing
) {
432
case
green
:
433
cons = iter -
VolMax
(
lastyellowiter
,
lastrediter
);
434
if
(parm.
printflag
& 4)
435
printf(
" G: Consecutive Gs = %3d\n\n"
, cons);
436
if
(cons >= parm.
greentestinvl
&& lambda < 2.0) {
437
lastgreeniter
=
lastyellowiter
=
lastrediter
= iter;
438
lambdafactor = 2.0;
439
if
(parm.
printflag
& 2)
440
printf(
"\n ---- increasing lamda to %g ----\n\n"
,
441
lambda * lambdafactor);
442
}
443
break
;
444
445
case
yellow
:
446
cons = iter -
VolMax
(
lastgreeniter
,
lastrediter
);
447
if
(parm.
printflag
& 4)
448
printf(
" Y: Consecutive Ys = %3d\n\n"
, cons);
449
if
(cons >= parm.
yellowtestinvl
) {
450
lastgreeniter
=
lastyellowiter
=
lastrediter
= iter;
451
lambdafactor = 1.1;
452
if
(parm.
printflag
& 2)
453
printf(
"\n **** increasing lamda to %g *****\n\n"
,
454
lambda * lambdafactor);
455
}
456
break
;
457
458
case
red
:
459
cons = iter -
VolMax
(
lastgreeniter
,
lastyellowiter
);
460
if
(parm.
printflag
& 4)
461
printf(
" R: Consecutive Rs = %3d\n\n"
, cons);
462
if
(cons >= parm.
redtestinvl
&& lambda > eps) {
463
lastgreeniter
=
lastyellowiter
=
lastrediter
= iter;
464
lambdafactor = 0.67;
465
if
(parm.
printflag
& 2)
466
printf(
"\n **** decreasing lamda to %g *****\n\n"
,
467
lambda * lambdafactor);
468
}
469
break
;
470
}
471
return
lambdafactor;
472
}
473
474
inline
void
475
print
() {
476
printf(
"**** G= %i, Y= %i, R= %i ****\n"
,
ngs
,
nys
,
nrs
);
477
ngs
=
nrs
=
nys
= 0;
478
}
479
};
480
481
//############################################################################
482
/* alpha should be decreased if after some number of iterations the objective
483
has increased less that 1% */
484
class
VOL_alpha_factor
{
485
private
:
486
VOL_alpha_factor
(
const
VOL_alpha_factor
&);
487
VOL_alpha_factor
&
operator=
(
const
VOL_alpha_factor
&);
488
public
:
489
double
lastvalue
;
490
491
VOL_alpha_factor
() {
lastvalue
= -DBL_MAX;}
492
~VOL_alpha_factor
() {}
493
494
inline
double
factor
(
const
VOL_parms
& parm,
const
double
lcost,
495
const
double
alpha) {
496
if
(alpha < parm.
alphamin
)
497
return
1.0;
498
register
const
double
ll =
VolAbs
(lcost);
499
const
double
x = ll > 10 ? (lcost-
lastvalue
)/ll : (lcost-
lastvalue
);
500
lastvalue
= lcost;
501
return
(x <= 0.01) ? parm.
alphafactor
: 1.0;
502
}
503
};
504
505
//############################################################################
506
/* here we compute the norm of the conjugate direction -hh-, the norm of the
507
subgradient -norm-, the inner product between the subgradient and the
508
last conjugate direction -vh-, and the inner product between the new
509
conjugate direction and the subgradient */
510
class
VOL_vh
{
511
private
:
512
VOL_vh
(
const
VOL_vh
&);
513
VOL_vh
&
operator=
(
const
VOL_vh
&);
514
public
:
515
double
hh
;
516
double
norm
;
517
double
vh
;
518
double
asc
;
519
520
VOL_vh
(
const
double
alpha,
521
const
VOL_dvector
& dual_lb,
const
VOL_dvector
& dual_ub,
522
const
VOL_dvector
& v,
const
VOL_dvector
& vstar,
523
const
VOL_dvector
& u);
524
~VOL_vh
(){}
525
};
526
527
//############################################################################
528
/* here we compute different parameter to be printed. v2 is the square of
529
the norm of the subgradient. vu is the inner product between the dual
530
variables and the subgradient. vabs is the maximum absolute value of
531
the violations of pstar. asc is the inner product between the conjugate
532
direction and the subgradient */
533
class
VOL_indc
{
534
private
:
535
VOL_indc
(
const
VOL_indc
&);
536
VOL_indc
&
operator=
(
const
VOL_indc
&);
537
public
:
538
double
v2
;
539
double
vu
;
540
double
vabs
;
541
double
asc
;
542
543
public
:
544
VOL_indc
(
const
VOL_dvector
& dual_lb,
const
VOL_dvector
& dual_ub,
545
const
VOL_primal
& primal,
const
VOL_primal
& pstar,
546
const
VOL_dual
& dual);
547
~VOL_indc
() {}
548
};
549
550
//#############################################################################
551
558
class
VOL_user_hooks
{
559
public
:
560
virtual
~VOL_user_hooks
() {}
561
public
:
562
// for all hooks: return value of -1 means that volume should quit
567
virtual
int
compute_rc
(
const
VOL_dvector
& u,
VOL_dvector
& rc) = 0;
568
577
virtual
int
solve_subproblem
(
const
VOL_dvector
& dual,
const
VOL_dvector
& rc,
578
double
& lcost,
VOL_dvector
& x,
VOL_dvector
& v,
579
double
& pcost) = 0;
586
virtual
int
heuristics
(
const
VOL_problem
& p,
587
const
VOL_dvector
& x,
double
& heur_val) = 0;
588
};
589
590
//#############################################################################
591
600
class
VOL_problem
{
601
private
:
602
VOL_problem
(
const
VOL_problem
&);
603
VOL_problem
&
operator=
(
const
VOL_problem
&);
604
void
set_default_parm
();
605
// ############ INPUT fields ########################
606
public
:
610
VOL_problem
();
613
VOL_problem
(
const
char
*filename);
615
~VOL_problem
();
617
623
int
solve
(
VOL_user_hooks
& hooks,
const
bool
use_preset_dual =
false
);
625
626
private
:
630
double
alpha_
;
632
double
lambda_
;
633
// This union is here for padding (so that data members would be
634
// double-aligned on x86 CPU
635
union
{
637
int
iter_
;
638
double
__pad0
;
639
};
641
642
public
:
643
647
double
value
;
649
VOL_dvector
dsol
;
651
VOL_dvector
psol
;
653
VOL_dvector
viol
;
655
659
VOL_parms
parm
;
661
int
psize
;
663
int
dsize
;
666
VOL_dvector
dual_lb
;
669
VOL_dvector
dual_ub
;
671
672
public
:
676
int
iter
()
const
{
return
iter_
; }
678
double
alpha
()
const
{
return
alpha_
; }
680
double
lambda
()
const
{
return
lambda_
; }
682
683
private
:
687
void
read_params
(
const
char
* filename);
688
690
int
initialize
(
const
bool
use_preset_dual);
691
693
void
print_info
(
const
int
iter
,
694
const
VOL_primal
& primal,
const
VOL_primal
& pstar,
695
const
VOL_dual
& dual);
696
699
double
readjust_target
(
const
double
oldtarget,
const
double
lcost)
const
;
700
708
double
power_heur
(
const
VOL_primal
& primal,
const
VOL_primal
& pstar,
709
const
VOL_dual
& dual)
const
;
711
};
712
713
#endif
Generated on Tue Mar 1 2016 22:28:24 by
1.8.4