1 /************ CRanLip - universal multivariate random ************************
2 * variate generatior based on acceptance/rejection
3 *
4 * begin : April 30 2004
5 * version : 1.0
6 * copyright : (C) 2004 by Gleb Beliakov
7 * email : gleb@deakin.edu.au
8 *
9 *
10 * implementation of CRanLip class
11 *
12 *
13 * Copyright (C) 2004 Gleb Beliakov (gleb@deakin.edu.au)
14 *
15 * This program is free software; you can redistribute it and/or modify
16 * it under the terms of the GNU General Public License as published by
17 * the Free Software Foundation; either version 2 of the License, or (at
18 * your option) any later version.
19 *
20 * This program is distributed in the hope that it will be useful, but
21 * WITHOUT ANY WARRANTY; without even the implied warranty of
22 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
23 * General Public License for more details.
24 *
25 * You should have received a copy of the GNU General Public License
26 * along with this program; if not, write to the Free Software
27 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
28 */
29
30 #include "ranlip.h"
31
32
33 extern "C" int TheSeed;
34 extern ranlux_state_t RLSTATE;
35
36 //////////////////////////////////////////////////////////////////////
37 // Construction/Destruction
38 //////////////////////////////////////////////////////////////////////
39
CRanLip()40 CRanLip::CRanLip()
41 {
42 Computed=0;
43 SetUniformGenerator(ranlux_get_double_V);
44
45 Probabilities=NULL;
46 Dist=NULL;
47 m_boundLeft=NULL;
48 m_boundRight=NULL;
49 m_tempLeft=NULL;
50 m_tempRight=NULL;
51 m_tempint=NULL;
52 h=hfine=NULL;
53 Lipschitz=0;
54 TheSeed=10;
55 }
56
~CRanLip()57 CRanLip::~CRanLip() {FreeMem();}
58
FreeMem()59 void CRanLip::FreeMem()
60 {
61 if(Dist !=NULL) gsl_ran_discrete_free(Dist);
62 if(Probabilities!=NULL) free(Probabilities);
63 if(m_boundLeft!=NULL) free(m_boundLeft);
64 if(m_boundRight!=NULL) free(m_boundRight);
65 if(m_tempLeft!=NULL) free(m_tempLeft);
66 if(m_tempRight!=NULL) free(m_tempRight);
67 if(m_tempint!=NULL) free(m_tempint);
68 if(m_tempintfine!=NULL) free(m_tempintfine);
69 if(h!=NULL) free(h);
70 if(hfine!=NULL) free(hfine);
71 if(V!=NULL) free(V);
72
73 Computed=0;
74 Probabilities=NULL;
75 Dist=NULL;
76 m_boundLeft=NULL;
77 m_boundRight=NULL;
78 m_tempLeft=NULL;
79 m_tempRight=NULL;
80 m_tempint=m_tempintfine=NULL;
81 h=hfine=V=NULL;
82 }
83
84
85
Distribution(double * p)86 double CRanLip::Distribution(double* p) // dummy method
87 { return 1; }
88
Seed(int seed)89 void CRanLip::Seed(int seed)
90 {
91 TheSeed=seed;
92 ranlux_set_seed(seed);
93 count_total= count_errors=0;
94 }
95
Init(int dim,double * left,double * right)96 void CRanLip::Init(int dim, double* left, double* right)
97 {
98 m_boundLeft=(double*) malloc(dim*sizeof(double));
99 m_boundRight=(double*) malloc(dim*sizeof(double));
100 m_tempLeft=(double*) malloc(dim*sizeof(double));
101 m_tempRight=(double*) malloc(dim*sizeof(double));
102 V =(double*) malloc(dim*sizeof(double));
103 h =(double*) malloc(dim*sizeof(double));
104 hfine =(double*) malloc(dim*sizeof(double));
105
106 m_tempint=(int*) malloc(dim*sizeof(int));
107 m_tempintfine=(int*) malloc(dim*sizeof(int));
108
109 Dimension=dim;
110 int i;
111 for(i=0;i<Dimension;i++)
112 {
113 m_boundLeft[i]=left[i];
114 m_boundRight[i]=right[i];
115 }
116 Computed=0;
117 }
118
RandomVecUniform(double * p)119 void CRanLip::RandomVecUniform(double* p)
120 {
121 size_t j = m_chosenElement = gsl_ran_discrete(Dist);
122
123 int i;
124 for(i =0;i<Dimension;i++) V[i] = UniformRNumber();
125
126 // calculate the position i,j,k,... in the array, starting at 0
127 GetIJK(j);
128
129 // m_tempint contains the indices
130 for(i =0;i<Dimension;i++) p[i] = m_boundLeft[i] + h[i]*m_tempint[i];
131 for(i =0;i<Dimension;i++) p[i] += V[i]*h[i];
132
133 count_total++;
134 }
135
RandomVec(double * p)136 void CRanLip::RandomVec(double* p)
137 {
138 while(Computed) { // must be 1 to generate random variates
139 RandomVecUniform(p);
140 double Y = UniformRNumber();
141 Y *= Probabilities[m_chosenElement];
142
143 if(Distribution(p) > Probabilities[m_chosenElement]) //return;
144 { count_errors++;
145 // printf("fail %f %f\n",Distribution(p),Probabilities[m_chosenElement]);
146 return;}
147
148 if(Y <= Distribution(p)) {/*count_success++;*/ return;}
149 }
150 }
151
152
PrepareHatFunction(int num,int numfine,double Lip)153 void CRanLip::PrepareHatFunction(int num, int numfine, double Lip)
154 {
155 int i,j;
156 double d;
157
158 num=max_(num,1);
159 numfine=max_(numfine-1,1);
160
161 for(bits=1;bits<32;bits++) if((1<<bits) >= numfine+1) break;
162 numfine=(1<<bits) -1;
163 mask1= (1<<(bits)) - 1;
164
165 Lipschitz=max_(Lip,1e-10);
166
167
168 TotalElements=1;
169 num_partition=num; num_small_partition=numfine; num_small_partition_p1=num_small_partition+1;
170 for(i=0;i<Dimension;i++) TotalElements*=num;
171
172 totvals=1;
173 for(i=0;i<Dimension;i++) totvals*=num_small_partition_p1;
174
175 Probabilities=(double*) malloc(TotalElements*sizeof(double));
176
177 // working arrays
178 LipschitzH=(double*)malloc(Dimension*sizeof(double));
179 m_delta=(int*)malloc(Dimension*sizeof(int));
180 vals=(double*)malloc(totvals*sizeof(double));
181
182
183 Volume=1;
184 m_delta[Dimension-1]=1;
185
186 for(i=0;i<Dimension;i++)
187 {
188 h[i]=(m_boundRight[i]-m_boundLeft[i])/num;
189 hfine[i]=h[i]/numfine;
190 Volume *= h[i];
191
192 LipschitzH[i]=0.5*Lipschitz*hfine[i] *2; //will divide by 2 later
193
194 if(i>0) m_delta[Dimension-i-1] = m_delta[Dimension-i] * num_small_partition_p1;
195 }
196
197 // make a special case if numfine==2, not to re-compute many values, store them in a table
198 // modify ComputeArray by using cache
199 if(numfine<=1) {
200 cache=(double*)malloc(TotalElements*sizeof(double));
201 for(j=0;j<TotalElements;j++) {
202 GetIJK(j);
203 for(i=0;i<Dimension;i++)
204 m_tempLeft[i] = m_boundLeft[i] + h[i]*m_tempint[i]; // left top corner
205 cache[j]=Distribution(m_tempLeft);
206 }
207
208 for(j=0;j<TotalElements;j++) {
209 ComputeArrayCache(j);
210 d=ComputeMaxBin();
211 Probabilities[j]=d;
212 }
213 free(cache);
214 } // special case
215 else { // general case
216 for(j=0;j<TotalElements;j++) {
217 GetIJK(j);
218 for(i=0;i<Dimension;i++) {
219 m_tempLeft[i] = m_boundLeft[i] + h[i]*m_tempint[i]; // left top corner
220 // m_tempRight[i]=m_tempLeft[i]+h[i]; // bottom right corner
221 }
222 ComputeArray();
223 d=ComputeMaxBin();
224 Probabilities[j]=d;
225 }
226 }
227
228 free(vals);
229 free(LipschitzH);
230 free(m_delta);
231
232 for(j=0;j<TotalElements;j++) Probabilities[j] *= Volume;
233
234 Dist = gsl_ran_discrete_preproc(TotalElements,Probabilities);
235
236 // now get back to function values in probabilities
237 for(j=0;j<TotalElements;j++) Probabilities[j] /= Volume;
238
239 count_total= count_errors=0;
240 Computed=1;
241 }
242
243
PrepareHatFunctionAuto(int num,int numfine,double minLip)244 void CRanLip::PrepareHatFunctionAuto(int num, int numfine, double minLip)
245 {
246 int i,j;
247
248 num=max_(num,1);
249 numfine=max_(numfine-1,1);
250
251 for(bits=1;bits<32;bits++) if((1<<bits) >= numfine+1) break;
252 numfine=(1<<bits) -1;
253 mask1= (1<<(bits)) - 1;
254
255
256 TotalElements=1;
257 num_partition=num; num_small_partition=numfine; num_small_partition_p1=num_small_partition+1;
258 for(i=0;i<Dimension;i++) TotalElements*=num;
259
260 totvals=1;
261 for(i=0;i<Dimension;i++) totvals*=num_small_partition_p1;
262
263 Probabilities=(double*) malloc(TotalElements*sizeof(double));
264
265 // working arrays
266 LipschitzH=(double*)malloc(Dimension*sizeof(double));
267 m_delta=(int*)malloc(Dimension*sizeof(int));
268 vals=(double*)malloc(totvals*sizeof(double));
269
270
271 Volume=1;
272 m_delta[Dimension-1]=1;
273
274 for(i=0;i<Dimension;i++)
275 {
276 h[i]=(m_boundRight[i]-m_boundLeft[i])/num;
277 hfine[i]=h[i]/numfine;
278 Volume *= h[i];
279
280 if(i>0) m_delta[Dimension-i-1] = m_delta[Dimension-i] * num_small_partition_p1;
281 }
282
283 Lipschitz=0;
284 double d;
285
286 // make a special case if numfine==2, not to re-compute many values, store them in a table
287 // modify ComputeArray by using cache
288 if(numfine<=1) {
289 cache=(double*)malloc(TotalElements*sizeof(double));
290 for(j=0;j<TotalElements;j++) {
291 GetIJK(j);
292 for(i=0;i<Dimension;i++)
293 m_tempLeft[i] = m_boundLeft[i] + h[i]*m_tempint[i]; // left top corner
294 cache[j]=Distribution(m_tempLeft);
295 }
296
297 for(j=0;j<TotalElements;j++) {
298 ComputeArrayCache(j);
299 d= Dimension* max_( ComputeLipschitzBin(), minLip); // an extra factor
300 // printf("%f\n",d);
301 if(d>Lipschitz) Lipschitz=d;
302 for(i=0;i<Dimension;i++) LipschitzH[i]=0.5 * d * hfine[i] *2 ; //will divide by 2 later
303
304 d=ComputeMaxBin();
305 Probabilities[j]=d;
306 }
307 free(cache);
308 } // special case
309 else { // general case
310
311 for(j=0;j<TotalElements;j++) {
312 GetIJK(j);
313 for(i=0;i<Dimension;i++) {
314 m_tempLeft[i] = m_boundLeft[i] + h[i]*m_tempint[i]; // left top corner
315 }
316 ComputeArray();
317 d =Dimension * max_(ComputeLipschitzBin(),minLip);
318 if(d>Lipschitz) Lipschitz=d;
319 for(i=0;i<Dimension;i++) LipschitzH[i]=0.5 * d * hfine[i] * 2; //will divide by 2 later
320
321 d=ComputeMaxBin();
322 Probabilities[j]=d;
323 }
324 }
325 Lipschitz=max_(Lipschitz, 1e-10);
326
327 free(vals);
328 free(LipschitzH);
329 free(m_delta);
330
331 for(j=0;j<TotalElements;j++) Probabilities[j] *= Volume;
332
333 Dist = gsl_ran_discrete_preproc(TotalElements,Probabilities);
334
335 // now get back to function values in probabilities
336 for(j=0;j<TotalElements;j++) Probabilities[j] /= Volume;
337
338 count_total= count_errors=0;
339 Computed=1;
340 }
341
ComputeArray()342 void CRanLip::ComputeArray()
343 {
344 int i,j,k,d1=Dimension-1;
345 double v=hfine[d1];
346 for(j=0;j<totvals;j++)
347 {
348 GetIJKfineBin(j);
349 for(i=0;i<Dimension;i++)
350 V[i] = m_tempLeft[i] + hfine[i]*m_tempintfine[i];
351
352 vals[j]=Distribution(V);
353 for(k=1;k<num_small_partition_p1;k++)
354 { //innermost index, no need for GetIJKfineBin
355 j++;
356 V[d1] += v;
357 vals[j]=Distribution(V);
358 }
359
360 }
361 }
362
363
ComputeArrayCache(int J)364 void CRanLip::ComputeArrayCache(int J)
365 { // use previously computed values
366 int j,i,k;
367 GetIJK(J);
368 vals[0]=cache[J];
369
370 for(i=0;i<Dimension;i++)
371 m_tempLeft[i] = m_boundLeft[i] + h[i]*m_tempint[i]; // left top corner
372
373 for(j=1;j<totvals;j++)
374 {
375 GetIJKfineBin(j);
376 k=GetIndexfromIJK(m_tempintfine);
377
378 if(k<TotalElements)
379 vals[j]=cache[k];
380 else { // was not computed
381 for(i=0;i<Dimension;i++)
382 V[i] = m_tempLeft[i] + hfine[i]*m_tempintfine[i];
383 vals[j]=Distribution(V);
384 }
385 }
386 }
387
GetIndexfromIJK(int * IJK)388 int CRanLip::GetIndexfromIJK( int* IJK)
389 {
390 // I have m_tempint and m_tempintfine=IJK, so that the needed element has indices
391 // m_tempint[i]+IJK[i]
392 int k,i,j,n;
393 k=0;
394 n=1;
395 for(i=Dimension-1;i>=0;i--) {
396 j=m_tempint[i]+IJK[i];
397 k += n * j;
398 if(j>=num_partition) return TotalElements+1; // outside the computed array
399 n *= num_partition;
400 }
401 return k;
402 }
403
404
ComputeMaxBin()405 double CRanLip::ComputeMaxBin()
406 {
407 double d=-10e20;
408 double u,v, D;
409 int i,j;
410 for(j=0;j<totvals;j++)
411 {
412 u=vals[j];
413 GetIJKfineBin(j);
414
415 for(i=0;i<Dimension;i++) {
416 if(m_tempintfine[i] < num_small_partition)
417 {
418 v=vals[j+m_delta[i]];
419 D=(u+v)+LipschitzH[i];
420 d=max_(d,D);
421 }
422 }
423 }
424
425 return d/2;
426 }
427
428
ComputeLipschitzBin()429 double CRanLip::ComputeLipschitzBin()
430 {
431 double d=-10e20;
432 double u,v;
433 int i,j;
434 for(j=0;j<totvals;j++)
435 {
436 u=vals[j];
437 GetIJKfineBin(j);
438
439 for(i=0;i<Dimension;i++) {
440 if(m_tempintfine[i] < num_small_partition)
441 {
442 v=vals[j+m_delta[i]];
443 v=fabs(u-v)/hfine[i];
444 d=max_(d,v);
445 }
446 }
447 }
448
449 return d;
450 }
451
GetIJK(int b)452 void CRanLip::GetIJK(int b)
453 {
454 int i;
455 div_t t;
456 for(i=1;i<Dimension;i++) {
457 t=div(b,num_partition);
458 m_tempint[Dimension-i]=t.rem;
459 b=t.quot;
460 }
461 m_tempint[0]=b;
462 }
463
464
GetIJKfineBin(int b)465 void CRanLip::GetIJKfineBin(int b)
466 {
467 int i;
468 for(i=Dimension-1;i>0;i--) {
469 m_tempintfine[i] = b & mask1;//rem;
470 b=b>>bits;//quot;
471 }
472 m_tempintfine[0]=b;
473 }
474
475
476
SavePartition(char * fname)477 void CRanLip::SavePartition(char * fname)
478 {
479 if(!Computed) return;
480
481 int j;
482 FILE * fp=fopen(fname,"w");
483 fprintf(fp,"Dim,Elements,Volume %d %d %d %f\n",Dimension,TotalElements,num_partition, Volume);
484
485 for(j=0;j<Dimension;j++) {
486 fprintf(fp,"%f %f\n",m_boundLeft[j],m_boundRight[j]);
487 }
488
489 for(j=0;j<TotalElements;j++) {
490 fprintf(fp,"%f\n",Probabilities[j]);
491 }
492 fclose(fp);
493 }
494
LoadPartition(char * fname)495 void CRanLip::LoadPartition(char * fname)
496 {
497 FreeMem();
498
499
500 int i,j;
501 FILE * fp=fopen(fname,"r");
502 fscanf(fp,"Dim,Elements,Volume %d %d %d %lf\n",&Dimension,&TotalElements, &num_partition, &Volume);
503
504 double* t_boundLeft=(double*) malloc(Dimension*sizeof(double));
505 double* t_boundRight=(double*) malloc(Dimension*sizeof(double));
506
507 for(j=0;j<Dimension;j++) {
508 fscanf(fp,"%lf %lf\n",&t_boundLeft[j], &t_boundRight[j]);
509 }
510
511 Init(Dimension,t_boundLeft,t_boundRight);
512 free(t_boundLeft);
513 free(t_boundRight);
514
515 double r;
516 Probabilities=(double*) calloc(TotalElements,sizeof(double));
517 if(Probabilities==NULL) return;
518
519 for(j=0;j<TotalElements;j++) {
520 fscanf(fp,"%lf\n",&r);
521 Probabilities[j] = r;
522 }
523 fclose(fp);
524
525 for(i=0;i<Dimension;i++)
526 {
527 h[i]=(m_boundRight[i]-m_boundLeft[i])/num_partition;
528 }
529
530 // for the discrete generator
531 Dist = gsl_ran_discrete_preproc(TotalElements,Probabilities);
532
533 count_total= count_errors=0;
534 Computed=1;
535 }
536
537
538
539