1 /*****************************************************************************
2 * *
3 * unuran -- Universal Non-Uniform Random number generator *
4 * *
5 *****************************************************************************
6 * *
7 * FILE: arou.h *
8 * *
9 * TYPE: continuous univariate random variate *
10 * METHOD: ratio-of-uniforms with enveloping polygon *
11 * *
12 * DESCRIPTION: *
13 * Given PDF of a T-concave distribution; *
14 * produce a value x consistent with its density *
15 * *
16 *****************************************************************************
17 * *
18 * Copyright (c) 2000-2006 Wolfgang Hoermann and Josef Leydold *
19 * Department of Statistics and Mathematics, WU Wien, Austria *
20 * *
21 * This program is free software; you can redistribute it and/or modify *
22 * it under the terms of the GNU General Public License as published by *
23 * the Free Software Foundation; either version 2 of the License, or *
24 * (at your option) any later version. *
25 * *
26 * This program is distributed in the hope that it will be useful, *
27 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
28 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
29 * GNU General Public License for more details. *
30 * *
31 * You should have received a copy of the GNU General Public License *
32 * along with this program; if not, write to the *
33 * Free Software Foundation, Inc., *
34 * 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA *
35 * *
36 *****************************************************************************
37 * *
38 * REFERENCES: *
39 * [1] Leydold J. (2000): Automatic Sampling with the ratio-of-uniforms *
40 * method, ACM TOMS, forthcoming *
41 * *
42 * [2] Kinderman, A.J. and Monahan, F.J. (1977): Computer generation of *
43 * random variables using the ratio of uniform deviates, *
44 * ACM Trans. Math. Software 3(3), pp. 257--260. *
45 * *
46 *****************************************************************************
47 * *
48 * The ratio-of-uniforms method introduced in [2] is a flexible method that *
49 * is based on the following theorem: *
50 * *
51 * THEOREM: *
52 * Let X be a random variable with density function f(x) = g(x) / G, *
53 * where g(x) is a positive integrable function with support (x_0,x_1) *
54 * not necessarily finite and G = integral g(x) dx. *
55 * If (V,U) is uniformly distributed in *
56 * A = {(v,u): 0 < u <= sqrt(g(v/u)), x_0 < v/u < x_1}, *
57 * then X = V/U has probability density function f(x). *
58 * *
59 * Generating point (V,U) uniformly distributed in A is done by rejection *
60 * from an enveloping region, usually from the minimal bounding rectangle. *
61 * *
62 * The implemented algorithm uses the fact, that for many distribtions, *
63 * A is convex. Then we easily can construct an enveloping polygon by means *
64 * of tangent lines and a squeeze region by means of secants. *
65 * The resulting algorithm is very fast, since we can sample from the *
66 * squeeze region with immedate acceptance (thus only one uniform random *
67 * number is necessary). The region between envelope and squeeze consists *
68 * of triangles and can be made arbitrarily small (see [1] for details). *
69 * *
70 * Distributions with a convex set A are characterized by the following *
71 * theorem that shows a connection to transformed density rejection TDR. *
72 * *
73 * THEOREM: *
74 * A is convex if and only if g is T-concave with transformation *
75 * T(x) = -1/sqrt(x), i.e., -1/sqrt(g(x)) is a concave function. *
76 * *
77 *...........................................................................*
78 * *
79 * The region A is divided into segments, that contain the origin. Each *
80 * segments consist of an inner triangle (the squeeze region) and the *
81 * outer triangle (the region between envelope and squeeze). We have to *
82 * compute the areas of these triangles. *
83 * *
84 * To generate from the distribution we have to sample from a discrete *
85 * random variate with probability vector proportional to these areas to get *
86 * one of these triangles (inner or outer). We use indexed search (or guide *
87 * tables) to perform this task ([1], see also description of DIS). *
88 * When we have an inner triangle (squeeze), we reuse the uniform random *
89 * variate to get a point (v,u) uniformly distributed on the edge opposite *
90 * to the origin and return the ratio x = v/u (Thus generation from the *
91 * squeeze region is equivalent to the inversion method.) *
92 * When whe have an outer triangle, we have to sample a point (v,u) *
93 * uniformly distributed in this triangle. If u <= g(v/u) then return the *
94 * ratio x = v/u, otherwise reject. *
95 * *
96 *...........................................................................*
97 * *
98 * Algorithm AROU *
99 * *
100 * [Required] *
101 * PDF f(x), construction points c_1,...,c_n *
102 * *
103 * [Setup] *
104 * 1: Construct inner triangles S_i^s and outer triangles S_i^o *
105 * 2: Foreach triangle Do *
106 * 3: Compute respective areas A_i^s and A_i^o. *
107 * *
108 * [Generate] *
109 * 4: Generate I proportional to (A_1^s, A_1^o; A_2^s, A_2^o; ...). *
110 * 5: If inner triangle S_i^s Then *
111 * 6: Compute point (v,u) on edge (Reuse u.r.n. from 4). *
112 * 7: Return v/u. *
113 * 8: Else *
114 * 9: Generate point (V,U) uniformly distributed in outer triangle S_i^o.*
115 * 10: If u <= f(v/u) Return V/U. *
116 * 11: Goto 4. *
117 * *
118 *****************************************************************************/
119
120 /*---------------------------------------------------------------------------*/
121
122 #include <unur_source.h>
123 #include <distr/distr.h>
124 #include <distr/distr_source.h>
125 #include <distr/cont.h>
126 #include <urng/urng.h>
127 #include "unur_methods_source.h"
128 #include "x_gen_source.h"
129 #include "arou.h"
130 #include "arou_struct.h"
131
132 /*---------------------------------------------------------------------------*/
133 /* Variants */
134
135 #define AROU_VARFLAG_VERIFY 0x01u /* flag for verifying mode */
136 #define AROU_VARFLAG_USECENTER 0x02u /* flag whether center is used as cpoint or not */
137 #define AROU_VARFLAG_PEDANTIC 0x04u /* whether pedantic checking is used */
138 #define AROU_VARFLAG_USEDARS 0x10u /* whether DARS is used in setup or not */
139
140 /*---------------------------------------------------------------------------*/
141 /* Debugging flags */
142 /* bit 01 ... pameters and structure of generator (do not use here) */
143 /* bits 02-12 ... setup */
144 /* bits 13-24 ... adaptive steps */
145 /* bits 25-32 ... trace sampling */
146
147 #define AROU_DEBUG_SEGMENTS 0x00000010u /* print list of segments */
148 #define AROU_DEBUG_SPLIT 0x00010000u /* trace splitting of segments */
149 #define AROU_DEBUG_DARS 0x00020000u
150
151 /*---------------------------------------------------------------------------*/
152 /* Flags for logging set calls */
153
154 #define AROU_SET_CENTER 0x001u
155 #define AROU_SET_STP 0x002u
156 #define AROU_SET_N_STP 0x004u
157 #define AROU_SET_GUIDEFACTOR 0x010u
158 #define AROU_SET_MAX_SQHRATIO 0x020u
159 #define AROU_SET_MAX_SEGS 0x040u
160 #define AROU_SET_USE_DARS 0x100u
161 #define AROU_SET_DARS_FACTOR 0x200u
162
163 /*---------------------------------------------------------------------------*/
164
165 #define GENTYPE "AROU" /* type of generator */
166
167 /*---------------------------------------------------------------------------*/
168
169 static struct unur_gen *_unur_arou_init( struct unur_par *par );
170 /*---------------------------------------------------------------------------*/
171 /* Initialize new generator. */
172 /*---------------------------------------------------------------------------*/
173
174 static struct unur_gen *_unur_arou_create( struct unur_par *par );
175 /*---------------------------------------------------------------------------*/
176 /* create new (almost empty) generator object. */
177 /*---------------------------------------------------------------------------*/
178
179 static struct unur_gen *_unur_arou_clone( const struct unur_gen *gen );
180 /*---------------------------------------------------------------------------*/
181 /* copy (clone) generator object. */
182 /*---------------------------------------------------------------------------*/
183
184 static void _unur_arou_free( struct unur_gen *gen);
185 /*---------------------------------------------------------------------------*/
186 /* destroy generator object. */
187 /*---------------------------------------------------------------------------*/
188
189 static double _unur_arou_sample( struct unur_gen *gen );
190 static double _unur_arou_sample_check( struct unur_gen *gen );
191 /*---------------------------------------------------------------------------*/
192 /* sample from generator */
193 /*---------------------------------------------------------------------------*/
194
195 static int _unur_arou_get_starting_cpoints( struct unur_par *par, struct unur_gen *gen );
196 /*---------------------------------------------------------------------------*/
197 /* create list of construction points for starting segments. */
198 /* if user has not provided such points compute these by means of the */
199 /* "equi-angle rule". */
200 /*---------------------------------------------------------------------------*/
201
202 static int _unur_arou_get_starting_segments( struct unur_gen *gen );
203 /*---------------------------------------------------------------------------*/
204 /* compute segments from given starting construction points. */
205 /*---------------------------------------------------------------------------*/
206
207 static double _unur_arou_compute_x( double v, double u );
208 /*---------------------------------------------------------------------------*/
209 /* compute point x from (v,u) tuple. */
210 /*---------------------------------------------------------------------------*/
211
212 static int _unur_arou_run_dars( struct unur_gen *gen );
213 /*---------------------------------------------------------------------------*/
214 /* run derandomized adaptive rejection sampling. */
215 /*---------------------------------------------------------------------------*/
216
217 static struct unur_arou_segment *_unur_arou_segment_new( struct unur_gen *gen, double x, double fx );
218 /*---------------------------------------------------------------------------*/
219 /* make a new segment with left construction point x = v/u. */
220 /*---------------------------------------------------------------------------*/
221
222 static int _unur_arou_segment_parameter( struct unur_gen *gen, struct unur_arou_segment *seg );
223 /*---------------------------------------------------------------------------*/
224 /* compute all necessary data for segment. */
225 /* return: */
226 /* UNUR_SUCCESS ... on success */
227 /* UNUR_ERR_SILENT ... do not add this construction point */
228 /* UNUR_ERR_INF ... if area = INFINITY */
229 /* other error code ... error (PDF not T-concave) */
230 /*---------------------------------------------------------------------------*/
231
232 static int _unur_arou_segment_split( struct unur_gen *gen, struct unur_arou_segment *seg_old, double x, double fx );
233 /*---------------------------------------------------------------------------*/
234 /* split a segment at point (direction) x. */
235 /*---------------------------------------------------------------------------*/
236
237 static int _unur_arou_make_guide_table( struct unur_gen *gen );
238 /*---------------------------------------------------------------------------*/
239 /* make a guide table for indexed search. */
240 /*---------------------------------------------------------------------------*/
241
242 static double _unur_arou_segment_arcmean( struct unur_arou_segment *seg );
243 /*---------------------------------------------------------------------------*/
244 /* compute the "arcmean" of the two construction points of a segement. */
245 /*---------------------------------------------------------------------------*/
246
247 #ifdef UNUR_ENABLE_LOGGING
248 /*---------------------------------------------------------------------------*/
249 /* the following functions print debugging information on output stream, */
250 /* i.e., into the LOG file if not specified otherwise. */
251 /*---------------------------------------------------------------------------*/
252
253 static void _unur_arou_debug_init( const struct unur_par *par, const struct unur_gen *gen );
254 /*---------------------------------------------------------------------------*/
255 /* print after generator has been initialized has completed. */
256 /*---------------------------------------------------------------------------*/
257
258 static void _unur_arou_debug_dars_start( const struct unur_gen *gen );
259 /*---------------------------------------------------------------------------*/
260 /* print header before runniung derandomized adaptive rejection sampling. */
261 /*---------------------------------------------------------------------------*/
262
263 static void _unur_arou_debug_dars( const struct unur_gen *gen );
264 /*---------------------------------------------------------------------------*/
265 /* print after generator has run derandomized adaptive rejection sampling. */
266 /*---------------------------------------------------------------------------*/
267
268 static void _unur_arou_debug_free( const struct unur_gen *gen );
269 /*---------------------------------------------------------------------------*/
270 /* print before generater is destroyed. */
271 /*---------------------------------------------------------------------------*/
272
273 static void _unur_arou_debug_segments( const struct unur_gen *gen );
274 /*---------------------------------------------------------------------------*/
275 /* print data for segments. */
276 /*---------------------------------------------------------------------------*/
277
278 static void _unur_arou_debug_split_start( const struct unur_gen *gen,
279 const struct unur_arou_segment *seg,
280 double x, double fx );
281 static void _unur_arou_debug_split_stop( const struct unur_gen *gen,
282 const struct unur_arou_segment *seg_left,
283 const struct unur_arou_segment *seg_right );
284 /*---------------------------------------------------------------------------*/
285 /* print before and after a segment has been split (not / successfully). */
286 /*---------------------------------------------------------------------------*/
287 #endif
288
289 #ifdef UNUR_ENABLE_INFO
290 static void _unur_arou_info( struct unur_gen *gen, int help );
291 /*---------------------------------------------------------------------------*/
292 /* create info string. */
293 /*---------------------------------------------------------------------------*/
294 #endif
295
296 /*---------------------------------------------------------------------------*/
297 /* abbreviations */
298
299 #define DISTR_IN distr->data.cont /* data for distribution object */
300
301 #define PAR ((struct unur_arou_par*)par->datap) /* data for parameter object */
302 #define GEN ((struct unur_arou_gen*)gen->datap) /* data for generator object */
303 #define DISTR gen->distr->data.cont /* data for distribution in generator object */
304
305 #define BD_LEFT domain[0] /* left boundary of domain of distribution */
306 #define BD_RIGHT domain[1] /* right boundary of domain of distribution */
307
308 #define SAMPLE gen->sample.cont /* pointer to sampling routine */
309
310 #define PDF(x) _unur_cont_PDF((x),(gen->distr)) /* call to PDF */
311 #define dPDF(x) _unur_cont_dPDF((x),(gen->distr)) /* call to derivative of PDF */
312
313 /*---------------------------------------------------------------------------*/
314
315 #define _unur_arou_getSAMPLE(gen) \
316 ( ((gen)->variant & AROU_VARFLAG_VERIFY) \
317 ? _unur_arou_sample_check : _unur_arou_sample )
318
319 /*---------------------------------------------------------------------------*/
320
321 /*****************************************************************************/
322 /** Public: User Interface (API) **/
323 /*****************************************************************************/
324
325 struct unur_par *
unur_arou_new(const struct unur_distr * distr)326 unur_arou_new( const struct unur_distr *distr )
327 /*----------------------------------------------------------------------*/
328 /* get default parameters */
329 /* */
330 /* parameters: */
331 /* distr ... pointer to distribution object */
332 /* */
333 /* return: */
334 /* default parameters (pointer to structure) */
335 /* */
336 /* error: */
337 /* return NULL */
338 /*----------------------------------------------------------------------*/
339 {
340 struct unur_par *par;
341
342 /* check arguments */
343 _unur_check_NULL(GENTYPE,distr,NULL);
344
345 /* check distribution */
346 if (distr->type != UNUR_DISTR_CONT) {
347 _unur_error(GENTYPE,UNUR_ERR_DISTR_INVALID,""); return NULL; }
348 COOKIE_CHECK(distr,CK_DISTR_CONT,NULL);
349
350 if (DISTR_IN.pdf == NULL) {
351 _unur_error(GENTYPE,UNUR_ERR_DISTR_REQUIRED,"PDF"); return NULL;
352 }
353 if (DISTR_IN.dpdf == NULL) {
354 _unur_error(GENTYPE,UNUR_ERR_DISTR_REQUIRED,"derivative of PDF"); return NULL; }
355
356 /* allocate structure */
357 par = _unur_par_new( sizeof(struct unur_arou_par) );
358 COOKIE_SET(par,CK_AROU_PAR);
359
360 /* copy input */
361 par->distr = distr; /* pointer to distribution object */
362
363 /* set default values */
364 PAR->guide_factor = 2.; /* size of guide table / number of intervals */
365 PAR->darsfactor = 0.99; /* factor for (derandomized) ARS.
366 do not add a new construction point in a segment,
367 where abiguous region is too small, i.e. if
368 area / (|S^e\S^s|/number of segments) < darsfactor */
369
370 PAR->starting_cpoints = NULL; /* pointer to array of starting points */
371 PAR->n_starting_cpoints = 30; /* number of starting points */
372 PAR->max_segs = 100; /* maximum number of segments */
373 PAR->max_ratio = 0.99; /* do not add construction points if
374 ratio r_n = |P^s| / |P^e| > max_ratio */
375
376 par->method = UNUR_METH_AROU; /* method */
377 par->variant = ( AROU_VARFLAG_USECENTER |
378 AROU_VARFLAG_USEDARS ); /* default variant */
379 par->set = 0u; /* inidicate default parameters */
380 par->urng = unur_get_default_urng(); /* use default urng */
381 par->urng_aux = par->urng; /* no special auxilliary URNG */
382
383 par->debug = _unur_default_debugflag; /* set default debugging flags */
384
385 /* routine for starting generator */
386 par->init = _unur_arou_init;
387
388 return par;
389
390 } /* end of unur_arou_new() */
391
392 /*****************************************************************************/
393
394 int
unur_arou_set_usedars(struct unur_par * par,int usedars)395 unur_arou_set_usedars( struct unur_par *par, int usedars )
396 /*----------------------------------------------------------------------*/
397 /* set flag for using DARS (derandomized adaptive rejection sampling). */
398 /* */
399 /* parameters: */
400 /* par ... pointer to parameter for building generator object */
401 /* usedars ... 0 = do not use, !0 = use DARS */
402 /* */
403 /* return: */
404 /* UNUR_SUCCESS ... on success */
405 /* error code ... on error */
406 /* */
407 /* comment: */
408 /* using not using DARS is the default */
409 /*----------------------------------------------------------------------*/
410 {
411 /* check arguments */
412 _unur_check_NULL( GENTYPE, par, UNUR_ERR_NULL );
413
414 /* check input */
415 _unur_check_par_object( par, AROU );
416
417 /* we use a bit in variant */
418 par->variant = (usedars) ? (par->variant | AROU_VARFLAG_USEDARS) : (par->variant & (~AROU_VARFLAG_USEDARS));
419
420 /* changelog */
421 par->set |= AROU_SET_USE_DARS;
422
423 /* o.k. */
424 return UNUR_SUCCESS;
425
426 } /* end of unur_arou_set_usedars() */
427
428 /*---------------------------------------------------------------------------*/
429
430 int
unur_arou_set_darsfactor(struct unur_par * par,double factor)431 unur_arou_set_darsfactor( struct unur_par *par, double factor )
432 /*----------------------------------------------------------------------*/
433 /* set factor for derandomized adaptive rejection sampling */
434 /* */
435 /* parameters: */
436 /* par ... pointer to parameter for building generator object */
437 /* factor ... parameter for DARS */
438 /* */
439 /* return: */
440 /* UNUR_SUCCESS ... on success */
441 /* error code ... on error */
442 /*----------------------------------------------------------------------*/
443 {
444 /* check arguments */
445 _unur_check_NULL( GENTYPE, par, UNUR_ERR_NULL );
446
447 /* check input */
448 _unur_check_par_object( par, AROU );
449
450 /* check new parameter for generator */
451 if (factor < 0.) {
452 _unur_warning(GENTYPE,UNUR_ERR_PAR_SET,"DARS factor < 0");
453 return UNUR_ERR_PAR_SET;
454 }
455
456 /* store date */
457 PAR->darsfactor = factor;
458
459 /* changelog */
460 par->set |= AROU_SET_DARS_FACTOR;
461
462 return UNUR_SUCCESS;
463
464 } /* end of unur_arou_set_darsfactor() */
465
466 /*---------------------------------------------------------------------------*/
467
468 int
unur_arou_set_cpoints(struct unur_par * par,int n_stp,const double * stp)469 unur_arou_set_cpoints( struct unur_par *par, int n_stp, const double *stp )
470 /*----------------------------------------------------------------------*/
471 /* set construction points for envelope */
472 /* and/or its number for initialization */
473 /* */
474 /* parameters: */
475 /* par ... pointer to parameter for building generator object */
476 /* n_stp ... number of starting points */
477 /* stp ... pointer to array of starting points */
478 /* (NULL for changing only the number of default points) */
479 /* */
480 /* return: */
481 /* UNUR_SUCCESS ... on success */
482 /* error code ... on error */
483 /*----------------------------------------------------------------------*/
484 {
485 int i;
486
487 /* check arguments */
488 _unur_check_NULL( GENTYPE, par, UNUR_ERR_NULL );
489
490 /* check input */
491 _unur_check_par_object( par, AROU );
492
493 /* check starting construction points */
494 /* we always use the boundary points as additional starting points,
495 so we do not count these here! */
496 if (n_stp < 0 ) {
497 _unur_warning(GENTYPE,UNUR_ERR_PAR_SET,"number of starting points < 0");
498 return UNUR_ERR_PAR_SET;
499 }
500
501 if (stp)
502 /* starting points must be strictly monontonically increasing */
503 for( i=1; i<n_stp; i++ )
504 if (stp[i] <= stp[i-1]) {
505 _unur_warning(GENTYPE,UNUR_ERR_PAR_SET,
506 "starting points not strictly monotonically increasing");
507 return UNUR_ERR_PAR_SET;
508 }
509
510 /* store date */
511 PAR->starting_cpoints = stp;
512 PAR->n_starting_cpoints = n_stp;
513
514 /* changelog */
515 par->set |= AROU_SET_N_STP | ((stp) ? AROU_SET_STP : 0);
516
517 return UNUR_SUCCESS;
518
519 } /* end of unur_arou_set_cpoints() */
520
521 /*---------------------------------------------------------------------------*/
522
523 int
unur_arou_set_guidefactor(struct unur_par * par,double factor)524 unur_arou_set_guidefactor( struct unur_par *par, double factor )
525 /*----------------------------------------------------------------------*/
526 /* set factor for relative size of guide table */
527 /* */
528 /* parameters: */
529 /* par ... pointer to parameter for building generator object */
530 /* factor ... relative size of table */
531 /* */
532 /* return: */
533 /* UNUR_SUCCESS ... on success */
534 /* error code ... on error */
535 /*----------------------------------------------------------------------*/
536 {
537 /* check arguments */
538 _unur_check_NULL( GENTYPE, par, UNUR_ERR_NULL );
539
540 /* check input */
541 _unur_check_par_object( par, AROU );
542
543 /* check new parameter for generator */
544 if (factor < 0) {
545 _unur_warning(GENTYPE,UNUR_ERR_PAR_SET,"guide table size < 0");
546 return UNUR_ERR_PAR_SET;
547 }
548
549 /* store date */
550 PAR->guide_factor = factor;
551
552 /* changelog */
553 par->set |= AROU_SET_GUIDEFACTOR;
554
555 return UNUR_SUCCESS;
556
557 } /* end of unur_arou_set_guidefactor() */
558
559 /*---------------------------------------------------------------------------*/
560
561 int
unur_arou_set_max_sqhratio(struct unur_par * par,double max_ratio)562 unur_arou_set_max_sqhratio( struct unur_par *par, double max_ratio )
563 /*----------------------------------------------------------------------*/
564 /* set bound for ratio A(squeeze) / A(hat) */
565 /* */
566 /* parameters: */
567 /* par ... pointer to parameter for building generator object */
568 /* max_ratio ... upper bound for ratio to add a new construction point*/
569 /* */
570 /* return: */
571 /* UNUR_SUCCESS ... on success */
572 /* error code ... on error */
573 /*----------------------------------------------------------------------*/
574 {
575 /* check arguments */
576 _unur_check_NULL( GENTYPE, par, UNUR_ERR_NULL );
577
578 /* check input */
579 _unur_check_par_object( par, AROU );
580
581 /* check new parameter for generator */
582 if (max_ratio < 0. || max_ratio > 1. ) {
583 _unur_warning(GENTYPE,UNUR_ERR_PAR_SET,"ratio A(squeeze)/A(hat) not in [0,1]");
584 return UNUR_ERR_PAR_SET;
585 }
586
587 /* store date */
588 PAR->max_ratio = max_ratio;
589
590 /* changelog */
591 par->set |= AROU_SET_MAX_SQHRATIO;
592
593 return UNUR_SUCCESS;
594
595 } /* end of unur_arou_set_max_sqhratio() */
596
597 /*---------------------------------------------------------------------------*/
598
599 double
unur_arou_get_sqhratio(const struct unur_gen * gen)600 unur_arou_get_sqhratio( const struct unur_gen *gen )
601 /*----------------------------------------------------------------------*/
602 /* get ratio A(squeeze) / A(hat) */
603 /* */
604 /* parameters: */
605 /* gen ... pointer to generator object */
606 /* */
607 /* return: */
608 /* ratio ... on success */
609 /* INFINITY ... on error */
610 /*----------------------------------------------------------------------*/
611 {
612 /* check input */
613 _unur_check_NULL( GENTYPE, gen, INFINITY );
614 _unur_check_gen_object( gen, AROU, INFINITY );
615
616 return (GEN->Asqueeze / GEN->Atotal);
617
618 } /* end of unur_arou_get_sqhratio() */
619
620 /*---------------------------------------------------------------------------*/
621
622 double
unur_arou_get_hatarea(const struct unur_gen * gen)623 unur_arou_get_hatarea( const struct unur_gen *gen )
624 /*----------------------------------------------------------------------*/
625 /* get area below hat */
626 /* */
627 /* parameters: */
628 /* gen ... pointer to generator object */
629 /* */
630 /* return: */
631 /* area ... on success */
632 /* INFINITY ... on error */
633 /*----------------------------------------------------------------------*/
634 {
635 /* check input */
636 _unur_check_NULL( GENTYPE, gen, INFINITY );
637 _unur_check_gen_object( gen, AROU, INFINITY );
638
639 return GEN->Atotal;
640
641 } /* end of unur_arou_get_hatarea() */
642
643 /*---------------------------------------------------------------------------*/
644
645 double
unur_arou_get_squeezearea(const struct unur_gen * gen)646 unur_arou_get_squeezearea( const struct unur_gen *gen )
647 /*----------------------------------------------------------------------*/
648 /* get area below squeeze */
649 /* */
650 /* parameters: */
651 /* gen ... pointer to generator object */
652 /* */
653 /* return: */
654 /* area ... on success */
655 /* INFINITY ... on error */
656 /*----------------------------------------------------------------------*/
657 {
658 /* check input */
659 _unur_check_NULL( GENTYPE, gen, INFINITY );
660 _unur_check_gen_object( gen, AROU, INFINITY );
661
662 return GEN->Asqueeze;
663
664 } /* end of unur_arou_get_squeezearea() */
665
666 /*---------------------------------------------------------------------------*/
667
668 int
unur_arou_set_max_segments(struct unur_par * par,int max_segs)669 unur_arou_set_max_segments( struct unur_par *par, int max_segs )
670 /*----------------------------------------------------------------------*/
671 /* set maximum number of segments */
672 /* */
673 /* parameters: */
674 /* par ... pointer to parameter for building generator object */
675 /* max_segs ... maximum number of segments */
676 /* */
677 /* return: */
678 /* UNUR_SUCCESS ... on success */
679 /* error code ... on error */
680 /*----------------------------------------------------------------------*/
681 {
682 /* check arguments */
683 _unur_check_NULL( GENTYPE, par, UNUR_ERR_NULL );
684
685 /* check input */
686 _unur_check_par_object( par, AROU );
687
688 /* check new parameter for generator */
689 if (max_segs < 1 ) {
690 _unur_warning(GENTYPE,UNUR_ERR_PAR_SET,"maximum number of segments < 1");
691 return UNUR_ERR_PAR_SET;
692 }
693
694 /* store date */
695 PAR->max_segs = max_segs;
696
697 /* changelog */
698 par->set |= AROU_SET_MAX_SEGS;
699
700 return UNUR_SUCCESS;
701
702 } /* end of unur_arou_set_max_segments() */
703
704 /*---------------------------------------------------------------------------*/
705
706 int
unur_arou_set_usecenter(struct unur_par * par,int usecenter)707 unur_arou_set_usecenter( struct unur_par *par, int usecenter )
708 /*----------------------------------------------------------------------*/
709 /* set flag for using center as construction point */
710 /* */
711 /* parameters: */
712 /* par ... pointer to parameter for building generator object */
713 /* usecenter ... 0 = do not use, !0 = use */
714 /* */
715 /* return: */
716 /* UNUR_SUCCESS ... on success */
717 /* error code ... on error */
718 /* */
719 /* comment: */
720 /* using center as construction point is the default */
721 /*----------------------------------------------------------------------*/
722 {
723 /* check arguments */
724 _unur_check_NULL( GENTYPE, par, UNUR_ERR_NULL );
725
726 /* check input */
727 _unur_check_par_object( par, AROU );
728
729 /* we use a bit in variant */
730 par->variant = (usecenter) ? (par->variant | AROU_VARFLAG_USECENTER) : (par->variant & (~AROU_VARFLAG_USECENTER));
731
732 /* o.k. */
733 return UNUR_SUCCESS;
734
735 } /* end of unur_arou_set_usecenter() */
736
737 /*---------------------------------------------------------------------------*/
738
739 int
unur_arou_set_verify(struct unur_par * par,int verify)740 unur_arou_set_verify( struct unur_par *par, int verify )
741 /*----------------------------------------------------------------------*/
742 /* turn verifying of algorithm while sampling on/off */
743 /* */
744 /* parameters: */
745 /* par ... pointer to parameter for building generator object */
746 /* verify ... 0 = no verifying, !0 = verifying */
747 /* */
748 /* return: */
749 /* UNUR_SUCCESS ... on success */
750 /* error code ... on error */
751 /* */
752 /* comment: */
753 /* no verifying is the default */
754 /*----------------------------------------------------------------------*/
755 {
756 /* check arguments */
757 _unur_check_NULL( GENTYPE, par, UNUR_ERR_NULL );
758
759 /* check input */
760 _unur_check_par_object( par, AROU );
761
762 /* we use a bit in variant */
763 par->variant = (verify) ? (par->variant | AROU_VARFLAG_VERIFY) : (par->variant & (~AROU_VARFLAG_VERIFY));
764
765 /* o.k. */
766 return UNUR_SUCCESS;
767
768 } /* end of unur_arou_set_verify() */
769
770 /*---------------------------------------------------------------------------*/
771
772 int
unur_arou_chg_verify(struct unur_gen * gen,int verify)773 unur_arou_chg_verify( struct unur_gen *gen, int verify )
774 /*----------------------------------------------------------------------*/
775 /* turn verifying of algorithm while sampling on/off */
776 /* */
777 /* parameters: */
778 /* gen ... pointer to generator object */
779 /* verify ... 0 = no verifying, !0 = verifying */
780 /* */
781 /* return: */
782 /* UNUR_SUCCESS ... on success */
783 /* error code ... on error */
784 /* */
785 /* comment: */
786 /* no verifying is the default */
787 /*----------------------------------------------------------------------*/
788 {
789 /* check input */
790 _unur_check_NULL( GENTYPE, gen, UNUR_ERR_NULL );
791 _unur_check_gen_object( gen, AROU, UNUR_ERR_GEN_INVALID );
792
793 /* we must not change this switch when sampling has been disabled by
794 using a pointer to the error producing routine */
795 if (SAMPLE == _unur_sample_cont_error)
796 return UNUR_FAILURE;
797
798 if (verify)
799 /* turn verify mode on */
800 gen->variant |= AROU_VARFLAG_VERIFY;
801 else
802 /* turn verify mode off */
803 gen->variant &= ~AROU_VARFLAG_VERIFY;
804
805 SAMPLE = _unur_arou_getSAMPLE(gen);
806
807 /* o.k. */
808 return UNUR_SUCCESS;
809
810 } /* end of unur_arou_chg_verify() */
811
812 /*---------------------------------------------------------------------------*/
813
814 int
unur_arou_set_pedantic(struct unur_par * par,int pedantic)815 unur_arou_set_pedantic( struct unur_par *par, int pedantic )
816 /*----------------------------------------------------------------------*/
817 /* turn pedantic mode on/off */
818 /* */
819 /* parameters: */
820 /* par ... pointer to parameter for building generator object */
821 /* pedantic ... 0 = no pedantic mode, !0 = use pedantic mode */
822 /* */
823 /* return: */
824 /* UNUR_SUCCESS ... on success */
825 /* error code ... on error */
826 /* */
827 /* comment: */
828 /* pedantic is the default */
829 /*----------------------------------------------------------------------*/
830 {
831 /* check arguments */
832 _unur_check_NULL( GENTYPE, par, UNUR_ERR_NULL );
833
834 /* check input */
835 _unur_check_par_object( par, AROU );
836
837 /* we use a bit in variant */
838 par->variant = (pedantic) ? (par->variant | AROU_VARFLAG_PEDANTIC) : (par->variant & (~AROU_VARFLAG_PEDANTIC));
839
840 /* o.k. */
841 return UNUR_SUCCESS;
842
843 } /* end of unur_arou_set_pedantic() */
844
845
846 /*****************************************************************************/
847 /** Private **/
848 /*****************************************************************************/
849
850 struct unur_gen *
_unur_arou_init(struct unur_par * par)851 _unur_arou_init( struct unur_par *par )
852 /*----------------------------------------------------------------------*/
853 /* initialize new generator */
854 /* */
855 /* parameters: */
856 /* par ... pointer to paramters for building generator object */
857 /* */
858 /* return: */
859 /* pointer to generator object */
860 /* */
861 /* error: */
862 /* return NULL */
863 /*----------------------------------------------------------------------*/
864 {
865 struct unur_gen *gen;
866 int i,k;
867
868 /* check arguments */
869 CHECK_NULL(par,NULL);
870
871 /* check input */
872 if ( par->method != UNUR_METH_AROU ) {
873 _unur_error(GENTYPE,UNUR_ERR_PAR_INVALID,"");
874 return NULL; }
875 COOKIE_CHECK(par,CK_AROU_PAR,NULL);
876
877 /* create a new empty generator object */
878 gen = _unur_arou_create(par);
879 if (!gen) { _unur_par_free(par); return NULL; }
880
881 /* get starting points */
882 if (_unur_arou_get_starting_cpoints(par,gen)!=UNUR_SUCCESS ) {
883 #ifdef UNUR_ENABLE_LOGGING
884 if (gen->debug) _unur_arou_debug_init(par,gen);
885 #endif
886 _unur_par_free(par); _unur_arou_free(gen);
887 return NULL;
888 }
889
890 /* compute segments for given starting points */
891 if ( _unur_arou_get_starting_segments(gen)!=UNUR_SUCCESS ) {
892 _unur_error(gen->genid,UNUR_ERR_GEN_CONDITION,"PDF not T-concave");
893 #ifdef UNUR_ENABLE_LOGGING
894 if (gen->debug) _unur_arou_debug_init(par,gen);
895 #endif
896 _unur_par_free(par); _unur_arou_free(gen);
897 return NULL;
898 }
899
900 /* we have to update the maximal number of segments,
901 if the user wants more starting points. */
902 if (GEN->n_segs > GEN->max_segs) {
903 _unur_warning(gen->genid,UNUR_ERR_GEN_DATA,"maximal number of segments too small. increase.");
904 GEN->max_segs = GEN->n_segs;
905 }
906
907 if (gen->variant & AROU_VARFLAG_USEDARS) {
908 /* run derandomized adaptive rejection sampling (DARS) */
909
910 #ifdef UNUR_ENABLE_LOGGING
911 if (gen->debug & AROU_DEBUG_DARS) {
912 /* make initial guide table (only necessary for writing debug info) */
913 _unur_arou_make_guide_table(gen);
914 /* write info into LOG file */
915 _unur_arou_debug_init(par,gen);
916 _unur_arou_debug_dars_start(gen);
917 }
918 #endif
919
920 for (i=0; i<3; i++) {
921 /* we make several tries */
922
923 /* run DARS */
924 if ( _unur_arou_run_dars(gen)!=UNUR_SUCCESS ) {
925 _unur_par_free(par); _unur_arou_free(gen);
926 return NULL;
927 }
928
929 /* make initial guide table */
930 _unur_arou_make_guide_table(gen);
931
932 /* check if DARS was completed */
933 if (GEN->n_segs < GEN->max_segs) {
934 /* ran ARS instead */
935 for (k=0; k<5; k++)
936 _unur_sample_cont(gen);
937 }
938 else
939 break;
940 }
941
942 #ifdef UNUR_ENABLE_LOGGING
943 /* write info into LOG file */
944 if (gen->debug) {
945 if (gen->debug & AROU_DEBUG_DARS)
946 _unur_arou_debug_dars(gen);
947 else
948 _unur_arou_debug_init(par,gen);
949 }
950 #endif
951 }
952
953 else { /* do not run DARS */
954 /* make initial guide table */
955 _unur_arou_make_guide_table(gen);
956
957 #ifdef UNUR_ENABLE_LOGGING
958 /* write info into LOG file */
959 if (gen->debug) _unur_arou_debug_init(par,gen);
960 #endif
961
962 }
963
964 /* free parameters */
965 _unur_par_free(par);
966
967 /* is there any envelope at all ? */
968 if (GEN->Atotal <= 0. || !_unur_isfinite(GEN->Atotal)) {
969 _unur_error(gen->genid,UNUR_ERR_GEN_DATA,"bad construction points");
970 _unur_arou_free(gen);
971 return NULL;
972 }
973
974 /* o.k. */
975 return gen;
976
977 } /* end of _unur_arou_init() */
978
979 /*---------------------------------------------------------------------------*/
980
981 struct unur_gen *
_unur_arou_create(struct unur_par * par)982 _unur_arou_create( struct unur_par *par )
983 /*----------------------------------------------------------------------*/
984 /* allocate memory for generator */
985 /* */
986 /* parameters: */
987 /* par ... pointer to parameter for building generator object */
988 /* */
989 /* return: */
990 /* pointer to (empty) generator object with default settings */
991 /* */
992 /* error: */
993 /* return NULL */
994 /*----------------------------------------------------------------------*/
995 {
996 struct unur_gen *gen;
997
998 /* check arguments */
999 CHECK_NULL(par,NULL); COOKIE_CHECK(par,CK_AROU_PAR,NULL);
1000
1001 /* create new generic generator object */
1002 gen = _unur_generic_create( par, sizeof(struct unur_arou_gen) );
1003
1004 /* magic cookies */
1005 COOKIE_SET(gen,CK_AROU_GEN);
1006
1007 /* set generator identifier */
1008 gen->genid = _unur_set_genid(GENTYPE);
1009
1010 /* routines for sampling and destroying generator */
1011 SAMPLE = _unur_arou_getSAMPLE(gen);
1012 gen->destroy = _unur_arou_free;
1013 gen->clone = _unur_arou_clone;
1014
1015 /* set all pointers to NULL */
1016 GEN->seg = NULL;
1017 GEN->n_segs = 0;
1018 GEN->guide = NULL;
1019 GEN->guide_size = 0;
1020 GEN->Atotal = 0.;
1021 GEN->Asqueeze = 0.;
1022
1023 /* copy some parameters into generator object */
1024 GEN->guide_factor = PAR->guide_factor; /* relative size of guide tables */
1025
1026 /* bounds for adding construction points */
1027 GEN->max_segs = PAR->max_segs; /* maximum number of segments */
1028 #ifdef UNUR_ENABLE_INFO
1029 GEN->max_segs_info = PAR->max_segs; /* ... for info string */
1030 #endif
1031 GEN->max_ratio = PAR->max_ratio;
1032 GEN->darsfactor = PAR->darsfactor;
1033
1034 /* get center */
1035 if ( (gen->distr->set & UNUR_DISTR_SET_CENTER) ||
1036 (gen->distr->set & UNUR_DISTR_SET_MODE) ) {
1037 GEN->center = unur_distr_cont_get_center(gen->distr);
1038 /* center must be in domain */
1039 GEN->center = _unur_max(GEN->center,DISTR.BD_LEFT);
1040 GEN->center = _unur_min(GEN->center,DISTR.BD_RIGHT);
1041 gen->set |= AROU_SET_CENTER;
1042 }
1043 else {
1044 GEN->center = 0.;
1045 /* we cannot use the center as construction point */
1046 gen->variant &= ~AROU_VARFLAG_USECENTER;
1047 }
1048
1049 #ifdef UNUR_ENABLE_INFO
1050 /* set function for creating info string */
1051 gen->info = _unur_arou_info;
1052 #endif
1053
1054 /* return pointer to (almost empty) generator object */
1055 return(gen);
1056
1057 } /* end of _unur_arou_create() */
1058
1059 /*---------------------------------------------------------------------------*/
1060
1061 struct unur_gen *
_unur_arou_clone(const struct unur_gen * gen)1062 _unur_arou_clone( const struct unur_gen *gen )
1063 /*----------------------------------------------------------------------*/
1064 /* copy (clone) generator object */
1065 /* */
1066 /* parameters: */
1067 /* gen ... pointer to generator object */
1068 /* */
1069 /* return: */
1070 /* pointer to clone of generator object */
1071 /* */
1072 /* error: */
1073 /* return NULL */
1074 /*----------------------------------------------------------------------*/
1075 {
1076 #define CLONE ((struct unur_arou_gen*)clone->datap)
1077
1078 struct unur_gen *clone;
1079 struct unur_arou_segment *seg,*next, *clone_seg, *clone_prev;
1080
1081 /* check arguments */
1082 CHECK_NULL(gen,NULL); COOKIE_CHECK(gen,CK_AROU_GEN,NULL);
1083
1084 /* create generic clone */
1085 clone = _unur_generic_clone( gen, GENTYPE );
1086
1087 /* copy linked list of segments */
1088 clone_seg = NULL;
1089 clone_prev = NULL;
1090 for (seg = GEN->seg; seg != NULL; seg = next) {
1091 /* copy segment */
1092 clone_seg = _unur_xmalloc( sizeof(struct unur_arou_segment) );
1093 memcpy( clone_seg, seg, sizeof(struct unur_arou_segment) );
1094 if (clone_prev == NULL) {
1095 /* starting point of linked list */
1096 CLONE->seg = clone_seg;
1097 }
1098 else {
1099 /* insert into linked list */
1100 clone_prev->next = clone_seg;
1101 clone_prev->rtp = clone_seg->ltp;
1102 clone_prev->drtp = clone_seg->dltp;
1103 }
1104 /* next step */
1105 next = seg->next;
1106 clone_prev = clone_seg;
1107 }
1108 /* terminate linked list */
1109 if (clone_seg) clone_seg->next = NULL;
1110
1111 /* make new guide table */
1112 CLONE->guide = NULL;
1113 _unur_arou_make_guide_table(clone);
1114
1115 /* finished clone */
1116 return clone;
1117
1118 #undef CLONE
1119 } /* end of _unur_arou_clone() */
1120
1121 /*---------------------------------------------------------------------------*/
1122
1123 void
_unur_arou_free(struct unur_gen * gen)1124 _unur_arou_free( struct unur_gen *gen )
1125 /*----------------------------------------------------------------------*/
1126 /* deallocate generator object */
1127 /* */
1128 /* parameters: */
1129 /* gen ... pointer to generator object */
1130 /*----------------------------------------------------------------------*/
1131 {
1132 /* check arguments */
1133 if( !gen ) /* nothing to do */
1134 return;
1135
1136 /* check input */
1137 if ( gen->method != UNUR_METH_AROU ) {
1138 _unur_warning(gen->genid,UNUR_ERR_GEN_INVALID,"");
1139 return; }
1140 COOKIE_CHECK(gen,CK_AROU_GEN,RETURN_VOID);
1141
1142 /* we cannot use this generator object any more */
1143 SAMPLE = NULL; /* make sure to show up a programming error */
1144
1145 /* write info into LOG file */
1146 #ifdef UNUR_ENABLE_LOGGING
1147 if (gen->debug) _unur_arou_debug_free(gen);
1148 #endif
1149
1150 /* free linked list of segments */
1151 {
1152 struct unur_arou_segment *seg,*next;
1153 for (seg = GEN->seg; seg != NULL; seg = next) {
1154 next = seg->next;
1155 free(seg);
1156 }
1157 }
1158
1159 /* free other memory not stored in list */
1160 if (GEN->guide) free(GEN->guide);
1161
1162 _unur_generic_free(gen);
1163
1164 } /* end of _unur_arou_free() */
1165
1166 /*****************************************************************************/
1167
1168 double
_unur_arou_sample(struct unur_gen * gen)1169 _unur_arou_sample( struct unur_gen *gen )
1170 /*----------------------------------------------------------------------*/
1171 /* sample from generator */
1172 /* */
1173 /* parameters: */
1174 /* gen ... pointer to generator object */
1175 /* */
1176 /* return: */
1177 /* double (sample from random variate) */
1178 /* */
1179 /* error: */
1180 /* return INFINITY */
1181 /*----------------------------------------------------------------------*/
1182 {
1183 /** TODO: check uniform random number: u != 0 and u != 1 ?? **/
1184
1185 UNUR_URNG *urng; /* pointer to uniform random number generator */
1186 struct unur_arou_segment *seg;
1187 int result_split;
1188 double R,R1,R2,R3,tmp,x,fx,u;
1189
1190 /* check arguments */
1191 CHECK_NULL(gen,INFINITY); COOKIE_CHECK(gen,CK_AROU_GEN,INFINITY);
1192
1193 /* main URNG */
1194 urng = gen->urng;
1195
1196 while (1) {
1197
1198 /* sample from U(0,1) */
1199 R = _unur_call_urng(urng);
1200
1201 /* look up in guide table and search for segment */
1202 seg = GEN->guide[(int) (R * GEN->guide_size)];
1203 R *= GEN->Atotal;
1204 while (seg->Acum < R) {
1205 seg = seg->next;
1206 }
1207 COOKIE_CHECK(seg,CK_AROU_SEG,INFINITY);
1208
1209 /* reuse of uniform random number */
1210 R = seg->Acum - R;
1211
1212 /* inside or outside squeeze */
1213 if (R < seg->Ain) {
1214 /* inside */
1215 /* reuse of random number. */
1216 /* We can avoid R = (seg->Ain - R) / seg->Ain */
1217 return( ( seg->Ain * seg->rtp[0] + R * (seg->ltp[0] - seg->rtp[0]) ) /
1218 ( seg->Ain * seg->rtp[1] + R * (seg->ltp[1] - seg->rtp[1]) ) );
1219 }
1220
1221 else {
1222 /* outside */
1223
1224 /* from now on we use the auxilliary generator
1225 (it can be the same as the main generator) */
1226 urng = gen->urng_aux;
1227
1228 /* three uniform random numbers with R1 + R2 + R3 = 1 */
1229 R1 = (R - seg->Ain) / seg->Aout; /* reuse of random number (good ?? ) */
1230 R2 = _unur_call_urng(urng);
1231 if (R1>R2) { tmp = R1; R1=R2; R2=tmp; } /* swap */
1232 R3 = 1.-R2;
1233 R2 -= R1;
1234
1235 /* point (v,u) and ratio x = v/u */
1236 u = seg->ltp[1]*R1 + seg->rtp[1]*R2 + seg->mid[1]*R3;
1237 x = (seg->ltp[0]*R1 + seg->rtp[0]*R2 + seg->mid[0]*R3) / u;
1238
1239 /* density at x */
1240 fx = PDF(x);
1241
1242 /* being outside the squeeze is bad. improve the situation! */
1243 if (GEN->n_segs < GEN->max_segs) {
1244 if (GEN->max_ratio * GEN->Atotal > GEN->Asqueeze) {
1245 result_split = _unur_arou_segment_split(gen,seg,x,fx);
1246 if ( !(result_split == UNUR_SUCCESS || result_split == UNUR_ERR_SILENT) ) {
1247 /* condition for PDF is violated! */
1248 _unur_error(gen->genid,UNUR_ERR_GEN_CONDITION,"");
1249 if (gen->variant & AROU_VARFLAG_PEDANTIC) {
1250 /* replace sampling routine by dummy routine that just returns INFINITY */
1251 SAMPLE = _unur_sample_cont_error;
1252 return INFINITY;
1253 }
1254 }
1255 else {
1256 /* splitting successful --> update guide table */
1257 _unur_arou_make_guide_table(gen);
1258 }
1259 }
1260 else
1261 /* no more construction points (avoid too many second if statements above) */
1262 GEN->max_segs = GEN->n_segs;
1263 }
1264
1265 /* if inside region of acceptance, return ratio x */
1266 if (u*u <= fx)
1267 return x;
1268 }
1269 }
1270 } /* end of _unur_arou_sample() */
1271
1272 /*---------------------------------------------------------------------------*/
1273
1274 double
_unur_arou_sample_check(struct unur_gen * gen)1275 _unur_arou_sample_check( struct unur_gen *gen )
1276 /*----------------------------------------------------------------------*/
1277 /* sample from generator and verify that method can be used */
1278 /* */
1279 /* parameters: */
1280 /* gen ... pointer to generator object */
1281 /* */
1282 /* return: */
1283 /* double (sample from random variate) */
1284 /* */
1285 /* error: */
1286 /* return INFINITY */
1287 /*----------------------------------------------------------------------*/
1288 {
1289 /** TODO: check uniform random number: u != 0 and u != 1 ?? **/
1290
1291 UNUR_URNG *urng; /* pointer to uniform random number generator */
1292 struct unur_arou_segment *seg;
1293 int result_split;
1294 double R,R1,R2,R3,tmp,x,fx,u,sqx,a;
1295
1296 /* check arguments */
1297 CHECK_NULL(gen,INFINITY); COOKIE_CHECK(gen,CK_AROU_GEN,INFINITY);
1298
1299 /* main URNG */
1300 urng = gen->urng;
1301
1302 while (1) {
1303
1304 /* sample from U(0,1) */
1305 R = _unur_call_urng(urng);
1306
1307 /* look up in guide table and search for segment */
1308 seg = GEN->guide[(int) (R * GEN->guide_size)];
1309 R *= GEN->Atotal;
1310 while (seg->Acum < R) {
1311 seg = seg->next;
1312 }
1313 COOKIE_CHECK(seg,CK_AROU_SEG,INFINITY);
1314
1315 /* reuse of uniform random number */
1316 R = seg->Acum - R;
1317
1318 /* inside or outside squeeze */
1319 if (R < seg->Ain) {
1320 /* inside */
1321 /* reuse of random number. */
1322 /* We can avoid R = (seg->Ain - R) / seg->Ain */
1323 x = ( ( seg->Ain * seg->rtp[0] + R * (seg->ltp[0] - seg->rtp[0]) ) /
1324 ( seg->Ain * seg->rtp[1] + R * (seg->ltp[1] - seg->rtp[1]) ) );
1325
1326 /* density at x */
1327 fx = PDF(x);
1328
1329 /* compute value of squeeze at x, i.e., we have to solve
1330 a*ltp[0] + (1-a)*rtp[0] == a*ltp[1] + (1-a)*rtp[1] */
1331 a = ( (seg->rtp[0] - x * seg->rtp[1]) /
1332 (seg->rtp[0] - seg->ltp[0] + x * (seg->ltp[1] - seg->rtp[1])) );
1333 sqx = a * seg->ltp[1] + (1.-a) * seg->rtp[1];
1334
1335 /* test for T-concavity */
1336 if (sqx*sqx > fx * (1.+UNUR_EPSILON))
1337 _unur_error(gen->genid,UNUR_ERR_GEN_CONDITION,"PDF not T-concave.");
1338
1339 return x;
1340 }
1341
1342 else {
1343 /* outside */
1344
1345 /* from now on we use the auxilliary generator
1346 (it can be the same as the main generator) */
1347 urng = gen->urng_aux;
1348
1349 /* three uniform random numbers with R1 + R2 + R3 = 1 */
1350 R1 = (R - seg->Ain) / seg->Aout; /* reuse of random number (good ?? ) */
1351 R2 = _unur_call_urng(urng);
1352 if (R1>R2) { tmp = R1; R1=R2; R2=tmp; } /* swap */
1353 R3 = 1.-R2;
1354 R2 -= R1;
1355
1356 /* point (v,u) and ratio x = v/u */
1357 u = seg->ltp[1]*R1 + seg->rtp[1]*R2 + seg->mid[1]*R3;
1358 x = (seg->ltp[0]*R1 + seg->rtp[0]*R2 + seg->mid[0]*R3) / u;
1359
1360 /* density at x */
1361 fx = PDF(x);
1362
1363 /* compute value of squeeze at x, i.e., we have to solve
1364 a*ltp[0] + (1-a)*rtp[0] == a*ltp[1] + (1-a)*rtp[1] */
1365 a = ( (seg->rtp[0] - x * seg->rtp[1]) /
1366 (seg->rtp[0] - seg->ltp[0] + x * (seg->ltp[1] - seg->rtp[1])) );
1367 sqx = a * seg->ltp[1] + (1.-a) * seg->rtp[1];
1368
1369 /* test for T-concavity */
1370 if (sqx*sqx > fx)
1371 _unur_error(gen->genid,UNUR_ERR_GEN_CONDITION,"PDF not T-concave.");
1372
1373 /* being outside the squeeze is bad. improve the situation! */
1374 if (GEN->n_segs < GEN->max_segs) {
1375 if (GEN->max_ratio * GEN->Atotal > GEN->Asqueeze) {
1376 result_split = _unur_arou_segment_split(gen,seg,x,fx);
1377 if ( !(result_split == UNUR_SUCCESS || result_split == UNUR_ERR_SILENT) ) {
1378 /* condition for PDF is violated! */
1379 _unur_error(gen->genid,UNUR_ERR_GEN_CONDITION,"");
1380 if (gen->variant & AROU_VARFLAG_PEDANTIC) {
1381 /* replace sampling routine by dummy routine that just returns INFINITY */
1382 SAMPLE = _unur_sample_cont_error;
1383 return INFINITY;
1384 }
1385 }
1386 else {
1387 /* splitting successful --> update guide table */
1388 _unur_arou_make_guide_table(gen);
1389 }
1390 }
1391 else
1392 /* no more construction points (avoid to many second if statement above */
1393 GEN->max_segs = GEN->n_segs;
1394 }
1395
1396 /* if inside region of acceptance, return ratio x */
1397 if (u*u <= fx)
1398 return x;
1399 }
1400 }
1401
1402 } /* end of _unur_arou_sample_check() */
1403
1404
1405 /*****************************************************************************/
1406 /** Auxilliary Routines **/
1407 /*****************************************************************************/
1408
1409 int
_unur_arou_get_starting_cpoints(struct unur_par * par,struct unur_gen * gen)1410 _unur_arou_get_starting_cpoints( struct unur_par *par, struct unur_gen *gen )
1411 /*----------------------------------------------------------------------*/
1412 /* list of construction points for starting segments. */
1413 /* if not provided as arguments compute these */
1414 /* by means of the "equi-angle rule". */
1415 /* */
1416 /* parameters: */
1417 /* par ... pointer to parameter for building generator object */
1418 /* gen ... pointer to generator object */
1419 /* */
1420 /* return: */
1421 /* UNUR_SUCCESS ... on success */
1422 /* error code ... on error */
1423 /*----------------------------------------------------------------------*/
1424 {
1425 struct unur_arou_segment *seg, *seg_new;
1426 double left_angle, right_angle, diff_angle, angle;
1427 double x, x_last, fx, fx_last;
1428 int i, use_center, is_center, is_increasing;
1429
1430 /* check arguments */
1431 CHECK_NULL(par,UNUR_ERR_NULL); COOKIE_CHECK(par,CK_AROU_PAR,UNUR_ERR_COOKIE);
1432 CHECK_NULL(gen,UNUR_ERR_NULL); COOKIE_CHECK(gen,CK_AROU_GEN,UNUR_ERR_COOKIE);
1433
1434 /* initialize boolean */
1435 is_center = FALSE;
1436
1437 /* use center as construction point ? */
1438 use_center = (gen->variant & AROU_VARFLAG_USECENTER) ? TRUE : FALSE;
1439
1440 /* reset counter of segments */
1441 GEN->n_segs = 0;
1442
1443 /* prepare for computing construction points */
1444 if (!PAR->starting_cpoints) {
1445 /* move center into x = 0 */
1446 /* angles of boundary of domain */
1447 left_angle = ( DISTR.BD_LEFT <= -INFINITY ) ? -M_PI/2. : atan(DISTR.BD_LEFT - GEN->center);
1448 right_angle = ( DISTR.BD_RIGHT >= INFINITY ) ? M_PI/2. : atan(DISTR.BD_RIGHT - GEN->center);
1449 /* we use equal distances between the angles of the cpoints */
1450 /* and the boundary points */
1451 diff_angle = (right_angle-left_angle) / (PAR->n_starting_cpoints + 1);
1452 angle = left_angle;
1453 }
1454 else
1455 diff_angle = angle = 0.; /* we do not need these variables in this case */
1456
1457 /* the left boundary point */
1458 x = x_last = DISTR.BD_LEFT;
1459 fx = fx_last = (x <= -INFINITY) ? 0. : PDF(x);
1460 seg = GEN->seg = _unur_arou_segment_new( gen, x, fx );
1461 if (seg == NULL) return UNUR_ERR_GEN_CONDITION; /* case of error */
1462 is_increasing = 1; /* assume PDF(x) is increasing for the first construction points */
1463
1464 /* now all the other points */
1465 for( i=0; i<=PAR->n_starting_cpoints; i++ ) {
1466
1467 /* starting point */
1468 if (i < PAR->n_starting_cpoints) {
1469 if (PAR->starting_cpoints) {
1470 /* construction points provided by user */
1471 x = PAR->starting_cpoints[i];
1472 /* check starting point */
1473 if (x <= DISTR.BD_LEFT || x >= DISTR.BD_RIGHT) {
1474 _unur_warning(gen->genid,UNUR_ERR_GEN_DATA,"starting point out of domain");
1475 continue;
1476 }
1477 if (x<=x_last) {
1478 _unur_warning(gen->genid,UNUR_ERR_GEN_DATA,"starting points not increasing -> skip");
1479 continue;
1480 }
1481 }
1482 else {
1483 /* compute construction points by means of "equidistance" rule */
1484 angle += diff_angle;
1485 x = tan( angle ) + GEN->center;
1486 }
1487 }
1488 else {
1489 /* the very last segment. it is rather a "virtual" segment to store
1490 the right vertex of the last segment, i.e., the right boundary point. */
1491 x = DISTR.BD_RIGHT;
1492 }
1493
1494 /* insert center ? */
1495 if (use_center && x >= GEN->center) {
1496 use_center = FALSE; /* we use the center only once (of course) */
1497 is_center = TRUE; /* the next construction point is the center */
1498 if (x>GEN->center) {
1499 x = GEN->center; /* use the center now ... */
1500 --i; /* and push the orignal starting point back on stack */
1501 if (!PAR->starting_cpoints)
1502 angle -= diff_angle; /* we have to compute the starting point in this case */
1503 }
1504 /* else: x == GEN->center --> nothing to do */
1505 }
1506 else
1507 is_center = FALSE;
1508
1509 /* value of PDF at starting point */
1510 fx = (x >= INFINITY) ? 0. : PDF(x);
1511
1512 /* check value of PDF at starting point */
1513 if (!is_increasing && fx > fx_last) {
1514 _unur_error(gen->genid,UNUR_ERR_GEN_CONDITION,"PDF not unimodal");
1515 return UNUR_ERR_GEN_CONDITION;
1516 }
1517
1518 if (fx <= 0. && fx_last <= 0.) {
1519 /* we do not need two such points */
1520 if (is_increasing) {
1521 /* PDF is still increasing, i.e., constant 0 til now */
1522 if (i<PAR->n_starting_cpoints) {
1523 /* and it is not the right boundary.
1524 otherwise the PDF is constant 0 on all construction points.
1525 then we need both boundary points. */
1526 /* we only have to change tangent line v/u = x,
1527 everything else remains unchanged */
1528 seg->dltp[0] = -1.;
1529 seg->dltp[1] = x;
1530 /* seg->dltp[0] = -1; seg->dltp[2] = 0.; not changed */
1531 x_last = x;
1532 continue; /* next construction point */
1533 }
1534 }
1535 else
1536 /* there should be no more points with PDF(x) > 0 */
1537 break;
1538 }
1539
1540 /* need a new segment */
1541 seg_new = _unur_arou_segment_new( gen, x, fx );
1542 if (seg_new == NULL) {
1543 /* case of error */
1544 seg->next = NULL; /* derminate list (for listing in debugging mode) */
1545 return UNUR_ERR_GEN_CONDITION;
1546 }
1547
1548 /* append to linked list */
1549 seg->next =seg_new;
1550 seg->rtp = seg_new->ltp;
1551 seg->drtp = seg_new->dltp;
1552
1553 /* next step */
1554 seg = seg_new;
1555
1556 /* PDF still increasing ? */
1557 if (is_increasing && fx < fx_last)
1558 is_increasing = 0;
1559
1560 /* store last computed values */
1561 x_last = x;
1562 fx_last = fx;
1563 }
1564
1565 /* we have left the loop with the right boundary of the support of PDF
1566 make shure that we will never use seg for sampling. */
1567 seg->Ain = seg->Aout = 0.;
1568 seg->Acum = INFINITY;
1569 seg->next = NULL; /* terminate list */
1570 --(GEN->n_segs); /* we do not count this segment */
1571
1572 /* o.k. */
1573 return UNUR_SUCCESS;
1574
1575 } /* end of _unur_arou_get_starting_cpoints() */
1576
1577 /*****************************************************************************/
1578
1579 int
_unur_arou_get_starting_segments(struct unur_gen * gen)1580 _unur_arou_get_starting_segments( struct unur_gen *gen )
1581 /*----------------------------------------------------------------------*/
1582 /* compute segments for starting points */
1583 /* */
1584 /* parameters: */
1585 /* gen ... pointer to generator object */
1586 /* */
1587 /* return: */
1588 /* UNUR_SUCCESS ... on success */
1589 /* error code ... on error */
1590 /*----------------------------------------------------------------------*/
1591 {
1592 #define MAX_IT (1000) /* maximal number of iterations to avoid
1593 infinite loop in case of numerical errors */
1594
1595 struct unur_arou_segment *seg, *seg_new, *seg_tmp;
1596 double x,fx; /* construction point, value of PDF at x */
1597 int n_it = 0; /* counter for number of iterations */
1598
1599 /* check arguments */
1600 CHECK_NULL(gen,UNUR_ERR_NULL); COOKIE_CHECK(gen,CK_AROU_GEN,UNUR_ERR_COOKIE);
1601
1602 /* compute paramters for all segments */
1603 for( seg=GEN->seg; seg->next != NULL; ) {
1604
1605 /* compute parameters for semgent */
1606 switch (_unur_arou_segment_parameter(gen, seg)) {
1607 case UNUR_SUCCESS: /* computation of parameters for segment successful */
1608 /* skip to next segment. */
1609 seg = seg->next;
1610 continue;
1611 case UNUR_ERR_SILENT: /* construction points too close */
1612 /* we have to remove this last segment from list */
1613 /* (the last construction point in the list is a boundary point.
1614 thus we might change the domain of the distribution.
1615 however, we only cut off a piece that is beyond the precesion
1616 of the floating point arithmetic.) */
1617 if (seg->next != NULL) {
1618 seg_tmp = seg->next;
1619 seg->next = seg->next->next;
1620 seg->rtp = seg->next->ltp;
1621 seg->drtp = seg->next->dltp;
1622 free(seg_tmp);
1623 --(GEN->n_segs);
1624 }
1625 else { /* seg->next==NULL */
1626 /* last (virtuel) interval in list.
1627 make shure that we will never use this segment */
1628 seg->Ain = seg->Aout = 0.;
1629 seg->Acum = INFINITY;
1630 }
1631 continue;
1632 case UNUR_ERR_INF: /* segment unbounded */
1633 /* split segment */
1634 break;
1635 default: /* PDF not T-concave */
1636 return UNUR_ERR_GEN_CONDITION;
1637 }
1638
1639 /* next iteration step */
1640 ++n_it;
1641
1642 if (n_it > MAX_IT) {
1643 /* maximal number of iterations exceeded */
1644 /* assume PDF not T-concave */
1645 return UNUR_ERR_GEN_CONDITION;
1646 }
1647
1648 /* area in segment infinite. insert new construction point. */
1649 x = _unur_arou_segment_arcmean(seg); /* use mean point in segment */
1650
1651 /* value of PDF at x */
1652 fx = PDF(x);
1653 if (fx < 0.) {
1654 _unur_error(gen->genid,UNUR_ERR_GEN_DATA,"PDF < 0");
1655 return UNUR_ERR_GEN_DATA;
1656 }
1657
1658 /* add a new segment, but check if we had to used too many segments */
1659 if (GEN->n_segs >= GEN->max_segs) {
1660 /* we do not want to create too many segments */
1661 _unur_error(gen->genid,UNUR_ERR_GEN_CONDITION,"cannot create bounded envelope!");
1662 return UNUR_ERR_GEN_CONDITION;
1663 }
1664 seg_new = _unur_arou_segment_new( gen, x, fx );
1665 if (seg_new == NULL) return UNUR_ERR_GEN_CONDITION; /* case of error */
1666
1667 /* insert into linked list */
1668 seg_new->next = seg->next;
1669 seg->next = seg_new;
1670
1671 /* right vertices */
1672 seg_new->rtp = seg->rtp;
1673 seg_new->drtp = seg->drtp;
1674 seg->rtp = seg_new->ltp;
1675 seg->drtp = seg_new->dltp;
1676
1677 }
1678
1679 /* o.k. */
1680 return UNUR_SUCCESS;
1681
1682 } /* end of _unur_arou_get_starting_segments() */
1683
1684 /*****************************************************************************/
1685
1686 struct unur_arou_segment *
_unur_arou_segment_new(struct unur_gen * gen,double x,double fx)1687 _unur_arou_segment_new( struct unur_gen *gen, double x, double fx )
1688 /*----------------------------------------------------------------------*/
1689 /* get new segment and compute left construction point at x. */
1690 /* */
1691 /* parameters: */
1692 /* gen ... pointer to generator object */
1693 /* x ... left point of new segment */
1694 /* fx ... value of PDF at x */
1695 /* */
1696 /* return: */
1697 /* pointer to new segment */
1698 /* */
1699 /* error: */
1700 /* return NULL */
1701 /*----------------------------------------------------------------------*/
1702 {
1703 struct unur_arou_segment *seg;
1704 double u,v,dfx;
1705
1706 /* check arguments */
1707 CHECK_NULL(gen,NULL); COOKIE_CHECK(gen,CK_AROU_GEN,NULL);
1708
1709 /* first check fx */
1710 if (fx<0.) {
1711 _unur_error(gen->genid,UNUR_ERR_GEN_DATA,"PDF(x) < 0.");
1712 return NULL;
1713 }
1714 if (_unur_FP_is_infinity(fx)) {
1715 /* over flow */
1716 _unur_error(gen->genid,UNUR_ERR_GEN_DATA,"PDF(x) overflow");
1717 return NULL;
1718 }
1719
1720 /* we need a new segment */
1721 seg = _unur_xmalloc( sizeof(struct unur_arou_segment) );
1722 seg->next = NULL; /* add eol marker */
1723 ++(GEN->n_segs); /* increment counter for segments */
1724 COOKIE_SET(seg,CK_AROU_SEG);
1725
1726 /* initialize some entries in seg */
1727 seg->Ain = seg->Aout = seg->Acum = 0.;
1728 seg->mid[0] = seg->mid[1] = 0.;
1729
1730 /* make left construction point in segment */
1731
1732 /* case: x out of support */
1733 if ( _unur_iszero(fx) ) {
1734 seg->ltp[0] = 0.; /* vertex == origin */
1735 seg->ltp[1] = 0.;
1736 if (x <= -INFINITY || x >= INFINITY ) {
1737 /* tangent line == line u = 0 (i.e., v-axis) */
1738 seg->dltp[0] = 0.; /* dv */
1739 seg->dltp[1] = 1.; /* du */
1740 seg->dltp[2] = 0.; /* v * dv + u * du */
1741 }
1742 else {
1743 /* tangent line == line v/u = x */
1744 seg->dltp[0] = -1.; /* dv */
1745 seg->dltp[1] = x; /* du */
1746 seg->dltp[2] = 0.; /* v * dv + u * du */
1747 }
1748 return seg;
1749 }
1750
1751 /* case: x in support */
1752 /* boundary point */
1753 u = sqrt( fx );
1754 v = x * u;
1755 seg->ltp[0] = v;
1756 seg->ltp[1] = u;
1757
1758 /* tangent line at tp */
1759
1760 /* compute derivative of PDF at tp x */
1761 dfx = dPDF(x);
1762
1763 /* subcase: derivative bounded
1764 use derivative for tangent line */
1765 if ( dfx > -INFINITY && dfx < INFINITY ) {
1766 seg->dltp[0] = -dfx / u; /* dv */ /** TODO: possible overflow **/
1767 seg->dltp[1] = 2 * u + dfx * x / u; /* du */ /** TODO: possible overflow **/
1768 seg->dltp[2] = seg->dltp[0] * v + seg->dltp[1] * u;
1769 return seg;
1770 }
1771
1772 /* subcase: derivative unbounded.
1773 use straight line through origin and vertex */
1774 seg->dltp[0] = -u; /* dv */
1775 seg->dltp[1] = v; /* du */
1776 seg->dltp[2] = 0.;
1777
1778 return seg;
1779
1780 } /* end of _unur_arou_segment_new() */
1781
1782 /*****************************************************************************/
1783
1784 /* maximal distance of intersection point from origin
1785 compared to distance of construction points to origin */
1786 #define MAX_NORM_OF_INTERSECTION_POINT 1.e6
1787
1788 /*---------------------------------------------------------------------------*/
1789
1790 int
_unur_arou_segment_parameter(struct unur_gen * gen,struct unur_arou_segment * seg)1791 _unur_arou_segment_parameter( struct unur_gen *gen, struct unur_arou_segment *seg )
1792 /*----------------------------------------------------------------------*/
1793 /* compute all parameters for a segment. */
1794 /* */
1795 /* parameters: */
1796 /* gen ... pointer to generator object */
1797 /* seg ... pointer to segment */
1798 /* */
1799 /* return: */
1800 /* UNUR_SUCCESS ... on success */
1801 /* UNUR_ERR_SILENT ... do not add this construction point */
1802 /* UNUR_ERR_INF ... if area = INFINITY */
1803 /* other error code ... error (PDF not T-concave) */
1804 /*----------------------------------------------------------------------*/
1805 {
1806 double coeff_det, cramer_det[2];
1807 double norm_vertex; /* sum of 1-norms of vertices */
1808 double diff_tangents; /* difference between coefficients of tangents */
1809 double det_bound; /* bound for determinant for Cramer's rule */
1810 double tmp_a, tmp_b;
1811
1812 /* check arguments */
1813 CHECK_NULL(gen,UNUR_ERR_NULL); COOKIE_CHECK(gen,CK_AROU_GEN,UNUR_ERR_COOKIE);
1814 CHECK_NULL(seg,UNUR_ERR_NULL); COOKIE_CHECK(seg,CK_AROU_SEG,UNUR_ERR_COOKIE);
1815
1816 /* sum of 1-norms of vertices */
1817 norm_vertex = fabs(seg->ltp[0]) + fabs(seg->ltp[1]) + fabs(seg->rtp[0]) + fabs(seg->rtp[1]);
1818
1819 /* area inside the squeeze */
1820 seg->Ain = (seg->ltp[1] * seg->rtp[0] - seg->ltp[0] * seg->rtp[1]) / 2.;
1821
1822 /* due to our ordering of construction points, seg->Ain must be >= 0 ! */
1823 if( seg->Ain < 0. ) {
1824 /* this could happen by round-off errors when the segment has angle extremely
1825 close to 0. Check this. */
1826 if (fabs(seg->Ain) < 1.e-8 * norm_vertex) {
1827 seg->Ain = seg->Aout = 0.;
1828 }
1829 else {
1830 /* This should not happen:
1831 non-ascending ordering of construction points */
1832 _unur_error(gen->genid,UNUR_ERR_SHOULD_NOT_HAPPEN,"");
1833 }
1834 return UNUR_ERR_SILENT;
1835 }
1836
1837 /* we use Cramer's rule to compute intersection of tangent lines.
1838 (well, we could save one multiplication otherwise) */
1839 coeff_det = seg->dltp[0] * seg->drtp[1] - seg->dltp[1] * seg->drtp[0];
1840 cramer_det[0] = seg->dltp[2] * seg->drtp[1] - seg->dltp[1] * seg->drtp[2];
1841 cramer_det[1] = seg->dltp[0] * seg->drtp[2] - seg->dltp[2] * seg->drtp[0];
1842
1843 /* we there are two possibilities for singular coefficent matrix:
1844
1845 either the outer triangle is unlimited.
1846 then the two tangents are distinct but parallel and
1847 the corresponding linear equation has no solution, i.e.
1848 coeff_det == 0 but cramer_det[0] != 0 or cramer_det[1] != 0.
1849
1850 or the outer triangle degenerates to a line segment and has area 0.
1851 then the two tangents are equal and
1852 the corresponding linear equation has a nonunique solution, i.e.
1853 coeff_det == cramer_det[0] == cramer_det[1] == 0.
1854 */
1855
1856 /* we to not allow that the outer triangles becomes too large.
1857 so if the 1-norm of intersection point is too large compared
1858 to norm_vertex we assume that this triangle is unbounded.
1859 we thus avoid numerical errors.
1860 (we use the 1-norm here since it much easier to handle.)
1861
1862 However, this might also happen due to roundoff errors,
1863 when the real position is extremely close to the secant.
1864 But at least we are on the save side.
1865 We only make an exception when coeff_det == 0, since otherwise there
1866 might be some problems with distributions like U(1,2).
1867 */
1868 det_bound = fabs(coeff_det) * norm_vertex * MAX_NORM_OF_INTERSECTION_POINT;
1869
1870 /* difference between coefficents of tangents */
1871 diff_tangents = ( fabs(seg->dltp[0] - seg->drtp[0]) + fabs(seg->dltp[1] - seg->drtp[1])
1872 + fabs(seg->dltp[2] - seg->drtp[2]) );
1873
1874 /* case: intersection point exists and is unique */
1875 if (!_unur_iszero(coeff_det) && !_unur_iszero(diff_tangents)) {
1876
1877 /* first check whether we can compute the intersection point */
1878 if ( fabs(cramer_det[0]) > det_bound || fabs(cramer_det[1]) > det_bound ) {
1879 /* case: triangle is assumed to be unbounded */
1880 /* _unur_warning(gen->genid,UNUR_ERR_GENERIC,"outer triangle assumed unbounded 1"); */
1881 seg->Aout = INFINITY;
1882 return UNUR_ERR_INF;
1883 }
1884
1885 /* compute intersection point */
1886 seg->mid[0] = cramer_det[0] / coeff_det;
1887 seg->mid[1] = cramer_det[1] / coeff_det;
1888
1889 /* area outside the squeeze */
1890 seg->Aout = ( (seg->ltp[0] - seg->mid[0]) * (seg->rtp[1] - seg->mid[1])
1891 - (seg->ltp[1] - seg->mid[1]) * (seg->rtp[0] - seg->mid[0])) / 2.;
1892
1893 /* due to our ordering of construction points, seg->Aout must be >= 0
1894 for a regular triangle.
1895 Thus if seg->Aout < 0, then the intersection point of tangents is on
1896 the WRONG side of the secant through vertices of segment,
1897 i.e. the "triangle" outside the squeeze region is unbounded.
1898
1899 However this might also happen due to roundoff errors when
1900 the real position is extremely close to the secant.
1901 We can distinguish between these two case by means of the u-coordinate
1902 of the intersection point. If it is on the wrong side of the secant,
1903 then we have seg->mid[1] < 0 (Otherwise we have either a round-off error
1904 or the PDF is not T-concave.) */
1905 if( seg->mid[1] < 0. ) {
1906 /* _unur_warning(gen->genid,UNUR_ERR_GENERIC,"outer triangle unbounded 2"); */
1907 seg->Aout = INFINITY;
1908 return UNUR_ERR_INF;
1909 }
1910
1911 /* at last check result.
1912 we must have:
1913 (*) seg->Aout > 0
1914 (*) intersection point right of left construction point
1915 (*) intersection point left of right construction point */
1916 if ( seg->Aout > 0. ) {
1917 tmp_a = seg->mid[0] * seg->ltp[1];
1918 tmp_b = seg->ltp[0] * seg->mid[1];
1919 if ( ! _unur_FP_less(tmp_a, tmp_b) ) {
1920 tmp_a = seg->mid[0] * seg->rtp[1];
1921 tmp_b = seg->rtp[0] * seg->mid[1];
1922 if ( ! _unur_FP_greater(tmp_a, tmp_b) )
1923 /* everything o.k. */
1924 return UNUR_SUCCESS;
1925 }
1926 }
1927
1928 /* there are two cases why the above check failed:
1929 (1) the PDF is not T-concave
1930 (2) small roundoff errors.
1931 */
1932
1933 if ( !_unur_iszero(seg->ltp[1]) && !_unur_iszero(seg->rtp[1]) ) {
1934 tmp_a = seg->ltp[0] * seg->rtp[1];
1935 tmp_b = seg->rtp[0] * seg->ltp[1];
1936 if ( _unur_FP_equal(tmp_a, tmp_b) ) {
1937 /* we assume that left construction point = right construction point */
1938 /* (construction points are too close) */
1939 seg->Ain = seg->Aout = 0.;
1940 return UNUR_ERR_SILENT;
1941 }
1942 }
1943
1944 /* _unur_warning(gen->genid,UNUR_ERR_GEN_CONDITION,"PDF not T-concave"); */
1945
1946 /* The outer area is unbounded (this might happen if the given PDF
1947 is not T-concave or due to round off errors.
1948 We assume round off errors. If the PDF is not T-concave we
1949 will exceed any bound for the number of segments. */
1950
1951 /* However due to round off errors, Aout might have become < 0
1952 when the boundary of region is (almost) a straight line
1953 and the area outside squeeze should be (almost) zero. */
1954
1955 if (!(fabs(seg->Aout) < fabs(seg->Ain) * UNUR_EPSILON)) {
1956 seg->Aout = INFINITY;
1957 return UNUR_ERR_INF;
1958 }
1959 /* else we assume round-off errors and set Aout = 0.
1960 i.e. go to the remaining case below */
1961 }
1962
1963 /* remaining case: triangle degenerates to a line segment, i.e.
1964 intersection point exists but is not unique */
1965
1966 /* boundary of region is (almost) a straight line
1967 and area outside squeeze is (almost) zero.
1968 use middle point as intersection point and
1969 set area outside the squeeze to 0.
1970 */
1971 /* _unur_warning(gen->genid,UNUR_ERR_GENERIC,"outer triangle is line"); */
1972 seg->mid[0] = 0.5 * (seg->ltp[0] + seg->rtp[0]);
1973 seg->mid[1] = 0.5 * (seg->ltp[1] + seg->rtp[1]);
1974 seg->Aout = 0.;
1975
1976 /* now it should be o.k. */
1977 return UNUR_SUCCESS;
1978
1979 } /* end if _unur_arou_segment_parameter() */
1980
1981 /*---------------------------------------------------------------------------*/
1982
1983 #undef MAX_NORM_INTERSECTION
1984
1985 /*****************************************************************************/
1986
1987 int
_unur_arou_segment_split(struct unur_gen * gen,struct unur_arou_segment * seg_oldl,double x,double fx)1988 _unur_arou_segment_split( struct unur_gen *gen, struct unur_arou_segment *seg_oldl, double x, double fx )
1989 /*----------------------------------------------------------------------*/
1990 /* insert new segment */
1991 /* old segment -> left hand side */
1992 /* new segment -> right hand side */
1993 /* */
1994 /* parameters: */
1995 /* gen ... pointer to generator object */
1996 /* seg_oldl ... pointer to segment */
1997 /* x ... left point of new segment */
1998 /* fx ... value of PDF at x */
1999 /* */
2000 /* return: */
2001 /* UNUR_SUCCESS ... on success */
2002 /* UNUR_ERR_SILENT ... if no intervals are splitted */
2003 /* other error codes ... on error */
2004 /* 1 ... if successful */
2005 /* 0 ... if */
2006 /* -1 ... error */
2007 /*----------------------------------------------------------------------*/
2008 {
2009 struct unur_arou_segment *seg_newr; /* pointer to newly created segment */
2010 struct unur_arou_segment seg_bak; /* space for saving data of segment */
2011 double Adiff;
2012
2013 /* check arguments */
2014 CHECK_NULL(gen,UNUR_ERR_NULL); COOKIE_CHECK(gen,CK_AROU_GEN,UNUR_ERR_COOKIE);
2015 CHECK_NULL(seg_oldl,UNUR_ERR_NULL); COOKIE_CHECK(seg_oldl,CK_AROU_SEG,UNUR_ERR_COOKIE);
2016
2017 #ifdef UNUR_ENABLE_LOGGING
2018 /* write info into LOG file */
2019 if (gen->debug & AROU_DEBUG_SPLIT)
2020 _unur_arou_debug_split_start( gen,seg_oldl,x,fx );
2021 #endif
2022
2023 /* we only add a new construction point, if the relative area is large enough */
2024 if (GEN->n_segs * seg_oldl->Aout / (GEN->Atotal - GEN->Asqueeze) < GEN->darsfactor )
2025 return UNUR_SUCCESS;
2026
2027 /* check for data error */
2028 if (fx < 0.) {
2029 _unur_error(gen->genid,UNUR_ERR_GEN_DATA,"PDF(x) < 0.!");
2030 return UNUR_ERR_GEN_DATA;
2031 }
2032
2033 /* backup data */
2034 memcpy(&seg_bak, seg_oldl, sizeof(struct unur_arou_segment));
2035
2036 /* PDF at x is 0. */
2037 if (fx <= 0.) {
2038 if (seg_oldl->rtp[1] <= 0. && seg_oldl->rtp[0] <= 0. ) {
2039 /* just chop off the right part of segment */
2040 /* we only have to change tangent line v/u = x
2041 at the right hand vertex */
2042 seg_oldl->drtp[1] = x; /* du */
2043 }
2044 else if (seg_oldl->ltp[1] <= 0. && seg_oldl->ltp[0] <= 0. ) {
2045 /* just chop off the left part of segment */
2046 /* we only have to change tangent line v/u = x
2047 at the left hand vertex */
2048 seg_oldl->dltp[1] = x; /* du */
2049 }
2050 else {
2051 _unur_error(gen->genid,UNUR_ERR_SHOULD_NOT_HAPPEN,"");
2052 return UNUR_ERR_SHOULD_NOT_HAPPEN;
2053 }
2054
2055 /* parameters of new segment */
2056 if( _unur_arou_segment_parameter(gen,seg_oldl)!=UNUR_SUCCESS ) {
2057 /* PDF not T-concave or area in segment not bounded */
2058 _unur_warning(gen->genid,UNUR_ERR_GEN_DATA,"Cannot chop segment at given point");
2059 /* error --> restore the old segment */
2060 memcpy(seg_oldl, &seg_bak, sizeof(struct unur_arou_segment));
2061 return UNUR_ERR_SILENT;
2062 }
2063
2064 /* for _unur_arou_debug_split_stop only */
2065 seg_newr = seg_oldl;
2066 }
2067
2068 else { /* fx > 0 */
2069
2070 /* need new segment */
2071 seg_newr = _unur_arou_segment_new(gen,x,fx);
2072 if (seg_newr == NULL) return UNUR_ERR_GEN_DATA; /* case of error */
2073
2074 /* link into list */
2075 seg_newr->next = seg_oldl->next;
2076 seg_oldl->next = seg_newr;
2077
2078 /* right vertices */
2079 seg_newr->rtp = seg_oldl->rtp;
2080 seg_newr->drtp = seg_oldl->drtp;
2081 seg_oldl->rtp = seg_newr->ltp;
2082 seg_oldl->drtp = seg_newr->dltp;
2083
2084 /* parameters of new segments */
2085 if( _unur_arou_segment_parameter(gen,seg_oldl)!=UNUR_SUCCESS ||
2086 /* PDF not T-concave or */
2087 _unur_arou_segment_parameter(gen,seg_newr)!=UNUR_SUCCESS
2088 /* area in segment not bounded */ ) {
2089
2090 /* new construction point not suitable --> do not add */
2091 _unur_warning(gen->genid,UNUR_ERR_GEN_DATA,"Cannot split segment at given point.");
2092 #ifdef UNUR_ENABLE_LOGGING
2093 /* write info into LOG file */
2094 if (gen->debug & AROU_DEBUG_SPLIT)
2095 _unur_arou_debug_split_stop( gen,seg_oldl,seg_newr );
2096 #endif
2097
2098 /* restore the old segment */
2099 memcpy(seg_oldl, &seg_bak, sizeof(struct unur_arou_segment));
2100
2101 /* decrement counter for segments and free unused segment */
2102 if (seg_newr) {
2103 --(GEN->n_segs);
2104 free( seg_newr );
2105 }
2106
2107 return UNUR_ERR_SILENT;
2108 }
2109 }
2110
2111 /* successful */
2112
2113 /* update total area below hat and squeeze */
2114 Adiff = - seg_bak.Ain + seg_oldl->Ain + ((seg_newr!=seg_oldl) ? seg_newr->Ain : 0. );
2115 GEN->Asqueeze += Adiff;
2116 Adiff += - seg_bak.Aout + seg_oldl->Aout + ((seg_newr!=seg_oldl) ? seg_newr->Aout : 0. );
2117 GEN->Atotal += Adiff;
2118
2119 #ifdef UNUR_ENABLE_LOGGING
2120 /* write info into LOG file */
2121 if (gen->debug & AROU_DEBUG_SPLIT)
2122 _unur_arou_debug_split_stop( gen,seg_oldl,seg_newr );
2123 #endif
2124
2125 /* o.k. */
2126 return UNUR_SUCCESS;
2127
2128 } /* end of _unur_arou_segment_split() */
2129
2130
2131 /*****************************************************************************/
2132
2133 double
_unur_arou_compute_x(double v,double u)2134 _unur_arou_compute_x( double v, double u )
2135 /*----------------------------------------------------------------------*/
2136 /* compute point x from (v,u) tuple. */
2137 /* */
2138 /* parameters: */
2139 /* v ... numerator */
2140 /* u ... denominator */
2141 /* */
2142 /* return: */
2143 /* point x of (x,y) tuple */
2144 /* */
2145 /* remark: */
2146 /* if (v,u)=(0,0) then x is set to INFINITY */
2147 /*----------------------------------------------------------------------*/
2148 {
2149 if (!_unur_iszero(u)) return v/u;
2150 else if (v<0.) return -INFINITY;
2151 else return INFINITY;
2152 } /* end of _unur_arou_compute_x() */
2153
2154 /*---------------------------------------------------------------------------*/
2155
2156 int
_unur_arou_run_dars(struct unur_gen * gen)2157 _unur_arou_run_dars( struct unur_gen *gen )
2158 /*----------------------------------------------------------------------*/
2159 /* run derandomized adaptive rejection sampling. */
2160 /* */
2161 /* parameters: */
2162 /* gen ... pointer to generator object */
2163 /* */
2164 /* return: */
2165 /* UNUR_SUCCESS ... on success */
2166 /* error code ... on error */
2167 /*----------------------------------------------------------------------*/
2168 {
2169 struct unur_arou_segment *seg, *seg_next;
2170 double Atot, Asqueezetot; /* total area below hat and squeeze, resp. */
2171 double Alimit; /* threshhold value for splitting interval */
2172 int n_splitted = 1; /* count splitted intervals */
2173 int splitted; /* result of splitting routine */
2174 double xl, xr; /* boundary of interval */
2175 double xsp, fxsp; /* splitting point in interval */
2176
2177 /* check arguments */
2178 CHECK_NULL(gen,UNUR_ERR_NULL); COOKIE_CHECK(gen,CK_AROU_GEN,UNUR_ERR_COOKIE);
2179
2180 /* there is no need to run DARS when the DARS factor is INFINITY */
2181 if (_unur_FP_is_infinity(GEN->darsfactor))
2182 return UNUR_SUCCESS;
2183
2184 /* first we need the total areas below hat and squeeze.
2185 (This is only necessary, when _unur_arou_make_guide_table() has not been
2186 called!) */
2187 Atot = 0.; /* area below hat */
2188 Asqueezetot = 0.; /* area below squeeze */
2189 for (seg = GEN->seg; seg != NULL; seg = seg->next ) {
2190 COOKIE_CHECK(seg,CK_AROU_SEG,UNUR_ERR_COOKIE);
2191 Asqueezetot += seg->Ain;
2192 Atot += seg->Ain + seg->Aout;
2193 }
2194 GEN->Atotal = Atot;
2195 GEN->Asqueeze = Asqueezetot;
2196
2197 /* now split intervals */
2198 while ( (GEN->max_ratio * GEN->Atotal > GEN->Asqueeze) &&
2199 (GEN->n_segs < GEN->max_segs) ) {
2200
2201 /* compute threshhold value. every interval with area between
2202 hat and squeeze greater than this value will be splitted. */
2203 if (GEN->n_segs > 1)
2204 Alimit = GEN->darsfactor * ( (GEN->Atotal - GEN->Asqueeze) / GEN->n_segs );
2205 else
2206 /* we split every interval if there are only one interval */
2207 Alimit = 0.;
2208
2209 /* reset counter for splitted intervals */
2210 n_splitted = 0;
2211
2212 /* for all intervals do ... */
2213 for (seg = GEN->seg; seg->next != NULL; seg = seg->next ) {
2214 COOKIE_CHECK(seg,CK_AROU_SEG,UNUR_ERR_COOKIE);
2215
2216 /* do not exceed the maximum number of intervals */
2217 if (GEN->n_segs >= GEN->max_segs)
2218 break;
2219
2220 /* we skip over all intervals where the area between hat and
2221 squeeze does not exceed the threshhold value. */
2222 if (seg->Aout <= Alimit)
2223 continue; /* goto next interval */
2224
2225 /* store pointer to next interval */
2226 seg_next = seg->next;
2227
2228 /* boundary of interval */
2229 xl = _unur_arou_compute_x(seg->ltp[0],seg->ltp[1]);
2230 xr = _unur_arou_compute_x(seg->rtp[0],seg->rtp[1]);
2231 if (xl>xr) xl = -INFINITY; /* (0,0) is mapped to INF instead of -INF */
2232
2233 /* However ... */
2234 if ( _unur_FP_is_minus_infinity(xl)
2235 && _unur_FP_same(seg->dltp[0],-1.) && _unur_iszero(seg->dltp[2]) )
2236 /* boundary of domain given by tangent */
2237 xl = seg->dltp[1];
2238
2239 if ( _unur_FP_is_infinity(xr)
2240 && _unur_FP_same(seg->drtp[0],-1.) && _unur_iszero(seg->drtp[2]) )
2241 /* boundary of domain given by tangent */
2242 xr = seg->drtp[1];
2243
2244 /* get splitting point (arc-mean rule) */
2245 xsp = _unur_arcmean(xl,xr);
2246
2247 /* value of PDF at splitting point */
2248 fxsp = PDF(xsp);
2249
2250 /* now split interval at given point */
2251 splitted = _unur_arou_segment_split(gen, seg, xsp, fxsp);
2252
2253 if (splitted == UNUR_SUCCESS) {
2254 /* splitting successful */
2255 ++n_splitted;
2256
2257 /* If we have chopped off the left or right tail,
2258 we can simply continue with seg->next.
2259 Otherwise, if we have split the segment, then
2260 seg->next points to the newly created segment
2261 and there is no need to split this in this loop;
2262 thus we skip seg->next and continue we seg->next->next.
2263 we can distinguish between these two cases by comparing
2264 the stored pointer seg_next with seg->next. */
2265 if (seg->next != seg_next)
2266 seg = seg->next;
2267 }
2268 else if (splitted != UNUR_ERR_SILENT) {
2269 /* some serious error occurred */
2270 /* condition for PDF is violated! */
2271 _unur_error(gen->genid,UNUR_ERR_GEN_CONDITION,"");
2272 return UNUR_ERR_GEN_CONDITION;
2273 }
2274 /* else: could not split construction points: too close (?) */
2275 }
2276
2277 if (n_splitted == 0) {
2278 /* we are not successful in splitting any inteval.
2279 abort to avoid endless loop */
2280 _unur_warning(gen->genid,UNUR_ERR_GENERIC,"DARS aborted: no intervals could be splitted.");
2281 break;
2282 }
2283 }
2284
2285 /* ratio between squeeze and hat o.k. ? */
2286 if ( GEN->max_ratio * GEN->Atotal > GEN->Asqueeze ) {
2287 if ( GEN->n_segs >= GEN->max_segs )
2288 _unur_warning(gen->genid,UNUR_ERR_GENERIC,"DARS aborted: maximum number of intervals exceeded.");
2289 _unur_warning(gen->genid,UNUR_ERR_GENERIC,"hat/squeeze ratio too small.");
2290 }
2291 else {
2292 /* no more construction points */
2293 GEN->max_segs = GEN->n_segs;
2294 }
2295
2296 /* o.k. */
2297 return UNUR_SUCCESS;
2298
2299 } /* end of _unur_arou_run_dars() */
2300
2301 /*****************************************************************************/
2302
2303 int
_unur_arou_make_guide_table(struct unur_gen * gen)2304 _unur_arou_make_guide_table( struct unur_gen *gen )
2305 /*----------------------------------------------------------------------*/
2306 /* make a guide table for indexed search */
2307 /* */
2308 /* parameters: */
2309 /* gen ... pointer to generator object */
2310 /* */
2311 /* return: */
2312 /* UNUR_SUCCESS ... on success */
2313 /* error code ... on error */
2314 /*----------------------------------------------------------------------*/
2315 {
2316 struct unur_arou_segment *seg;
2317 double Acum, Aincum, Astep;
2318 int max_guide_size;
2319 int j;
2320
2321 /* check arguments */
2322 CHECK_NULL(gen,UNUR_ERR_NULL); COOKIE_CHECK(gen,CK_AROU_GEN,UNUR_ERR_COOKIE);
2323
2324 /* allocate blocks for guide table (if necessary).
2325 (we allocate blocks for maximal guide table.) */
2326 if (!GEN->guide) {
2327 max_guide_size = (GEN->guide_factor > 0.) ? ((int)(GEN->max_segs * GEN->guide_factor)) : 1;
2328 if (max_guide_size <= 0) max_guide_size = 1; /* protect against overflow */
2329 GEN->guide = _unur_xmalloc( max_guide_size * sizeof(struct unur_arou_segment*) );
2330 }
2331
2332 /* first we need cumulated areas in segments */
2333 Acum = 0.; /* area in enveloping polygon */
2334 Aincum = 0.; /* area in squeeze */
2335 for (seg = GEN->seg; seg != NULL; seg = seg->next ) {
2336 COOKIE_CHECK(seg,CK_AROU_SEG,UNUR_ERR_COOKIE);
2337 Acum += seg->Ain + seg->Aout;
2338 Aincum += seg->Ain;
2339 seg->Acum = Acum;
2340 }
2341
2342 /* total area below hat */
2343 GEN->Atotal = Acum;
2344 GEN->Asqueeze = Aincum;
2345
2346 /* actual size of guide table */
2347 GEN->guide_size = (int)(GEN->n_segs * GEN->guide_factor);
2348 /* we do not vary the relative size of the guide table,
2349 since it has very little influence on speed */
2350
2351 /* make table (use variant 2; see dgt.c) */
2352 Astep = GEN->Atotal / GEN->guide_size;
2353 Acum=0.;
2354 for( j=0, seg=GEN->seg; j < GEN->guide_size; j++ ) {
2355 COOKIE_CHECK(seg,CK_AROU_SEG,UNUR_ERR_COOKIE);
2356 while( seg->Acum < Acum )
2357 if( seg->next != NULL ) /* skip to next segment if it exists */
2358 seg = seg->next;
2359 else {
2360 _unur_warning(gen->genid,UNUR_ERR_ROUNDOFF,"guide table");
2361 break;
2362 }
2363 GEN->guide[j] = seg;
2364 Acum += Astep;
2365 }
2366
2367 /* if there has been an round off error, we have to complete the guide table */
2368 for( ; j<GEN->guide_size ;j++ )
2369 GEN->guide[j] = seg;
2370
2371 return UNUR_SUCCESS;
2372 } /* end of _unur_arou_make_guide_table() */
2373
2374 /*****************************************************************************/
2375
2376 double
_unur_arou_segment_arcmean(struct unur_arou_segment * seg)2377 _unur_arou_segment_arcmean( struct unur_arou_segment *seg )
2378 /*----------------------------------------------------------------------*/
2379 /* compute "arctan mean" of two numbers expressed as v/u, u>=0 */
2380 /* */
2381 /* parameters: */
2382 /* seg ... pointer to segment */
2383 /* */
2384 /* return: */
2385 /* mean */
2386 /* */
2387 /* error: */
2388 /* return INFINITY */
2389 /* */
2390 /* comment: */
2391 /* "arctan mean" = tan(0.5*(arctan(x0)+arctan(x1))) */
2392 /* */
2393 /* a combination of arithmetical mean (for x0 and x1 close to 0) */
2394 /* and the harmonic mean (for |x0| and |x1| large). */
2395 /*----------------------------------------------------------------------*/
2396 {
2397 double xl, xr;
2398
2399 /* check arguments */
2400 CHECK_NULL(seg,INFINITY); COOKIE_CHECK(seg,CK_AROU_SEG,INFINITY);
2401
2402 /* if u != 0 ... x is stored in tp (= v/u) */
2403 /* else ... x is stored in tangent dltp */
2404 xl = (seg->ltp[1] > 0.) ? (seg->ltp[0] / seg->ltp[1]) :
2405 ( _unur_iszero(seg->dltp[0]) ? -INFINITY : (seg->dltp[1]) );
2406
2407 xr = (seg->rtp[1] > 0.) ? (seg->rtp[0] / seg->rtp[1]) :
2408 ( _unur_iszero(seg->drtp[0]) ? INFINITY : (seg->drtp[1]) );
2409
2410 return _unur_arcmean(xl,xr);
2411
2412 } /* end of _unur_arou_segment_arcmean() */
2413
2414 /*---------------------------------------------------------------------------*/
2415
2416 /*****************************************************************************/
2417 /** Debugging utilities **/
2418 /*****************************************************************************/
2419
2420 /*---------------------------------------------------------------------------*/
2421 #ifdef UNUR_ENABLE_LOGGING
2422 /*---------------------------------------------------------------------------*/
2423
2424 void
_unur_arou_debug_init(const struct unur_par * par,const struct unur_gen * gen)2425 _unur_arou_debug_init( const struct unur_par *par, const struct unur_gen *gen )
2426 /*----------------------------------------------------------------------*/
2427 /* write info about generator after setup into LOG file */
2428 /* */
2429 /* parameters: */
2430 /* par ... pointer to parameter for building generator object */
2431 /* gen ... pointer to generator object */
2432 /*----------------------------------------------------------------------*/
2433 {
2434 FILE *LOG;
2435 int i;
2436
2437 /* check arguments */
2438 CHECK_NULL(gen,RETURN_VOID); COOKIE_CHECK(gen,CK_AROU_GEN,RETURN_VOID);
2439 CHECK_NULL(par,RETURN_VOID); COOKIE_CHECK(par,CK_AROU_PAR,RETURN_VOID);
2440
2441 LOG = unur_get_stream();
2442
2443 fprintf(LOG,"%s:\n",gen->genid);
2444 fprintf(LOG,"%s: type = continuous univariate random variates\n",gen->genid);
2445 fprintf(LOG,"%s: method = ratio-of-uniforms method with enveloping polygon\n",gen->genid);
2446 fprintf(LOG,"%s:\n",gen->genid);
2447
2448 _unur_distr_cont_debug( gen->distr, gen->genid );
2449
2450 fprintf(LOG,"%s: sampling routine = _unur_arou_sample",gen->genid);
2451 if (gen->variant & AROU_VARFLAG_VERIFY)
2452 fprintf(LOG,"_check()\n");
2453 else
2454 fprintf(LOG,"()\n");
2455 fprintf(LOG,"%s:\n",gen->genid);
2456
2457 fprintf(LOG,"%s: center = %g",gen->genid,GEN->center);
2458 _unur_print_if_default(gen,AROU_SET_CENTER);
2459 if (gen->variant & AROU_VARFLAG_USECENTER)
2460 fprintf(LOG,"\n%s: use center as construction point",gen->genid);
2461 fprintf(LOG,"\n%s:\n",gen->genid);
2462
2463 fprintf(LOG,"%s: maximum number of segments = %d",gen->genid,GEN->max_segs);
2464 _unur_print_if_default(gen,AROU_SET_MAX_SEGS);
2465 fprintf(LOG,"\n%s: bound for ratio Asqueeze / Atotal = %g%%",gen->genid,GEN->max_ratio*100.);
2466 _unur_print_if_default(gen,AROU_SET_MAX_SQHRATIO);
2467 fprintf(LOG,"\n%s:\n",gen->genid);
2468
2469 if (gen->variant & AROU_VARFLAG_USEDARS) {
2470 fprintf(LOG,"%s: Derandomized ARS enabled ",gen->genid);
2471 _unur_print_if_default(gen,AROU_SET_USE_DARS);
2472 fprintf(LOG,"\n%s:\tDARS factor = %g",gen->genid,GEN->darsfactor);
2473 _unur_print_if_default(gen,AROU_SET_DARS_FACTOR);
2474 }
2475 else {
2476 fprintf(LOG,"%s: Derandomized ARS disabled ",gen->genid);
2477 _unur_print_if_default(gen,AROU_SET_USE_DARS);
2478 }
2479 fprintf(LOG,"\n%s:\n",gen->genid);
2480
2481 fprintf(LOG,"%s: sampling from list of segments: indexed search (guide table method)\n",gen->genid);
2482 fprintf(LOG,"%s: relative guide table size = %g%%",gen->genid,100.*GEN->guide_factor);
2483 _unur_print_if_default(gen,AROU_SET_GUIDEFACTOR);
2484 fprintf(LOG,"\n%s:\n",gen->genid);
2485
2486 fprintf(LOG,"%s: number of starting points = %d",gen->genid,PAR->n_starting_cpoints);
2487 _unur_print_if_default(gen,AROU_SET_N_STP);
2488 fprintf(LOG,"\n%s: starting points:",gen->genid);
2489 if (gen->set & AROU_SET_STP)
2490 for (i=0; i<PAR->n_starting_cpoints; i++) {
2491 if (i%5==0) fprintf(LOG,"\n%s:\t",gen->genid);
2492 fprintf(LOG," %#g,",PAR->starting_cpoints[i]);
2493 }
2494 else
2495 fprintf(LOG," use \"equidistribution\" rule [default]");
2496 fprintf(LOG,"\n%s:\n",gen->genid);
2497
2498 _unur_arou_debug_segments(gen);
2499
2500 fprintf(LOG,"%s: INIT completed **********************\n",gen->genid);
2501 fprintf(LOG,"%s:\n",gen->genid);
2502
2503 fflush(LOG);
2504
2505 } /* end of _unur_arou_debug_init() */
2506
2507 /*****************************************************************************/
2508
2509 void
_unur_arou_debug_dars_start(const struct unur_gen * gen)2510 _unur_arou_debug_dars_start( const struct unur_gen *gen )
2511 /*----------------------------------------------------------------------*/
2512 /* print header before runniung DARS into LOG file */
2513 /* */
2514 /* parameters: */
2515 /* gen ... pointer to generator object */
2516 /*----------------------------------------------------------------------*/
2517 {
2518 FILE *LOG;
2519
2520 /* check arguments */
2521 CHECK_NULL(gen,RETURN_VOID); COOKIE_CHECK(gen,CK_AROU_GEN,RETURN_VOID);
2522
2523 LOG = unur_get_stream();
2524
2525 fprintf(LOG,"%s: DARS started **********************\n",gen->genid);
2526 fprintf(LOG,"%s:\n",gen->genid);
2527 fprintf(LOG,"%s: DARS factor = %g",gen->genid,GEN->darsfactor);
2528 _unur_print_if_default(gen,AROU_SET_DARS_FACTOR);
2529 fprintf(LOG,"\n%s:\n",gen->genid);
2530
2531 fflush(LOG);
2532 } /* end of _unur_arou_debug_dars_start() */
2533
2534 /*---------------------------------------------------------------------------*/
2535
2536 void
_unur_arou_debug_dars(const struct unur_gen * gen)2537 _unur_arou_debug_dars( const struct unur_gen *gen )
2538 /*----------------------------------------------------------------------*/
2539 /* print infor after generator has run DARS into LOG file */
2540 /* */
2541 /* parameters: */
2542 /* gen ... pointer to generator object */
2543 /*----------------------------------------------------------------------*/
2544 {
2545 FILE *LOG;
2546
2547 /* check arguments */
2548 CHECK_NULL(gen,RETURN_VOID); COOKIE_CHECK(gen,CK_AROU_GEN,RETURN_VOID);
2549
2550 LOG = unur_get_stream();
2551
2552 fprintf(LOG,"%s:\n",gen->genid);
2553 fprintf(LOG,"%s: DARS finished **********************\n",gen->genid);
2554 fprintf(LOG,"%s:\n",gen->genid);
2555 _unur_arou_debug_segments(gen);
2556 fprintf(LOG,"%s:\n",gen->genid);
2557 fprintf(LOG,"%s: DARS completed **********************\n",gen->genid);
2558 fprintf(LOG,"%s:\n",gen->genid);
2559
2560 fflush(LOG);
2561 } /* end of _unur_arou_debug_dars() */
2562
2563 /*****************************************************************************/
2564
2565 void
_unur_arou_debug_free(const struct unur_gen * gen)2566 _unur_arou_debug_free( const struct unur_gen *gen )
2567 /*----------------------------------------------------------------------*/
2568 /* write info about generator before destroying into LOG file */
2569 /* */
2570 /* parameters: */
2571 /* gen ... pointer to generator object */
2572 /*----------------------------------------------------------------------*/
2573 {
2574 FILE *LOG;
2575
2576 /* check arguments */
2577 CHECK_NULL(gen,RETURN_VOID); COOKIE_CHECK(gen,CK_AROU_GEN,RETURN_VOID);
2578
2579 LOG = unur_get_stream();
2580
2581 fprintf(LOG,"%s:\n",gen->genid);
2582 fprintf(LOG,"%s: GENERATOR destroyed **********************\n",gen->genid);
2583 fprintf(LOG,"%s:\n",gen->genid);
2584 _unur_arou_debug_segments(gen);
2585 fprintf(LOG,"%s:\n",gen->genid);
2586
2587 fflush(LOG);
2588
2589 } /* end of _unur_arou_debug_free() */
2590
2591 /*****************************************************************************/
2592
2593 void
_unur_arou_debug_segments(const struct unur_gen * gen)2594 _unur_arou_debug_segments( const struct unur_gen *gen )
2595 /*----------------------------------------------------------------------*/
2596 /* write list of segments into LOGfile */
2597 /* */
2598 /* parameters: */
2599 /* gen ... pointer to generator object */
2600 /*----------------------------------------------------------------------*/
2601 {
2602 FILE *LOG;
2603 struct unur_arou_segment *seg;
2604 double sAin, sAout, Atotal;
2605 int i;
2606
2607 /* check arguments */
2608 CHECK_NULL(gen,RETURN_VOID); COOKIE_CHECK(gen,CK_AROU_GEN,RETURN_VOID);
2609
2610 LOG = unur_get_stream();
2611
2612 fprintf(LOG,"%s:Segments: %d\n",gen->genid,GEN->n_segs);
2613 if ((gen->debug & AROU_DEBUG_SEGMENTS) && GEN->seg != NULL) {
2614 fprintf(LOG,"%s: Nr.\t left touching point\t\t intersection point\t\t tangent at left touching point\n",gen->genid);
2615 for (seg = GEN->seg, i=0; seg->next!=NULL; seg=seg->next, i++) {
2616 COOKIE_CHECK(seg,CK_AROU_SEG,RETURN_VOID);
2617 fprintf(LOG,"%s:[%3d]: (%-12.6g,%-12.6g) (%-12.6g,%-12.6g) (%-12.6g,%-12.6g,%-12.6g)\n", gen->genid, i,
2618 seg->ltp[0],seg->ltp[1],
2619 seg->mid[0],seg->mid[1],
2620 seg->dltp[0],seg->dltp[1],seg->dltp[2]);
2621 }
2622 COOKIE_CHECK(seg,CK_AROU_SEG,RETURN_VOID);
2623 fprintf(LOG,"%s:[...]: (%-12.6g,%-12.6g)\n", gen->genid,seg->ltp[0],seg->ltp[1]);
2624 }
2625 fprintf(LOG,"%s:\n",gen->genid);
2626
2627 if (GEN->Atotal <= 0.) {
2628 fprintf(LOG,"%s: Construction of enveloping polygon not successful\n",gen->genid);
2629 fprintf(LOG,"%s: Areas may be meaningless !!!!!!!!!!!!!!!!!!!!!!!!\n",gen->genid);
2630 fprintf(LOG,"%s:\n",gen->genid);
2631 Atotal = -1.; /* to avoid floating point exceptions */
2632 }
2633 else {
2634 Atotal = GEN->Atotal;
2635 }
2636
2637 /* print and sum areas inside and outside of squeeze */
2638 if ((gen->debug & AROU_DEBUG_SEGMENTS) && GEN->seg != NULL) {
2639 fprintf(LOG,"%s:Areas in segments:\n",gen->genid);
2640 fprintf(LOG,"%s: Nr.\t inside squeeze\t\t outside squeeze\t total segment\t\tcumulated\n",gen->genid);
2641 sAin = sAout = 0.;
2642 for (seg = GEN->seg, i=0; seg->next!=NULL; seg=seg->next, i++) {
2643 COOKIE_CHECK(seg,CK_AROU_SEG,RETURN_VOID);
2644 sAin += seg->Ain;
2645 sAout += seg->Aout;
2646 fprintf(LOG,"%s:[%3d]: %-12.6g(%6.3f%%) | %-12.6g(%6.3f%%) | %-12.6g(%6.3f%%) | %-12.6g(%6.3f%%)\n",
2647 gen->genid,i,
2648 seg->Ain, seg->Ain * 100. / Atotal,
2649 seg->Aout, seg->Aout * 100. / Atotal,
2650 seg->Ain + seg->Aout, (seg->Ain + seg->Aout) * 100. / Atotal,
2651 seg->Acum, seg->Acum * 100. / Atotal);
2652 }
2653 fprintf(LOG,"%s:\t---------- --------- | ---------- --------- | ---------- --------- +\n",gen->genid);
2654 fprintf(LOG,"%s: Sum : %-12.6g(%6.3f%%) | %-12.6g(%6.3f%%) | %-11.6g(%6.3f%%)\n",
2655 gen->genid,
2656 sAin, sAin * 100./Atotal,
2657 sAout, sAout * 100./Atotal,
2658 sAin + sAout, (sAin + sAout) * 100./Atotal);
2659 fprintf(LOG,"%s:\n",gen->genid);
2660 }
2661
2662 /* summary of areas */
2663 fprintf(LOG,"%s: A(squeeze) = %-12.6g (%6.3f%%)\n",gen->genid,
2664 GEN->Asqueeze, GEN->Asqueeze * 100./Atotal);
2665 fprintf(LOG,"%s: A(hat\\squeeze) = %-12.6g (%6.3f%%)\n",gen->genid,
2666 GEN->Atotal - GEN->Asqueeze, (Atotal - GEN->Asqueeze) * 100./Atotal);
2667 fprintf(LOG,"%s: A(total) = %-12.6g\n",gen->genid, Atotal);
2668
2669 fprintf(LOG,"%s:\n",gen->genid);
2670
2671 } /* end of _unur_arou_debug_segments() */
2672
2673 /*****************************************************************************/
2674
2675 void
_unur_arou_debug_split_start(const struct unur_gen * gen,const struct unur_arou_segment * seg,double x,double fx)2676 _unur_arou_debug_split_start( const struct unur_gen *gen,
2677 const struct unur_arou_segment *seg,
2678 double x, double fx )
2679 /*----------------------------------------------------------------------*/
2680 /* write info about splitting segment */
2681 /* */
2682 /* parameters: */
2683 /* gen ... pointer to generator object */
2684 /* seg ... pointer to segment */
2685 /* x ... split at this point */
2686 /* fx ... value of PDF at x */
2687 /*----------------------------------------------------------------------*/
2688 {
2689 FILE *LOG;
2690
2691 /* check arguments */
2692 CHECK_NULL(gen,RETURN_VOID); COOKIE_CHECK(gen,CK_AROU_GEN,RETURN_VOID);
2693 CHECK_NULL(seg,RETURN_VOID); COOKIE_CHECK(seg,CK_AROU_SEG,RETURN_VOID);
2694
2695 LOG = unur_get_stream();
2696
2697 fprintf(LOG,"%s: split segment at x = %g \t\tf(x) = %g\n",gen->genid,x,fx);
2698 fprintf(LOG,"%s: old segment:\n",gen->genid);
2699
2700 fprintf(LOG,"%s: left construction point = (%-12.6g,%-12.6g)\t x = v/u = %-12.6g\tf(x) = %-12.6g\n",
2701 gen->genid, seg->ltp[0], seg->ltp[1], seg->ltp[0]/seg->ltp[1], sqrt(seg->ltp[1]) );
2702
2703 fprintf(LOG,"%s: intersection point = (%-12.6g,%-12.6g)\t x = v/u = %-12.6g\n",
2704 gen->genid, seg->mid[0], seg->mid[1], seg->mid[0]/seg->mid[1]);
2705
2706 fprintf(LOG,"%s: right construction point = (%-12.6g,%-12.6g)\t x = v/u = %-12.6g\tf(x) = %-12.6g\n",
2707 gen->genid, seg->rtp[0], seg->rtp[1], seg->rtp[0]/seg->rtp[1], sqrt(seg->rtp[1]) );
2708
2709 fprintf(LOG,"%s: A(squeeze) = %-12.6g\t(%6.3f%%)\n",gen->genid,
2710 seg->Ain, seg->Ain * 100./GEN->Atotal);
2711 fprintf(LOG,"%s: A(hat\\squeeze) = %-12.6g\t(%6.3f%%)\n",gen->genid,
2712 seg->Aout, seg->Aout * 100./GEN->Atotal);
2713 fprintf(LOG,"%s: A(hat) = %-12.6g\t(%6.3f%%)\n",gen->genid,
2714 (seg->Ain + seg->Aout), (seg->Ain +seg->Aout) * 100./GEN->Atotal);
2715
2716 fflush(LOG);
2717
2718 } /* end of _unur_arou_debug_split_start() */
2719
2720 /*****************************************************************************/
2721
2722 void
_unur_arou_debug_split_stop(const struct unur_gen * gen,const struct unur_arou_segment * seg_left,const struct unur_arou_segment * seg_right)2723 _unur_arou_debug_split_stop( const struct unur_gen *gen,
2724 const struct unur_arou_segment *seg_left,
2725 const struct unur_arou_segment *seg_right )
2726 /*----------------------------------------------------------------------*/
2727 /* write info about new splitted segments */
2728 /* */
2729 /* parameters: */
2730 /* gen ... pointer to generator object */
2731 /* iv_left ... pointer to new left hand segment */
2732 /* iv_right ... pointer to new right hand segment */
2733 /*----------------------------------------------------------------------*/
2734 {
2735 FILE *LOG;
2736 int chopped;
2737
2738 /* check arguments */
2739 CHECK_NULL(gen,RETURN_VOID); COOKIE_CHECK(gen,CK_AROU_GEN,RETURN_VOID);
2740 CHECK_NULL(seg_left,RETURN_VOID); COOKIE_CHECK(seg_left,CK_AROU_SEG,RETURN_VOID);
2741 CHECK_NULL(seg_right,RETURN_VOID); COOKIE_CHECK(seg_right,CK_AROU_SEG,RETURN_VOID);
2742
2743 LOG = unur_get_stream();
2744
2745 /* whether segment was split or just chopped */
2746 chopped = (seg_left==seg_right) ? 1 : 0;
2747
2748 if (chopped)
2749 fprintf(LOG,"%s: new segment (chopped):\n",gen->genid);
2750 else
2751 fprintf(LOG,"%s: new segments:\n",gen->genid);
2752
2753 fprintf(LOG,"%s: left construction point = (%-12.6g,%-12.6g)\t x = v/u = %-12.6g\tf(x) = %-12.6g\n",
2754 gen->genid, seg_left->ltp[0], seg_left->ltp[1], seg_left->ltp[0]/seg_left->ltp[1], sqrt(seg_left->ltp[1]) );
2755
2756 fprintf(LOG,"%s: intersection point = (%-12.6g,%-12.6g)\t x = v/u = %-12.6g\n",
2757 gen->genid, seg_left->mid[0], seg_left->mid[1], seg_left->mid[0]/seg_left->mid[1] );
2758
2759 if (chopped) {
2760 fprintf(LOG,"%s: right construction point = (%-12.6g,%-12.6g)\t x = v/u = %-12.6g\tf(x) = %-12.6g\n",
2761 gen->genid, seg_left->rtp[0], seg_left->rtp[1], seg_left->rtp[0]/seg_left->rtp[1], sqrt(seg_left->rtp[1]) );
2762 }
2763 else {
2764 fprintf(LOG,"%s: middle construction point = (%-12.6g,%-12.6g)\t x = v/u = %-12.6g\tf(x) = %-12.6g\n",
2765 gen->genid, seg_left->rtp[0], seg_left->rtp[1], seg_left->rtp[0]/seg_left->rtp[1], sqrt(seg_left->rtp[1]) );
2766
2767 fprintf(LOG,"%s: intersection point = (%-12.6g,%-12.6g)\t x = v/u = %-12.6g\n",
2768 gen->genid, seg_right->mid[0], seg_right->mid[1], seg_right->mid[0]/seg_right->mid[1] );
2769
2770 fprintf(LOG,"%s: right construction point = (%-12.6g,%-12.6g)\t x = v/u = %-12.6g\tf(x) = %-12.6g\n",
2771 gen->genid, seg_right->rtp[0], seg_right->rtp[1], seg_right->rtp[0]/seg_right->rtp[1], sqrt(seg_right->rtp[1]) );
2772 }
2773
2774 if (!chopped) {
2775 fprintf(LOG,"%s: left segment:\n",gen->genid);
2776 fprintf(LOG,"%s: A(squeeze) = %-12.6g\t(%6.3f%%)\n",gen->genid,
2777 seg_left->Ain, seg_left->Ain * 100./GEN->Atotal);
2778 fprintf(LOG,"%s: A(hat\\squeeze) = %-12.6g\t(%6.3f%%)\n",gen->genid,
2779 seg_left->Aout, seg_left->Aout * 100./GEN->Atotal);
2780 fprintf(LOG,"%s: A(hat) = %-12.6g\t(%6.3f%%)\n",gen->genid,
2781 (seg_left->Ain + seg_left->Aout), (seg_left->Ain +seg_left->Aout) * 100./GEN->Atotal);
2782
2783 fprintf(LOG,"%s: right segment:\n",gen->genid);
2784 fprintf(LOG,"%s: A(squeeze) = %-12.6g\t(%6.3f%%)\n",gen->genid,
2785 seg_right->Ain, seg_right->Ain * 100./GEN->Atotal);
2786 fprintf(LOG,"%s: A(hat\\squeeze) = %-12.6g\t(%6.3f%%)\n",gen->genid,
2787 seg_right->Aout, seg_right->Aout * 100./GEN->Atotal);
2788 fprintf(LOG,"%s: A(hat) = %-12.6g\t(%6.3f%%)\n",gen->genid,
2789 (seg_right->Ain + seg_right->Aout), (seg_right->Ain +seg_right->Aout) * 100./GEN->Atotal);
2790 }
2791
2792 fprintf(LOG,"%s: total areas:\n",gen->genid);
2793 fprintf(LOG,"%s: A(squeeze) = %-12.6g\t(%6.3f%%)\n",gen->genid,
2794 GEN->Asqueeze, GEN->Asqueeze * 100./GEN->Atotal);
2795 fprintf(LOG,"%s: A(hat\\squeeze) = %-12.6g\t(%6.3f%%)\n",gen->genid,
2796 GEN->Atotal - GEN->Asqueeze, (GEN->Atotal - GEN->Asqueeze) * 100./GEN->Atotal);
2797 fprintf(LOG,"%s: A(total) = %-12.6g\n",gen->genid, GEN->Atotal);
2798
2799 fprintf(LOG,"%s:\n",gen->genid);
2800
2801 fflush(LOG);
2802
2803 } /* end of _unur_arou_debug_split_stop() */
2804
2805 /*---------------------------------------------------------------------------*/
2806 #endif /* end UNUR_ENABLE_LOGGING */
2807 /*---------------------------------------------------------------------------*/
2808
2809 /*---------------------------------------------------------------------------*/
2810 #ifdef UNUR_ENABLE_INFO
2811 /*---------------------------------------------------------------------------*/
2812
2813 void
_unur_arou_info(struct unur_gen * gen,int help)2814 _unur_arou_info( struct unur_gen *gen, int help )
2815 /*----------------------------------------------------------------------*/
2816 /* create character string that contains information about the */
2817 /* given generator object. */
2818 /* */
2819 /* parameters: */
2820 /* gen ... pointer to generator object */
2821 /* help ... whether to print additional comments */
2822 /*----------------------------------------------------------------------*/
2823 {
2824 struct unur_string *info = gen->infostr;
2825 struct unur_distr *distr = gen->distr;
2826
2827 /* generator ID */
2828 _unur_string_append(info,"generator ID: %s\n\n", gen->genid);
2829
2830 /* distribution */
2831 _unur_string_append(info,"distribution:\n");
2832 _unur_distr_info_typename(gen);
2833 _unur_string_append(info," functions = PDF dPDF\n");
2834 _unur_string_append(info," domain = (%g, %g)\n", DISTR.domain[0],DISTR.domain[1]);
2835 _unur_string_append(info," center = %g", unur_distr_cont_get_center(distr));
2836 if ( !(distr->set & UNUR_DISTR_SET_CENTER) ) {
2837 if ( distr->set & UNUR_DISTR_SET_MODE )
2838 _unur_string_append(info," [= mode]\n");
2839 else
2840 _unur_string_append(info," [default]\n");
2841 }
2842 else {
2843 _unur_string_append(info,"\n");
2844 }
2845
2846
2847 if (help) {
2848 if ( !(distr->set & (UNUR_DISTR_SET_CENTER | UNUR_DISTR_SET_MODE )) )
2849 _unur_string_append(info,"\n[ Hint: %s ]\n",
2850 "You may provide a point near the mode as \"center\".");
2851 }
2852 _unur_string_append(info,"\n");
2853
2854 /* method */
2855 _unur_string_append(info,"method: AROU (Automatic Ratio-Of-Uniforms)\n");
2856 _unur_string_append(info,"\n");
2857
2858 /* performance */
2859 _unur_string_append(info,"performance characteristics:\n");
2860 _unur_string_append(info," area(hat) = %g\n", GEN->Atotal);
2861
2862 _unur_string_append(info," rejection constant ");
2863 if (distr->set & UNUR_DISTR_SET_PDFAREA)
2864 _unur_string_append(info,"= %g\n", GEN->Atotal/(0.5*DISTR.area));
2865 else
2866 _unur_string_append(info,"<= %g\n", GEN->Atotal/GEN->Asqueeze);
2867
2868 _unur_string_append(info," area ratio squeeze/hat = %g\n",
2869 GEN->Asqueeze/GEN->Atotal);
2870
2871 _unur_string_append(info," # segments = %d\n", GEN->n_segs);
2872 _unur_string_append(info,"\n");
2873
2874 /* parameters */
2875 if (help) {
2876 _unur_string_append(info,"parameters:\n");
2877
2878 _unur_string_append(info," max_sqhratio = %g %s\n", GEN->max_ratio,
2879 (gen->set & AROU_SET_MAX_SQHRATIO) ? "" : "[default]");
2880 _unur_string_append(info," max_segments = %d %s\n", GEN->max_segs_info,
2881 (gen->set & AROU_SET_MAX_SEGS) ? "" : "[default]");
2882
2883 if (gen->variant & AROU_VARFLAG_VERIFY)
2884 _unur_string_append(info," verify = on\n");
2885
2886 if (gen->variant & AROU_VARFLAG_PEDANTIC)
2887 _unur_string_append(info," pedantic = on\n");
2888
2889 /* Not displayed:
2890 int unur_arou_set_usedars( UNUR_PAR *parameters, int usedars );
2891 int unur_arou_set_darsfactor( UNUR_PAR *parameters, double factor );
2892 int unur_arou_set_cpoints( UNUR_PAR *parameters, int n_stp, const double *stp );
2893 int unur_arou_set_usecenter( UNUR_PAR *parameters, int usecenter );
2894 int unur_arou_set_guidefactor( UNUR_PAR *parameters, double factor );
2895 */
2896 _unur_string_append(info,"\n");
2897 }
2898
2899 /* Hints */
2900 if (help) {
2901 if ( !(gen->set & AROU_SET_MAX_SQHRATIO) )
2902 _unur_string_append(info,"[ Hint: %s ]\n",
2903 "You can set \"max_sqhratio\" closer to 1 to decrease rejection constant." );
2904 if (GEN->Asqueeze/GEN->Atotal < GEN->max_ratio)
2905 _unur_string_append(info,"[ Hint: %s ]\n",
2906 "You should increase \"max_segments\" to obtain the desired rejection constant." );
2907 _unur_string_append(info,"\n");
2908 }
2909
2910 } /* end of _unur_arou_info() */
2911
2912 /*---------------------------------------------------------------------------*/
2913 #endif /* end UNUR_ENABLE_INFO */
2914 /*---------------------------------------------------------------------------*/
2915