1 /*****************************************************************************
2 * *
3 * UNURAN -- Universal Non-Uniform Random number generator *
4 * *
5 *****************************************************************************
6 * *
7 * FILE: srou.c *
8 * *
9 * TYPE: continuous univariate random variate *
10 * METHOD: simple universal method (ratio-of-uniforms method) *
11 * *
12 * DESCRIPTION: *
13 * Given PDF and mode of a T_{-1/2}-concave distribution *
14 * produce a value x consistent with its density *
15 * *
16 * REQUIRED: *
17 * pointer to the density function *
18 * mode of the density *
19 * area below PDF *
20 * OPTIONAL: *
21 * CDF at mode *
22 * *
23 *****************************************************************************
24 * *
25 * Copyright (c) 2000-2006 Wolfgang Hoermann and Josef Leydold *
26 * Department of Statistics and Mathematics, WU Wien, Austria *
27 * *
28 * This program is free software; you can redistribute it and/or modify *
29 * it under the terms of the GNU General Public License as published by *
30 * the Free Software Foundation; either version 2 of the License, or *
31 * (at your option) any later version. *
32 * *
33 * This program is distributed in the hope that it will be useful, *
34 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
35 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
36 * GNU General Public License for more details. *
37 * *
38 * You should have received a copy of the GNU General Public License *
39 * along with this program; if not, write to the *
40 * Free Software Foundation, Inc., *
41 * 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA *
42 * *
43 *****************************************************************************
44 * *
45 * REFERENCES: *
46 * [1] Leydold J. (2001): A simple universal generator for continuous and *
47 * discrete univariate T-concave distributions, *
48 * ACM Trans. Math. Software 27(1), pp. 66--82. *
49 * *
50 * [2] Kinderman, A.J. and Monahan, F.J. (1977): Computer generation of *
51 * random variables using the ratio of uniform deviates, *
52 * ACM Trans. Math. Software 3(3), pp. 257--260. *
53 * *
54 * [3] Leydold J. (2002): Short universal generators via generalized *
55 * ratio-of-uniforms method, Preprint. *
56 * *
57 * [4] Wakefield J. C., Gelfand A. E., and Smith A. F. M. (1991): *
58 * Efficient generation of random variates via the ratio-of-uniforms *
59 * method, Statist. Comput. 1(2), pp. 129--133. *
60 * *
61 *****************************************************************************
62 * *
63 * The ratio-of-uniforms method introduced in [2] is a flexible method that *
64 * is based on the following theorem: *
65 * *
66 * THEOREM: *
67 * Let X be a random variable with density function f(x) = g(x) / G, *
68 * where g(x) is a positive integrable function with support (x_0,x_1) *
69 * not necessarily finite and G = integral g(x) dx. *
70 * If (V,U) is uniformly distributed in *
71 * A = {(v,u): 0 < u <= sqrt(g(v/u)), x_0 < v/u < x_1}, *
72 * then X = V/U has probability density function f(x). *
73 * *
74 * Generating point (V,U) uniformly distributed in A is done by rejection *
75 * from an enveloping region, usually from the minimal bounding rectangle. *
76 * *
77 * The implemented algorithm uses the fact, that for many distribtions, *
78 * A is convex. Then we easily can construct an enveloping rectangle. *
79 * Define *
80 * R = {(v,u): v_l <= v <= v_r, 0 <= u <= u_m}, *
81 * Q = {(v,u): -v_m <= v <= v_m, 0 <= u <= u_m}, *
82 * where *
83 * u_m = sqrt(f(mode)), v_m = (\int f dx) / u_m *
84 * v_l = -F(\mode) v_m, v_r = (1-F(mode)) v_m *
85 * Then *
86 * A subset R subset Q *
87 * *
88 * Thus we can use R to generate whenever the CDF F(mode) at the mode *
89 * is known, and Q otherwise. *
90 * Notice, that the rection constant is 2 in the first case and 4 and the *
91 * latter. *
92 * *
93 * If F(mode) it known, it is even possible to get an universal squeeze *
94 * (see [1] for details). However its usage is only recommended when *
95 * the PDF is (very) expensive. *
96 * *
97 * When F(mode) is not known the mirror principle can be used. i.e., make *
98 * an enveloping rectangle for f(x)+f(-x). It reduces the rejection constant *
99 * to 2 * sqrt(2) at the expense of more evaluations of the PDF. *
100 * Its usage is only recommended when the generation time for the underlying *
101 * uniform prng is extreme large. *
102 * *
103 * Distributions with a convex set A are characterized by the following *
104 * theorem that shows a connection to transformed density rejection TDR. *
105 * *
106 * THEOREM: *
107 * A is convex if and only if g is T-concave with transformation *
108 * T(x) = -1/sqrt(x), i.e., -1/sqrt(g(x)) is a concave function. *
109 * *
110 *****************************************************************************
111 * *
112 * The theory for r!=1 is described in [3]. *
113 * *
114 *****************************************************************************/
115
116 /*---------------------------------------------------------------------------*/
117
118 #include <unur_source.h>
119 #include <distr/distr.h>
120 #include <distr/distr_source.h>
121 #include <distr/cont.h>
122 #include <urng/urng.h>
123 #include "unur_methods_source.h"
124 #include "x_gen_source.h"
125 #include "srou.h"
126 #include "srou_struct.h"
127
128 #ifdef UNUR_ENABLE_INFO
129 # include <tests/unuran_tests.h>
130 #endif
131
132 /*---------------------------------------------------------------------------*/
133 /* Variants: */
134
135 #define SROU_VARFLAG_VERIFY 0x002u /* run verify mode */
136 #define SROU_VARFLAG_SQUEEZE 0x004u /* use universal squeeze if possible */
137 #define SROU_VARFLAG_MIRROR 0x008u /* use mirror principle */
138
139 /*---------------------------------------------------------------------------*/
140 /* Debugging flags */
141 /* bit 01 ... pameters and structure of generator (do not use here) */
142 /* bits 02-12 ... setup */
143 /* bits 13-24 ... adaptive steps */
144 /* bits 25-32 ... trace sampling */
145
146 #define SROU_DEBUG_REINIT 0x00000010u /* print parameters after reinit */
147
148 /*---------------------------------------------------------------------------*/
149 /* Flags for logging set calls */
150
151 #define SROU_SET_R 0x001u /* parameter r for power transform. */
152 #define SROU_SET_CDFMODE 0x002u /* CDF at mode is known */
153 #define SROU_SET_PDFMODE 0x004u /* PDF at mode is set */
154
155 /*---------------------------------------------------------------------------*/
156
157 #define GENTYPE "SROU" /* type of generator */
158
159 /*---------------------------------------------------------------------------*/
160
161 static struct unur_gen *_unur_srou_init( struct unur_par *par );
162 /*---------------------------------------------------------------------------*/
163 /* Initialize new generator. */
164 /*---------------------------------------------------------------------------*/
165
166 static int _unur_srou_reinit( struct unur_gen *gen );
167 /*---------------------------------------------------------------------------*/
168 /* Reinitialize generator. */
169 /*---------------------------------------------------------------------------*/
170
171 static struct unur_gen *_unur_srou_create( struct unur_par *par );
172 /*---------------------------------------------------------------------------*/
173 /* create new (almost empty) generator object. */
174 /*---------------------------------------------------------------------------*/
175
176 static int _unur_srou_check_par( struct unur_gen *gen );
177 /*---------------------------------------------------------------------------*/
178 /* Check parameters of given distribution and method */
179 /*---------------------------------------------------------------------------*/
180
181 static struct unur_gen *_unur_srou_clone( const struct unur_gen *gen );
182 /*---------------------------------------------------------------------------*/
183 /* copy (clone) generator object. */
184 /*---------------------------------------------------------------------------*/
185
186 static void _unur_srou_free( struct unur_gen *gen);
187 /*---------------------------------------------------------------------------*/
188 /* destroy generator object. */
189 /*---------------------------------------------------------------------------*/
190
191 static double _unur_srou_sample( struct unur_gen *gen );
192 static double _unur_srou_sample_mirror( struct unur_gen *gen );
193 static double _unur_srou_sample_check( struct unur_gen *gen );
194 /*---------------------------------------------------------------------------*/
195 /* sample from generator (case r=1). */
196 /*---------------------------------------------------------------------------*/
197
198 static double _unur_gsrou_sample( struct unur_gen *gen );
199 static double _unur_gsrou_sample_check( struct unur_gen *gen );
200 /*---------------------------------------------------------------------------*/
201 /* sample from generator (case r>1). */
202 /*---------------------------------------------------------------------------*/
203
204 static int _unur_srou_rectangle( struct unur_gen *gen );
205 /*---------------------------------------------------------------------------*/
206 /* compute universal bounding rectangle (case r=1). */
207 /*---------------------------------------------------------------------------*/
208
209 static int _unur_gsrou_envelope( struct unur_gen *gen );
210 /*---------------------------------------------------------------------------*/
211 /* compute parameters for universal bounding envelope (case r>1). */
212 /*---------------------------------------------------------------------------*/
213
214 #ifdef UNUR_ENABLE_LOGGING
215 /*---------------------------------------------------------------------------*/
216 /* the following functions print debugging information on output stream, */
217 /* i.e., into the LOG file if not specified otherwise. */
218 /*---------------------------------------------------------------------------*/
219
220 static void _unur_srou_debug_init( const struct unur_gen *gen, int is_reinit );
221 /*---------------------------------------------------------------------------*/
222 /* print after generator has been initialized has completed. */
223 /*---------------------------------------------------------------------------*/
224 #endif
225
226 #ifdef UNUR_ENABLE_INFO
227 static void _unur_srou_info( struct unur_gen *gen, int help );
228 /*---------------------------------------------------------------------------*/
229 /* create info string. */
230 /*---------------------------------------------------------------------------*/
231 #endif
232
233 /*---------------------------------------------------------------------------*/
234 /* abbreviations */
235
236 #define DISTR_IN distr->data.cont /* data for distribution object */
237
238 #define PAR ((struct unur_srou_par*)par->datap) /* data for parameter object */
239 #define GEN ((struct unur_srou_gen*)gen->datap) /* data for generator object */
240 #define DISTR gen->distr->data.cont /* data for distribution in generator object */
241
242 #define BD_LEFT domain[0] /* left boundary of domain of distribution */
243 #define BD_RIGHT domain[1] /* right boundary of domain of distribution */
244
245 #define SAMPLE gen->sample.cont /* pointer to sampling routine */
246
247 #define PDF(x) _unur_cont_PDF((x),(gen->distr)) /* call to PDF */
248
249 /*---------------------------------------------------------------------------*/
250
251 static UNUR_SAMPLING_ROUTINE_CONT *
_unur_srou_getSAMPLE(struct unur_gen * gen)252 _unur_srou_getSAMPLE( struct unur_gen *gen )
253 {
254 if (gen->variant & SROU_VARFLAG_VERIFY)
255 return (gen->set & SROU_SET_R) ? _unur_gsrou_sample_check : _unur_srou_sample_check;
256 else {
257 if (gen->set & SROU_SET_R)
258 return _unur_gsrou_sample;
259 else
260 return (gen->variant & SROU_VARFLAG_MIRROR) ? _unur_srou_sample_mirror : _unur_srou_sample;
261 }
262 } /* end of _unur_srou_getSAMPLE() */
263
264 /*---------------------------------------------------------------------------*/
265 /* constants */
266
267 #define SQRT2 (M_SQRT2)
268
269 /*---------------------------------------------------------------------------*/
270
271 /*****************************************************************************/
272 /** Public: User Interface (API) **/
273 /*****************************************************************************/
274
275 struct unur_par *
unur_srou_new(const struct unur_distr * distr)276 unur_srou_new( const struct unur_distr *distr )
277 /*----------------------------------------------------------------------*/
278 /* get default parameters */
279 /* */
280 /* parameters: */
281 /* distr ... pointer to distribution object */
282 /* */
283 /* return: */
284 /* default parameters (pointer to structure) */
285 /* */
286 /* error: */
287 /* return NULL */
288 /*----------------------------------------------------------------------*/
289 {
290 struct unur_par *par;
291
292 /* check arguments */
293 _unur_check_NULL( GENTYPE,distr,NULL );
294
295 /* check distribution */
296 if (distr->type != UNUR_DISTR_CONT) {
297 _unur_error(GENTYPE,UNUR_ERR_DISTR_INVALID,""); return NULL; }
298 COOKIE_CHECK(distr,CK_DISTR_CONT,NULL);
299
300 if (DISTR_IN.pdf == NULL) {
301 _unur_error(GENTYPE,UNUR_ERR_DISTR_REQUIRED,"PDF");
302 return NULL;
303 }
304
305 /* allocate structure */
306 par = _unur_par_new( sizeof(struct unur_srou_par) );
307 COOKIE_SET(par,CK_SROU_PAR);
308
309 /* copy input */
310 par->distr = distr; /* pointer to distribution object */
311
312 /* set default values */
313 PAR->r = 1.; /* parameter for power transformation */
314 PAR->Fmode = -1.; /* CDF at mode (unknown yet ) */
315 PAR->um = -1.; /* (square) root of PDF at mode (unknown) */
316
317 par->method = UNUR_METH_SROU; /* method and default variant */
318 par->variant = 0u; /* default variant */
319 par->set = 0u; /* inidicate default parameters */
320 par->urng = unur_get_default_urng(); /* use default urng */
321 par->urng_aux = NULL; /* no auxilliary URNG required */
322
323 par->debug = _unur_default_debugflag; /* set default debugging flags */
324
325 /* routine for starting generator */
326 par->init = _unur_srou_init;
327
328 return par;
329
330 } /* end of unur_srou_new() */
331
332 /*****************************************************************************/
333
334 int
unur_srou_set_r(struct unur_par * par,double r)335 unur_srou_set_r( struct unur_par *par, double r )
336 /*----------------------------------------------------------------------*/
337 /* set parameter r for power transformation */
338 /* */
339 /* parameters: */
340 /* par ... pointer to parameter for building generator object */
341 /* r ... parameter r */
342 /* */
343 /* return: */
344 /* UNUR_SUCCESS ... on success */
345 /* error code ... on error */
346 /*----------------------------------------------------------------------*/
347 {
348 /* check arguments */
349 _unur_check_NULL( GENTYPE, par, UNUR_ERR_NULL );
350 _unur_check_par_object( par, SROU );
351
352 /* check new parameter for generator */
353 if (r < 1.) {
354 _unur_warning(GENTYPE,UNUR_ERR_PAR_SET,"r < 1");
355 return UNUR_ERR_PAR_SET;
356 }
357
358 if (_unur_isone(r)) {
359 /* simple version, same as R is not set */
360 PAR->r = r;
361 par->set &= ~SROU_SET_R;
362 }
363 else {
364 /* for computational reasons r should be at least 1.01 */
365 if (r<1.01) r = 1.01;
366 PAR->r = r;
367 par->set |= SROU_SET_R;
368 }
369
370 /* we have to reset the marker for the PDF at the mode */
371 par->set &= ~SROU_SET_PDFMODE;
372
373 return UNUR_SUCCESS;
374
375 } /* end of unur_srou_set_r() */
376
377 /*---------------------------------------------------------------------------*/
378
379 int
unur_srou_set_cdfatmode(struct unur_par * par,double Fmode)380 unur_srou_set_cdfatmode( struct unur_par *par, double Fmode )
381 /*----------------------------------------------------------------------*/
382 /* set value of cdf at mode */
383 /* */
384 /* parameters: */
385 /* par ... pointer to parameter for building generator object */
386 /* Fmode ... cdf at mode */
387 /* */
388 /* return: */
389 /* UNUR_SUCCESS ... on success */
390 /* error code ... on error */
391 /*----------------------------------------------------------------------*/
392 {
393 /* check arguments */
394 _unur_check_NULL( GENTYPE, par, UNUR_ERR_NULL );
395 _unur_check_par_object( par, SROU );
396
397 /* check new parameter for generator */
398 if (Fmode < 0. || Fmode > 1.) {
399 _unur_warning(GENTYPE,UNUR_ERR_PAR_SET,"CDF(mode)");
400 return UNUR_ERR_PAR_SET;
401 }
402
403 /* store date */
404 PAR->Fmode = Fmode;
405
406 /* changelog */
407 par->set |= SROU_SET_CDFMODE;
408
409 return UNUR_SUCCESS;
410
411 } /* end of unur_srou_set_cdfatmode() */
412
413 /*---------------------------------------------------------------------------*/
414
415 int
unur_srou_set_pdfatmode(UNUR_PAR * par,double fmode)416 unur_srou_set_pdfatmode( UNUR_PAR *par, double fmode )
417 /*----------------------------------------------------------------------*/
418 /* Set pdf at mode. if set the PDF at the mode is never changed. */
419 /* */
420 /* parameters: */
421 /* par ... pointer to parameter for building generator object */
422 /* fmode ... PDF at mode */
423 /* */
424 /* return: */
425 /* UNUR_SUCCESS ... on success */
426 /* error code ... on error */
427 /*----------------------------------------------------------------------*/
428 {
429 /* check arguments */
430 _unur_check_NULL( GENTYPE, par, UNUR_ERR_NULL );
431 _unur_check_par_object( par, SROU );
432
433 /* check new parameter for generator */
434 if (fmode <= 0.) {
435 _unur_warning(GENTYPE,UNUR_ERR_PAR_SET,"PDF(mode)");
436 return UNUR_ERR_PAR_SET;
437 }
438 if (!_unur_isfinite(fmode)) {
439 _unur_warning(GENTYPE,UNUR_ERR_PAR_SET,"PDF(mode) overflow");
440 return UNUR_ERR_PAR_SET;
441 }
442
443 /* store date ((square) root of fmode) */
444 PAR->um = (par->set & SROU_SET_R) ? pow(fmode,1./(PAR->r+1.)) : sqrt(fmode);
445
446 /* changelog */
447 par->set |= SROU_SET_PDFMODE;
448
449 return UNUR_SUCCESS;
450
451 } /* end of unur_srou_set_pdfatmode() */
452
453 /*---------------------------------------------------------------------------*/
454
455 int
unur_srou_set_verify(struct unur_par * par,int verify)456 unur_srou_set_verify( struct unur_par *par, int verify )
457 /*----------------------------------------------------------------------*/
458 /* turn verifying of algorithm while sampling on/off */
459 /* */
460 /* parameters: */
461 /* par ... pointer to parameter for building generator object */
462 /* verify ... 0 = no verifying, !0 = verifying */
463 /* */
464 /* return: */
465 /* UNUR_SUCCESS ... on success */
466 /* error code ... on error */
467 /* */
468 /* comment: */
469 /* no verifying is the default */
470 /*----------------------------------------------------------------------*/
471 {
472 /* check arguments */
473 _unur_check_NULL( GENTYPE, par, UNUR_ERR_NULL );
474 _unur_check_par_object( par, SROU );
475
476 /* we use a bit in variant */
477 par->variant = (verify) ? (par->variant | SROU_VARFLAG_VERIFY) : (par->variant & (~SROU_VARFLAG_VERIFY));
478
479 /* o.k. */
480 return UNUR_SUCCESS;
481
482 } /* end of unur_srou_set_verify() */
483
484 /*---------------------------------------------------------------------------*/
485
486 int
unur_srou_chg_verify(struct unur_gen * gen,int verify)487 unur_srou_chg_verify( struct unur_gen *gen, int verify )
488 /*----------------------------------------------------------------------*/
489 /* turn verifying of algorithm while sampling on/off */
490 /* */
491 /* parameters: */
492 /* gen ... pointer to generator object */
493 /* verify ... 0 = no verifying, !0 = verifying */
494 /* */
495 /* return: */
496 /* UNUR_SUCCESS ... on success */
497 /* error code ... on error */
498 /* */
499 /* comment: */
500 /* no verifying is the default */
501 /*----------------------------------------------------------------------*/
502 {
503 /* check input */
504 _unur_check_NULL( GENTYPE, gen, UNUR_ERR_NULL );
505 _unur_check_gen_object( gen, SROU, UNUR_ERR_GEN_INVALID );
506
507 /* we must not change this switch when sampling has been disabled by
508 using a pointer to the error producing routine */
509 if (SAMPLE == _unur_sample_cont_error)
510 return UNUR_FAILURE;
511
512 if (verify)
513 /* turn verify mode on */
514 gen->variant |= SROU_VARFLAG_VERIFY;
515 else
516 /* turn verify mode off */
517 gen->variant &= ~SROU_VARFLAG_VERIFY;
518
519 SAMPLE = _unur_srou_getSAMPLE(gen);
520
521 /* o.k. */
522 return UNUR_SUCCESS;
523
524 } /* end of unur_srou_chg_verify() */
525
526 /*---------------------------------------------------------------------------*/
527
528 int
unur_srou_set_usesqueeze(struct unur_par * par,int usesqueeze)529 unur_srou_set_usesqueeze( struct unur_par *par, int usesqueeze )
530 /*----------------------------------------------------------------------*/
531 /* set flag for using universal squeeze (default: off) */
532 /* */
533 /* parameters: */
534 /* par ... pointer to parameter for building generator object */
535 /* usesqueeze ... 0 = no squeeze, !0 = use squeeze */
536 /* */
537 /* return: */
538 /* UNUR_SUCCESS ... on success */
539 /* error code ... on error */
540 /* */
541 /* comment: */
542 /* no squeeze is the default */
543 /*----------------------------------------------------------------------*/
544 {
545 /* check arguments */
546 _unur_check_NULL( GENTYPE, par, UNUR_ERR_NULL );
547 _unur_check_par_object( par, SROU );
548
549 /* we use a bit in variant */
550 par->variant = (usesqueeze) ? (par->variant | SROU_VARFLAG_SQUEEZE) : (par->variant & (~SROU_VARFLAG_SQUEEZE));
551
552 /* o.k. */
553 return UNUR_SUCCESS;
554
555 } /* end of unur_srou_set_usesqueeze() */
556
557 /*---------------------------------------------------------------------------*/
558
559 int
unur_srou_set_usemirror(struct unur_par * par,int usemirror)560 unur_srou_set_usemirror( struct unur_par *par, int usemirror )
561 /*----------------------------------------------------------------------*/
562 /* set flag for using mirror principle (default: off) */
563 /* */
564 /* parameters: */
565 /* par ... pointer to parameter for building generator object */
566 /* usemirror ... 0 = no mirror princ., !0 = use mirror principle */
567 /* */
568 /* return: */
569 /* UNUR_SUCCESS ... on success */
570 /* error code ... on error */
571 /* */
572 /* comment: */
573 /* do not use mirror principle is the default */
574 /*----------------------------------------------------------------------*/
575 {
576 /* check arguments */
577 _unur_check_NULL( GENTYPE, par, UNUR_ERR_NULL );
578 _unur_check_par_object( par, SROU );
579
580 /* we use a bit in variant */
581 par->variant = (usemirror) ? (par->variant | SROU_VARFLAG_MIRROR) : (par->variant & (~SROU_VARFLAG_MIRROR));
582
583 /* o.k. */
584 return UNUR_SUCCESS;
585
586 } /* end of unur_srou_set_usemirror() */
587
588 /*****************************************************************************/
589
590 int
unur_srou_chg_cdfatmode(struct unur_gen * gen,double Fmode)591 unur_srou_chg_cdfatmode( struct unur_gen *gen, double Fmode )
592 /*----------------------------------------------------------------------*/
593 /* change value of cdf at mode */
594 /* */
595 /* parameters: */
596 /* gen ... pointer to generator object */
597 /* Fmode ... cdf at mode */
598 /* */
599 /* return: */
600 /* UNUR_SUCCESS ... on success */
601 /* error code ... on error */
602 /*----------------------------------------------------------------------*/
603 {
604 /* check arguments */
605 _unur_check_NULL( GENTYPE, gen, UNUR_ERR_NULL );
606 _unur_check_gen_object( gen, SROU, UNUR_ERR_GEN_INVALID );
607
608 /* check new parameter for generator */
609 if (Fmode < 0. || Fmode > 1.) {
610 _unur_warning(gen->genid,UNUR_ERR_PAR_SET,"CDF(mode)");
611 return UNUR_ERR_PAR_SET;
612 }
613
614 /* copy parameters */
615 GEN->Fmode = Fmode;
616
617 /* changelog */
618 gen->set |= SROU_SET_CDFMODE;
619
620 /* o.k. */
621 return UNUR_SUCCESS;
622 } /* end of unur_srou_chg_cdfatmode() */
623
624 /*---------------------------------------------------------------------------*/
625
626 int
unur_srou_chg_pdfatmode(struct unur_gen * gen,double fmode)627 unur_srou_chg_pdfatmode( struct unur_gen *gen, double fmode )
628 /*----------------------------------------------------------------------*/
629 /* change value of PDF at mode */
630 /* */
631 /* parameters: */
632 /* gen ... pointer to generator object */
633 /* fmode ... PDF at mode */
634 /* */
635 /* return: */
636 /* UNUR_SUCCESS ... on success */
637 /* error code ... on error */
638 /*----------------------------------------------------------------------*/
639 {
640 /* check arguments */
641 _unur_check_NULL( GENTYPE, gen, UNUR_ERR_NULL );
642 _unur_check_gen_object( gen, SROU, UNUR_ERR_GEN_INVALID );
643
644 /* check new parameter for generator */
645 if (fmode <= 0.) {
646 _unur_warning(gen->genid,UNUR_ERR_PAR_SET,"PDF(mode)");
647 return UNUR_ERR_PAR_SET;
648 }
649 if (!_unur_isfinite(fmode)) {
650 _unur_warning(GENTYPE,UNUR_ERR_PAR_SET,"PDF(mode) overflow");
651 return UNUR_ERR_PAR_SET;
652 }
653
654 /* store date ((square) root of fmode) */
655 GEN->um = (gen->set & SROU_SET_R) ? pow(fmode,1./(GEN->r+1.)) : sqrt(fmode);
656
657 /* changelog */
658 gen->set |= SROU_SET_PDFMODE;
659
660 /* o.k. */
661 return UNUR_SUCCESS;
662 } /* end of unur_srou_chg_pdfatmode() */
663
664
665 /*****************************************************************************/
666 /** Private **/
667 /*****************************************************************************/
668
669 struct unur_gen *
_unur_srou_init(struct unur_par * par)670 _unur_srou_init( struct unur_par *par )
671 /*----------------------------------------------------------------------*/
672 /* initialize new generator */
673 /* */
674 /* parameters: */
675 /* params ... pointer to paramters for building generator object */
676 /* */
677 /* return: */
678 /* pointer to generator object */
679 /* */
680 /* error: */
681 /* return NULL */
682 /*----------------------------------------------------------------------*/
683 {
684 struct unur_gen *gen;
685 int rcode;
686
687 /* check arguments */
688 CHECK_NULL(par,NULL);
689
690 /* check input */
691 if ( par->method != UNUR_METH_SROU ) {
692 _unur_error(GENTYPE,UNUR_ERR_PAR_INVALID,"");
693 return NULL; }
694 COOKIE_CHECK(par,CK_SROU_PAR,NULL);
695
696 /* there are no squeezes and no mirrir principle when r is changed */
697 if (par->set & SROU_SET_R) {
698 par->variant &= ~SROU_VARFLAG_MIRROR;
699 par->variant &= ~SROU_VARFLAG_SQUEEZE;
700 }
701
702 if (par->set & SROU_SET_CDFMODE)
703 /* cdf at mode known -->
704 thus it does not make sense to use the mirror principle */
705 par->variant &= ~SROU_VARFLAG_MIRROR;
706 else
707 /* cdf at mode unknown -->
708 thus we cannot use universal squeeze */
709 par->variant &= ~SROU_VARFLAG_SQUEEZE;
710
711 /* create a new empty generator object */
712 gen = _unur_srou_create(par);
713
714 /* free parameters */
715 _unur_par_free(par);
716
717 if (!gen) return NULL;
718
719 /* check parameters */
720 if (_unur_srou_check_par(gen) != UNUR_SUCCESS) {
721 _unur_srou_free(gen); return NULL;
722 }
723
724 /* compute universal bounding envelope */
725 if (gen->set & SROU_SET_R)
726 rcode = _unur_gsrou_envelope( gen );
727 else
728 rcode = _unur_srou_rectangle( gen );
729
730 if (rcode!=UNUR_SUCCESS) {
731 /* error */
732 _unur_srou_free(gen); return NULL;
733 }
734
735 #ifdef UNUR_ENABLE_LOGGING
736 /* write info into LOG file */
737 if (gen->debug) _unur_srou_debug_init(gen, FALSE);
738 #endif
739
740 return gen;
741
742 } /* end of _unur_srou_init() */
743
744 /*---------------------------------------------------------------------------*/
745
746 int
_unur_srou_reinit(struct unur_gen * gen)747 _unur_srou_reinit( struct unur_gen *gen )
748 /*----------------------------------------------------------------------*/
749 /* re-initialize (existing) generator. */
750 /* */
751 /* parameters: */
752 /* gen ... pointer to generator object */
753 /* */
754 /* return: */
755 /* UNUR_SUCCESS ... on success */
756 /* error code ... on error */
757 /*----------------------------------------------------------------------*/
758 {
759 int rcode;
760
761 /* check parameters */
762 if ( (rcode = _unur_srou_check_par(gen)) != UNUR_SUCCESS)
763 return rcode;
764
765 /* compute universal bounding envelope */
766 if (gen->set & SROU_SET_R)
767 rcode = _unur_gsrou_envelope( gen );
768 else
769 rcode = _unur_srou_rectangle( gen );
770
771 /* (re)set sampling routine */
772 SAMPLE = _unur_srou_getSAMPLE(gen);
773
774 #ifdef UNUR_ENABLE_LOGGING
775 /* write info into LOG file */
776 if (gen->debug & SROU_DEBUG_REINIT) _unur_srou_debug_init(gen,TRUE);
777 #endif
778
779 return rcode;
780 } /* end of _unur_srou_reinit() */
781
782 /*---------------------------------------------------------------------------*/
783
784 struct unur_gen *
_unur_srou_create(struct unur_par * par)785 _unur_srou_create( struct unur_par *par )
786 /*----------------------------------------------------------------------*/
787 /* allocate memory for generator */
788 /* */
789 /* parameters: */
790 /* par ... pointer to parameter for building generator object */
791 /* */
792 /* return: */
793 /* pointer to (empty) generator object with default settings */
794 /* */
795 /* error: */
796 /* return NULL */
797 /*----------------------------------------------------------------------*/
798 {
799 struct unur_gen *gen;
800
801 /* check arguments */
802 CHECK_NULL(par,NULL); COOKIE_CHECK(par,CK_SROU_PAR,NULL);
803
804 /* create new generic generator object */
805 gen = _unur_generic_create( par, sizeof(struct unur_srou_gen) );
806
807 /* magic cookies */
808 COOKIE_SET(gen,CK_SROU_GEN);
809
810 /* set generator identifier */
811 gen->genid = _unur_set_genid(GENTYPE);
812
813 /* routines for sampling and destroying generator */
814 SAMPLE = _unur_srou_getSAMPLE(gen);
815 gen->destroy = _unur_srou_free;
816 gen->clone = _unur_srou_clone;
817 gen->reinit = _unur_srou_reinit;
818
819 /* copy some parameters into generator object */
820 GEN->r = PAR->r; /* parameter for power transformation */
821 GEN->Fmode = PAR->Fmode; /* CDF at mode */
822 GEN->um = PAR->um; /* square root of PDF at mode */
823
824 /* initialize variables */
825 GEN->vl = GEN->vr = 0.;
826 GEN->xl = GEN->xr = 0.;
827 GEN->p = 0.;
828 GEN->a = GEN->b = 0.;
829 GEN->log_ab = 0.;
830
831 #ifdef UNUR_ENABLE_INFO
832 /* set function for creating info string */
833 gen->info = _unur_srou_info;
834 #endif
835
836 /* return pointer to (almost empty) generator object */
837 return gen;
838
839 } /* end of _unur_srou_create() */
840
841 /*---------------------------------------------------------------------------*/
842
843 int
_unur_srou_check_par(struct unur_gen * gen)844 _unur_srou_check_par( struct unur_gen *gen )
845 /*----------------------------------------------------------------------*/
846 /* check parameters of given distribution and method */
847 /* */
848 /* parameters: */
849 /* gen ... pointer to generator object */
850 /* */
851 /* return: */
852 /* UNUR_SUCCESS ... on success */
853 /* error code ... on error */
854 /*----------------------------------------------------------------------*/
855 {
856 /* check for required data: mode */
857 if (!(gen->distr->set & UNUR_DISTR_SET_MODE)) {
858 _unur_warning(GENTYPE,UNUR_ERR_DISTR_REQUIRED,"mode: try finding it (numerically)");
859 if (unur_distr_cont_upd_mode(gen->distr)!=UNUR_SUCCESS) {
860 _unur_error(GENTYPE,UNUR_ERR_DISTR_REQUIRED,"mode");
861 return UNUR_ERR_DISTR_REQUIRED;
862 }
863 }
864
865 /* check for required data: area */
866 if (!(gen->distr->set & UNUR_DISTR_SET_PDFAREA)) {
867 if (unur_distr_cont_upd_pdfarea(gen->distr)!=UNUR_SUCCESS) {
868 _unur_error(GENTYPE,UNUR_ERR_DISTR_REQUIRED,"area below PDF");
869 return UNUR_ERR_DISTR_REQUIRED;
870 }
871 }
872
873 /* mode must be in domain */
874 if ( (DISTR.mode < DISTR.BD_LEFT) ||
875 (DISTR.mode > DISTR.BD_RIGHT) ) {
876 /* there is something wrong.
877 assume: user has change domain without changing mode.
878 but then, she probably has not updated area and is to large */
879 _unur_warning(GENTYPE,UNUR_ERR_GEN_DATA,"area and/or CDF at mode");
880 DISTR.mode = _unur_max(DISTR.mode,DISTR.BD_LEFT);
881 DISTR.mode = _unur_min(DISTR.mode,DISTR.BD_RIGHT);
882 }
883
884 return UNUR_SUCCESS;
885 } /* end of _unur_srou_check_par() */
886
887 /*---------------------------------------------------------------------------*/
888
889 struct unur_gen *
_unur_srou_clone(const struct unur_gen * gen)890 _unur_srou_clone( const struct unur_gen *gen )
891 /*----------------------------------------------------------------------*/
892 /* copy (clone) generator object */
893 /* */
894 /* parameters: */
895 /* gen ... pointer to generator object */
896 /* */
897 /* return: */
898 /* pointer to clone of generator object */
899 /* */
900 /* error: */
901 /* return NULL */
902 /*----------------------------------------------------------------------*/
903 {
904 #define CLONE ((struct unur_srou_gen*)clone->datap)
905
906 struct unur_gen *clone;
907
908 /* check arguments */
909 CHECK_NULL(gen,NULL); COOKIE_CHECK(gen,CK_SROU_GEN,NULL);
910
911 /* create generic clone */
912 clone = _unur_generic_clone( gen, GENTYPE );
913
914 return clone;
915
916 #undef CLONE
917 } /* end of _unur_srou_clone() */
918
919 /*---------------------------------------------------------------------------*/
920
921 void
_unur_srou_free(struct unur_gen * gen)922 _unur_srou_free( struct unur_gen *gen )
923 /*----------------------------------------------------------------------*/
924 /* deallocate generator object */
925 /* */
926 /* parameters: */
927 /* gen ... pointer to generator object */
928 /*----------------------------------------------------------------------*/
929 {
930 /* check arguments */
931 if( !gen ) /* nothing to do */
932 return;
933
934 /* check input */
935 if ( gen->method != UNUR_METH_SROU ) {
936 _unur_warning(gen->genid,UNUR_ERR_GEN_INVALID,"");
937 return; }
938 COOKIE_CHECK(gen,CK_SROU_GEN,RETURN_VOID);
939
940 /* we cannot use this generator object any more */
941 SAMPLE = NULL; /* make sure to show up a programming error */
942
943 /* free memory */
944 _unur_generic_free(gen);
945
946 } /* end of _unur_srou_free() */
947
948 /*****************************************************************************/
949
950 double
_unur_srou_sample(struct unur_gen * gen)951 _unur_srou_sample( struct unur_gen *gen )
952 /*----------------------------------------------------------------------*/
953 /* sample from generator */
954 /* */
955 /* parameters: */
956 /* gen ... pointer to generator object */
957 /* */
958 /* return: */
959 /* double (sample from random variate) */
960 /* */
961 /* error: */
962 /* return INFINITY */
963 /*----------------------------------------------------------------------*/
964 {
965 double U,V,X,x,xx;
966
967 /* check arguments */
968 CHECK_NULL(gen,INFINITY); COOKIE_CHECK(gen,CK_SROU_GEN,INFINITY);
969
970 while (1) {
971 /* generate point uniformly on rectangle */
972 while ( _unur_iszero(U = _unur_call_urng(gen->urng)) );
973 U *= GEN->um;
974 V = GEN->vl + _unur_call_urng(gen->urng) * (GEN->vr - GEN->vl);
975
976 /* ratio */
977 X = V/U;
978
979 /* compute x */
980 x = X + DISTR.mode;
981
982 /* inside domain ? */
983 if ( (x < DISTR.BD_LEFT) || (x > DISTR.BD_RIGHT) )
984 continue;
985
986 /* evaluate squeeze */
987 if ( (gen->variant & SROU_VARFLAG_SQUEEZE) &&
988 (X >= GEN->xl) &&
989 (X <= GEN->xr ) &&
990 (U < GEN->um) ) {
991 xx = V / (GEN->um - U);
992 if ( (xx >= GEN->xl) && (xx <= GEN->xr ) )
993 return x;
994 }
995
996 /* accept or reject */
997 if (U*U <= PDF(x))
998 return x;
999 }
1000
1001 } /* end of _unur_srou_sample() */
1002
1003 /*---------------------------------------------------------------------------*/
1004
1005 double
_unur_srou_sample_mirror(struct unur_gen * gen)1006 _unur_srou_sample_mirror( struct unur_gen *gen )
1007 /*----------------------------------------------------------------------*/
1008 /* sample from generator, use mirror principle */
1009 /* */
1010 /* parameters: */
1011 /* gen ... pointer to generator object */
1012 /* */
1013 /* return: */
1014 /* double (sample from random variate) */
1015 /* */
1016 /* error: */
1017 /* return INFINITY */
1018 /*----------------------------------------------------------------------*/
1019 {
1020 double U,V,X,x,fx,fnx,uu;
1021
1022 /* check arguments */
1023 CHECK_NULL(gen,INFINITY); COOKIE_CHECK(gen,CK_SROU_GEN,INFINITY);
1024
1025 while (1) {
1026 /* generate point uniformly on rectangle */
1027 while ( _unur_iszero(U = _unur_call_urng(gen->urng)) );
1028 U *= GEN->um * SQRT2;
1029 V = 2. * (_unur_call_urng(gen->urng) - 0.5) * GEN->vr;
1030 /* vr = vm when the CDF at the mode is not known */
1031
1032 /* ratio */
1033 X = V/U;
1034
1035 /* evaluate PDF */
1036 x = X + DISTR.mode;
1037 fx = (x < DISTR.BD_LEFT || x > DISTR.BD_RIGHT) ? 0. : PDF(x);
1038 uu = U * U;
1039
1040 /* accept or reject */
1041 if (uu <= fx)
1042 return x;
1043
1044 /* try mirrored PDF */
1045 x = -X + DISTR.mode;
1046 fnx = (x < DISTR.BD_LEFT || x > DISTR.BD_RIGHT) ? 0. : PDF(x);
1047 if (uu <= fx + fnx)
1048 return x;
1049 }
1050
1051 } /* end of _unur_srou_sample_mirror() */
1052
1053 /*---------------------------------------------------------------------------*/
1054
1055 double
_unur_srou_sample_check(struct unur_gen * gen)1056 _unur_srou_sample_check( struct unur_gen *gen )
1057 /*----------------------------------------------------------------------*/
1058 /* sample from generator and verify that method can be used */
1059 /* */
1060 /* parameters: */
1061 /* gen ... pointer to generator object */
1062 /* */
1063 /* return: */
1064 /* double (sample from random variate) */
1065 /* */
1066 /* error: */
1067 /* return INFINITY */
1068 /*----------------------------------------------------------------------*/
1069 {
1070 double U,uu,V,X,x,nx,fx,sfx,fnx,xfx,xfnx,xx;
1071
1072 /* check arguments */
1073 CHECK_NULL(gen,INFINITY); COOKIE_CHECK(gen,CK_SROU_GEN,INFINITY);
1074
1075 if (gen->variant & SROU_VARFLAG_MIRROR) {
1076 /* use mirror principle */
1077
1078 while (1) {
1079 /* generate point uniformly on rectangle */
1080 while ( _unur_iszero(U = _unur_call_urng(gen->urng)) );
1081 U *= GEN->um * SQRT2;
1082 V = 2. * (_unur_call_urng(gen->urng) - 0.5) * GEN->vr;
1083 /* vr = vm when the CDF at the mode is not known */
1084
1085 /* ratio */
1086 X = V/U;
1087
1088 /* x values */
1089 x = X + DISTR.mode;
1090 nx = -X + DISTR.mode;
1091
1092 /* evaluate PDF */
1093 fx = (x < DISTR.BD_LEFT || x > DISTR.BD_RIGHT) ? 0. : PDF(x);
1094 fnx = (nx < DISTR.BD_LEFT || nx > DISTR.BD_RIGHT) ? 0. : PDF(nx);
1095 uu = U * U;
1096
1097 /* check hat */
1098 xfx = (x - DISTR.mode) * sqrt(fx);
1099 xfnx = (nx - DISTR.mode) * sqrt(fnx);
1100
1101 if ( ((2.+4.*DBL_EPSILON) * GEN->um*GEN->um < fx + fnx) /* avoid roundoff error with FP registers */
1102 || (xfx < (1.+UNUR_EPSILON) * GEN->vl)
1103 || (xfx > (1.+UNUR_EPSILON) * GEN->vr)
1104 || (xfnx < (1.+UNUR_EPSILON) * GEN->vl)
1105 || (xfnx > (1.+UNUR_EPSILON) * GEN->vr) )
1106 _unur_error(gen->genid,UNUR_ERR_GEN_CONDITION,"PDF(x) > hat(x)");
1107
1108 /* accept or reject */
1109 if (uu <= fx)
1110 return x;
1111
1112 /* try mirrored PDF */
1113 if (uu <= fx + fnx)
1114 return nx;
1115 }
1116 }
1117
1118 else { /* do not use mirror principle */
1119
1120 while (1) {
1121 /* generate point uniformly on rectangle */
1122 while ( _unur_iszero(U = _unur_call_urng(gen->urng)) );
1123 U *= GEN->um;
1124 V = GEN->vl + _unur_call_urng(gen->urng) * (GEN->vr - GEN->vl);
1125
1126 /* ratio */
1127 X = V/U;
1128
1129 /* compute x */
1130 x = X + DISTR.mode;
1131
1132 /* inside domain ? */
1133 if ( (x < DISTR.BD_LEFT) || (x > DISTR.BD_RIGHT) )
1134 continue;
1135
1136 /* evaluate PDF */
1137 fx = PDF(x);
1138
1139 /* the point on the boundary of the region of acceptance
1140 in direction X = V/U has the coordinates
1141 ( X * sqrt(fx), sqrt(fx) ). */
1142 sfx = sqrt(fx);
1143 xfx = X * sfx;
1144
1145 /* check hat */
1146 if ( ( sfx > (1.+DBL_EPSILON) * GEN->um ) /* avoid roundoff error with FP registers */
1147 || (xfx < (1.+UNUR_EPSILON) * GEN->vl)
1148 || (xfx > (1.+UNUR_EPSILON) * GEN->vr) )
1149 _unur_error(gen->genid,UNUR_ERR_GEN_CONDITION,"PDF(x) > hat(x)");
1150
1151 /* evaluate and check squeeze */
1152 if ( (gen->variant & SROU_VARFLAG_SQUEEZE) &&
1153 (X >= GEN->xl) &&
1154 (X <= GEN->xr ) &&
1155 (U < GEN->um) ) {
1156
1157 /* check squeeze */
1158 xx = xfx / (GEN->um - sfx);
1159 if ( (xx > (1.-UNUR_EPSILON) * GEN->xl) &&
1160 (xx < (1.-UNUR_EPSILON) * GEN->xr) )
1161 _unur_error(gen->genid,UNUR_ERR_GEN_CONDITION,"PDF(x) < squeeze(x)");
1162
1163 /* squeeze acceptance */
1164 xx = V / (GEN->um - U);
1165 if ( (xx >= GEN->xl) && (xx <= GEN->xr ) )
1166 return x;
1167 }
1168
1169 /* accept or reject */
1170 if (U*U <= PDF(x))
1171 return x;
1172 }
1173 }
1174
1175 } /* end of _unur_srou_sample_check() */
1176
1177 /*****************************************************************************/
1178
1179 double
_unur_gsrou_sample(struct unur_gen * gen)1180 _unur_gsrou_sample( struct unur_gen *gen )
1181 /*----------------------------------------------------------------------*/
1182 /* sample from generator */
1183 /* */
1184 /* parameters: */
1185 /* gen ... pointer to generator object */
1186 /* */
1187 /* return: */
1188 /* double (sample from random variate) */
1189 /* */
1190 /* error: */
1191 /* return INFINITY */
1192 /*----------------------------------------------------------------------*/
1193 {
1194 double U,Ur,V,W,X,Z;
1195
1196 /* check arguments */
1197 CHECK_NULL(gen,INFINITY); COOKIE_CHECK(gen,CK_SROU_GEN,INFINITY);
1198
1199 while (1) {
1200 W = GEN->log_ab *_unur_call_urng(gen->urng);
1201 Z = GEN->vl + _unur_call_urng(gen->urng) * (GEN->vr - GEN->vl);
1202 U = (exp(-W)-1.) * GEN->a/GEN->b;
1203 V = -Z/(GEN->a + GEN->b*U);
1204 U *= GEN->um;
1205 Ur = pow(U,GEN->r);
1206 X = V/Ur + DISTR.mode;
1207
1208 /* inside domain ? */
1209 if ( (X < DISTR.BD_LEFT) || (X > DISTR.BD_RIGHT) )
1210 continue;
1211
1212 /* accept or reject */
1213 if (Ur*U <= PDF(X))
1214 return X;
1215 }
1216
1217 } /* end of _unur_gsrou_sample() */
1218
1219 /*---------------------------------------------------------------------------*/
1220
1221 double
_unur_gsrou_sample_check(struct unur_gen * gen)1222 _unur_gsrou_sample_check( struct unur_gen *gen )
1223 /*----------------------------------------------------------------------*/
1224 /* sample from generator and verify that method can be used */
1225 /* */
1226 /* parameters: */
1227 /* gen ... pointer to generator object */
1228 /* */
1229 /* return: */
1230 /* double (sample from random variate) */
1231 /* */
1232 /* error: */
1233 /* return INFINITY */
1234 /*----------------------------------------------------------------------*/
1235 {
1236 double U,Ur,V,W,X,x,Z;
1237 double fx,uf,vf,vhl,vhr;
1238
1239 /* check arguments */
1240 CHECK_NULL(gen,INFINITY); COOKIE_CHECK(gen,CK_SROU_GEN,INFINITY);
1241
1242 while (1) {
1243
1244 W = GEN->log_ab *_unur_call_urng(gen->urng);
1245 Z = GEN->vl + _unur_call_urng(gen->urng) * (GEN->vr - GEN->vl);
1246 U = (exp(-W)-1.) * GEN->a/GEN->b;
1247 V = -Z/(GEN->a + GEN->b*U);
1248 U *= GEN->um;
1249 Ur = pow(U,GEN->r);
1250 X = V/Ur;
1251
1252 /* compute x */
1253 x = X + DISTR.mode;
1254 /* inside domain ? */
1255 if ( (x < DISTR.BD_LEFT) || (x > DISTR.BD_RIGHT) )
1256 continue;
1257
1258 /* evaluate density */
1259 fx = PDF(x);
1260
1261 /* the point on the boundary of the region of acceptance
1262 in direction X = V/U^r has the coordinates
1263 ( (x-mode) * (fx)^(r/(r+1)), sqrt[r+1](fx) ). */
1264 uf = pow(fx,1./(GEN->r+1));
1265 vf = X * pow(fx,GEN->r/(GEN->r+1.));
1266
1267 /* the corresponding point on boundary of the enveloping region */
1268 /* with same u-coordinate. */
1269 vhl = - GEN->vl /(GEN->a + GEN->b*(uf/GEN->um));
1270 vhr = - GEN->vr /(GEN->a + GEN->b*(uf/GEN->um));
1271
1272 /* check hat */
1273 if ( ( uf > (1.+DBL_EPSILON) * GEN->um ) /* avoid roundoff error with FP registers */
1274 || (vf < (1.+UNUR_EPSILON) * vhl )
1275 || (vf > (1.+UNUR_EPSILON) * vhr ) )
1276 {
1277 _unur_error(gen->genid,UNUR_ERR_GEN_CONDITION,"PDF(x) > hat(x)");
1278 }
1279
1280 /* accept or reject */
1281 if (Ur*U <= fx)
1282 return x;
1283 }
1284
1285 } /* end of _unur_gsrou_sample_check() */
1286
1287 /*****************************************************************************/
1288 /** Auxilliary Routines **/
1289 /*****************************************************************************/
1290
1291 int
_unur_srou_rectangle(struct unur_gen * gen)1292 _unur_srou_rectangle( struct unur_gen *gen )
1293 /*----------------------------------------------------------------------*/
1294 /* compute universal bounding rectangle */
1295 /* */
1296 /* parameters: */
1297 /* gen ... pointer to generator object */
1298 /* */
1299 /* return: */
1300 /* UNUR_SUCCESS ... on success */
1301 /* error code ... on error */
1302 /*----------------------------------------------------------------------*/
1303 {
1304 double vm, fm; /* width of rectangle, PDF at mode */
1305
1306 /* check arguments */
1307 CHECK_NULL( gen, UNUR_ERR_NULL );
1308 COOKIE_CHECK( gen,CK_SROU_GEN, UNUR_ERR_COOKIE );
1309
1310 /* compute PDF at mode (if not given by user) */
1311 if (!(gen->set & SROU_SET_PDFMODE)) {
1312 fm = PDF(DISTR.mode);
1313 /* fm must be positive */
1314 if (fm <= 0.) {
1315 _unur_error(gen->genid,UNUR_ERR_GEN_DATA,"PDF(mode) <= 0.");
1316 return UNUR_ERR_GEN_DATA;
1317 }
1318 if (!_unur_isfinite(fm)) {
1319 _unur_warning(GENTYPE,UNUR_ERR_PAR_SET,"PDF(mode) overflow");
1320 return UNUR_ERR_PAR_SET;
1321 }
1322 GEN->um = sqrt(fm); /* height of rectangle */
1323 }
1324
1325 /* width of rectangle */
1326 vm = DISTR.area / GEN->um;
1327
1328 if (gen->set & SROU_SET_CDFMODE) {
1329 /* cdf at mode known */
1330 GEN->vl = -GEN->Fmode * vm;
1331 GEN->vr = vm + GEN->vl;
1332 GEN->xl = GEN->vl/GEN->um;
1333 GEN->xr = GEN->vr/GEN->um;
1334 }
1335 else {
1336 /* cdf at mode unknown */
1337 GEN->vl = -vm;
1338 GEN->vr = vm;
1339 GEN->xl = GEN->vl/GEN->um;
1340 GEN->xr = GEN->vr/GEN->um;
1341 /* we cannot use universal squeeze */
1342 gen->variant &= ~SROU_VARFLAG_SQUEEZE;
1343 }
1344
1345 /* o.k. */
1346 return UNUR_SUCCESS;
1347
1348 } /* end of _unur_srou_rectangle() */
1349
1350 /*---------------------------------------------------------------------------*/
1351
1352 int
_unur_gsrou_envelope(struct unur_gen * gen)1353 _unur_gsrou_envelope( struct unur_gen *gen )
1354 /*----------------------------------------------------------------------*/
1355 /* compute parameters for universal bounding envelope. */
1356 /* */
1357 /* parameters: */
1358 /* gen ... pointer to generator object */
1359 /* */
1360 /* return: */
1361 /* UNUR_SUCCESS ... on success */
1362 /* error code ... on error */
1363 /*----------------------------------------------------------------------*/
1364 {
1365 double fm; /* PDF at mode */
1366 double vm; /* maximal width of envelope */
1367 double pr; /* p^r */
1368
1369 double p; /* short cuts */
1370 double r = GEN->r;
1371
1372 /* check arguments */
1373 CHECK_NULL( gen, UNUR_ERR_NULL );
1374 COOKIE_CHECK( gen, CK_SROU_GEN, UNUR_ERR_COOKIE );
1375
1376 if (!(gen->set & SROU_SET_PDFMODE)) {
1377 /* compute PDF at mode */
1378 fm = PDF(DISTR.mode);
1379 /* fm must be positive */
1380 if (fm <= 0.) {
1381 _unur_error(gen->genid,UNUR_ERR_GEN_DATA,"PDF(mode) <= 0.");
1382 return UNUR_ERR_GEN_DATA;
1383 }
1384 if (!_unur_isfinite(fm)) {
1385 _unur_warning(GENTYPE,UNUR_ERR_PAR_SET,"PDF(mode) overflow");
1386 return UNUR_ERR_PAR_SET;
1387 }
1388 GEN->um = pow(fm,1./(r+1.)); /* height of envelope */
1389 }
1390
1391 /* maximal width of envelope */
1392 vm = DISTR.area / (GEN->r*GEN->um);
1393
1394 if (gen->set & SROU_SET_CDFMODE) {
1395 /* cdf at mode known */
1396 GEN->vl = -GEN->Fmode * vm;
1397 GEN->vr = vm + GEN->vl;
1398 }
1399 else {
1400 /* cdf at mode unknown */
1401 GEN->vl = -vm;
1402 GEN->vr = vm;
1403 }
1404
1405 /* construction point for bounding curve */
1406 GEN->p = p = 1. - 2.187/pow(r + 5 - 1.28/r, 0.9460 );
1407 pr = pow(p,r);
1408
1409 /* parameters for bounding envelope */
1410 GEN->b = (1. - r * pr/p + (r-1.)*pr) / ((pr-1.)*(pr-1));
1411 GEN->a = -(p-1.)/(pr-1.) - p * GEN->b;
1412 GEN->log_ab = log(GEN->a/(GEN->a+GEN->b));
1413
1414 /* o.k. */
1415 return UNUR_SUCCESS;
1416
1417 } /* end of _unur_gsrou_envelope() */
1418
1419 /*****************************************************************************/
1420 /** Debugging utilities **/
1421 /*****************************************************************************/
1422
1423 /*---------------------------------------------------------------------------*/
1424 #ifdef UNUR_ENABLE_LOGGING
1425 /*---------------------------------------------------------------------------*/
1426
1427 void
_unur_srou_debug_init(const struct unur_gen * gen,int is_reinit)1428 _unur_srou_debug_init( const struct unur_gen *gen, int is_reinit )
1429 /*----------------------------------------------------------------------*/
1430 /* write info about generator into LOG file */
1431 /* */
1432 /* parameters: */
1433 /* gen ... pointer to generator object */
1434 /* is_reinit ... if TRUE the generator has been reinitialized */
1435 /*----------------------------------------------------------------------*/
1436 {
1437 FILE *LOG;
1438
1439 /* check arguments */
1440 CHECK_NULL(gen,RETURN_VOID); COOKIE_CHECK(gen,CK_SROU_GEN,RETURN_VOID);
1441
1442 LOG = unur_get_stream();
1443
1444 fprintf(LOG,"%s:\n",gen->genid);
1445 if (!is_reinit) {
1446 fprintf(LOG,"%s: type = continuous univariate random variates\n",gen->genid);
1447 fprintf(LOG,"%s: method = srou (simple universal ratio-of-uniforms)\n",gen->genid);
1448 }
1449 else
1450 fprintf(LOG,"%s: reinit!\n",gen->genid);
1451 fprintf(LOG,"%s:\n",gen->genid);
1452
1453 _unur_distr_cont_debug( gen->distr, gen->genid );
1454
1455 if (gen->set & SROU_SET_R) {
1456 fprintf(LOG,"%s: Generalized version: r = %g\n",gen->genid,GEN->r);
1457 fprintf(LOG,"%s:\n",gen->genid);
1458
1459 fprintf(LOG,"%s: sampling routine = _unur_gsrou_sample",gen->genid);
1460 if (gen->variant & SROU_VARFLAG_VERIFY)
1461 fprintf(LOG,"_check");
1462 fprintf(LOG,"()\n%s:\n",gen->genid);
1463 }
1464 else {
1465 fprintf(LOG,"%s: Simple version (r = 1) [default]\n",gen->genid);
1466 fprintf(LOG,"%s:\n",gen->genid);
1467
1468 fprintf(LOG,"%s: sampling routine = _unur_srou_sample",gen->genid);
1469 if (gen->variant & SROU_VARFLAG_VERIFY)
1470 fprintf(LOG,"_check");
1471 else if (gen->variant & SROU_VARFLAG_MIRROR)
1472 fprintf(LOG,"_mirror");
1473 fprintf(LOG,"()\n%s:\n",gen->genid);
1474 }
1475
1476 if (gen->set & SROU_SET_CDFMODE)
1477 fprintf(LOG,"%s: F(mode) = %g\n",gen->genid,GEN->Fmode);
1478 else
1479 fprintf(LOG,"%s: F(mode) unknown\n",gen->genid);
1480
1481 if (gen->variant & SROU_VARFLAG_SQUEEZE)
1482 fprintf(LOG,"%s: use universal squeeze\n",gen->genid);
1483 else
1484 fprintf(LOG,"%s: no (universal) squeeze\n",gen->genid);
1485
1486 if (gen->variant & SROU_VARFLAG_MIRROR)
1487 fprintf(LOG,"%s: use mirror principle\n",gen->genid);
1488
1489 fprintf(LOG,"%s:\n",gen->genid);
1490
1491 if (gen->set & SROU_SET_R) {
1492 fprintf(LOG,"%s: Enveloping region:\n",gen->genid);
1493 fprintf(LOG,"%s: um = %g\n",gen->genid,GEN->um);
1494 fprintf(LOG,"%s: vl = %g\n",gen->genid,GEN->vl);
1495 fprintf(LOG,"%s: vr = %g\n",gen->genid,GEN->vr);
1496 fprintf(LOG,"%s: p = %g\n",gen->genid,GEN->p);
1497 fprintf(LOG,"%s: a = %g\n",gen->genid,GEN->a);
1498 fprintf(LOG,"%s: b = %g\n",gen->genid,GEN->b);
1499 }
1500 else {
1501 fprintf(LOG,"%s: Rectangle:\n",gen->genid);
1502 fprintf(LOG,"%s: left upper point = (%g,%g)\n",gen->genid,GEN->vl,GEN->um);
1503 fprintf(LOG,"%s: right upper point = (%g,%g)\n",gen->genid,GEN->vr,GEN->um);
1504 }
1505
1506 fprintf(LOG,"%s:\n",gen->genid);
1507
1508 } /* end of _unur_srou_debug_init() */
1509
1510 /*---------------------------------------------------------------------------*/
1511 #endif /* end UNUR_ENABLE_LOGGING */
1512 /*---------------------------------------------------------------------------*/
1513
1514
1515 /*---------------------------------------------------------------------------*/
1516 #ifdef UNUR_ENABLE_INFO
1517 /*---------------------------------------------------------------------------*/
1518
1519 void
_unur_srou_info(struct unur_gen * gen,int help)1520 _unur_srou_info( struct unur_gen *gen, int help )
1521 /*----------------------------------------------------------------------*/
1522 /* create character string that contains information about the */
1523 /* given generator object. */
1524 /* */
1525 /* parameters: */
1526 /* gen ... pointer to generator object */
1527 /* help ... whether to print additional comments */
1528 /*----------------------------------------------------------------------*/
1529 {
1530 struct unur_string *info = gen->infostr;
1531 struct unur_distr *distr = gen->distr;
1532 int samplesize = 10000;
1533 double h_area, rc;
1534
1535 /* generator ID */
1536 _unur_string_append(info,"generator ID: %s\n\n", gen->genid);
1537
1538 /* distribution */
1539 _unur_string_append(info,"distribution:\n");
1540 _unur_distr_info_typename(gen);
1541 _unur_string_append(info," functions = PDF\n");
1542 _unur_string_append(info," domain = (%g, %g)\n", DISTR.domain[0],DISTR.domain[1]);
1543 _unur_string_append(info," mode = %g %s\n", DISTR.mode,
1544 (distr->set & UNUR_DISTR_SET_MODE_APPROX) ? "[numeric.]" : "");
1545 _unur_string_append(info," area(PDF) = %g\n", DISTR.area);
1546 if (gen->set & SROU_SET_CDFMODE)
1547 _unur_string_append(info," F(mode) = %g\n", GEN->Fmode);
1548 else
1549 _unur_string_append(info," F(mode) = [unknown]\n");
1550
1551 if (help) {
1552 if ( distr->set & UNUR_DISTR_SET_MODE_APPROX )
1553 _unur_string_append(info,"\n[ Hint: %s ]\n",
1554 "You may provide the \"mode\"");
1555 }
1556 _unur_string_append(info,"\n");
1557
1558 /* method */
1559 _unur_string_append(info,"method: SROU (Simple Ratio-Of-Uniforms)\n");
1560 _unur_string_append(info," r = %g %s\n", GEN->r,
1561 (gen->set & SROU_SET_R) ? "[generalized version]" : "");
1562 if (gen->set & SROU_SET_CDFMODE)
1563 _unur_string_append(info," use CDF at mode\n");
1564 if (gen->variant & SROU_VARFLAG_SQUEEZE)
1565 _unur_string_append(info," use squeeze\n");
1566 if (gen->variant & SROU_VARFLAG_MIRROR)
1567 _unur_string_append(info," use mirror principle\n");
1568 _unur_string_append(info,"\n");
1569
1570 /* performance */
1571 _unur_string_append(info,"performance characteristics:\n");
1572 if (gen->set & SROU_SET_R) {
1573 rc = unur_test_count_urn(gen,samplesize,0,NULL)/(2.*samplesize);
1574 _unur_string_append(info," enveloping rectangle = (%g,%g) x (%g,%g)\n",
1575 GEN->vl,GEN->vr, 0.,GEN->um);
1576 _unur_string_append(info," rejection constant = %.2f [approx.]\n", rc);
1577 }
1578 else {
1579 _unur_string_append(info," bounding rectangle = (%g,%g) x (%g,%g)\n",
1580 GEN->vl,GEN->vr, 0.,GEN->um);
1581 h_area = (GEN->vr - GEN->vl) * GEN->um;
1582 _unur_string_append(info," area(hat) = %g\n", h_area);
1583 if (gen->set & SROU_SET_CDFMODE)
1584 rc = 2.;
1585 else
1586 rc = (gen->variant & SROU_VARFLAG_MIRROR) ? 2.829 : 4.;
1587
1588 _unur_string_append(info," rejection constant = %g\n", rc);
1589 }
1590 _unur_string_append(info,"\n");
1591
1592 /* parameters */
1593 if (help) {
1594 _unur_string_append(info,"parameters:\n");
1595 _unur_string_append(info," r = %g %s\n", GEN->r,
1596 (gen->set & SROU_SET_R) ? "" : "[default]");
1597 if (gen->set & SROU_SET_CDFMODE)
1598 _unur_string_append(info," cdfatmode = %g\n", GEN->Fmode);
1599 else
1600 _unur_string_append(info," cdfatmode = [not set]\n");
1601
1602 if (gen->variant & SROU_VARFLAG_SQUEEZE)
1603 _unur_string_append(info," usesqueeze\n");
1604
1605 if (gen->variant & SROU_VARFLAG_MIRROR)
1606 _unur_string_append(info," usemirror\n");
1607
1608 if (gen->variant & SROU_VARFLAG_VERIFY)
1609 _unur_string_append(info," verify = on\n");
1610
1611 _unur_string_append(info,"\n");
1612
1613 /* Not displayed:
1614 int unur_srou_set_pdfatmode( UNUR_PAR *parameters, double fmode );
1615 */
1616 }
1617
1618
1619 /* Hints */
1620 if (help) {
1621 if ( !(gen->set & SROU_SET_CDFMODE))
1622 _unur_string_append(info,"[ Hint: %s ]\n",
1623 "You can set \"cdfatmode\" to reduce the rejection constant.");
1624 _unur_string_append(info,"\n");
1625 }
1626
1627 } /* end of _unur_srou_info() */
1628
1629 /*---------------------------------------------------------------------------*/
1630 #endif /* end UNUR_ENABLE_INFO */
1631 /*---------------------------------------------------------------------------*/
1632