1 /*****************************************************************************
2 * *
3 * UNURAN -- Universal Non-Uniform Random number generator *
4 * *
5 *****************************************************************************
6 * *
7 * FILE: dsrou.c *
8 * *
9 * TYPE: discrete univariate random variate *
10 * METHOD: discrete simple universal method (ratio-of-uniforms method) *
11 * *
12 * DESCRIPTION: *
13 * Given PMF and mode of a T_{-1/2}-concave distribution *
14 * produce a value X consistent with its PMF *
15 * *
16 * REQUIRED: *
17 * pointer to the mass function *
18 * mode of distribution *
19 * sum over PMF *
20 * *
21 * OPTIONAL: *
22 * CDF at mode *
23 * *
24 *****************************************************************************
25 * *
26 * Copyright (c) 2000-2006 Wolfgang Hoermann and Josef Leydold *
27 * Department of Statistics and Mathematics, WU Wien, Austria *
28 * *
29 * This program is free software; you can redistribute it and/or modify *
30 * it under the terms of the GNU General Public License as published by *
31 * the Free Software Foundation; either version 2 of the License, or *
32 * (at your option) any later version. *
33 * *
34 * This program is distributed in the hope that it will be useful, *
35 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
36 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
37 * GNU General Public License for more details. *
38 * *
39 * You should have received a copy of the GNU General Public License *
40 * along with this program; if not, write to the *
41 * Free Software Foundation, Inc., *
42 * 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA *
43 * *
44 *****************************************************************************
45 * *
46 * REFERENCES: *
47 * [1] Leydold J. (2001): A simple universal generator for continuous and *
48 * discrete univariate T-concave distributions, *
49 * ACM Trans. Math. Software 27(1), pp. 66--82. *
50 * *
51 * [2] Kinderman, A.J. and Monahan, F.J. (1977): Computer generation of *
52 * random variables using the ratio of uniform deviates, *
53 * ACM Trans. Math. Software 3(3), pp. 257--260. *
54 * *
55 *****************************************************************************
56 * *
57 * The ratio-of-uniforms method introduced in [2] is a flexible method that *
58 * is based on the following theorem: *
59 * *
60 * THEOREM: *
61 * Let X be a random variable with density function f(x) = g(x) / G, *
62 * where g(x) is a positive integrable function with support (x_0,x_1) *
63 * not necessarily finite and G = integral g(x) dx. *
64 * If (V,U) is uniformly distributed in *
65 * A = {(v,u): 0 < u <= sqrt(g(v/u)), x_0 < v/u < x_1}, *
66 * then X = V/U has probability density function f(x). *
67 * *
68 * Generating point (V,U) uniformly distributed in A is done by rejection *
69 * from an enveloping region, usually from the minimal bounding rectangle. *
70 * *
71 * For discrete random variates the continuous PDF *
72 * PDF(x) = PMF(floor(x)) *
73 * is used. *
74 * *
75 * The implemented algorithm uses the fact, that for many distributions, *
76 * the polygon having the "spikes" of A as its verticeds is convex. *
77 * Then we can find the follow bounding rectangles: *
78 * (For simplicity we have assumed that the sum over the PMF is 1) *
79 * *
80 * Define *
81 * R = {(v,u): -1/u_l <= v <= 0, 0 <= u <= u_l} \cup *
82 * {(v,u): 0 <= v <= 1/u_r, 0 <= u <= u_r} *
83 * Q = {(v,u): v_l <= v <= 0, 0 <= u <= u_l} \cup *
84 * {(v,u): 0 <= v <= v_r, 0 <= u <= u_r} *
85 * where *
86 * u_l = sqrt(PMF(mode-1)), u_r = sqrt(PMF(mode)), *
87 * v_l = -F(mode-1)/u_l, v_r = (1-F(mode-1))/u_r *
88 * Then *
89 * A subset R subset Q *
90 * *
91 * Thus we can use R to generate whenever the CDF F(mode-1) at the mode *
92 * is known, and Q otherwise. *
93 * Notice, that the rection constant is 2 in the first case and 4 and the *
94 * latter. *
95 * *
96 * Distributions with a convex set A are characterized by the following *
97 * theorem that shows a connection to transformed density rejection TDR. *
98 * *
99 * THEOREM: *
100 * A is convex if and only if g is T-concave with transformation *
101 * T(x) = -1/sqrt(x), i.e., -1/sqrt(g(x)) is a concave function. *
102 * *
103 *****************************************************************************/
104
105 /*---------------------------------------------------------------------------*/
106
107 #include <unur_source.h>
108 #include <distr/distr.h>
109 #include <distr/distr_source.h>
110 #include <distr/discr.h>
111 #include <urng/urng.h>
112 #include "unur_methods_source.h"
113 #include "x_gen_source.h"
114 #include "dsrou.h"
115 #include "dsrou_struct.h"
116
117 /*---------------------------------------------------------------------------*/
118 /* Variants: */
119
120 #define DSROU_VARFLAG_VERIFY 0x002u /* run verify mode */
121
122 /*---------------------------------------------------------------------------*/
123 /* Debugging flags */
124 /* bit 01 ... pameters and structure of generator (do not use here) */
125 /* bits 02-12 ... setup */
126 /* bits 13-24 ... adaptive steps */
127 /* bits 25-32 ... trace sampling */
128
129 #define DSROU_DEBUG_REINIT 0x00000010u /* print parameters after reinit */
130
131 /*---------------------------------------------------------------------------*/
132 /* Flags for logging set calls */
133
134 #define DSROU_SET_CDFMODE 0x001u /* CDF at mode is known */
135
136 /*---------------------------------------------------------------------------*/
137
138 #define GENTYPE "DSROU" /* type of generator */
139
140 /*---------------------------------------------------------------------------*/
141
142 static struct unur_gen *_unur_dsrou_init( struct unur_par *par );
143 /*---------------------------------------------------------------------------*/
144 /* Initialize new generator. */
145 /*---------------------------------------------------------------------------*/
146
147 static int _unur_dsrou_reinit( struct unur_gen *gen );
148 /*---------------------------------------------------------------------------*/
149 /* Reinitialize generator. */
150 /*---------------------------------------------------------------------------*/
151
152 static struct unur_gen *_unur_dsrou_create( struct unur_par *par );
153 /*---------------------------------------------------------------------------*/
154 /* create new (almost empty) generator object. */
155 /*---------------------------------------------------------------------------*/
156
157 static int _unur_dsrou_check_par( struct unur_gen *gen );
158 /*---------------------------------------------------------------------------*/
159 /* Check parameters of given distribution and method */
160 /*---------------------------------------------------------------------------*/
161
162 static struct unur_gen *_unur_dsrou_clone( const struct unur_gen *gen );
163 /*---------------------------------------------------------------------------*/
164 /* copy (clone) generator object. */
165 /*---------------------------------------------------------------------------*/
166
167 static void _unur_dsrou_free( struct unur_gen *gen);
168 /*---------------------------------------------------------------------------*/
169 /* destroy generator object. */
170 /*---------------------------------------------------------------------------*/
171
172 static int _unur_dsrou_sample( struct unur_gen *gen );
173 static int _unur_dsrou_sample_check( struct unur_gen *gen );
174 /*---------------------------------------------------------------------------*/
175 /* sample from generator */
176 /*---------------------------------------------------------------------------*/
177
178 static int _unur_dsrou_rectangle( struct unur_gen *gen );
179 /*---------------------------------------------------------------------------*/
180 /* compute universal bounding rectangle. */
181 /*---------------------------------------------------------------------------*/
182
183 #ifdef UNUR_ENABLE_LOGGING
184 /*---------------------------------------------------------------------------*/
185 /* the following functions print debugging information on output stream, */
186 /* i.e., into the LOG file if not specified otherwise. */
187 /*---------------------------------------------------------------------------*/
188 static void _unur_dsrou_debug_init( const struct unur_gen *gen, int is_reinit );
189
190 /*---------------------------------------------------------------------------*/
191 /* print after generator has been initialized has completed. */
192 /*---------------------------------------------------------------------------*/
193 #endif
194
195 #ifdef UNUR_ENABLE_INFO
196 static void _unur_dsrou_info( struct unur_gen *gen, int help );
197 /*---------------------------------------------------------------------------*/
198 /* create info string. */
199 /*---------------------------------------------------------------------------*/
200 #endif
201
202 /*---------------------------------------------------------------------------*/
203 /* abbreviations */
204
205 #define DISTR_IN distr->data.discr /* data for distribution object */
206
207 #define PAR ((struct unur_dsrou_par*)par->datap) /* data for parameter object */
208 #define GEN ((struct unur_dsrou_gen*)gen->datap) /* data for generator object */
209 #define DISTR gen->distr->data.discr /* data for distribution in generator object */
210
211 #define BD_LEFT domain[0] /* left boundary of domain of distribution */
212 #define BD_RIGHT domain[1] /* right boundary of domain of distribution */
213
214 #define SAMPLE gen->sample.discr /* pointer to sampling routine */
215
216 #define PMF(x) _unur_discr_PMF((x),(gen->distr)) /* call to PMF */
217
218 /*---------------------------------------------------------------------------*/
219
220 #define _unur_dsrou_getSAMPLE(gen) \
221 ( ((gen)->variant & DSROU_VARFLAG_VERIFY) \
222 ? _unur_dsrou_sample_check : _unur_dsrou_sample )
223
224 /*---------------------------------------------------------------------------*/
225 /* constants */
226
227 #define SQRT2 (M_SQRT2)
228
229 /*---------------------------------------------------------------------------*/
230
231 /*****************************************************************************/
232 /** Public: User Interface (API) **/
233 /*****************************************************************************/
234
235 struct unur_par *
unur_dsrou_new(const struct unur_distr * distr)236 unur_dsrou_new( const struct unur_distr *distr )
237 /*----------------------------------------------------------------------*/
238 /* get default parameters */
239 /* */
240 /* parameters: */
241 /* distr ... pointer to distribution object */
242 /* */
243 /* return: */
244 /* default parameters (pointer to structure) */
245 /* */
246 /* error: */
247 /* return NULL */
248 /*----------------------------------------------------------------------*/
249 {
250 struct unur_par *par;
251
252 /* check arguments */
253 _unur_check_NULL( GENTYPE,distr,NULL );
254
255 /* check distribution */
256 if (distr->type != UNUR_DISTR_DISCR) {
257 _unur_error(GENTYPE,UNUR_ERR_DISTR_INVALID,""); return NULL; }
258 COOKIE_CHECK(distr,CK_DISTR_DISCR,NULL);
259
260 if (DISTR_IN.pmf == NULL) {
261 _unur_error(GENTYPE,UNUR_ERR_DISTR_REQUIRED,"PMF");
262 return NULL;
263 }
264
265 /* allocate structure */
266 par = _unur_par_new( sizeof(struct unur_dsrou_par) );
267 COOKIE_SET(par,CK_DSROU_PAR);
268
269 /* copy input */
270 par->distr = distr; /* pointer to distribution object */
271
272 /* set default values */
273 PAR->Fmode = -1.; /* CDF at mode (unknown yet) */
274
275 par->method = UNUR_METH_DSROU; /* method and default variant */
276 par->variant = 0u; /* default variant */
277 par->set = 0u; /* inidicate default parameters */
278 par->urng = unur_get_default_urng(); /* use default urng */
279 par->urng_aux = NULL; /* no auxilliary URNG required */
280
281 par->debug = _unur_default_debugflag; /* set default debugging flags */
282
283 /* routine for starting generator */
284 par->init = _unur_dsrou_init;
285
286 return par;
287
288 } /* end of unur_dsrou_new() */
289
290 /*****************************************************************************/
291
292 int
unur_dsrou_set_cdfatmode(struct unur_par * par,double Fmode)293 unur_dsrou_set_cdfatmode( struct unur_par *par, double Fmode )
294 /*----------------------------------------------------------------------*/
295 /* set value of cdf at mode */
296 /* */
297 /* parameters: */
298 /* par ... pointer to parameter for building generator object */
299 /* Fmode ... CDF at mode */
300 /* */
301 /* return: */
302 /* UNUR_SUCCESS ... on success */
303 /* error code ... on error */
304 /*----------------------------------------------------------------------*/
305 {
306 /* check arguments */
307 _unur_check_NULL( GENTYPE, par, UNUR_ERR_NULL );
308 _unur_check_par_object( par, DSROU );
309
310 /* check new parameter for generator */
311 if (Fmode < 0. || Fmode > 1.) {
312 _unur_warning(GENTYPE,UNUR_ERR_PAR_SET,"CDF(mode)");
313 return UNUR_ERR_PAR_SET;
314 }
315
316 /* store date */
317 PAR->Fmode = Fmode;
318
319 /* changelog */
320 par->set |= DSROU_SET_CDFMODE;
321
322 return UNUR_SUCCESS;
323
324 } /* end of unur_dsrou_set_cdfatmode() */
325
326 /*---------------------------------------------------------------------------*/
327
328 int
unur_dsrou_set_verify(struct unur_par * par,int verify)329 unur_dsrou_set_verify( struct unur_par *par, int verify )
330 /*----------------------------------------------------------------------*/
331 /* turn verifying of algorithm while sampling on/off */
332 /* */
333 /* parameters: */
334 /* par ... pointer to parameter for building generator object */
335 /* verify ... 0 = no verifying, !0 = verifying */
336 /* */
337 /* return: */
338 /* UNUR_SUCCESS ... on success */
339 /* error code ... on error */
340 /* */
341 /* comment: */
342 /* no verifying is the default */
343 /*----------------------------------------------------------------------*/
344 {
345 /* check arguments */
346 _unur_check_NULL( GENTYPE, par, UNUR_ERR_NULL );
347 _unur_check_par_object( par, DSROU );
348
349 /* we use a bit in variant */
350 par->variant = (verify) ? (par->variant | DSROU_VARFLAG_VERIFY) : (par->variant & (~DSROU_VARFLAG_VERIFY));
351
352 /* o.k. */
353 return UNUR_ERR_NULL;
354 } /* end of unur_dsrou_set_verify() */
355
356 /*---------------------------------------------------------------------------*/
357
358 int
unur_dsrou_chg_verify(struct unur_gen * gen,int verify)359 unur_dsrou_chg_verify( struct unur_gen *gen, int verify )
360 /*----------------------------------------------------------------------*/
361 /* turn verifying of algorithm while sampling on/off */
362 /* */
363 /* parameters: */
364 /* gen ... pointer to generator object */
365 /* verify ... 0 = no verifying, !0 = verifying */
366 /* */
367 /* return: */
368 /* UNUR_SUCCESS ... on success */
369 /* error code ... on error */
370 /* */
371 /* comment: */
372 /* no verifying is the default */
373 /*----------------------------------------------------------------------*/
374 {
375 /* check input */
376 _unur_check_NULL( GENTYPE, gen, UNUR_ERR_NULL );
377 _unur_check_gen_object( gen, DSROU, UNUR_ERR_GEN_INVALID );
378
379 /* we must not change this switch when sampling has been disabled by
380 using a pointer to the error producing routine */
381 if (SAMPLE == _unur_sample_discr_error)
382 return UNUR_FAILURE;
383
384 if (verify)
385 /* turn verify mode on */
386 gen->variant |= DSROU_VARFLAG_VERIFY;
387 else
388 /* turn verify mode off */
389 gen->variant &= ~DSROU_VARFLAG_VERIFY;
390
391 SAMPLE = _unur_dsrou_getSAMPLE(gen);
392
393 /* o.k. */
394 return UNUR_SUCCESS;
395
396 } /* end of unur_dsrou_chg_verify() */
397
398 /*****************************************************************************/
399
400 int
unur_dsrou_chg_cdfatmode(struct unur_gen * gen,double Fmode)401 unur_dsrou_chg_cdfatmode( struct unur_gen *gen, double Fmode )
402 /*----------------------------------------------------------------------*/
403 /* change value of cdf at mode */
404 /* */
405 /* parameters: */
406 /* gen ... pointer to generator object */
407 /* Fmode ... CDF at mode */
408 /* */
409 /* return: */
410 /* UNUR_SUCCESS ... on success */
411 /* error code ... on error */
412 /*----------------------------------------------------------------------*/
413 {
414 /* check arguments */
415 _unur_check_NULL( GENTYPE, gen, UNUR_ERR_NULL );
416 _unur_check_gen_object( gen, DSROU, UNUR_ERR_GEN_INVALID );
417
418 /* check new parameter for generator */
419 if (Fmode < 0. || Fmode > 1.) {
420 _unur_warning(gen->genid,UNUR_ERR_PAR_SET,"CDF(mode)");
421 return UNUR_ERR_PAR_SET;
422 }
423
424 /* copy parameters */
425 GEN->Fmode = Fmode;
426
427 /* changelog */
428 gen->set |= DSROU_SET_CDFMODE;
429
430 /* o.k. */
431 return UNUR_SUCCESS;
432 } /* end of unur_dsrou_chg_cdfatmode() */
433
434
435 /*****************************************************************************/
436 /** Private **/
437 /*****************************************************************************/
438
439 struct unur_gen *
_unur_dsrou_init(struct unur_par * par)440 _unur_dsrou_init( struct unur_par *par )
441 /*----------------------------------------------------------------------*/
442 /* initialize new generator */
443 /* */
444 /* parameters: */
445 /* params ... pointer to paramters for building generator object */
446 /* */
447 /* return: */
448 /* pointer to generator object */
449 /* */
450 /* error: */
451 /* return NULL */
452 /*----------------------------------------------------------------------*/
453 {
454 struct unur_gen *gen;
455
456 /* check arguments */
457 CHECK_NULL(par,NULL);
458
459 /* check input */
460 if ( par->method != UNUR_METH_DSROU ) {
461 _unur_error(GENTYPE,UNUR_ERR_PAR_INVALID,"");
462 return NULL; }
463 COOKIE_CHECK(par,CK_DSROU_PAR,NULL);
464
465 /* create a new empty generator object */
466 gen = _unur_dsrou_create(par);
467
468 /* free parameters */
469 _unur_par_free(par);
470
471 if (!gen) return NULL;
472
473 /* check parameters */
474 if (_unur_dsrou_check_par(gen) != UNUR_SUCCESS) {
475 _unur_dsrou_free(gen); return NULL;
476 }
477
478 /* compute universal bounding rectangle */
479 if ( _unur_dsrou_rectangle(gen)!=UNUR_SUCCESS ) {
480 _unur_dsrou_free(gen); return NULL;
481 }
482
483 #ifdef UNUR_ENABLE_LOGGING
484 /* write info into LOG file */
485 if (gen->debug) _unur_dsrou_debug_init(gen, FALSE);
486 #endif
487
488 return gen;
489
490 } /* end of _unur_dsrou_init() */
491
492 /*---------------------------------------------------------------------------*/
493
494 int
_unur_dsrou_reinit(struct unur_gen * gen)495 _unur_dsrou_reinit( struct unur_gen *gen )
496 /*----------------------------------------------------------------------*/
497 /* re-initialize (existing) generator. */
498 /* */
499 /* parameters: */
500 /* gen ... pointer to generator object */
501 /* */
502 /* return: */
503 /* UNUR_SUCCESS ... on success */
504 /* error code ... on error */
505 /*----------------------------------------------------------------------*/
506 {
507 int result;
508
509 /* check parameters */
510 if ( (result = _unur_dsrou_check_par(gen)) != UNUR_SUCCESS)
511 return result;
512
513 /* compute universal bounding rectangle */
514 if ( (result = _unur_dsrou_rectangle(gen)) != UNUR_SUCCESS)
515 return result;
516
517 /* (re)set sampling routine */
518 SAMPLE = _unur_dsrou_getSAMPLE(gen);
519
520 #ifdef UNUR_ENABLE_LOGGING
521 /* write info into LOG file */
522 if (gen->debug & DSROU_DEBUG_REINIT) _unur_dsrou_debug_init(gen,TRUE);
523 #endif
524
525 return UNUR_SUCCESS;
526 } /* end of _unur_dsrou_reinit() */
527
528 /*---------------------------------------------------------------------------*/
529
530 struct unur_gen *
_unur_dsrou_create(struct unur_par * par)531 _unur_dsrou_create( struct unur_par *par )
532 /*----------------------------------------------------------------------*/
533 /* allocate memory for generator */
534 /* */
535 /* parameters: */
536 /* par ... pointer to parameter for building generator object */
537 /* */
538 /* return: */
539 /* pointer to (empty) generator object with default settings */
540 /* */
541 /* error: */
542 /* return NULL */
543 /*----------------------------------------------------------------------*/
544 {
545 struct unur_gen *gen;
546
547 /* check arguments */
548 CHECK_NULL(par,NULL); COOKIE_CHECK(par,CK_DSROU_PAR,NULL);
549
550 /* create new generic generator object */
551 gen = _unur_generic_create( par, sizeof(struct unur_dsrou_gen) );
552
553 /* magic cookies */
554 COOKIE_SET(gen,CK_DSROU_GEN);
555
556 /* set generator identifier */
557 gen->genid = _unur_set_genid(GENTYPE);
558
559 /* routines for sampling and destroying generator */
560 SAMPLE = _unur_dsrou_getSAMPLE(gen);
561 gen->destroy = _unur_dsrou_free;
562 gen->clone = _unur_dsrou_clone;
563 gen->reinit = _unur_dsrou_reinit;
564
565 /* copy some parameters into generator object */
566 GEN->Fmode = PAR->Fmode; /* CDF at mode */
567
568 #ifdef UNUR_ENABLE_INFO
569 /* set function for creating info string */
570 gen->info = _unur_dsrou_info;
571 #endif
572
573 /* return pointer to (almost empty) generator object */
574 return gen;
575
576 } /* end of _unur_dsrou_create() */
577
578 /*---------------------------------------------------------------------------*/
579
580 int
_unur_dsrou_check_par(struct unur_gen * gen)581 _unur_dsrou_check_par( struct unur_gen *gen )
582 /*----------------------------------------------------------------------*/
583 /* check parameters of given distribution and method */
584 /* */
585 /* parameters: */
586 /* gen ... pointer to generator object */
587 /* */
588 /* return: */
589 /* UNUR_SUCCESS ... on success */
590 /* error code ... on error */
591 /*----------------------------------------------------------------------*/
592 {
593 /* check for required data: mode */
594 if (!(gen->distr->set & UNUR_DISTR_SET_MODE)) {
595 _unur_warning(GENTYPE,UNUR_ERR_DISTR_REQUIRED,"mode: try finding it (numerically)");
596 if (unur_distr_discr_upd_mode(gen->distr)!=UNUR_SUCCESS) {
597 _unur_error(GENTYPE,UNUR_ERR_DISTR_REQUIRED,"mode");
598 return UNUR_ERR_DISTR_REQUIRED;
599 }
600 }
601
602 /* check for required data: sum over PMF */
603 if (!(gen->distr->set & UNUR_DISTR_SET_PMFSUM))
604 if (unur_distr_discr_upd_pmfsum(gen->distr)!=UNUR_SUCCESS) {
605 _unur_error(GENTYPE,UNUR_ERR_DISTR_REQUIRED,"sum over PMF");
606 return UNUR_ERR_DISTR_REQUIRED;
607 }
608
609 /* mode must be in domain */
610 if ( (DISTR.mode < DISTR.BD_LEFT) ||
611 (DISTR.mode > DISTR.BD_RIGHT) ) {
612 /* there is something wrong.
613 assume: user has change domain without changing mode.
614 but then, she probably has not updated area and is to large */
615 _unur_warning(GENTYPE,UNUR_ERR_GEN_DATA,"area and/or CDF at mode");
616 DISTR.mode = _unur_max(DISTR.mode,DISTR.BD_LEFT);
617 DISTR.mode = _unur_min(DISTR.mode,DISTR.BD_RIGHT);
618 }
619
620 return UNUR_SUCCESS;
621 } /* end of _unur_dsrou_check_par() */
622
623 /*---------------------------------------------------------------------------*/
624
625 struct unur_gen *
_unur_dsrou_clone(const struct unur_gen * gen)626 _unur_dsrou_clone( const struct unur_gen *gen )
627 /*----------------------------------------------------------------------*/
628 /* copy (clone) generator object */
629 /* */
630 /* parameters: */
631 /* gen ... pointer to generator object */
632 /* */
633 /* return: */
634 /* pointer to clone of generator object */
635 /* */
636 /* error: */
637 /* return NULL */
638 /*----------------------------------------------------------------------*/
639 {
640 #define CLONE ((struct unur_dsrou_gen*)clone->datap)
641
642 struct unur_gen *clone;
643
644 /* check arguments */
645 CHECK_NULL(gen,NULL); COOKIE_CHECK(gen,CK_DSROU_GEN,NULL);
646
647 /* create generic clone */
648 clone = _unur_generic_clone( gen, GENTYPE );
649
650 return clone;
651
652 #undef CLONE
653 } /* end of _unur_dsrou_clone() */
654
655 /*---------------------------------------------------------------------------*/
656
657 void
_unur_dsrou_free(struct unur_gen * gen)658 _unur_dsrou_free( struct unur_gen *gen )
659 /*----------------------------------------------------------------------*/
660 /* deallocate generator object */
661 /* */
662 /* parameters: */
663 /* gen ... pointer to generator object */
664 /*----------------------------------------------------------------------*/
665 {
666 /* check arguments */
667 if( !gen ) /* nothing to do */
668 return;
669
670 /* check input */
671 if ( gen->method != UNUR_METH_DSROU ) {
672 _unur_warning(gen->genid,UNUR_ERR_GEN_INVALID,"");
673 return; }
674 COOKIE_CHECK(gen,CK_DSROU_GEN,RETURN_VOID);
675
676 /* we cannot use this generator object any more */
677 SAMPLE = NULL; /* make sure to show up a programming error */
678
679 /* free memory */
680 _unur_generic_free(gen);
681
682 } /* end of _unur_dsrou_free() */
683
684 /*****************************************************************************/
685
686 int
_unur_dsrou_sample(struct unur_gen * gen)687 _unur_dsrou_sample( struct unur_gen *gen )
688 /*----------------------------------------------------------------------*/
689 /* sample from generator */
690 /* */
691 /* parameters: */
692 /* gen ... pointer to generator object */
693 /* */
694 /* return: */
695 /* int (sample from random variate) */
696 /* */
697 /* error: */
698 /* return INT_MAX */
699 /*----------------------------------------------------------------------*/
700 {
701 double U,V,X;
702 int I;
703
704 /* check arguments */
705 CHECK_NULL(gen,INT_MAX); COOKIE_CHECK(gen,CK_DSROU_GEN,INT_MAX);
706
707 while (1) {
708 /* generate point uniformly in union of rectangles */
709 V = GEN->al + _unur_call_urng(gen->urng) * (GEN->ar - GEN->al);
710 V /= (V<0.) ? GEN->ul : GEN->ur; /* if ul==0. then al==0. and thus V>=0. */
711
712 while ( _unur_iszero(U = _unur_call_urng(gen->urng)));
713 U *= (V<0.) ? GEN->ul : GEN->ur;
714
715 /* ratio */
716 X = floor(V/U) + DISTR.mode;
717
718 /* inside domain ? */
719 if ( (X < DISTR.BD_LEFT) || (X > DISTR.BD_RIGHT) )
720 continue;
721
722 /* convert to int */
723 I = (int) X;
724
725 /* accept or reject */
726 if (U*U <= PMF(I))
727 return I;
728 }
729 } /* end of _unur_dsrou_sample() */
730
731 /*---------------------------------------------------------------------------*/
732
733 int
_unur_dsrou_sample_check(struct unur_gen * gen)734 _unur_dsrou_sample_check( struct unur_gen *gen )
735 /*----------------------------------------------------------------------*/
736 /* sample from generator and verify that method can be used */
737 /* */
738 /* parameters: */
739 /* gen ... pointer to generator object */
740 /* */
741 /* return: */
742 /* int (sample from random variate) */
743 /* */
744 /* error: */
745 /* return INT_MAX */
746 /*----------------------------------------------------------------------*/
747 {
748 double U,V,pI,VI,X;
749 double um2, vl, vr;
750 int I;
751
752 /* check arguments */
753 CHECK_NULL(gen,INT_MAX); COOKIE_CHECK(gen,CK_DSROU_GEN,INT_MAX);
754
755 while (1) {
756 /* generate point uniformly in union of rectangles */
757 V = GEN->al + _unur_call_urng(gen->urng) * (GEN->ar - GEN->al);
758 V /= (V<0.) ? GEN->ul : GEN->ur;
759
760 while ( _unur_iszero(U = _unur_call_urng(gen->urng)));
761 U *= (V<0.) ? GEN->ul : GEN->ur;
762
763 /* ratios */
764 X = floor(V/U) + DISTR.mode;
765
766 /* inside domain ? */
767 if ( (X < DISTR.BD_LEFT) || (X > DISTR.BD_RIGHT) )
768 continue;
769
770 /* convert to int */
771 I = (int) X;
772
773 /* values of PMF and v-coordinate of point */
774 pI = PMF(I);
775 VI = V/U * sqrt(pI);
776
777 /* values of boundary of rectangle */
778 /* (avoid roundoff error with FP registers) */
779 um2 = (2.+4.*DBL_EPSILON) * ((V<0) ? GEN->ul*GEN->ul : GEN->ur*GEN->ur);
780 vl = (GEN->ul>0.) ? (1.+UNUR_EPSILON) * GEN->al/GEN->ul : 0.;
781 vr = (1.+UNUR_EPSILON) * GEN->ar/GEN->ur;
782
783 /* check hat */
784 if ( pI > um2 || VI < vl || VI > vr ) {
785 /* printf("pI = %g < %g VI = %g < %g < %g\n", */
786 /* pI, ur2, vl, VI, vr); */
787 _unur_error(gen->genid,UNUR_ERR_GEN_CONDITION,"PMF(x) > hat(x)");
788 }
789
790 /* accept or reject */
791 if (U*U <= pI)
792 return I;
793 }
794 } /* end of _unur_dsrou_sample_check() */
795
796 /*****************************************************************************/
797 /** Auxilliary Routines **/
798 /*****************************************************************************/
799
800 int
_unur_dsrou_rectangle(struct unur_gen * gen)801 _unur_dsrou_rectangle( struct unur_gen *gen )
802 /*----------------------------------------------------------------------*/
803 /* compute universal bounding rectangle */
804 /* */
805 /* parameters: */
806 /* gen ... pointer to generator object */
807 /* */
808 /* return: */
809 /* UNUR_SUCCESS ... on success */
810 /* error code ... on error */
811 /*----------------------------------------------------------------------*/
812 {
813 double pm, pbm; /* PMF at mode and mode-1 */
814
815 /* check arguments */
816 CHECK_NULL( gen, UNUR_ERR_NULL );
817 COOKIE_CHECK( gen,CK_DSROU_GEN, UNUR_ERR_COOKIE );
818
819 /* compute PMF at mode and mode-1 */
820 pm = PMF(DISTR.mode);
821 pbm = (DISTR.mode-1 < DISTR.BD_LEFT) ? 0. : PMF(DISTR.mode-1);
822
823 /* pm and pbm must be positive */
824 if (pm <= 0. || pbm < 0.) {
825 _unur_error(gen->genid,UNUR_ERR_GEN_DATA,"PMF(mode) <= 0.");
826 return UNUR_ERR_GEN_DATA;
827 }
828
829 /* heights of rectangles */
830 GEN->ul = sqrt(pbm);
831 GEN->ur = sqrt(pm);
832
833 /* areas of rectangle */
834 if (_unur_iszero(GEN->ul)) {
835 /* PMF monotonically decreasing */
836 GEN->al = 0.;
837 GEN->ar = DISTR.sum;
838 }
839 else if (gen->set & DSROU_SET_CDFMODE) {
840 /* CDF at mode known */
841 GEN->al = -(GEN->Fmode * DISTR.sum)+pm;
842 GEN->ar = DISTR.sum + GEN->al;
843 }
844 else {
845 GEN->al = -(DISTR.sum - pm);
846 GEN->ar = DISTR.sum;
847 }
848
849 /* o.k. */
850 return UNUR_SUCCESS;
851
852 } /* end of _unur_dsrou_rectangle() */
853
854 /*****************************************************************************/
855 /** Debugging utilities **/
856 /*****************************************************************************/
857
858 /*---------------------------------------------------------------------------*/
859 #ifdef UNUR_ENABLE_LOGGING
860 /*---------------------------------------------------------------------------*/
861
862 static void
_unur_dsrou_debug_init(const struct unur_gen * gen,int is_reinit)863 _unur_dsrou_debug_init( const struct unur_gen *gen, int is_reinit )
864 /*----------------------------------------------------------------------*/
865 /* write info about generator into LOG file */
866 /* */
867 /* parameters: */
868 /* gen ... pointer to generator object */
869 /* is_reinit ... if TRUE the generator has been reinitialized */
870 /*----------------------------------------------------------------------*/
871 {
872 FILE *LOG;
873
874 /* check arguments */
875 CHECK_NULL(gen,RETURN_VOID); COOKIE_CHECK(gen,CK_DSROU_GEN,RETURN_VOID);
876
877 LOG = unur_get_stream();
878
879 fprintf(LOG,"%s:\n",gen->genid);
880 if (!is_reinit) {
881 fprintf(LOG,"%s: type = discrete univariate random variates\n",gen->genid);
882 fprintf(LOG,"%s: method = dsrou (discrete simple universal ratio-of-uniforms)\n",gen->genid);
883 }
884 else
885 fprintf(LOG,"%s: reinit!\n",gen->genid);
886 fprintf(LOG,"%s:\n",gen->genid);
887
888 _unur_distr_discr_debug( gen->distr, gen->genid, FALSE );
889
890 fprintf(LOG,"%s: sampling routine = _unur_dsrou_sample",gen->genid);
891 if (gen->variant & DSROU_VARFLAG_VERIFY)
892 fprintf(LOG,"_check");
893 fprintf(LOG,"()\n%s:\n",gen->genid);
894
895 if (gen->set & DSROU_SET_CDFMODE)
896 fprintf(LOG,"%s: CDF(mode) = %g\n",gen->genid,GEN->Fmode);
897 else
898 fprintf(LOG,"%s: CDF(mode) unknown\n",gen->genid);
899
900 fprintf(LOG,"%s: no (universal) squeeze\n",gen->genid);
901 fprintf(LOG,"%s: no mirror principle\n",gen->genid);
902 fprintf(LOG,"%s:\n",gen->genid);
903
904 fprintf(LOG,"%s: Rectangles:\n",gen->genid);
905 if (GEN->ul > 0.)
906 fprintf(LOG,"%s: left upper point = (%g,%g) \tarea = %g (%5.2f%%)\n",
907 gen->genid,GEN->al/GEN->ul,GEN->ul,fabs(GEN->al),100.*fabs(GEN->al)/(-GEN->al+GEN->ar));
908 else
909 fprintf(LOG,"%s: left upper point = (0,0) \tarea = 0 (0.00%%)\n",gen->genid);
910
911 fprintf(LOG,"%s: right upper point = (%g,%g) \tarea = %g (%5.2f%%)\n",
912 gen->genid,GEN->ar/GEN->ur,GEN->ur,GEN->ar,100.*GEN->ar/(-GEN->al+GEN->ar));
913
914 fprintf(LOG,"%s:\n",gen->genid);
915
916 fflush(LOG);
917
918 } /* end of _unur_dsrou_debug_init() */
919
920 /*---------------------------------------------------------------------------*/
921 #endif /* end UNUR_ENABLE_LOGGING */
922 /*---------------------------------------------------------------------------*/
923
924
925 /*---------------------------------------------------------------------------*/
926 #ifdef UNUR_ENABLE_INFO
927 /*---------------------------------------------------------------------------*/
928
929 void
_unur_dsrou_info(struct unur_gen * gen,int help)930 _unur_dsrou_info( struct unur_gen *gen, int help )
931 /*----------------------------------------------------------------------*/
932 /* create character string that contains information about the */
933 /* given generator object. */
934 /* */
935 /* parameters: */
936 /* gen ... pointer to generator object */
937 /* help ... whether to print additional comments */
938 /*----------------------------------------------------------------------*/
939 {
940 struct unur_string *info = gen->infostr;
941 struct unur_distr *distr = gen->distr;
942
943 /* generator ID */
944 _unur_string_append(info,"generator ID: %s\n\n", gen->genid);
945
946 /* distribution */
947 _unur_string_append(info,"distribution:\n");
948 _unur_distr_info_typename(gen);
949 _unur_string_append(info," functions = PMF\n");
950 _unur_string_append(info," domain = (%d, %d)\n", DISTR.domain[0],DISTR.domain[1]);
951 _unur_string_append(info," mode = %d %s\n", DISTR.mode,
952 (distr->set & UNUR_DISTR_SET_MODE_APPROX) ? "[numeric.]" : "");
953 _unur_string_append(info," sum(PMF) = %g\n", DISTR.sum);
954 if (gen->set & DSROU_SET_CDFMODE)
955 _unur_string_append(info," F(mode) = %g\n", GEN->Fmode);
956 else
957 _unur_string_append(info," F(mode) = [unknown]\n");
958 _unur_string_append(info,"\n");
959
960 if (help) {
961 if ( distr->set & UNUR_DISTR_SET_MODE_APPROX )
962 _unur_string_append(info,"[ Hint: %s ]\n",
963 "You may provide the \"mode\"");
964 _unur_string_append(info,"\n");
965 }
966
967 /* method */
968 _unur_string_append(info,"method: DSROU (Discrete Simple Ratio-Of-Uniforms)\n");
969 _unur_string_append(info,"\n");
970
971 /* performance */
972 _unur_string_append(info,"performance characteristics:\n");
973
974 _unur_string_append(info," enveloping rectangle = (%g,%g) x (%g,%g) [left]\n",
975 (GEN->ul > 0.)?GEN->al/GEN->ul:0., 0.,
976 0., (GEN->ul > 0.)?GEN->ul:0.);
977 _unur_string_append(info," (%g,%g) x (%g,%g) [right]\n",
978 0.,GEN->ar/GEN->ur, 0., GEN->ur);
979
980 _unur_string_append(info," area(hat) = %g + %g = %g\n",
981 fabs(GEN->al), GEN->ar, -GEN->al+GEN->ar);
982
983 _unur_string_append(info," rejection constant = %g\n",
984 2. * (-GEN->al+GEN->ar) / DISTR.sum);
985 _unur_string_append(info,"\n");
986
987 /* parameters */
988 if (help) {
989 _unur_string_append(info,"parameters:\n");
990 if (gen->set & DSROU_SET_CDFMODE)
991 _unur_string_append(info," cdfatmode = %g\n", GEN->Fmode);
992 else
993 _unur_string_append(info," cdfatmode = [not set]\n");
994
995 if (gen->variant & DSROU_VARFLAG_VERIFY)
996 _unur_string_append(info," verify = on\n");
997
998 _unur_string_append(info,"\n");
999 }
1000
1001 /* Hints */
1002 if (help) {
1003 if ( !(gen->set & DSROU_SET_CDFMODE))
1004 _unur_string_append(info,"[ Hint: %s ]\n",
1005 "You can set \"cdfatmode\" to reduce the rejection constant.");
1006 _unur_string_append(info,"\n");
1007 }
1008
1009 } /* end of _unur_dsrou_info() */
1010
1011 /*---------------------------------------------------------------------------*/
1012 #endif /* end UNUR_ENABLE_INFO */
1013 /*---------------------------------------------------------------------------*/
1014