1 /*
2 * Name: chooseblock.c
3 * Author: Pietro Belotti
4 * Purpose: select subset of violated inequalities for the update direction
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 <unistd.h>
12 #include <stdlib.h>
13 #include <math.h>
14 #include <stdio.h>
15 #include <string.h>
16
17 #ifdef RTR_MPI
18 #include <mpi.h>
19 #endif
20
21 #include "sparse.h"
22 #include "init.h"
23 #include "rtr.h"
24 #include "move.h"
25 #include "locsrch.h"
26 #include "linopt.h"
27
28 /*
29 * block choice routine
30 */
31
choose_block(sparseLP * lp,int * block,char * sat,int nSat,int cardinality,double * b_Ax,double sumviol)32 void choose_block (sparseLP *lp /* LP data */
33 ,int *block /* arrays with row indices to use (last is -1) */
34 ,char *sat /* fulfilled constraints */
35 ,int nSat /* number of fulfilled constraints */
36 ,int cardinality /* cardinality of the block */
37 ,double *b_Ax /* violation vector */
38 ,double sumviol /* sum of all violations */
39 ) {
40
41 register int j=0, k;
42
43 const int r = lp->rk;
44 char* chosen = lp->chosen;
45
46 /*
47 double cum_weight = 0.0;
48 double* cum_weights = lp->cum_weights;
49 */
50
51 if (r == nSat) return;
52
53 if (!(lp->dblrand))
54
55 while (j<cardinality) {
56
57 if (!(sat [k = floor (drand48() * r)]) &&
58 !(chosen [k])) {
59
60 chosen [k] = 1;
61 block [j++] = k;
62 }
63 }
64 else {
65
66 double relaxed_divider = sumviol;
67 int cnt = 0;
68
69 while (j<cardinality) {
70
71 if (!(sat [k = floor (drand48() * r)]) &&
72 !(chosen [k]) &&
73 (drand48 () < exp (- b_Ax [k] / relaxed_divider))) {
74
75 chosen [k] = 1;
76 block [j++] = k;
77 }
78
79 if (cnt++ > r) {
80
81 relaxed_divider *= 10;
82 cnt = 0;
83 }
84 }
85 }
86
87 /*
88
89 for (i = 0, j = 0; i < r; ++i) {
90 if (!sat[i]) {
91 cum_weight += exp(-b_Ax[i]/sumviol);
92 }
93 cum_weights[i] = cum_weight;
94 }
95
96 for (i = 5 * cardinality - 1, j = 0; j < cardinality && i >= 0; --i) {
97
98 const double key = drand48() * cum_weight;
99 const int select = dsearch (key, cum_weights, r);
100 if (chosen[select] == 0) {
101
102 block[j] = select;
103 chosen[select] = 1;
104 ++j;
105 }
106 }
107 */
108
109 block [j--] = -1;
110
111 for (; j >= 0; --j) chosen [block [j]] = 0;
112 }
113
114
115 /*
116
117 inline int cum_weight_comp2(const void* w0, const void* w1)
118 {
119 const double diff = (*(double*)w0) - (*(double*)w1);
120 return diff < 0 ? -1 : ( diff > 0 ? 1 : 0);
121 }
122
123 */
124
125 /*
126 * TO DO: make it faster!
127 */
128
129 /*
130 static inline int lookup (const double key,
131 const double *v,
132 const int length) {
133
134 register int i;
135
136 for (i=0; i<length; i++) if (v [i] >= key) return i;
137
138 return length;
139 }
140 */
141
142 /*
143 * interval binary search
144 */
145
146 /*
147 inline int dsearch (double key, register double *v, int length) {
148
149 register double result = 0;
150 register int cursor;
151
152 cursor = length /= 2;
153
154 while (length != 0) {
155
156 result = key - v [cursor];
157
158 if (result > EPSILON) cursor += length/2;
159 else if (result < - EPSILON) cursor -= length/2;
160 else return cursor;
161
162 if ((length>1) && (length % 2)) length++;
163 length /= 2;
164 }
165
166 if (result < - EPSILON) return cursor - 1;
167 else return cursor;
168 }
169 */
170
171 /*
172 * block choice routine (vector linear)
173 */
174
175 /*
176 void choose_block2 (sparseLP *lp, // LP data
177 int *block, // arrays with row indices to use (last is -1)
178 char *sat, // fulfilled constraints
179 int nSat, // number of fulfilled constraints
180 int cardinality, // cardinality of the block
181 double *b_Ax, // violation vector
182 double sumviol // sum of all violations
183 ) {
184
185 register int i, j, rk = lp -> rk;
186
187 double *random = NULL;
188 double cum_rand = 0;
189
190 double *cum_weights = lp -> cum_weights;
191 double cum_weight = 0;
192
193 char *violated = lp -> chosen;
194
195 //
196 // generate cuWei
197 //
198
199 b_Ax += rk;
200
201 for (j=0, i=rk; i>0; --i, --b_Ax)
202
203 if (!(*sat++)) {
204 cum_weight += exp (- *b_Ax / sumviol);
205 *cum_weights++ = cum_weight;
206 *violated++ = i;
207 ++j;
208 }
209
210 cum_weights -= j;
211 violated -= j;
212 sat -= rk;
213
214 //
215 // generate y = a set of doubles (cumulatively)
216 //
217
218 for (i=cardinality; i>0; --i) {
219
220 cum_rand += drand48 ();
221 *random++ = cum_rand;
222 }
223
224 random -= cardinality;
225
226 //
227 // scale them down to cuWei [n]
228 //
229
230 for (i=cardinality; i>0; --i) {
231
232 *random++ /= cum_rand;
233 }
234
235 random -= cardinality;
236
237 // repeat 10 times
238 // for i=1..card
239 // linear search until right cuWei [k] found -> output k, card--
240
241
242 int i, j;
243 double cum_weight = 0.0;
244 const int r = lp->rk;
245 const int nViol = r - nSat;
246
247 if (nViol == 0)
248 return;
249
250 for (i = 0, j = 0; i < r; ++i) {
251 if (!sat[i]) {
252 cum_weight += exp(-b_Ax[i]/sumviol);
253 }
254 cum_weights[i] = cum_weight;
255 }
256
257 for (i = 5 * cardinality - 1, j = 0; j < cardinality && i >= 0; --i) {
258
259 const double key = drand48() * cum_weight;
260 const int select = dsearch (key, cum_weights, r);
261 if (chosen[select] == 0) {
262
263 block[j] = select;
264 chosen[select] = 1;
265 ++j;
266 }
267 }
268
269 block [j--] = -1;
270
271 for (; j >= 0; --j) {
272 chosen[block[j]] = 0;
273 }
274 }
275 */
276