1 /*
2 * Name: rtr.c
3 * Author: Pietro Belotti
4 * Purpose: find a maximum feasible subsystem of an infeasible LP
5 *
6 * This code is published under the Eclipse Public License (EPL).
7 * See http://www.eclipse.org/legal/epl-v10.html
8 */
9
10 #include <unistd.h>
11 #include <stdlib.h>
12 #include <math.h>
13 #include <stdio.h>
14 #include <string.h>
15
16 #include "sparse.h"
17 #include "init.h"
18 #include "rtr.h"
19 #include "move.h"
20 #include "locsrch.h"
21 #include "linopt.h"
22 #include "chooseblock.h"
23 #include "isfeas.h"
24
25 #define MAX_NIMPROV 2000
26
27 #define MIN_TEMP 0.0001
28 #define MAX_TEMP 1e-4
29
30 #define INIT_MOMENTUM 40
31
32 extern int interrupt;
33
34 /*
35 * Randomized Thermal Relaxation algorithm
36 */
37
rtr(sparseLP * lp,char * sol)38 int rtr (sparseLP *lp, /* description of infeasible LP */
39 char *sol /* contains feasible subsystem */
40 ) {
41
42 char lincool = lp -> lincool,
43 invcool = lp -> invcool;
44
45 int maxIter = lp -> nIter;
46
47 double alpha = lp -> alpha,
48 beta = lp -> beta,
49 gammaRate = lp -> gammaRate,
50 muRate = lp -> muRate;
51
52 int nSatd = 0,
53 nSatd_loc = 0,
54 best_nSatd = 0,
55 nImprov = 0,
56 nIter;
57
58 int dnLoc, /* improvement with variable local search */
59 dnAms; /* improvement with rtr */
60
61 int rtrMomentum = INIT_MOMENTUM;
62
63 char whichmove = USE_RTR;
64
65 double *x; /* current solution */
66 double *b_Ax; /* current ineq. violation (>0 if violating) */
67 char *satd; /* satd [i] = 1 if i-th ineq. is satisfied, 0 otherwise */
68 int *block; /* set of inequalities used for computing x (k+1) */
69
70 double temperature = 0, gamma = 1, sumViol, sumViol_loc, mu = 1;
71
72 double time = 0;
73
74 if (lp -> my_id == 0) time = CoinCpuTime();
75
76 if (!(x =(double*) malloc ((lp->c0) * sizeof(double)))) fprintf (stderr,"Error: malloc failed.\n");
77 if (!(b_Ax =(double*) malloc ((lp->rk) * sizeof(double)))) fprintf (stderr,"Error: malloc failed.\n");
78 if (!(satd =(char *) malloc ((lp->rk) * sizeof(char)))) fprintf (stderr,"Error: malloc failed.\n");
79 if (!(block=(int *) malloc ((lp->rk+1) * sizeof(int)))) fprintf (stderr,"Error: malloc failed.\n");
80
81 /*
82 * Random initialization of all variables
83 */
84
85 if (lp -> my_id == 0) {
86
87 printf (" #iter |mfs| Temp. time |bl.| %%covg.\n");
88 printf ("==========================================================\n");
89 }
90
91 for (nIter = 0;
92
93 ((lp->timelimit < 0) || (CoinCpuTime() - time <= lp -> timelimit)) &&
94 (nIter < maxIter) &&
95 !interrupt &&
96 (nSatd < lp->r0);
97
98 nIter++) {
99
100 /*
101 * (Re)start from randomized solution
102 */
103
104 if ((!(nIter % (lp->restFreq)) || /* time to restart */
105 (nImprov > MAX_NIMPROV) || /* too many non-improving iterations */
106 (temperature*gamma<MIN_TEMP)) /* temperature too low */
107 && (nSatd <= best_nSatd)) { /* aspiration criteria: do not restart */
108 /* if good solution found */
109
110 nImprov = 0;
111
112 /* if ((anneal *= 2.0) > 0.4) anneal = 0.4; */
113
114 init_x (lp, x);
115
116 /*
117 * satd & b_Ax initialized using current x values.
118 */
119
120 nSatd_loc = init_sat (lp, satd, b_Ax, x, &sumViol_loc);
121
122 #ifdef RTR_MPI
123
124 //
125 // local sum of violation -- to be dealt with
126 //
127
128 // MPI_Reduce (&sumViol_loc, &sum_Viol, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
129 sumViol = sumViol_loc;
130 #else
131 sumViol = sumViol_loc;
132 #endif
133
134 /*
135 * TO DO: Should alpha be computed locally?
136 * Should we use sumViol_loc instead?
137 */
138
139 /*
140 * CAUTION! mu decreases in non improving iterations, but if
141 * other rules are used the alpha needs be recomputed globally
142 */
143
144 temperature = alpha * sumViol / (lp -> rk - nSatd_loc);
145
146 if (temperature > MAX_TEMP)
147 temperature = MAX_TEMP;
148
149 gamma = 1;
150 mu = 0.3;
151 }
152
153 /* if ((maxIter<1000) || !(nIter % (maxIter/1000))) */
154 /* fprintf (stderr, "\r%.1f%%", (nIter+1) * 100.0 / maxIter); */
155
156 #ifdef RTR_MPI
157 MPI_Allreduce (&nSatd_loc, &nSatd, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
158 #else
159 nSatd = nSatd_loc;
160 #endif
161
162 /*
163 * Found a better point?
164 */
165
166 if (nSatd > best_nSatd) {
167
168 int i, feas = 1;
169
170 int feask = (isFeas (lp, satd, x, &nSatd_loc)) == 3;
171
172 #ifdef RTR_MPI
173 MPI_Reduce (&feask, &feas, 1, MPI_INT, MPI_LAND, 0, MPI_COMM_WORLD);
174 #else
175 feas = feask;
176 #endif
177
178 if (lp -> my_id == 0) {
179
180 printf ("%8d%c %9d %11.2f %8.2f %9d %7.3f %c\n",
181 nIter, ((whichmove == USE_LOCSRCH)? '+' : ' '),
182 nSatd, temperature * gamma,
183 CoinCpuTime () - time,
184 mymax (1, (int) (mu * (lp->r0 - nSatd))),
185 (int) (1000.0 * (100.0 * nSatd / lp -> r0)) / 1000.0,
186 feas? '*' : '!');
187
188 fflush (stdout);
189
190 /*
191 * test actual feasibility of the subsystem found
192 */
193
194 /* fprintf (stderr, "\n\r%.1f%%", (nIter+1)*100.0/maxIter); fflush (stderr); */
195 }
196
197 best_nSatd = nSatd;
198
199 /*
200 * save current solution
201 */
202
203 /* #pragma disjoint */
204
205 for (i=lp->rk-1; i>=0; --i) *sol++ = *satd++;
206 satd -= lp->rk;
207 sol -= lp->rk;
208
209 nImprov = 0;
210
211 /*
212 * this is a good point, increase block percentage
213 */
214
215 if ((mu *= 1.1) > 1) mu = 1;
216
217 rtrMomentum = INIT_MOMENTUM;
218 }
219
220 else {
221 /* no improvement: */
222 if (rtrMomentum > 0) --rtrMomentum; /* - decrease importance of rtr vs. variable local search */
223 mu /= muRate; /* - shrink block fraction */
224 nImprov++; /* - increase no. non-improving iterations */
225 }
226
227 /*
228 * Main loop
229 */
230
231 if (whichmove == USE_RTR) { /* primal move (pure rtr) */
232
233 /*
234 * select a set of inequalities to be used in computing x(k+1).
235 */
236
237 choose_block (lp, block, satd, nSatd_loc,
238 mymax (1, (int) (mu * (lp -> rk - nSatd_loc))),
239 b_Ax, sumViol);
240
241 /*
242 * calc. x(k+1), update b_Ax, satd
243 */
244
245 nSatd_loc += (dnAms = move (lp, x, b_Ax, block, satd, temperature * gamma, &sumViol));
246
247 if ((dnAms <= 0) && (lp -> locsea) && (rtrMomentum == 0)) {
248
249 /*
250 * too many non-improving iterations with rtr, next iteration
251 * try variable local search
252 */
253
254 whichmove = USE_LOCSRCH;
255 rtrMomentum = INIT_MOMENTUM;
256 }
257 }
258
259 else { /* dual move: variable local search */
260
261 nSatd += (dnLoc = locsrch (lp, x, b_Ax, satd, &sumViol));
262
263 if (dnLoc <= 0) whichmove = USE_RTR; /* no improvement with variable local search */
264 /* choose rtr as next method */
265 }
266
267 /*
268 * End of main loop
269 */
270
271 if (lincool) gamma -= 1.0 / maxIter; /* update gamma linearly (temperature scaling factor) */
272 else if (invcool) gamma *= (nIter+1) / (nIter+2); /* update gamma inverse-linearly */
273 else gamma /= gammaRate; /* geometrically */
274
275 /* update temperature: */
276 /* convex combination (beta) of previous temperature and */
277 /* current total violation averaged on violated constraints */
278
279 if (nSatd_loc < lp -> rk)
280 temperature = beta * temperature +
281 (1 - beta) * alpha * sumViol / (lp -> rk - nSatd_loc);
282 }
283
284 if (lp -> my_id == 0) {
285
286 if (interrupt) printf ("User interrupt\n");
287 else if (nSatd >= lp->r0) printf ("All ineqs satisfied\n");
288 else printf ("Completed\n");
289 }
290
291 if (lp -> my_id == 0) printf ("Total time: %.2f\n", CoinCpuTime() - time);
292
293 /*
294 * check actual feasibility
295 */
296
297 switch (isFeas (lp, satd, x, &nSatd_loc)) { /* 3 indicates the mfs is actually feasible */
298
299 case -1: fprintf (stderr, "\rBounds Violated\n"); break;
300 case 0: fprintf (stderr, "\rError: mfs with no sense\n"); break;
301 case 1: fprintf (stderr, "\rWarning: mfs includes unsatisfied ineqs\n"); break;
302 case 2: fprintf (stderr, "\rWarning: mfs does not include satisfied ineqs\n"); break;
303 case 3: if (nSatd == lp->r0) fprintf (stderr, "\rProblem Feasible\n"); break;
304 default: fprintf (stderr, "\rCase out of bounds\n");
305 }
306
307 /*
308 * free static malloc'ed pointers in subroutines
309 */
310
311 move (NULL, NULL, NULL, NULL, NULL, 0, NULL);
312 locsrch (NULL, NULL, NULL, NULL, NULL);
313
314 free (b_Ax);
315 free (satd);
316 free (block);
317 free (x);
318
319 return best_nSatd;
320 }
321