1 /*****************************************************************************
2 * *
3 * UNURAN -- Universal Non-Uniform Random number generator *
4 * *
5 *****************************************************************************
6 * *
7 * FILE: dgt.c *
8 * *
9 * TYPE: discrete univariate random variate *
10 * METHOD: guide table (indexed search) *
11 * *
12 * DESCRIPTION: *
13 * Given N discrete events with different probabilities P[k] *
14 * produce a value k consistent with its probability. *
15 * *
16 * REQUIRED: pointer to probability vector *
17 * *
18 *****************************************************************************
19 * *
20 * Copyright (c) 2000-2009 Wolfgang Hoermann and Josef Leydold *
21 * Department of Statistics and Mathematics, WU Wien, Austria *
22 * *
23 * This program is free software; you can redistribute it and/or modify *
24 * it under the terms of the GNU General Public License as published by *
25 * the Free Software Foundation; either version 2 of the License, or *
26 * (at your option) any later version. *
27 * *
28 * This program is distributed in the hope that it will be useful, *
29 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
30 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
31 * GNU General Public License for more details. *
32 * *
33 * You should have received a copy of the GNU General Public License *
34 * along with this program; if not, write to the *
35 * Free Software Foundation, Inc., *
36 * 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA *
37 * *
38 *****************************************************************************
39 * *
40 * REFERENCES: *
41 * [1] Chen, H. C. and Asau, Y. (1974): On generating random variates *
42 * from an empirical distribution, AIIE Trans. 6, pp. 163-166 *
43 * *
44 * SUMMARY: *
45 * [2] Devroye, L. (1986): Non-Uniform Random Variate Generation, New-York *
46 * *
47 *****************************************************************************
48 * *
49 * This method a varient of the inversion method, i.e. *
50 * *
51 * (1) Generate a random number U ~ U(0,1). *
52 * (2) Find largest integer I such that F(I) = P(X<=I) <= U. *
53 * *
54 * Step (2) is the crucial step. Using sequential search requires O(N) *
55 * comparisons. Indexed search however uses a guide table to jump to some *
56 * I' <= I near I. In a preprossing step (0,1) is partitioned into N *
57 * equal intervals and P(X<=I) <= (k-1)/N is solved for all k=1,...,N, *
58 * and the solutions are stored in a table (guide table). Setup this *
59 * table can be done in O(N). [2] has shown that the expected number of *
60 * of comparisons is at most 2. *
61 * *
62 * In the current implementation we use a variant that allows a different *
63 * size for the guide table. Bigger guide tables reduce the expected *
64 * number of comparisions but needs more memory and need more setup time. *
65 * On the other hand, the guide table can be made arbitrarily small to *
66 * save memory and setup time. Indeed, for size = 1 we have sequential *
67 * search again that requires no preprocessing. *
68 * *
69 *****************************************************************************
70 * *
71 * VARIANTS: *
72 * *
73 * We have three variants for the setup procedure: *
74 * *
75 * 1 ... compute (k-1)/N for each k=1,...,N *
76 * 2 ... compute (k-1)/N by summing up 1/N *
77 * 0 ... use 1 for large N and 2 for small N (default) *
78 * *
79 * variant 2 is faster but is more sensitive to roundoff errors when *
80 * N is large. *
81 * *
82 *****************************************************************************/
83
84 /*---------------------------------------------------------------------------*/
85
86 #include <unur_source.h>
87 #include <distr/distr.h>
88 #include <distr/distr_source.h>
89 #include <distr/discr.h>
90 #include <urng/urng.h>
91 #include "unur_methods_source.h"
92 #include "x_gen_source.h"
93 #include "dgt.h"
94 #include "dgt_struct.h"
95
96 /*---------------------------------------------------------------------------*/
97 /* Variants */
98
99 #define DGT_VARFLAG_DIV 0x01u /* compute guide table by division n/k */
100 #define DGT_VARFLAG_ADD 0x02u /* compute guide table by adding */
101
102 #define DGT_VAR_THRESHOLD 1000 /* above this value: use variant 1, else 2 */
103
104 /*---------------------------------------------------------------------------*/
105 /* Debugging flags */
106 /* bit 01 ... pameters and structure of generator (do not use here) */
107 /* bits 02-12 ... setup */
108 /* bits 13-24 ... adaptive steps */
109 /* bits 25-32 ... trace sampling */
110
111 #define DGT_DEBUG_REINIT 0x00000010u /* print parameters after reinit */
112 #define DGT_DEBUG_PRINTVECTOR 0x00000100u
113 #define DGT_DEBUG_TABLE 0x00000200u
114
115 /*---------------------------------------------------------------------------*/
116 /* Flags for logging set calls */
117
118 #define DGT_SET_GUIDEFACTOR 0x010u
119 #define DGT_SET_VARIANT 0x020u
120
121 /*---------------------------------------------------------------------------*/
122
123 #define GENTYPE "DGT" /* type of generator */
124
125 /*---------------------------------------------------------------------------*/
126
127 static struct unur_gen *_unur_dgt_init( struct unur_par *par );
128 /*---------------------------------------------------------------------------*/
129 /* Initialize new generator. */
130 /*---------------------------------------------------------------------------*/
131
132 static int _unur_dgt_reinit( struct unur_gen *gen );
133 /*---------------------------------------------------------------------------*/
134 /* Reinitialize generator. */
135 /*---------------------------------------------------------------------------*/
136
137 static struct unur_gen *_unur_dgt_create( struct unur_par *par );
138 /*---------------------------------------------------------------------------*/
139 /* create new (almost empty) generator object. */
140 /*---------------------------------------------------------------------------*/
141
142 static int _unur_dgt_check_par( struct unur_gen *gen );
143 /*---------------------------------------------------------------------------*/
144 /* Check parameters of given distribution and method */
145 /*---------------------------------------------------------------------------*/
146
147 static struct unur_gen *_unur_dgt_clone( const struct unur_gen *gen );
148 /*---------------------------------------------------------------------------*/
149 /* copy (clone) generator object. */
150 /*---------------------------------------------------------------------------*/
151
152 static void _unur_dgt_free( struct unur_gen *gen);
153 /*---------------------------------------------------------------------------*/
154 /* destroy generator object. */
155 /*---------------------------------------------------------------------------*/
156
157 static int _unur_dgt_sample( struct unur_gen *gen );
158 /*---------------------------------------------------------------------------*/
159 /* sample from generator */
160 /*---------------------------------------------------------------------------*/
161
162 static int _unur_dgt_create_tables( struct unur_gen *gen );
163 /*---------------------------------------------------------------------------*/
164 /* create (allocate) tables */
165 /*---------------------------------------------------------------------------*/
166
167 static int _unur_dgt_make_guidetable( struct unur_gen *gen );
168 /*---------------------------------------------------------------------------*/
169 /* create table for indexed search */
170 /*---------------------------------------------------------------------------*/
171
172 #ifdef UNUR_ENABLE_LOGGING
173 /*---------------------------------------------------------------------------*/
174 /* the following functions print debugging information on output stream, */
175 /* i.e., into the LOG file if not specified otherwise. */
176 /*---------------------------------------------------------------------------*/
177
178 static void _unur_dgt_debug_init( struct unur_gen *gen );
179 /*---------------------------------------------------------------------------*/
180 /* print after generator has been initialized has completed. */
181 /*---------------------------------------------------------------------------*/
182
183 static void _unur_dgt_debug_table( struct unur_gen *gen );
184 /*---------------------------------------------------------------------------*/
185 /* print data for guide table. */
186 /*---------------------------------------------------------------------------*/
187 #endif
188
189 #ifdef UNUR_ENABLE_INFO
190 static void _unur_dgt_info( struct unur_gen *gen, int help );
191 /*---------------------------------------------------------------------------*/
192 /* create info string. */
193 /*---------------------------------------------------------------------------*/
194 #endif
195
196 /*---------------------------------------------------------------------------*/
197 /* abbreviations */
198
199 #define DISTR_IN distr->data.discr /* data for distribution object */
200
201 #define PAR ((struct unur_dgt_par*)par->datap) /* data for parameter object */
202 #define GEN ((struct unur_dgt_gen*)gen->datap) /* data for generator object */
203 #define DISTR gen->distr->data.discr /* data for distribution in generator object */
204
205 #define SAMPLE gen->sample.discr /* pointer to sampling routine */
206
207 /*---------------------------------------------------------------------------*/
208
209 #define _unur_dgt_getSAMPLE(gen) (_unur_dgt_sample)
210
211 /*---------------------------------------------------------------------------*/
212
213 /*****************************************************************************/
214 /** Public: User Interface (API) **/
215 /*****************************************************************************/
216
217 struct unur_par *
unur_dgt_new(const struct unur_distr * distr)218 unur_dgt_new( const struct unur_distr *distr )
219 /*----------------------------------------------------------------------*/
220 /* get default parameters */
221 /* */
222 /* parameters: */
223 /* distr ... pointer to distribution object */
224 /* */
225 /* return: */
226 /* default parameters (pointer to structure) */
227 /* */
228 /* error: */
229 /* return NULL */
230 /*----------------------------------------------------------------------*/
231 {
232 struct unur_par *par;
233
234 /* check arguments */
235 _unur_check_NULL( GENTYPE,distr,NULL );
236
237 /* check distribution */
238 if (distr->type != UNUR_DISTR_DISCR) {
239 _unur_error(GENTYPE,UNUR_ERR_DISTR_INVALID,""); return NULL; }
240 COOKIE_CHECK(distr,CK_DISTR_DISCR,NULL);
241
242 if (DISTR_IN.pv == NULL) {
243 /* There is no PV try to compute it. */
244 if ( DISTR_IN.pmf
245 && ( (((unsigned)DISTR_IN.domain[1] - (unsigned)DISTR_IN.domain[0]) < UNUR_MAX_AUTO_PV)
246 || ( (distr->set & UNUR_DISTR_SET_PMFSUM) && DISTR_IN.domain[0] > INT_MIN ) ) ) {
247 /* However this requires a PMF and either a bounded domain */
248 /* or the sum over the PMF. */
249 _unur_warning(GENTYPE,UNUR_ERR_DISTR_REQUIRED,"PV. Try to compute it.");
250 }
251 else {
252 _unur_error(GENTYPE,UNUR_ERR_DISTR_REQUIRED,"PV"); return NULL;
253 }
254 }
255
256 /* allocate structure */
257 par = _unur_par_new( sizeof(struct unur_dgt_par) );
258 COOKIE_SET(par,CK_DGT_PAR);
259
260 /* copy input */
261 par->distr = distr; /* pointer to distribution object */
262
263 /* set default values */
264 PAR->guide_factor = 1.; /* use same size for guide table */
265
266 par->method = UNUR_METH_DGT; /* method */
267 par->variant = 0u; /* default variant */
268 par->set = 0u; /* inidicate default parameters */
269 par->urng = unur_get_default_urng(); /* use default urng */
270 par->urng_aux = NULL; /* no auxilliary URNG required */
271
272 par->debug = _unur_default_debugflag; /* set default debugging flags */
273
274 /* routine for starting generator */
275 par->init = _unur_dgt_init;
276
277 return par;
278
279 } /* end of unur_dgt_new() */
280
281 /*---------------------------------------------------------------------------*/
282
283 int
unur_dgt_set_variant(struct unur_par * par,unsigned variant)284 unur_dgt_set_variant( struct unur_par *par, unsigned variant )
285 /*----------------------------------------------------------------------*/
286 /* set variant of method */
287 /* */
288 /* parameters: */
289 /* par ... pointer to parameter for building generator object */
290 /* variant ... indicator for variant */
291 /* */
292 /* return: */
293 /* UNUR_SUCCESS ... on success */
294 /* error code ... on error */
295 /*----------------------------------------------------------------------*/
296 {
297 /* check arguments */
298 _unur_check_NULL( GENTYPE, par, UNUR_ERR_NULL );
299 _unur_check_par_object( par, DGT );
300
301 /* check new parameter for generator */
302 if (variant != DGT_VARFLAG_ADD && variant != DGT_VARFLAG_DIV) {
303 _unur_warning(GENTYPE,UNUR_ERR_PAR_VARIANT,"");
304 return UNUR_ERR_PAR_VARIANT;
305 }
306
307 /* changelog */
308 par->set |= DGT_SET_VARIANT;
309
310 par->variant = variant;
311
312 return UNUR_SUCCESS;
313 } /* end of unur_dgt_set_variant() */
314
315 /*---------------------------------------------------------------------------*/
316
317 int
unur_dgt_set_guidefactor(struct unur_par * par,double factor)318 unur_dgt_set_guidefactor( struct unur_par *par, double factor )
319 /*----------------------------------------------------------------------*/
320 /* set factor for relative size of guide table */
321 /* */
322 /* parameters: */
323 /* par ... pointer to parameter for building generator object */
324 /* factor ... relative size of table */
325 /* */
326 /* return: */
327 /* UNUR_SUCCESS ... on success */
328 /* error code ... on error */
329 /*----------------------------------------------------------------------*/
330 {
331 /* check arguments */
332 _unur_check_NULL( GENTYPE, par, UNUR_ERR_NULL );
333 _unur_check_par_object( par, DGT );
334
335 /* check new parameter for generator */
336 if (factor < 0) {
337 _unur_warning(GENTYPE,UNUR_ERR_PAR_SET,"relative table size < 0");
338 return UNUR_ERR_PAR_SET;
339 }
340
341 /* store date */
342 PAR->guide_factor = factor;
343
344 /* changelog */
345 par->set |= DGT_SET_GUIDEFACTOR;
346
347 return UNUR_SUCCESS;
348
349 } /* end of unur_dgt_set_guidefactor() */
350
351
352 /*****************************************************************************/
353 /** Private **/
354 /*****************************************************************************/
355
356 struct unur_gen *
_unur_dgt_init(struct unur_par * par)357 _unur_dgt_init( struct unur_par *par )
358 /*----------------------------------------------------------------------*/
359 /* initialize new generator */
360 /* */
361 /* parameters: */
362 /* par ... pointer to parameter for building generator object */
363 /* */
364 /* return: */
365 /* pointer to generator object */
366 /* */
367 /* error: */
368 /* return NULL */
369 /*----------------------------------------------------------------------*/
370 {
371 struct unur_gen *gen; /* pointer to generator object */
372
373 /* check arguments */
374 CHECK_NULL(par,NULL);
375
376 /* check input */
377 if ( par->method != UNUR_METH_DGT ) {
378 _unur_error(GENTYPE,UNUR_ERR_PAR_INVALID,"");
379 return NULL; }
380 COOKIE_CHECK(par,CK_DGT_PAR,NULL);
381
382 /* create a new empty generator object */
383 gen = _unur_dgt_create(par);
384 _unur_par_free(par);
385 if (!gen) return NULL;
386
387 /* check parameters */
388 if ( _unur_dgt_check_par(gen) != UNUR_SUCCESS ) {
389 _unur_dgt_free(gen); return NULL;
390 }
391
392 /* compute guide table */
393 if ( (_unur_dgt_create_tables(gen) != UNUR_SUCCESS) ||
394 (_unur_dgt_make_guidetable(gen) != UNUR_SUCCESS) ) {
395 _unur_dgt_free(gen); return NULL;
396 }
397
398 #ifdef UNUR_ENABLE_LOGGING
399 /* write info into LOG file */
400 if (gen->debug) _unur_dgt_debug_init(gen);
401 #endif
402
403 return gen;
404 } /* end of _unur_dgt_init() */
405
406 /*---------------------------------------------------------------------------*/
407
408 int
_unur_dgt_reinit(struct unur_gen * gen)409 _unur_dgt_reinit( struct unur_gen *gen )
410 /*----------------------------------------------------------------------*/
411 /* re-initialize (existing) generator. */
412 /* */
413 /* parameters: */
414 /* gen ... pointer to generator object */
415 /* */
416 /* return: */
417 /* UNUR_SUCCESS ... on success */
418 /* error code ... on error */
419 /*----------------------------------------------------------------------*/
420 {
421 int rcode;
422
423 /* check parameters */
424 if ( (rcode = _unur_dgt_check_par(gen)) != UNUR_SUCCESS)
425 return rcode;
426
427 /* compute table */
428 if ( ((rcode = _unur_dgt_create_tables(gen)) != UNUR_SUCCESS) ||
429 ((rcode = _unur_dgt_make_guidetable(gen)) != UNUR_SUCCESS) ) {
430 return rcode;
431 }
432
433 /* (re)set sampling routine */
434 SAMPLE = _unur_dgt_getSAMPLE(gen);
435
436 #ifdef UNUR_ENABLE_LOGGING
437 /* write info into LOG file */
438 if (gen->debug & DGT_DEBUG_REINIT) _unur_dgt_debug_init(gen);
439 #endif
440
441 return UNUR_SUCCESS;
442 } /* end of _unur_dgt_reinit() */
443
444 /*---------------------------------------------------------------------------*/
445
446 struct unur_gen *
_unur_dgt_create(struct unur_par * par)447 _unur_dgt_create( struct unur_par *par )
448 /*----------------------------------------------------------------------*/
449 /* allocate memory for generator */
450 /* */
451 /* parameters: */
452 /* par ... pointer to parameter for building generator object */
453 /* */
454 /* return: */
455 /* pointer to (empty) generator object with default settings */
456 /* */
457 /* error: */
458 /* return NULL */
459 /*----------------------------------------------------------------------*/
460 {
461 struct unur_gen *gen; /* pointer to generator object */
462
463 /* check arguments */
464 CHECK_NULL(par,NULL); COOKIE_CHECK(par,CK_DGT_PAR,NULL);
465
466 /* create new generic generator object */
467 gen = _unur_generic_create( par, sizeof(struct unur_dgt_gen) );
468
469 /* magic cookies */
470 COOKIE_SET(gen,CK_DGT_GEN);
471
472 /* set generator identifier */
473 gen->genid = _unur_set_genid(GENTYPE);
474
475 /* routines for sampling and destroying generator */
476 SAMPLE = _unur_dgt_getSAMPLE(gen);
477 gen->destroy = _unur_dgt_free;
478 gen->clone = _unur_dgt_clone;
479 gen->reinit = _unur_dgt_reinit;
480
481 /* copy some parameters into generator object */
482 GEN->guide_factor = PAR->guide_factor;
483
484 /* set all pointers to NULL */
485 GEN->cumpv = NULL;
486 GEN->guide_table = NULL;
487
488 #ifdef UNUR_ENABLE_INFO
489 /* set function for creating info string */
490 gen->info = _unur_dgt_info;
491 #endif
492
493 /* return pointer to (almost empty) generator object */
494 return gen;
495
496 } /* end of _unur_dgt_create() */
497
498 /*---------------------------------------------------------------------------*/
499
500 int
_unur_dgt_check_par(struct unur_gen * gen)501 _unur_dgt_check_par( struct unur_gen *gen )
502 /*----------------------------------------------------------------------*/
503 /* check parameters of given distribution and method */
504 /* */
505 /* parameters: */
506 /* gen ... pointer to generator object */
507 /* */
508 /* return: */
509 /* UNUR_SUCCESS ... on success */
510 /* error code ... on error */
511 /*----------------------------------------------------------------------*/
512 {
513 /* we need a PV */
514 if (DISTR.pv == NULL) {
515 /* try to compute PV */
516 if (unur_distr_discr_make_pv( gen->distr ) <= 0) {
517 /* not successful */
518 _unur_error(GENTYPE,UNUR_ERR_DISTR_REQUIRED,"PV");
519 return UNUR_ERR_DISTR_REQUIRED;
520 }
521 }
522
523 /* default variant? */
524 if (gen->variant == 0) /* default variant */
525 gen->variant = (DISTR.n_pv > DGT_VAR_THRESHOLD)
526 ? DGT_VARFLAG_DIV : DGT_VARFLAG_ADD;
527
528 return UNUR_SUCCESS;
529 } /* end of _unur_dgt_check_par() */
530
531 /*---------------------------------------------------------------------------*/
532
533 struct unur_gen *
_unur_dgt_clone(const struct unur_gen * gen)534 _unur_dgt_clone( const struct unur_gen *gen )
535 /*----------------------------------------------------------------------*/
536 /* copy (clone) generator object */
537 /* */
538 /* parameters: */
539 /* gen ... pointer to generator object */
540 /* */
541 /* return: */
542 /* pointer to clone of generator object */
543 /* */
544 /* error: */
545 /* return NULL */
546 /*----------------------------------------------------------------------*/
547 {
548 #define CLONE ((struct unur_dgt_gen*)clone->datap)
549
550 struct unur_gen *clone;
551
552 /* check arguments */
553 CHECK_NULL(gen,NULL); COOKIE_CHECK(gen,CK_DGT_GEN,NULL);
554
555 /* create generic clone */
556 clone = _unur_generic_clone( gen, GENTYPE );
557
558 /* copy data for distribution */
559 CLONE->cumpv = _unur_xmalloc( DISTR.n_pv * sizeof(double) );
560 memcpy( CLONE->cumpv, GEN->cumpv, DISTR.n_pv * sizeof(double) );
561 CLONE->guide_table = _unur_xmalloc( GEN->guide_size * sizeof(int) );
562 memcpy( CLONE->guide_table, GEN->guide_table, GEN->guide_size * sizeof(int) );
563
564 return clone;
565
566 #undef CLONE
567 } /* end of _unur_dgt_clone() */
568
569 /*---------------------------------------------------------------------------*/
570
571 void
_unur_dgt_free(struct unur_gen * gen)572 _unur_dgt_free( struct unur_gen *gen )
573 /*----------------------------------------------------------------------*/
574 /* deallocate generator object */
575 /* */
576 /* parameters: */
577 /* gen ... pointer to generator object */
578 /*----------------------------------------------------------------------*/
579 {
580
581 /* check arguments */
582 if (!gen) /* nothing to do */
583 return;
584
585 /* check input */
586 if ( gen->method != UNUR_METH_DGT ) {
587 _unur_warning(gen->genid,UNUR_ERR_GEN_INVALID,"");
588 return; }
589 COOKIE_CHECK(gen,CK_DGT_GEN,RETURN_VOID);
590
591 /* we cannot use this generator object any more */
592 SAMPLE = NULL; /* make sure to show up a programming error */
593
594 /* free two auxiliary tables */
595 if (GEN->guide_table) free(GEN->guide_table);
596 if (GEN->cumpv) free(GEN->cumpv);
597
598 /* free memory */
599 _unur_generic_free(gen);
600
601 } /* end of _unur_dgt_free() */
602
603 /*****************************************************************************/
604
605 int
_unur_dgt_sample(struct unur_gen * gen)606 _unur_dgt_sample( struct unur_gen *gen )
607 /*----------------------------------------------------------------------*/
608 /* sample from generator */
609 /* */
610 /* parameters: */
611 /* gen ... pointer to generator object */
612 /* */
613 /* return: */
614 /* integer (sample from random variate) */
615 /* */
616 /* error: */
617 /* return INT_MAX */
618 /*----------------------------------------------------------------------*/
619 {
620 int j;
621 double u;
622
623 /* check arguments */
624 CHECK_NULL(gen,INT_MAX); COOKIE_CHECK(gen,CK_DGT_GEN,INT_MAX);
625
626 /* sample from U(0,1) */
627 u = _unur_call_urng(gen->urng);
628
629 /* look up in guide table ... */
630 j = GEN->guide_table[(int)(u * GEN->guide_size)];
631 /* ... and search */
632 u *= GEN->sum;
633 while (GEN->cumpv[j] < u) j++;
634
635 return (j + DISTR.domain[0]);
636
637 } /* end of _unur_dgt_sample() */
638
639
640 /*---------------------------------------------------------------------------*/
641
642 int
unur_dgt_eval_invcdf_recycle(const struct unur_gen * gen,double u,double * recycle)643 unur_dgt_eval_invcdf_recycle( const struct unur_gen *gen, double u, double *recycle )
644 /*----------------------------------------------------------------------*/
645 /* evaluate inverse CDF at u. */
646 /* */
647 /* parameters: */
648 /* gen ... pointer to generator object */
649 /* u ... argument for inverse CDF (0<=u<=1) */
650 /* recycle ... if not NULL then store recycled 'u' */
651 /* */
652 /* return: */
653 /* integer (inverse CDF) */
654 /* */
655 /* error: */
656 /* return INT_MAX */
657 /*----------------------------------------------------------------------*/
658 {
659 int j;
660
661 /* set default */
662 if (recycle) *recycle = 0.;
663
664 /* check arguments */
665 _unur_check_NULL( GENTYPE, gen, INT_MAX );
666 if ( gen->method != UNUR_METH_DGT ) {
667 _unur_error(gen->genid,UNUR_ERR_GEN_INVALID,"");
668 return INT_MAX;
669 }
670 COOKIE_CHECK(gen,CK_DGT_GEN,INT_MAX);
671
672 /* check range of u */
673 if ( ! (u>0. && u<1.)) {
674 if ( ! (u>=0. && u<=1.)) {
675 _unur_warning(gen->genid,UNUR_ERR_DOMAIN,"U not in [0,1]");
676 }
677 if (u<=0.) return DISTR.domain[0];
678 if (u>=1.) return DISTR.domain[1];
679 return INT_MAX; /* u == NaN */
680 }
681
682 /* look up in guide table ... */
683 j = GEN->guide_table[(int)(u * GEN->guide_size)];
684 /* ... and search */
685 u *= GEN->sum;
686 while (GEN->cumpv[j] < u) j++;
687
688 if (recycle) {
689 *recycle = 1. - (GEN->cumpv[j] - u) / DISTR.pv[j];
690 }
691
692 j+=DISTR.domain[0];
693
694 /* validate range */
695 if (j<DISTR.domain[0]) j = DISTR.domain[0];
696 if (j>DISTR.domain[1]) j = DISTR.domain[1];
697
698 return j;
699
700 } /* end of unur_dgt_eval_invcdf_recycle() */
701
702 /*---------------------------------------------------------------------------*/
703
704 int
unur_dgt_eval_invcdf(const struct unur_gen * gen,double u)705 unur_dgt_eval_invcdf( const struct unur_gen *gen, double u )
706 /*----------------------------------------------------------------------*/
707 /* evaluate inverse CDF at u. */
708 /* */
709 /* parameters: */
710 /* gen ... pointer to generator object */
711 /* u ... argument for inverse CDF (0<=u<=1) */
712 /* */
713 /* return: */
714 /* integer (inverse CDF) */
715 /* */
716 /* error: */
717 /* return INT_MAX */
718 /*----------------------------------------------------------------------*/
719 {
720 return unur_dgt_eval_invcdf_recycle(gen,u,NULL);
721 } /* end of unur_dgt_eval_invcdf() */
722
723 /*****************************************************************************/
724 /** Auxilliary Routines **/
725 /*****************************************************************************/
726
727 int
_unur_dgt_create_tables(struct unur_gen * gen)728 _unur_dgt_create_tables( struct unur_gen *gen )
729 /*----------------------------------------------------------------------*/
730 /* create (allocate) tables */
731 /* */
732 /* parameters: */
733 /* gen ... pointer to generator object */
734 /* */
735 /* return: */
736 /* UNUR_SUCCESS ... on success */
737 /* error code ... on error */
738 /*----------------------------------------------------------------------*/
739 {
740 /* size of guide table */
741 GEN->guide_size = (int)( DISTR.n_pv * GEN->guide_factor);
742 if (GEN->guide_size <= 0)
743 /* do not use a guide table whenever params->guide_factor is 0 or less */
744 GEN->guide_size = 1;
745
746 /* allocation for cummulated probabilities */
747 GEN->cumpv = _unur_xrealloc( GEN->cumpv, DISTR.n_pv * sizeof(double) );
748
749 /* allocate memory for the guide table */
750 GEN->guide_table = _unur_xrealloc( GEN->guide_table, GEN->guide_size * sizeof(int) );
751
752 /* o.k. */
753 return UNUR_SUCCESS;
754 } /* end of _unur_dgt_create_tables() */
755
756 /*---------------------------------------------------------------------------*/
757
758 int
_unur_dgt_make_guidetable(struct unur_gen * gen)759 _unur_dgt_make_guidetable( struct unur_gen *gen )
760 /*----------------------------------------------------------------------*/
761 /* create guide table */
762 /* */
763 /* parameters: */
764 /* gen ... pointer to generator object */
765 /* */
766 /* return: */
767 /* UNUR_SUCCESS ... on success */
768 /* error code ... on error */
769 /*----------------------------------------------------------------------*/
770 {
771 double *pv; /* pointer to probability vector */
772 int n_pv; /* length of probability vector */
773 double pvh; /* aux variable for computing cumulated sums */
774 double gstep; /* step size when computing guide table */
775 int i,j;
776
777 /* probability vector */
778 pv = DISTR.pv;
779 n_pv = DISTR.n_pv;
780
781 /* computation of cumulated probabilities */
782 for( i=0, pvh=0.; i<n_pv; i++ ) {
783 GEN->cumpv[i] = ( pvh += pv[i] );
784 /* ... and check probability vector */
785 if (pv[i] < 0.) {
786 _unur_error(gen->genid,UNUR_ERR_GEN_DATA,"probability < 0");
787 return UNUR_ERR_GEN_DATA;
788 }
789 }
790 GEN->sum = GEN->cumpv[n_pv-1];
791
792 /* computation of guide-table */
793
794 if (gen->variant == DGT_VARFLAG_DIV) {
795 GEN->guide_table[0] = 0;
796 for( j=1, i=0; j<GEN->guide_size ;j++ ) {
797 while( GEN->cumpv[i]/GEN->sum < ((double)j)/GEN->guide_size )
798 i++;
799 if (i >= n_pv) {
800 _unur_warning(gen->genid,UNUR_ERR_ROUNDOFF,"guide table");
801 break;
802 }
803 GEN->guide_table[j]=i;
804 }
805 }
806
807 else { /* gen->variant == DGT_VARFLAG_ADD */
808 gstep = GEN->sum / GEN->guide_size;
809 pvh = 0.;
810 for( j=0, i=0; j<GEN->guide_size ;j++ ) {
811 while (GEN->cumpv[i] < pvh)
812 i++;
813 if (i >= n_pv) {
814 _unur_warning(gen->genid,UNUR_ERR_ROUNDOFF,"guide table");
815 break;
816 }
817 GEN->guide_table[j] = i;
818 pvh += gstep;
819 }
820 }
821
822 /* if there has been an round off error, we have to complete the guide table */
823 for( ; j<GEN->guide_size ;j++ )
824 GEN->guide_table[j] = n_pv - 1;
825 /* o.k. */
826 return UNUR_SUCCESS;
827
828 } /* end of _unur_dgt_make_guidetable() */
829
830
831 /*****************************************************************************/
832 /** Debugging utilities **/
833 /*****************************************************************************/
834
835 /*---------------------------------------------------------------------------*/
836 #ifdef UNUR_ENABLE_LOGGING
837 /*---------------------------------------------------------------------------*/
838
839 void
_unur_dgt_debug_init(struct unur_gen * gen)840 _unur_dgt_debug_init( struct unur_gen *gen )
841 /*----------------------------------------------------------------------*/
842 /* write info about generator into LOG file */
843 /* */
844 /* parameters: */
845 /* gen ... pointer to generator object */
846 /*----------------------------------------------------------------------*/
847 {
848 FILE *LOG;
849
850 /* check arguments */
851 CHECK_NULL(gen,RETURN_VOID); COOKIE_CHECK(gen,CK_DGT_GEN,RETURN_VOID);
852
853 LOG = unur_get_stream();
854
855 fprintf(LOG,"%s:\n",gen->genid);
856 fprintf(LOG,"%s: type = discrete univariate random variates\n",gen->genid);
857 fprintf(LOG,"%s: method = indexed search (guide table)\n",gen->genid);
858
859 fprintf(LOG,"%s: variant = %u ",gen->genid,gen->variant);
860 _unur_print_if_default(gen,DGT_SET_VARIANT);
861 fprintf(LOG,"\n%s:\n",gen->genid);
862
863 _unur_distr_discr_debug( gen->distr,gen->genid,(gen->debug & DGT_DEBUG_PRINTVECTOR));
864
865 fprintf(LOG,"%s: sampling routine = _unur_dgt_sample()\n",gen->genid);
866 fprintf(LOG,"%s:\n",gen->genid);
867
868 fprintf(LOG,"%s: length of probability vector = %d\n",gen->genid,DISTR.n_pv);
869 fprintf(LOG,"%s: length of guide table = %d (rel. = %g%%",
870 gen->genid,GEN->guide_size,100.*GEN->guide_factor);
871 _unur_print_if_default(gen,DGT_SET_GUIDEFACTOR);
872 if (GEN->guide_size == 1)
873 fprintf(LOG,") \t (-->sequential search");
874 fprintf(LOG,")\n%s:\n",gen->genid);
875
876 fprintf(LOG,"%s: sum over PMF (as computed) = %#-20.16g\n",gen->genid,GEN->sum);
877
878 if (gen->debug & DGT_DEBUG_TABLE)
879 _unur_dgt_debug_table(gen);
880
881 } /* end of _unur_dgt_debug_init() */
882
883 /*---------------------------------------------------------------------------*/
884
885 void
_unur_dgt_debug_table(struct unur_gen * gen)886 _unur_dgt_debug_table( struct unur_gen *gen )
887 /*----------------------------------------------------------------------*/
888 /* write guide table into LOG file */
889 /* */
890 /* parameters: */
891 /* gen ... pointer to generator object */
892 /*----------------------------------------------------------------------*/
893 {
894 FILE *LOG;
895 int i,j,m;
896 int n_asts;
897
898 /* check arguments */
899 CHECK_NULL(gen,RETURN_VOID); COOKIE_CHECK(gen,CK_DGT_GEN,RETURN_VOID);
900
901 LOG = unur_get_stream();
902
903 fprintf(LOG,"%s: guide table:\n", gen->genid);
904 fprintf(LOG,"%s:\n", gen->genid);
905 n_asts = 0;
906 for (i=0; i<GEN->guide_size; i++){
907 fprintf(LOG,"%s: [%5d] -> %5d ", gen->genid, i, GEN->guide_table[i]);
908 /* print row of asterisks */
909 if (i == GEN->guide_size-1)
910 j = GEN->guide_size - GEN->guide_table[i];
911 else
912 j = GEN->guide_table[i+1] - GEN->guide_table[i] + 1;
913 for (m=0; m<j && m<10; m++ ) {
914 fprintf(LOG," *");
915 ++n_asts;
916 }
917 /* too many asterisks print */
918 if (m<j){
919 n_asts += j-m;
920 fprintf(LOG," ... %d", j);
921 }
922 fprintf(LOG,"\n");
923 }
924
925 /* print expected number of comparisons */
926 fprintf(LOG,"%s:\n", gen->genid);
927 fprintf(LOG,"%s: expected number of comparisons = %g\n",gen->genid,
928 ((double)n_asts)/GEN->guide_size);
929 fprintf(LOG,"%s:\n", gen->genid);
930
931 fprintf(LOG,"%s:\n",gen->genid);
932
933 } /* end of _unur_dgt_debug_table() */
934
935 /*---------------------------------------------------------------------------*/
936 #endif /* end UNUR_ENABLE_LOGGING */
937 /*---------------------------------------------------------------------------*/
938
939
940 /*---------------------------------------------------------------------------*/
941 #ifdef UNUR_ENABLE_INFO
942 /*---------------------------------------------------------------------------*/
943
944 void
_unur_dgt_info(struct unur_gen * gen,int help)945 _unur_dgt_info( struct unur_gen *gen, int help )
946 /*----------------------------------------------------------------------*/
947 /* create character string that contains information about the */
948 /* given generator object. */
949 /* */
950 /* parameters: */
951 /* gen ... pointer to generator object */
952 /* help ... whether to print additional comments */
953 /*----------------------------------------------------------------------*/
954 {
955 struct unur_string *info = gen->infostr;
956
957 /* generator ID */
958 _unur_string_append(info,"generator ID: %s\n\n", gen->genid);
959
960 /* distribution */
961 _unur_string_append(info,"distribution:\n");
962 _unur_distr_info_typename(gen);
963 _unur_string_append(info," functions = PV [length=%d%s]\n",
964 DISTR.domain[1]-DISTR.domain[0]+1,
965 (DISTR.pmf==NULL) ? "" : ", created from PMF");
966 _unur_string_append(info," domain = (%d, %d)\n", DISTR.domain[0],DISTR.domain[1]);
967 _unur_string_append(info,"\n");
968
969 /* method */
970 _unur_string_append(info,"method: DGT (Guide Table)\n");
971 _unur_string_append(info,"\n");
972
973 /* performance */
974 _unur_string_append(info,"performance characteristics:\n");
975 _unur_string_append(info," E [#look-ups] = %g\n", 1+1./GEN->guide_factor);
976 _unur_string_append(info,"\n");
977
978 /* parameters */
979 if (help) {
980 _unur_string_append(info,"parameters:\n");
981 _unur_string_append(info," guidefactor = %g %s\n", GEN->guide_factor,
982 (gen->set & DGT_SET_GUIDEFACTOR) ? "" : "[default]");
983 if (gen->set & DGT_SET_VARIANT)
984 _unur_string_append(info," variant = %d\n", gen->variant);
985 _unur_string_append(info,"\n");
986 }
987
988 /* Hints */
989 /* if (help) { */
990 /* _unur_string_append(info,"\n"); */
991 /* } */
992
993 } /* end of _unur_dgt_info() */
994
995 /*---------------------------------------------------------------------------*/
996 #endif /* end UNUR_ENABLE_INFO */
997 /*---------------------------------------------------------------------------*/
998