1 /*
2 * Name: locsrch.c
3 * Author: Pietro Belotti
4 * Purpose: search for best point by fixing all variables except one
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
11 #include <stdlib.h>
12 #include <stdio.h>
13
14 #include "sparse.h"
15 #include "linopt.h"
16 #include "locsrch.h"
17
18 #define EPS 1e-6
19 #define ALPHA 0.1
20
compare_frontier(const void * one,const void * two)21 int compare_frontier (const void *one, const void *two) {
22
23 if (((const frontier *) one) -> thres < ((const frontier *) two) -> thres) return -1;
24 if (((const frontier *) one) -> thres > ((const frontier *) two) -> thres) return 1;
25 return 0;
26 }
27
compare_dvar(const void * one,const void * two)28 int compare_dvar (const void *one, const void *two) {
29
30 if (((const dvar *) one) -> nsi > ((const dvar *) two) -> nsi) return -1;
31 if (((const dvar *) one) -> nsi < ((const dvar *) two) -> nsi) return 1;
32 return 0;
33 }
34
locsrch2(sparseLP * lp,double * b_Ax)35 int locsrch2 (sparseLP *lp,
36 /*double *x, */
37 double *b_Ax
38 /*double *dx */
39 /* char *satd, */
40 /* double *sumViol */
41 ) {
42
43 /*
44 * for each var {
45 *
46 * 1) get lambdas, pluses, minuses;
47 *
48 * 2a) sort them according to [alpha * pluses + (1-alpha) *
49 * pluses/(pluses+minuses)] select J:|J|=lp->blkcard, first |J|
50 * in ordering
51 *
52 * or
53 *
54 * 2b) select randomly with probability dependent on [alpha *
55 * pluses + (1-alpha) * pluses/(pluses+minuses)] or uniformly
56 *
57 * 3) call one_opt on convex combination
58 *
59 * }
60 */
61
62
63 int r;
64 int c;
65
66 int i, j;
67
68 static double *impr = NULL;
69 static double *avg = NULL;
70
71 double **vc = NULL;
72 int **vp = NULL;
73 int *vl = NULL;
74
75 if (!lp) {
76
77 if (impr) free (impr);
78 if (avg) free (avg);
79
80 return 0;
81 }
82
83 r = lp -> rk;
84 c = lp -> c0;
85
86 if (!impr) {
87
88 impr = (double *) malloc (c * sizeof (double));
89 avg = (double *) malloc (c * sizeof (double));
90 }
91
92 /*
93 * Roughly check variables that can have good improvement
94 */
95
96 for (i=0; i<c; i++) {
97
98 int up=0, down=0, bad=0;
99
100 register int *pp = *vp++;
101 register double *pc = *vc++;
102
103 int l = *vl++;
104 int maximp;
105
106 for (j=l; j>0; j--, pp++, pc++) {
107
108 register double viol = b_Ax [*pp];
109
110 if (viol < - EPSILON) { /* constraint fulfilled */
111 if (*pc > EPSILON) bad++; /* would violate it */
112 else {if (*pc < - EPSILON) down++;}}
113 else if (viol > EPSILON) {
114 if (*pc < - EPSILON) bad++;
115 else {if (*pc > EPSILON) up++;}}
116 }
117
118 maximp = mymax (up, down);
119
120 /*
121 * proportion is important, but the total is more
122 */
123
124 *impr++ = ALPHA * (maximp / l) + (1-ALPHA) * (maximp / r);
125
126 /*
127
128 intersect at
129
130 lb < x+d < ub
131
132 e_j > 0
133 b - A (x + d e_j) = b - Ax - d A e_j
134 b-Ax >0
135 b-Ax - d a_j <0 -> a_j > 0
136
137 b-Ax <0
138 b-Ax - d a_j >0 -> a_j < 0
139
140 up++
141
142 lb < x + d e_j < ub
143
144 (lb - x) < d < (ub - x)
145
146 e_j < 0
147 b - A (x + d e_j) = b - Ax - d A e_j
148
149 */
150
151 }
152
153 impr -= c;
154
155 vc -= c;
156 vp -= c;
157
158 #ifdef RTR_MPI
159 // MPI_Allreduce (impr, imprs, k, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
160 #endif
161 /*
162 for (i=0; i<c; i++) {
163 if (drand48 () < *impr++) *dx++ =
164 }
165 */
166
167 return 0;
168 }
169
locsrch(sparseLP * lp,double * x,double * b_Ax,char * satd,double * sumViol)170 int locsrch (sparseLP *lp,
171 double *x,
172 double *b_Ax,
173 char *satd,
174 double *sumViol
175 ) {
176
177 static frontier *cuts = NULL; /* threshold */
178 static dvar *deltas = NULL; /* delta vars */
179 static double *dx = NULL; /* delta x */
180
181 dvar *deltas_save;
182
183 int i,j,k, max, argmax, ncuts, cursat, ndelta;
184
185 double th, Delta;
186
187 double *coe;
188 int *pos;
189
190 double **vc;
191 int **vp;
192
193 int *vl;
194 double *lb;
195 double *ub;
196
197 int r, c;
198
199 double *rhs;
200
201 double *pV;
202 char *pS;
203
204 double step;
205 int dSatd;
206
207 if (!lp) {
208 if (cuts) {
209 free (deltas);
210 free (cuts);
211 free (dx);
212 }
213 return 0;
214 }
215
216 if (!cuts) {
217
218 /*
219 * non-skeleton implementation, rk or r0?
220 */
221
222 deltas = (dvar *) malloc ( (lp -> c0) * sizeof (dvar));
223 cuts = (frontier *) malloc (((lp -> rk)+1) * sizeof (frontier));
224 dx = (double *) malloc ( (lp -> c0) * sizeof (double));
225 }
226
227 vc = lp -> vc;
228 vp = lp -> vp;
229
230 vl = lp -> vl;
231 lb = lp -> lb;
232 ub = lp -> ub;
233
234 r = lp -> rk;
235 c = lp -> c0;
236
237 rhs = lp -> rhs;
238
239 for (i=ndelta=0; i<c; i++, x++, vl++, vp++, vc++, ub++, lb++) if (*vl) {
240
241 coe = *vc;
242 pos = *vp;
243
244 /* printf ("x%d: ", i);*/
245
246 for (ncuts=k=cursat=0, j=*vl; j>0; j--, pos++, coe++) {
247
248 if (((th = b_Ax [*pos] / *coe + *x) > *lb + EPS) && (th < *ub - EPS)) {
249
250 cuts -> thres = th;
251
252 if (*coe > EPS) {cuts -> sense = 'G'; if (*x > th + EPS) cursat ++;}
253 else {cuts -> sense = 'L'; k++; if (*x < th - EPS) cursat ++;}
254
255 /* printf ("(%d) %c=%.4f ", *pos, (cuts -> sense=='G')?'>':'<', th);*/
256 cuts++; ncuts++;
257 }
258 }
259
260 /* printf ("\n");*/
261
262 cuts -> thres = *ub;
263 cuts -> sense = 'L'; /* cut */
264
265 cuts -= ncuts;
266
267 qsort ((void *) cuts, (size_t) ncuts, (size_t) sizeof (frontier), compare_frontier);
268
269 max = k;
270 argmax = 0;
271
272 for (j=0; j<ncuts; j++, cuts++)
273
274 if (cuts -> sense == 'L') k--;
275 else if (++k>max) {max = k; argmax = j+1;}
276
277 cuts -= ncuts;
278
279 if (cursat < max) {
280
281 deltas -> nsi = max - cursat;
282
283 if (argmax==0) deltas -> delta = (*lb + cuts -> thres) / 2 - *x;
284 else deltas -> delta = (cuts [argmax].thres + cuts [argmax-1].thres) / 2 - *x;
285
286 deltas -> index = i;
287 deltas++; ndelta++;
288 }
289 }
290
291 deltas -= ndelta;
292
293 vl -= c;
294 lb -= c;
295 ub -= c;
296
297 x -= c;
298
299 if (!ndelta) return 0;
300
301 qsort ((void *) deltas, (size_t) ndelta, (size_t) sizeof (dvar), compare_dvar);
302
303 deltas_save = deltas;
304
305 for (j=0; (j<ndelta) && (deltas -> nsi); j++, deltas++) {
306 dx [deltas -> index] = deltas -> delta;
307 }
308
309 deltas = deltas_save;
310
311 /* Update satd, x, b_Ax */
312
313 if (lp -> blkcard == 1) {
314
315 step = deltas -> delta;
316 k = deltas -> index;
317
318 /*
319 printf ("x%d: %d (%.3f%+.3f)\n", deltas [0] . index,
320 deltas [0] . nsi, x [k],
321 deltas [0] . delta); fflush (stdout);
322 */
323
324 dSatd = 0;
325
326 pos = lp -> vp [k];
327 coe = lp -> vc [k];
328
329 if ((x [k] += step) > ub [k]) {x [k] = ub [k];}
330 else if ( x [k] < lb [k]) {x [k] = lb [k];}
331
332 for (j = vl [k]; j>0 ; j--) {
333
334 pV = b_Ax + *pos;
335 pS = satd + *pos++;
336 Delta = - step * *coe++;
337
338 if ((*pV) > 0) {
339
340 if ((Delta > 0) || ((*pV) > - Delta))
341 *sumViol += Delta;
342 else {
343 *sumViol -= *pV;
344 if (!*pS) { dSatd++; *pS = SATD; }
345 /* printf ("con %d now satd", *(pos-1));*/
346 }
347 }
348 else
349 if (Delta > - (*pV)) {
350 /* printf ("con %d now violated\n", *(pos-1));*/
351 *sumViol += (*pV + Delta);
352 if (*pS) { dSatd--; *pS = UNSATD; }
353 }
354 *pV += Delta;
355 }
356
357 return dSatd;
358 }
359 else {
360
361 for (j=0; j<c; j++) *dx++ = 0;
362
363 dx -= c;
364
365 deltas_save = deltas;
366
367 for (j=0; (j<ndelta) && (j<lp->blkcard) && (deltas -> nsi); j++, deltas++) {
368 /* printf ("[%5d]: %4d (%.2f)\n", deltas [j] . index, deltas [j] . nsi, deltas [j] . delta);
369 fflush (stdout);*/
370 dx [deltas -> index] = deltas -> nsi * deltas -> delta;
371 }
372
373 deltas = deltas_save;
374
375 /*
376 for (i=0; i<c;) {
377 printf ("%3d %.2f ", i, dx [i]);
378 if (!(++i % 20)) printf ("\n");
379 fflush (stdout);
380 }
381 printf ("\n");
382 */
383
384 if (0) { /* (lp->algorithm == RTR) { */
385 /*
386 dSatd = 0;
387 k = 0;
388
389 maxMod = 100.0;
390
391 for (i=0; i<c; i++, dx++, x++, ub++, lb++) if (*tv++) {
392
393 if (*dx < - EPSILON) {
394 if (*x < *lb + EPSILON) *dx = 0;
395 else if ((z = (*lb - *x) / *dx) < maxMod) maxMod = z;
396 }
397 else if (*dx > EPSILON) {
398 if (*x > *ub - EPSILON) *dx = 0;
399 else if ((z = (*ub - *x) / *dx) < maxMod) maxMod = z;
400 }
401 else *dx = 0;
402
403 // if (((*dx>0) && ((z = (*ub - *x) / *dx) < maxMod)) ||
404 // ((*dx<0) && ((z = (*lb - *x) / *dx) < maxMod))) maxMod = z;
405
406 }
407
408 tv -= c;
409
410 dx -= c;
411 x -= c;
412
413 ub -= c;
414 lb -= c;
415
416 if (maxMod <= 0.0) maxMod = 1.0;
417
418 for (i=0; i<c; i++, vl++, colC++, colP++, dx++, x++, ub++, lb++, tv++) {
419 if (*tv)
420 {
421 *tv = 0;
422
423 vp = *colP;
424 vc = *colC;
425
426 z = *x;
427
428 if ((*x += maxMod * *dx) > *ub) {*x = *ub;}
429 else if (*x < *lb) {*x = *lb;}
430
431 z -= *x;
432
433 for (j = *vl; j>0 ; j--) {
434
435 pV = b_Ax + *vp;
436 pS = sat + *vp++;
437 delta = z * *vc++;
438
439 if ((*pV) > 0)
440 if ((delta > 0) || ((*pV) > - delta))
441 *sum += delta;
442 else {
443 *sum -= *pV;
444 if (!*pS) { dSatd++; *pS = SATD; }
445 }
446 else
447 if (delta > - (*pV)) {
448 *sum += (*pV + delta);
449 if (*pS) { dSatd--; *pS = UNSATD; }
450 }
451 *pV += delta;
452 }
453 *dx = 0;
454 }
455 }
456 dx -= c;
457 tv -= c;
458
459 ub -= c;
460 lb -= c;
461 x -= c;
462 */
463 }
464 /* else return linOpt (lp, x, dx, b_Ax, satd, sumViol, 1/temp); */
465 }
466
467 return 0;
468 }
469