1 /*****************************************************************************
2 * *
3 * UNURAN -- Universal Non-Uniform Random number generator *
4 * *
5 *****************************************************************************
6 * *
7 * FILE: utdr.h *
8 * *
9 * TYPE: continuous univariate random variate *
10 * METHOD: transformed density rejection with three points of contact *
11 * T(x) = -1/sqrt(x) (T_c with c= -1/2) *
12 * *
13 * DESCRIPTION: *
14 * Given PDF and mode of a T-concave distribution *
15 * produce a value x consistent with its density *
16 * *
17 * REQUIRED: pointer to the density, mode of the density *
18 * *
19 * PARAMETERS: *
20 * double il, ir ... left and right boundary of domain *
21 * (default: +/- INFINITY) *
22 * double pdf_area ... area below PDF (need not be 1) *
23 * (default: 1.) *
24 * double *pdf_param ... parameters of PDF *
25 * (default: NULL) *
26 * int n_pdf_param ... number of parameters of PDF *
27 * (default: 0) *
28 * double c_factor ... constant for choosing the constr. points *
29 * (default: 0.664) *
30 * double delta_factor ... constant for approx. first derivative *
31 * (default: 0.00001) *
32 * *
33 *****************************************************************************
34 * *
35 * Copyright (c) 2000-2006 Wolfgang Hoermann and Josef Leydold *
36 * Department of Statistics and Mathematics, WU Wien, Austria *
37 * *
38 * This program is free software; you can redistribute it and/or modify *
39 * it under the terms of the GNU General Public License as published by *
40 * the Free Software Foundation; either version 2 of the License, or *
41 * (at your option) any later version. *
42 * *
43 * This program is distributed in the hope that it will be useful, *
44 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
45 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
46 * GNU General Public License for more details. *
47 * *
48 * You should have received a copy of the GNU General Public License *
49 * along with this program; if not, write to the *
50 * Free Software Foundation, Inc., *
51 * 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA *
52 * *
53 *****************************************************************************
54 * *
55 * REFERENCES: *
56 * [1] Hoermann W. (1995): A rejection technique for sampling from *
57 * T-concave distributions, ACM TOMS 21, p. 182-193 *
58 * (see Algorithm UTDR) *
59 * *
60 *****************************************************************************
61 * *
62 * ..... beschreibung .... *
63 * *
64 * *
65 *****************************************************************************/
66
67 /*---------------------------------------------------------------------------*/
68
69 #include <unur_source.h>
70 #include <distr/distr.h>
71 #include <distr/distr_source.h>
72 #include <distr/cont.h>
73 #include <urng/urng.h>
74 #include "unur_methods_source.h"
75 #include "x_gen_source.h"
76 #include "utdr.h"
77 #include "utdr_struct.h"
78
79 #ifdef UNUR_ENABLE_INFO
80 # include <tests/unuran_tests.h>
81 #endif
82
83 /*---------------------------------------------------------------------------*/
84 /* Variants */
85
86 #define UTDR_VARFLAG_VERIFY 0x01u /* flag for verifying mode */
87
88 /*---------------------------------------------------------------------------*/
89 /* Debugging flags */
90 /* bit 01 ... pameters and structure of generator (do not use here) */
91 /* bits 02-12 ... setup */
92 /* bits 13-24 ... adaptive steps */
93 /* bits 25-32 ... trace sampling */
94
95 #define UTDR_DEBUG_REINIT 0x00000010u /* print parameters after reinit */
96
97 /*---------------------------------------------------------------------------*/
98 /* Flags for logging set calls */
99
100 #define UTDR_SET_CPFACTOR 0x001u
101 #define UTDR_SET_DELTA 0x002u
102 #define UTDR_SET_PDFMODE 0x004u /* PDF at mode is set */
103
104 /*---------------------------------------------------------------------------*/
105
106 #define GENTYPE "UTDR" /* type of generator */
107
108 /*---------------------------------------------------------------------------*/
109
110 static struct unur_gen *_unur_utdr_init( struct unur_par *par );
111 /*---------------------------------------------------------------------------*/
112 /* Initialize new generator. */
113 /*---------------------------------------------------------------------------*/
114
115 static int _unur_utdr_reinit( struct unur_gen *gen );
116 /*---------------------------------------------------------------------------*/
117 /* Reinitialize generator. */
118 /*---------------------------------------------------------------------------*/
119
120 static struct unur_gen *_unur_utdr_create( struct unur_par *par );
121 /*---------------------------------------------------------------------------*/
122 /* create new (almost empty) generator object. */
123 /*---------------------------------------------------------------------------*/
124
125 static int _unur_utdr_check_par( struct unur_gen *gen );
126 /*---------------------------------------------------------------------------*/
127 /* Check parameters of given distribution and method */
128 /*---------------------------------------------------------------------------*/
129
130 static struct unur_gen *_unur_utdr_clone( const struct unur_gen *gen );
131 /*---------------------------------------------------------------------------*/
132 /* copy (clone) generator object. */
133 /*---------------------------------------------------------------------------*/
134
135 static void _unur_utdr_free( struct unur_gen *gen);
136 /*---------------------------------------------------------------------------*/
137 /* destroy generator object. */
138 /*---------------------------------------------------------------------------*/
139
140 static double _unur_utdr_sample( struct unur_gen *generator );
141 static double _unur_utdr_sample_check( struct unur_gen *generator );
142 /*---------------------------------------------------------------------------*/
143 /* sample from generator */
144 /*---------------------------------------------------------------------------*/
145
146 static int _unur_utdr_hat( struct unur_gen *gen );
147 /*---------------------------------------------------------------------------*/
148 /* compute hat and squeezes. */
149 /*---------------------------------------------------------------------------*/
150
151 #ifdef UNUR_ENABLE_LOGGING
152 /*---------------------------------------------------------------------------*/
153 /* the following functions print debugging information on output stream, */
154 /* i.e., into the LOG file if not specified otherwise. */
155 /*---------------------------------------------------------------------------*/
156
157 static void _unur_utdr_debug_init( const struct unur_gen *gen,
158 double ttly, double ttlys, double ttry, double ttrys,
159 double cfac, int setupok, double c );
160 /*---------------------------------------------------------------------------*/
161 /* print after generator has been initialized has completed. */
162 /*---------------------------------------------------------------------------*/
163 #endif
164
165 #ifdef UNUR_ENABLE_INFO
166 static void _unur_utdr_info( struct unur_gen *gen, int help );
167 /*---------------------------------------------------------------------------*/
168 /* create info string. */
169 /*---------------------------------------------------------------------------*/
170 #endif
171
172 /*---------------------------------------------------------------------------*/
173 /* abbreviations */
174
175 #define DISTR_IN distr->data.cont /* data for distribution object */
176
177 #define PAR ((struct unur_utdr_par*)par->datap) /* data for parameter object */
178 #define GEN ((struct unur_utdr_gen*)gen->datap) /* data for generator object */
179 #define DISTR gen->distr->data.cont /* data for distribution in generator object */
180
181 #define BD_LEFT domain[0] /* left boundary of domain of distribution */
182 #define BD_RIGHT domain[1] /* right boundary of domain of distribution */
183
184 #define SAMPLE gen->sample.cont /* pointer to sampling routine */
185
186 #define PDF(x) _unur_cont_PDF((x),(gen->distr)) /* call to PDF */
187
188 /*---------------------------------------------------------------------------*/
189
190 #define _unur_utdr_getSAMPLE(gen) \
191 ( ((gen)->variant & UTDR_VARFLAG_VERIFY) \
192 ? _unur_utdr_sample_check : _unur_utdr_sample )
193
194 /*---------------------------------------------------------------------------*/
195
196 /*****************************************************************************/
197 /** Public: User Interface (API) **/
198 /*****************************************************************************/
199
200 struct unur_par *
unur_utdr_new(const struct unur_distr * distr)201 unur_utdr_new( const struct unur_distr *distr )
202 /*----------------------------------------------------------------------*/
203 /* get default parameters */
204 /* */
205 /* parameters: */
206 /* distr ... pointer to distribution object */
207 /* */
208 /* return: */
209 /* default parameters (pointer to structure) */
210 /* */
211 /* error: */
212 /* return NULL */
213 /* */
214 /* comment: */
215 /* if the area below the PDF is not close to 1 it is necessary to */
216 /* set pdf_area to an approximate value of its area (+/- 30 % is ok). */
217 /*----------------------------------------------------------------------*/
218 {
219 struct unur_par *par;
220
221 /* check arguments */
222 _unur_check_NULL( GENTYPE,distr,NULL );
223
224 /* check distribution */
225 if (distr->type != UNUR_DISTR_CONT) {
226 _unur_error(GENTYPE,UNUR_ERR_DISTR_INVALID,""); return NULL; }
227 COOKIE_CHECK(distr,CK_DISTR_CONT,NULL);
228
229 if (DISTR_IN.pdf == NULL) {
230 _unur_error(GENTYPE,UNUR_ERR_DISTR_REQUIRED,"PDF");
231 return NULL;
232 }
233
234 /**TODOWH:
235 Kann man eigentlich einen bound fuer die area auch
236 im Distribution-objekt eingeben?????
237 so dass es klar ist, dass da nur ein bound ist?
238
239 JL: ich habe ich SSR die moeglichkeit eingebaut, dass nur eine
240 obere schranke eingegeben wurde. dass einzige problem dabei ist,
241 dass dann der squeeze nicht mehr verwendet werden kann.
242 man kann aber die verwendung des squeezes (der eigentlich nur bei
243 sehr teuren dichten etwas bringt, mittels
244 int unur_ssr_set_usesqueeze( UNUR_PAR *parameters, int usesqueeze );
245 abschalten.
246
247 waere das ein brauchbarer weg fuer Dich?
248
249 alternativ koennten wir auch ein unur_distr_set_approx_pdfarea()
250 einbauen. die unterscheidung zwischen exactem wert und naehrung kann
251 man mittels set-flag machen. (irgendwie ist mir bis jetzt nie die idee
252 gekommen, dass so zu machen.) ist das so besser?
253 **/
254
255 /* allocate structure */
256 par = _unur_par_new( sizeof(struct unur_utdr_par) );
257 COOKIE_SET(par,CK_UTDR_PAR);
258
259 /* copy input */
260 par->distr = distr; /* pointer to distribution object */
261
262 /* set default values */
263 PAR->c_factor = 0.664;
264 /* optimal value for the normal distribution, which is good for
265 all bell-shaped densities. The minimax approach for that
266 transformation has c_factor=2. */
267 PAR->delta_factor = 0.00001;
268 /* constant for choosing delta to replace the tangent.
269 default should not be changed if used doubles have at least
270 10 decimal digits precision. */
271
272 PAR->fm = -1.; /* PDF at mode (unknown) */
273 PAR->hm = -1.; /* square of PDF at mode (unknown) */
274
275 par->method = UNUR_METH_UTDR; /* method */
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_utdr_init;
285
286 return par;
287
288 } /* end of unur_utdr_new() */
289
290 /*****************************************************************************/
291
292 int
unur_utdr_set_pdfatmode(UNUR_PAR * par,double fmode)293 unur_utdr_set_pdfatmode( UNUR_PAR *par, double fmode )
294 /*----------------------------------------------------------------------*/
295 /* Set PDF at mode. if set, the PDF at the mode is never changed. */
296 /* */
297 /* parameters: */
298 /* par ... pointer to parameter for building generator object */
299 /* fmode ... PDF 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, UTDR );
309
310 /* check new parameter for generator */
311 if (fmode <= 0.) {
312 _unur_warning(GENTYPE,UNUR_ERR_PAR_SET,"PDF(mode)");
313 return UNUR_ERR_PAR_SET;
314 }
315
316 /* store date */
317 PAR->fm = fmode; /* PDF at mode */
318 PAR->hm = -1./sqrt(fmode); /* transformed PDF at mode */
319
320 /* changelog */
321 par->set |= UTDR_SET_PDFMODE;
322
323 return UNUR_SUCCESS;
324
325 } /* end of unur_utdr_set_pdfatmode() */
326
327 /*---------------------------------------------------------------------------*/
328
329 int
unur_utdr_set_cpfactor(struct unur_par * par,double cp_factor)330 unur_utdr_set_cpfactor( struct unur_par *par, double cp_factor )
331 /*----------------------------------------------------------------------*/
332 /* set factor for position of left and right construction point */
333 /* */
334 /* parameters: */
335 /* par ... pointer to parameter for building generator object */
336 /* cp_factor ... factor */
337 /* */
338 /* return: */
339 /* UNUR_SUCCESS ... on success */
340 /* error code ... on error */
341 /*----------------------------------------------------------------------*/
342 {
343 /* check arguments */
344 _unur_check_NULL( GENTYPE, par, UNUR_ERR_NULL );
345 _unur_check_par_object( par, UTDR );
346
347 /* check new parameter for generator */
348 /** TODO: welche werte fuer c sind zulaessig / sinnvoll ?
349 zulaessig ist jedes c>0, man koennte damit falsche Flaechenangaben kompensieren.
350 Wenn area genau bekannt ist, ist ein c > 2 (2 ist der minimax approach) so weit
351 ich weiss nie sinnvoll. Ich denke aber, das sollte man besser nicht prinzipiell
352 verbieten, hoechstens eine warnung.**/
353 if (cp_factor <= 0.) {
354 _unur_warning(GENTYPE,UNUR_ERR_PAR_SET,"cp-factor <= 0");
355 return UNUR_ERR_PAR_SET;
356 }
357
358 if (cp_factor > 2.1)
359 _unur_warning(GENTYPE,UNUR_ERR_PAR_SET,"cp-factor > 2 not recommended. skip");
360
361 /* store date */
362 PAR->c_factor = cp_factor;
363
364 /* changelog */
365 par->set |= UTDR_SET_CPFACTOR;
366
367 return UNUR_SUCCESS;
368
369 } /* end of unur_utdr_set_cpfactor() */
370
371 /*---------------------------------------------------------------------------*/
372
373 int
unur_utdr_set_deltafactor(struct unur_par * par,double delta)374 unur_utdr_set_deltafactor( struct unur_par *par, double delta )
375 /*----------------------------------------------------------------------*/
376 /* set factor for replacing tangents by secants */
377 /* (which then are move above the transformed density) */
378 /* */
379 /* parameters: */
380 /* par ... pointer to parameter for building generator object */
381 /* delta ... delta-factor */
382 /* */
383 /* return: */
384 /* UNUR_SUCCESS ... on success */
385 /* error code ... on error */
386 /*----------------------------------------------------------------------*/
387 {
388 /* check arguments */
389 _unur_check_NULL( GENTYPE, par, UNUR_ERR_NULL );
390 _unur_check_par_object( par, UTDR );
391
392 /* check new parameter for generator */
393 /** TODO: welche werte fuer delta sind zulaessig / sinnvoll ?
394 von einem rein numerischen Gesichtspunkt kann man sagen: Delta_factor
395 soll nicht kleiner sein als sqrt(epsilon) (mit epsilon meine ich die
396 kleinste Zahl die einen Unterschied zwischen 1 und 1+epsilon erkennen
397 laesst). Dann ist garantiert, dass wir durch Ausloeschung nicht mehr
398 als die haelfte der mantisse verlieren
399 Von Algorithmus ueberlegungen her sollte delta>=xl*deltafactor sicher
400 nicht groesser sein als ein Bruchteil (wohl hoechstens 1/10) zwischen
401 mode und Design-punkt, weil sonst die Anahmewahrscheinlichkeit schlechter
402 wird. Das ergibt sicher probleme fuer Dichten wo der mode sehr gross
403 und der skale-parameter im Vergleich dazu klein ist.
404 Ich habe den Algorithmus so geaendert, dass er, falls es dieses
405 Problem gibt, eine (roundoff-problem) Warnung ausspuckt, und als
406 Delta 1/100*|designpoint-mode| nimmt Empirisch habe ich dafuer
407 festgestellt, dass damit die Flaeche unterm Hut Kaum mehr als 1%
408 zunimmt und das ist ja nicht schlimm**/
409 if (delta <= 0.) {
410 _unur_warning(GENTYPE,UNUR_ERR_PAR_SET,"delta <= 0");
411 return UNUR_ERR_PAR_SET;
412 }
413 if (delta > 0.1) {
414 _unur_warning(GENTYPE,UNUR_ERR_PAR_SET,"delta must be small");
415 return UNUR_ERR_PAR_SET;
416 }
417
418 /* store date */
419 PAR->delta_factor = delta;
420
421 /* changelog */
422 par->set |= UTDR_SET_DELTA;
423
424 return UNUR_SUCCESS;
425
426 } /* end of unur_utdr_set_delta() */
427
428 /*---------------------------------------------------------------------------*/
429
430 int
unur_utdr_set_verify(struct unur_par * par,int verify)431 unur_utdr_set_verify( struct unur_par *par, int verify )
432 /*----------------------------------------------------------------------*/
433 /* turn verifying of algorithm while sampling on/off */
434 /* */
435 /* parameters: */
436 /* par ... pointer to parameter for building generator object */
437 /* verify ... 0 = no verifying, !0 = verifying */
438 /* */
439 /* return: */
440 /* UNUR_SUCCESS ... on success */
441 /* error code ... on error */
442 /* */
443 /* comment: */
444 /* no verifying is the default */
445 /*----------------------------------------------------------------------*/
446 {
447 /* check arguments */
448 _unur_check_NULL( GENTYPE, par, UNUR_ERR_NULL );
449 _unur_check_par_object( par, UTDR );
450
451 /* we use a bit in variant */
452 par->variant = (verify) ? (par->variant | UTDR_VARFLAG_VERIFY) : (par->variant & (~UTDR_VARFLAG_VERIFY));
453
454 /* o.k. */
455 return UNUR_SUCCESS;
456
457 } /* end of unur_utdr_set_verify() */
458
459 /*---------------------------------------------------------------------------*/
460
461 int
unur_utdr_chg_verify(struct unur_gen * gen,int verify)462 unur_utdr_chg_verify( struct unur_gen *gen, int verify )
463 /*----------------------------------------------------------------------*/
464 /* turn verifying of algorithm while sampling on/off */
465 /* */
466 /* parameters: */
467 /* gen ... pointer to generator object */
468 /* verify ... 0 = no verifying, !0 = verifying */
469 /* */
470 /* return: */
471 /* UNUR_SUCCESS ... on success */
472 /* error code ... on error */
473 /* */
474 /* comment: */
475 /* no verifying is the default */
476 /*----------------------------------------------------------------------*/
477 {
478 /* check input */
479 _unur_check_NULL( GENTYPE, gen, UNUR_ERR_NULL );
480 _unur_check_gen_object( gen, UTDR, UNUR_ERR_GEN_INVALID );
481
482 /* we must not change this switch when sampling has been disabled by
483 using a pointer to the error producing routine */
484 if (SAMPLE == _unur_sample_cont_error)
485 return UNUR_FAILURE;
486
487 if (verify)
488 /* turn verify mode on */
489 gen->variant |= UTDR_VARFLAG_VERIFY;
490 else
491 /* turn verify mode off */
492 gen->variant &= ~UTDR_VARFLAG_VERIFY;
493
494 SAMPLE = _unur_utdr_getSAMPLE(gen);
495
496 /* o.k. */
497 return UNUR_SUCCESS;
498
499 } /* end of unur_utdr_chg_verify() */
500
501 /*****************************************************************************/
502
503 int
unur_utdr_chg_pdfatmode(struct unur_gen * gen,double fmode)504 unur_utdr_chg_pdfatmode( struct unur_gen *gen, double fmode )
505 /*----------------------------------------------------------------------*/
506 /* change value of PDF at mode */
507 /* */
508 /* parameters: */
509 /* gen ... pointer to generator object */
510 /* fmode ... PDF at mode */
511 /* */
512 /* return: */
513 /* UNUR_SUCCESS ... on success */
514 /* error code ... on error */
515 /*----------------------------------------------------------------------*/
516 {
517 /* check arguments */
518 _unur_check_NULL( GENTYPE, gen, UNUR_ERR_NULL );
519 _unur_check_gen_object( gen, UTDR, UNUR_ERR_GEN_INVALID );
520
521 /* check new parameter for generator */
522 if (fmode <= 0.) {
523 _unur_warning(gen->genid,UNUR_ERR_PAR_SET,"PDF(mode)");
524 return UNUR_ERR_PAR_SET;
525 }
526
527 /* store date */
528 GEN->fm = fmode; /* PDF at mode */
529 GEN->hm = -1./sqrt(fmode); /* transformed PDF at mode */
530
531 /* changelog */
532 gen->set |= UTDR_SET_PDFMODE;
533
534 /* o.k. */
535 return UNUR_SUCCESS;
536 } /* end of unur_utdr_chg_pdfatmode() */
537
538
539 /*****************************************************************************/
540 /** Private **/
541 /*****************************************************************************/
542
543 struct unur_gen *
_unur_utdr_init(struct unur_par * par)544 _unur_utdr_init( struct unur_par *par )
545 /*----------------------------------------------------------------------*/
546 /* initialize new generator */
547 /* */
548 /* parameters: */
549 /* par ... pointer to paramters for building generator object */
550 /* */
551 /* return: */
552 /* pointer to generator object */
553 /* */
554 /* error: */
555 /* return NULL */
556 /*----------------------------------------------------------------------*/
557 {
558 struct unur_gen *gen;
559
560 /* check arguments */
561 _unur_check_NULL( GENTYPE,par,NULL );
562
563 /* check input */
564 if ( par->method != UNUR_METH_UTDR ) {
565 _unur_error(GENTYPE,UNUR_ERR_PAR_INVALID,"");
566 return NULL; }
567 COOKIE_CHECK(par,CK_UTDR_PAR,NULL);
568
569 /* create a new empty generator object */
570 gen = _unur_utdr_create(par);
571 _unur_par_free(par);
572 if (!gen) return NULL;
573
574 /* check parameters */
575 if (_unur_utdr_check_par(gen) != UNUR_SUCCESS) {
576 _unur_utdr_free(gen); return NULL;
577 }
578
579 /* create hat and squeeze (setup procedure) */
580 if ( _unur_utdr_hat(gen)!=UNUR_SUCCESS ) {
581 _unur_utdr_free(gen); return NULL;
582 }
583
584 /* hat successfully created */
585 return gen;
586
587 } /* end of _unur_utdr_init() */
588
589 /*---------------------------------------------------------------------------*/
590
591 int
_unur_utdr_reinit(struct unur_gen * gen)592 _unur_utdr_reinit( struct unur_gen *gen )
593 /*----------------------------------------------------------------------*/
594 /* re-initialize (existing) generator. */
595 /* */
596 /* parameters: */
597 /* gen ... pointer to generator object */
598 /* */
599 /* return: */
600 /* UNUR_SUCCESS ... on success */
601 /* error code ... on error */
602 /*----------------------------------------------------------------------*/
603 {
604 int rcode;
605
606 /* check parameters */
607 if ( (rcode = _unur_utdr_check_par(gen)) != UNUR_SUCCESS)
608 return rcode;
609
610 /* update left and right boundary for algorithm */
611 GEN->il = DISTR.BD_LEFT;
612 GEN->ir = DISTR.BD_RIGHT;
613
614 /* (re)set sampling routine */
615 SAMPLE = _unur_utdr_getSAMPLE(gen);
616
617 /* compute universal bounding rectangle */
618 return _unur_utdr_hat( gen );
619 } /* end of _unur_utdr_reinit() */
620
621 /*---------------------------------------------------------------------------*/
622
623 struct unur_gen *
_unur_utdr_create(struct unur_par * par)624 _unur_utdr_create( struct unur_par *par )
625 /*----------------------------------------------------------------------*/
626 /* allocate memory for generator */
627 /* */
628 /* parameters: */
629 /* par ... pointer to parameter for building generator object */
630 /* */
631 /* return: */
632 /* pointer to (empty) generator object with default settings */
633 /* */
634 /* error: */
635 /* return NULL */
636 /*----------------------------------------------------------------------*/
637 {
638 struct unur_gen *gen;
639
640 /* check arguments */
641 CHECK_NULL(par,NULL); COOKIE_CHECK(par,CK_UTDR_PAR,NULL);
642
643 /* create new generic generator object */
644 gen = _unur_generic_create( par, sizeof(struct unur_utdr_gen) );
645
646 /* magic cookies */
647 COOKIE_SET(gen,CK_UTDR_GEN);
648
649 /* set generator identifier */
650 gen->genid = _unur_set_genid(GENTYPE);
651
652 /* routines for sampling and destroying generator */
653 SAMPLE = _unur_utdr_getSAMPLE(gen);
654 gen->destroy = _unur_utdr_free;
655 gen->clone = _unur_utdr_clone;
656 gen->reinit = _unur_utdr_reinit;
657
658 /* copy some parameters into generator object */
659 GEN->il = DISTR.BD_LEFT; /* left boundary of domain */
660 GEN->ir = DISTR.BD_RIGHT; /* right boundary of domain */
661 GEN->fm = PAR->fm; /* PDF at mode */
662 GEN->hm = PAR->hm; /* square root of PDF at mode */
663 GEN->c_factor = PAR->c_factor;
664 GEN->delta_factor = PAR->delta_factor;
665
666 /* initialize parameters */
667 /** TODO !!! **/
668 /** ist das wirklich so noetig ?? **/
669 GEN->vollc = 0.;
670 GEN->volcompl = 0.;
671 GEN->voll = 0.;
672 GEN->al = 0.;
673 GEN->ar = 0.;
674 GEN->col = 0.;
675 GEN->cor = 0.;
676 GEN->sal = 0.;
677 GEN->sar = 0.;
678 GEN->bl = 0.;
679 GEN->br = 0.;
680 GEN->ttlx = 0.;
681 GEN->ttrx = 0.;
682 GEN->brblvolc = 0.;
683 GEN->drar = 0.;
684 GEN->dlal = 0.;
685 GEN->ooar2 = 0.;
686 GEN->ooal2 = 0.;
687 /* constants of the hat and for generation*/
688
689 #ifdef UNUR_ENABLE_INFO
690 /* set function for creating info string */
691 gen->info = _unur_utdr_info;
692 #endif
693
694 /* return pointer to (almost empty) generator object */
695 return gen;
696
697 } /* end of _unur_utdr_create() */
698
699 /*---------------------------------------------------------------------------*/
700
701 int
_unur_utdr_check_par(struct unur_gen * gen)702 _unur_utdr_check_par( struct unur_gen *gen )
703 /*----------------------------------------------------------------------*/
704 /* check parameters of given distribution and method */
705 /* */
706 /* parameters: */
707 /* gen ... pointer to generator object */
708 /* */
709 /* return: */
710 /* UNUR_SUCCESS ... on success */
711 /* error code ... on error */
712 /*----------------------------------------------------------------------*/
713 {
714
715 /* check for required data: mode */
716 if (!(gen->distr->set & UNUR_DISTR_SET_MODE)) {
717 _unur_warning(GENTYPE,UNUR_ERR_DISTR_REQUIRED,"mode: try finding it (numerically)");
718 if (unur_distr_cont_upd_mode(gen->distr)!=UNUR_SUCCESS) {
719 _unur_error(GENTYPE,UNUR_ERR_DISTR_REQUIRED,"mode");
720 return UNUR_ERR_DISTR_REQUIRED;
721 }
722 }
723
724 /* check for required data: area */
725 if (!(gen->distr->set & UNUR_DISTR_SET_PDFAREA)) {
726 if (unur_distr_cont_upd_pdfarea(gen->distr)!=UNUR_SUCCESS) {
727 _unur_error(GENTYPE,UNUR_ERR_DISTR_REQUIRED,"area below PDF");
728 return UNUR_ERR_DISTR_REQUIRED;
729 }
730 }
731
732 /* mode must be in domain */
733 if ( (DISTR.mode < DISTR.BD_LEFT) ||
734 (DISTR.mode > DISTR.BD_RIGHT) ) {
735 /* there is something wrong.
736 assume: user has change domain without changing mode.
737 but then, she probably has not updated area and is to large */
738 _unur_warning(GENTYPE,UNUR_ERR_GEN_DATA,"area and/or CDF at mode");
739 DISTR.mode = _unur_max(DISTR.mode,DISTR.BD_LEFT);
740 DISTR.mode = _unur_min(DISTR.mode,DISTR.BD_RIGHT);
741 }
742
743 return UNUR_SUCCESS;
744 } /* end of _unur_utdr_check_par() */
745
746 /*---------------------------------------------------------------------------*/
747
748 struct unur_gen *
_unur_utdr_clone(const struct unur_gen * gen)749 _unur_utdr_clone( const struct unur_gen *gen )
750 /*----------------------------------------------------------------------*/
751 /* copy (clone) generator object */
752 /* */
753 /* parameters: */
754 /* gen ... pointer to generator object */
755 /* */
756 /* return: */
757 /* pointer to clone of generator object */
758 /* */
759 /* error: */
760 /* return NULL */
761 /*----------------------------------------------------------------------*/
762 {
763 #define CLONE ((struct unur_utdr_gen*)clone->datap)
764
765 struct unur_gen *clone;
766
767 /* check arguments */
768 CHECK_NULL(gen,NULL); COOKIE_CHECK(gen,CK_UTDR_GEN,NULL);
769
770 /* create generic clone */
771 clone = _unur_generic_clone( gen, GENTYPE );
772
773 return clone;
774
775 #undef CLONE
776 } /* end of _unur_utdr_clone() */
777
778 /*---------------------------------------------------------------------------*/
779
780 void
_unur_utdr_free(struct unur_gen * gen)781 _unur_utdr_free( struct unur_gen *gen )
782 /*----------------------------------------------------------------------*/
783 /* deallocate generator object */
784 /* */
785 /* parameters: */
786 /* gen ... pointer to generator object */
787 /*----------------------------------------------------------------------*/
788 {
789
790 /* check arguments */
791 if( !gen ) /* nothing to do */
792 return;
793
794 /* check input */
795 if ( gen->method != UNUR_METH_UTDR ) {
796 _unur_warning(gen->genid,UNUR_ERR_GEN_INVALID,"");
797 return; }
798 COOKIE_CHECK(gen,CK_UTDR_GEN,RETURN_VOID);
799
800 /* we cannot use this generator object any more */
801 SAMPLE = NULL; /* make sure to show up a programming error */
802
803 /* free memory */
804 _unur_generic_free(gen);
805
806 } /* end of _unur_utdr_free() */
807
808 /*****************************************************************************/
809
810 double
_unur_utdr_sample(struct unur_gen * gen)811 _unur_utdr_sample( struct unur_gen *gen )
812 /*----------------------------------------------------------------------*/
813 /* sample from generator */
814 /* */
815 /* parameters: */
816 /* gen ... pointer to generator object */
817 /* */
818 /* return: */
819 /* double (sample from random variate) */
820 /* */
821 /* error: */
822 /* return INFINITY */
823 /*----------------------------------------------------------------------*/
824 {
825 double u,v,x,help,linx;
826
827 /* check arguments */
828 CHECK_NULL(gen,INFINITY); COOKIE_CHECK(gen,CK_UTDR_GEN,INFINITY);
829
830 while (1) {
831 /*2*/
832 u = _unur_call_urng(gen->urng) * GEN->volcompl;
833 /*2.1*/
834 if (u <= GEN->voll) {
835 u = GEN->voll-u; /*added to ensure inversion for the hat-generation*/
836 x = -GEN->dlal+GEN->ooal2/(u-GEN->col);
837 help = GEN->al*(u-GEN->col);
838 linx = help*help;
839 }
840 else {
841 if (u <= GEN->vollc) {
842 x = (u-GEN->voll) * GEN->brblvolc + GEN->bl;
843 linx = GEN->fm;
844 }
845 else {
846 x = - GEN->drar - GEN->ooar2 / (u-GEN->vollc - GEN->cor);
847 help = GEN->ar * (u-GEN->vollc - GEN->cor);
848 linx = help*help;
849 }
850 }
851 /*2.2*/
852 v = _unur_call_urng(gen->urng) * linx;
853 /*2.3*/
854 if (x<DISTR.mode) {
855 if (x >= GEN->ttlx) {
856 help = GEN->hm - (DISTR.mode - x) * GEN->sal;
857 if (v * help * help <= 1.) return x;
858 }
859 }
860 else {
861 if (x <= GEN->ttrx) {
862 help = GEN->hm - (DISTR.mode - x) * GEN->sar;
863 if (v * help * help <= 1.) return x;
864 }
865 }
866 if (v <= PDF(x)) return x;
867 }
868
869 } /* end of _unur_utdr_sample() */
870
871 /*---------------------------------------------------------------------------*/
872
873 double
_unur_utdr_sample_check(struct unur_gen * gen)874 _unur_utdr_sample_check( struct unur_gen *gen )
875 /*----------------------------------------------------------------------*/
876 /* sample from generator */
877 /* */
878 /* parameters: */
879 /* gen ... pointer to generator object */
880 /* */
881 /* return: */
882 /* double (sample from random variate) */
883 /* */
884 /* error: */
885 /* return INFINITY */
886 /*----------------------------------------------------------------------*/
887 {
888 double u,v,x,help,linx,pdfx,squeezex;
889
890 /* check arguments */
891 CHECK_NULL(gen,INFINITY); COOKIE_CHECK(gen,CK_UTDR_GEN,INFINITY);
892
893 while (1) {
894 /*2*/
895 u = _unur_call_urng(gen->urng) * GEN->volcompl;
896 /*2.1*/
897 if (u <= GEN->voll) {
898 u = GEN->voll-u; /* added to ensure inversion for the hat-generation */
899 x = -GEN->dlal+GEN->ooal2/(u-GEN->col);
900 help = GEN->al*(u-GEN->col);
901 linx = help*help;
902 }
903 else {
904 if (u <= GEN->vollc) {
905 x = (u-GEN->voll) * GEN->brblvolc + GEN->bl;
906 linx = GEN->fm;
907 }
908 else {
909 x = - GEN->drar - GEN->ooar2 / (u-GEN->vollc - GEN->cor);
910 help = GEN->ar * (u-GEN->vollc - GEN->cor);
911 linx = help*help;
912 }
913 }
914 /*2.2*/
915 v = _unur_call_urng(gen->urng) * linx;
916 /*2.3*/
917 squeezex=0.;
918 if (x<DISTR.mode) {
919 if (x >= GEN->ttlx) {
920 help = GEN->hm - (DISTR.mode - x) * GEN->sal;
921 squeezex=1./(help*help);
922 /* if (v * help * help <= 1.) return x;*/
923 }
924 }
925 else {
926 if (x <= GEN->ttrx) {
927 help = GEN->hm - (DISTR.mode - x) * GEN->sar;
928 squeezex=1./(help*help);
929 /* if (v * help * help <= 1.) return x; */
930 }
931 }
932 /*evaluate density-function*/
933 pdfx=PDF(x);
934
935 /* verify hat function */
936 if(_unur_FP_less(linx,pdfx))
937 { _unur_error(gen->genid,UNUR_ERR_GEN_CONDITION,"PDF(x) > hat(x)");
938 _unur_log_printf(gen->genid,__FILE__,__LINE__,"x %e PDF(x) %e hat(x) %e squeeze(x) %e", \
939 x,pdfx,linx,squeezex );
940 }
941 /* verify squeeze function */
942 if(_unur_FP_less(pdfx,squeezex))
943 { _unur_error(gen->genid,UNUR_ERR_GEN_CONDITION,"PDF(x) < squeeze(x)");
944 _unur_log_printf(gen->genid,__FILE__,__LINE__,"x %e PDF(x) %e hat(x) %e squeeze(x) %e", \
945 x,pdfx,linx,squeezex );
946 }
947 if (v <= PDF(x)) return x;
948 }
949
950 } /* end of _unur_utdr_sample_check() */
951
952 /*****************************************************************************/
953 /** Auxilliary Routines **/
954 /*****************************************************************************/
955
956 /** TODO gibts da schon eine UNURAN constante? ich hab sie nicht gefunden!! **/
957 #define SMALL_VAL 1.e-50
958 /** JL: was soll die leisten?
959 wenn sie gebraucht wird, wuerde ich sie mit einer entsprenden beschreibung
960 in unuran_config.h stellen.
961 **/
962
963 int
_unur_utdr_hat(struct unur_gen * gen)964 _unur_utdr_hat( struct unur_gen *gen )
965 /*----------------------------------------------------------------------*/
966 /* compute hat and squeeze */
967 /* */
968 /* parameters: */
969 /* gen ... pointer to generator object */
970 /* */
971 /* return: */
972 /* pointer to generator object */
973 /* */
974 /* return: */
975 /* UNUR_SUCCESS ... on success */
976 /* error code ... on error */
977 /*----------------------------------------------------------------------*/
978 {
979 double fm;
980
981 int setupok=1;
982 double c,cfac,volc,volr,ttly,ttlys,ttry,ttrys,dl,dr,delta,delta1,delta2,pdfx;
983
984 /* check arguments */
985 CHECK_NULL( gen, UNUR_ERR_NULL );
986 COOKIE_CHECK( gen, CK_UTDR_GEN, UNUR_ERR_COOKIE );
987
988 /* compute PDF at mode (if not given by user) */
989 if (!(gen->set & UTDR_SET_PDFMODE)) {
990 fm = PDF(DISTR.mode);
991 /* fm must be positive */
992 if (fm <= 0.) {
993 _unur_warning(gen->genid,UNUR_ERR_GEN_DATA,"PDF(mode) <= 0.");
994 return UNUR_ERR_GEN_DATA;
995 }
996 /* step 1.0 of algorithm UTDR */
997 GEN->fm = fm; /* PDF at mode */
998 GEN->hm = -1/sqrt(fm); /* transformed PDF at mode */
999 }
1000
1001 /** TODO: inititialisieren notwendig ?? **/
1002
1003 /** ich habe diese variablen initialisiert, weil ich auf einem
1004 alpha rechner eine flaoting exception bei der ersten verwendung
1005 von ttry oder ttrys bekommen habe.
1006 das ist keine IEEE 764 architektur.
1007 es stand dann ttry irgend ein bitmuster, dass keiner regulaeren
1008 zahl entsprochen hat und daher die FPE versucht hat
1009 (in eine if abfrage !!)
1010 das zeigt, das da ttry tatsaechlich nicht initialiesiert verwendet
1011 wurde. das sollte man sauberer programmieren.
1012 **/
1013
1014 ttry = 0.;
1015 ttrys = 0.;
1016 ttly = 0.;
1017 ttlys = 0.;
1018 dl = 0.;
1019 dr = 0.;
1020 volr = 0.;
1021
1022 /* start of the set-up procedure */
1023
1024 /* step 1.0 of algorithm UTDR */
1025 /* see above or in unur_utdr_set_pdfatmode() */
1026
1027 do {
1028
1029 /* 1.1 */
1030 cfac = (setupok) ? GEN->c_factor : 2.; /* gibt es hier nur zwei varianten ?? ja*/
1031 c = cfac * DISTR.area/GEN->fm;
1032 setupok=1;
1033
1034 GEN->ttlx = DISTR.mode - c;
1035 GEN->ttrx = DISTR.mode + c;
1036
1037 /* 1.2 */
1038 /** TODO: kann man das nicht loeschen ?? **/
1039 if (/*GEN->il > -INFINITY &&*/ GEN->ttlx < GEN->il) {
1040 /* this is the case of no left tail*/
1041 GEN->bl = GEN->il;
1042 GEN->al = 0.;
1043 GEN->voll = 0.;
1044 if (GEN->il < DISTR.mode) {
1045 /* if the left domain border is left of the mode we set ttlx
1046 only to use it for the squeeze*/
1047 GEN->ttlx = DISTR.mode + (GEN->il - DISTR.mode) * 0.6;
1048 pdfx=PDF(GEN->ttlx);
1049 if (pdfx > SMALL_VAL)
1050 GEN->sal = (GEN->hm + 1./sqrt(pdfx)) / (DISTR.mode - GEN->ttlx);
1051 else
1052 GEN->ttlx = DISTR.mode;
1053 /* pdfx is too small. We set ttlx=mode that no squeeze is used */
1054 }
1055 }
1056 else {
1057 ttlys = PDF(GEN->ttlx);
1058 if (ttlys < SMALL_VAL) {
1059 /* in this case we cut off the left tail*/
1060 GEN->il = GEN->ttlx;
1061 GEN->bl = GEN->il;
1062 GEN->al = 0.;
1063 GEN->voll = 0.;
1064 GEN->ttlx=DISTR.mode;
1065 /* pdfx is too small. We set ttlx=mode that no squeeze is used */
1066 }
1067 else {
1068 ttlys = -1./sqrt(ttlys);
1069 GEN->sal = (GEN->hm - ttlys) / (DISTR.mode - GEN->ttlx);
1070
1071 /* delta1> 0 as -ttlys>0 und sal >0; */
1072 delta2 = ( GEN->sal > 0. ) ? -ttlys/GEN->sal : -ttlys;
1073 delta1 = fabs(GEN->ttlx);
1074 delta = GEN->delta_factor * ((delta1<=delta2) ? delta2 : delta1);
1075 if (delta > c * 0.01) {
1076 delta = UNUR_SQRT_DBL_EPSILON * ((delta1<=delta2) ? delta2 : delta1);
1077 /* Using sqrt(DBL_EPSILON) as delta_factor guarantees that not more than
1078 half of the precision will be lost when computing the slope GEN->al */
1079 if (delta > c * 0.01) {
1080 delta = c * 0.01;
1081 /* delta is forced to be c * 0.01 although this can
1082 result in numerical inaccuracies when computing
1083 the slope GEN->al. Therefore a warning is
1084 issued */
1085 _unur_warning(gen->genid,UNUR_ERR_GEN_DATA,
1086 "Delta larger than c/100!!, perhaps you can use a mode closer to 0 to remove this problem?");
1087 }
1088 }
1089
1090 ttly = -1./sqrt(PDF(GEN->ttlx+delta));
1091 GEN->al = (ttly-ttlys)/delta;
1092
1093 if (GEN->al <= 0.)
1094 /* setupok==0 means that this setup is quitted and set-up is restarted
1095 with different value for cfac */
1096 setupok = 0;
1097 else {
1098 GEN->bl = GEN->ttlx + (GEN->hm - ttly)/GEN->al;
1099 dl = ttly - GEN->al * GEN->ttlx;
1100 GEN->voll = -1./(GEN->al * GEN->hm);
1101 GEN->col = GEN->voll;
1102 if (GEN->il > -INFINITY)
1103 GEN->voll += 1./(GEN->al * (GEN->al * GEN->il + dl));
1104 }
1105 }
1106 }
1107
1108 /* 1.3 */
1109 if(setupok) {
1110 if (/*GEN->ir < INFINITY &&*/ GEN->ttrx > GEN->ir) {
1111 /* this is the case of no right tail */
1112 GEN->br = GEN->ir;
1113 GEN->ar = 0.;
1114 volr = 0.;
1115 if (GEN->ir > DISTR.mode) {
1116 /* if the right domain border is right of the mode we set ttrx
1117 only to use it for the squeeze */
1118 GEN->ttrx = DISTR.mode + (GEN->ir - DISTR.mode) * 0.6;
1119 pdfx = PDF(GEN->ttrx);
1120 if (pdfx > SMALL_VAL)
1121 GEN->sar = (GEN->hm + 1./sqrt(PDF(GEN->ttrx))) / (DISTR.mode - GEN->ttrx);
1122 else
1123 GEN->ttrx = DISTR.mode;
1124 /* pdfx is too small. We set ttrx=mode that no squeeze is used */
1125 }
1126 }
1127 else {
1128 ttrys = PDF(GEN->ttrx);
1129 if (ttrys < SMALL_VAL){
1130 /* in this case we cut off the right tail */
1131 GEN->ir = GEN->ttrx;
1132 GEN->br = GEN->ir;
1133 GEN->ar = 0.;
1134 volr = 0.;
1135 GEN->ttrx = DISTR.mode;
1136 /* pdfx is too small. We set ttrx=mode that no squeeze is used */
1137 }
1138 else {
1139 ttrys= -1./sqrt(ttrys);
1140 /* see 1.2. for explanations */
1141 GEN->sar = (GEN->hm - ttrys) / (DISTR.mode - GEN->ttrx);
1142 /* delta is positive, da ttrys<0 und sar <0 */
1143 delta2 = (GEN->sar<0.) ? ttrys/GEN->sar : -ttrys;
1144 delta1 = fabs(GEN->ttrx);
1145 delta = GEN->delta_factor * ((delta1<=delta2) ? delta2 : delta1);
1146 if (delta > c*0.01) {
1147 delta = UNUR_SQRT_DBL_EPSILON * ((delta1<=delta2) ? delta2 : delta1);
1148 /* Using sqrt(DBL_EPSILON) as delta_factor guarantees that not more than
1149 have of the precision will be lost when computing the slope GEN->al */
1150 if (delta > c*0.01) {
1151 delta=c*0.01;
1152 /* delta is forced to be c*0.01 allthough this can result in numerical
1153 inaccuracies when computing the slope GEN->al. Therefore a warning is
1154 issued*/
1155 _unur_warning(gen->genid,UNUR_ERR_ROUNDOFF,
1156 "Delta larger than c/100!!, perhaps you can use a mode closer to 0 to remove this problem?");
1157 }
1158 }
1159
1160 ttry = -1./sqrt(PDF(GEN->ttrx-delta));
1161 GEN->ar = (ttrys - ttry)/delta;
1162 if (GEN->ar >= 0.)
1163 /* setupok==0 means that this setup is quitted and set-up is
1164 restarted with different value for cfac */
1165 setupok = 0;
1166 else {
1167 GEN->br = GEN->ttrx + (GEN->hm - ttry) / GEN->ar;
1168 dr = ttry - GEN->ar * GEN->ttrx;
1169 volr = 1./(GEN->ar * GEN->hm);
1170 GEN->cor = volr;
1171 if (GEN->ir<INFINITY)
1172 volr -= 1./(GEN->ar * (GEN->ar * GEN->ir + dr));
1173 }
1174 }
1175 }
1176 }
1177
1178 /* 1.4 */
1179 if(setupok) {
1180 volc = (GEN->br - GEN->bl) * GEN->fm;
1181 GEN->vollc = GEN->voll + volc;
1182 GEN->volcompl = GEN->vollc + volr;
1183 if (volc>0.)
1184 GEN->brblvolc = (GEN->br - GEN->bl)/volc;
1185 if (!_unur_iszero(GEN->ar)) {
1186 GEN->drar = dr/GEN->ar;
1187 GEN->ooar2 = 1./(GEN->ar*GEN->ar);
1188 }
1189 if (!_unur_iszero(GEN->al)) {
1190 GEN->dlal = dl/GEN->al;
1191 GEN->ooal2 = 1./(GEN->al*GEN->al);
1192 }
1193 }
1194
1195 #ifdef UNUR_ENABLE_LOGGING
1196 /* write info into LOG file */
1197 if (gen->debug) _unur_utdr_debug_init(gen,ttly,ttlys,ttry,ttrys,cfac,setupok,c);
1198 #endif
1199
1200 if (!_unur_isfsame(cfac,2.)) {
1201 if(setupok)
1202 if (GEN->volcompl > 4. * DISTR.area || GEN->volcompl < 0.5 * DISTR.area)
1203 setupok=0;
1204 }
1205 else {
1206 if (setupok==0 || GEN->volcompl > 8. * DISTR.area || GEN->volcompl < 0.5 * DISTR.area) {
1207 _unur_error(gen->genid,UNUR_ERR_GEN_DATA,"; Area below hat too large or zero!! possible reasons: PDF, mode or area below PDF wrong; density not T-concave\n");
1208 return 0;
1209 }
1210 }
1211
1212 } while (!setupok);
1213
1214 return UNUR_SUCCESS;
1215
1216 } /* end of _unur_utdr_hat() */
1217
1218 #undef SMALL_VAL
1219
1220 /*****************************************************************************/
1221 /** Debugging utilities **/
1222 /*****************************************************************************/
1223
1224 /*---------------------------------------------------------------------------*/
1225 #ifdef UNUR_ENABLE_LOGGING
1226 /*---------------------------------------------------------------------------*/
1227
1228 static void
_unur_utdr_debug_init(const struct unur_gen * gen,double ttly,double ttlys,double ttry,double ttrys,double cfac,int setupok,double c)1229 _unur_utdr_debug_init( const struct unur_gen *gen,
1230 double ttly, double ttlys, double ttry, double ttrys,
1231 double cfac, int setupok, double c )
1232 /*----------------------------------------------------------------------*/
1233 /* write info about generator into LOG file */
1234 /* */
1235 /* parameters: */
1236 /* gen ... pointer to generator object */
1237 /* */
1238 /*----------------------------------------------------------------------*/
1239 {
1240 FILE *LOG;
1241
1242 /* check arguments */
1243 CHECK_NULL(gen,RETURN_VOID); COOKIE_CHECK(gen,CK_UTDR_GEN,RETURN_VOID);
1244
1245 LOG = unur_get_stream();
1246
1247 fprintf(LOG,"%s:\n",gen->genid);
1248 fprintf(LOG,"%s: type = continuous univariate random variates\n",gen->genid);
1249 fprintf(LOG,"%s: method = utdr(transformed density rejection with 3 points of contact)\n",gen->genid);
1250 fprintf(LOG,"%s:\n",gen->genid);
1251
1252 _unur_distr_cont_debug( gen->distr, gen->genid );
1253
1254 fprintf(LOG,"%s: sampling routine = _unur_utdr_sample",gen->genid);
1255 if (gen->variant & UTDR_VARFLAG_VERIFY)
1256 fprintf(LOG,"_check()\n");
1257 else
1258 fprintf(LOG,"()\n");
1259 fprintf(LOG,"%s:\n",gen->genid);
1260
1261 fprintf(LOG,"%s: Data for hat and squeeze:\n",gen->genid);
1262 fprintf(LOG,"%s:\n",gen->genid);
1263 fprintf(LOG,"%s:\tc_factor=%e delta_factor=%e real c=%e\n",gen->genid,GEN->c_factor,GEN->delta_factor,c);
1264 fprintf(LOG,"%s:\ttlx=%e bl=%e mode=%e\n",gen->genid,GEN->ttlx,GEN->bl,DISTR.mode);
1265 fprintf(LOG,"%s:\tbr=%e trx=%e\n",gen->genid,GEN->br,GEN->ttrx);
1266 fprintf(LOG,"%s:\ttly=%e tlys=%e al=%e \n",gen->genid,ttly,ttlys,GEN->al);
1267 fprintf(LOG,"%s:\ttry=%e trys=%e ar=%e \n",gen->genid,ttry,ttrys,GEN->ar);
1268 fprintf(LOG,"%s:\tcfac=%e setupok=%d volcompl=%e pdf_area=%e\n",gen->genid,cfac,setupok,GEN->volcompl,DISTR.area);
1269 fprintf(LOG,"%s:\n",gen->genid);
1270 fprintf(LOG,"%s: INIT completed **********************\n",gen->genid);
1271
1272 } /* end of _unur_utdr_debug_init() */
1273
1274 /*---------------------------------------------------------------------------*/
1275 #endif /* end UNUR_ENABLE_LOGGING */
1276 /*---------------------------------------------------------------------------*/
1277
1278
1279 /*---------------------------------------------------------------------------*/
1280 #ifdef UNUR_ENABLE_INFO
1281 /*---------------------------------------------------------------------------*/
1282
1283 void
_unur_utdr_info(struct unur_gen * gen,int help)1284 _unur_utdr_info( struct unur_gen *gen, int help )
1285 /*----------------------------------------------------------------------*/
1286 /* create character string that contains information about the */
1287 /* given generator object. */
1288 /* */
1289 /* parameters: */
1290 /* gen ... pointer to generator object */
1291 /* help ... whether to print additional comments */
1292 /*----------------------------------------------------------------------*/
1293 {
1294 struct unur_string *info = gen->infostr;
1295 struct unur_distr *distr = gen->distr;
1296 int samplesize = 10000;
1297
1298 /* generator ID */
1299 _unur_string_append(info,"generator ID: %s\n\n", gen->genid);
1300
1301 /* distribution */
1302 _unur_string_append(info,"distribution:\n");
1303 _unur_distr_info_typename(gen);
1304 _unur_string_append(info," functions = PDF\n");
1305 _unur_string_append(info," domain = (%g, %g)\n", DISTR.domain[0],DISTR.domain[1]);
1306 _unur_string_append(info," mode = %g %s\n", unur_distr_cont_get_mode(distr),
1307 (distr->set & UNUR_DISTR_SET_MODE_APPROX) ? "[numeric.]" : "");
1308 _unur_string_append(info," area(PDF) = %g\n", DISTR.area);
1309 _unur_string_append(info,"\n");
1310
1311 /* method */
1312 _unur_string_append(info,"method: UTDR (Universal Transformed Density Rejection -- 3 point method)\n");
1313 _unur_string_append(info,"\n");
1314
1315 /* performance */
1316 _unur_string_append(info,"performance characteristics:\n");
1317 _unur_string_append(info," rejection constant = %.2f [approx.]\n",
1318 unur_test_count_urn(gen,samplesize,0,NULL)/(2.*samplesize));
1319 _unur_string_append(info,"\n");
1320
1321 /* parameters */
1322 if (help) {
1323 _unur_string_append(info,"parameters:\n");
1324 _unur_string_append(info," deltafactor = %g %s\n", GEN->delta_factor,
1325 (gen->set & UTDR_SET_DELTA) ? "" : "[default]");
1326
1327 if (gen->set & UTDR_SET_PDFMODE)
1328 _unur_string_append(info," pdfatmode = %g\n", GEN->fm);
1329
1330 if (gen->set & UTDR_SET_CPFACTOR)
1331 _unur_string_append(info," cpfactor = %g\n", GEN->c_factor);
1332
1333 if (gen->variant & UTDR_VARFLAG_VERIFY)
1334 _unur_string_append(info," verify = on\n");
1335
1336 _unur_string_append(info,"\n");
1337 }
1338
1339 /* Hints */
1340 /* if (help) { */
1341 /* _unur_string_append(info,"\n"); */
1342 /* } */
1343
1344 } /* end of _unur_utdr_info() */
1345
1346 /*---------------------------------------------------------------------------*/
1347 #endif /* end UNUR_ENABLE_INFO */
1348 /*---------------------------------------------------------------------------*/
1349