1 /*****************************************************************************
2 * *
3 * UNURAN -- Universal Non-Uniform Random number generator *
4 * *
5 *****************************************************************************
6 * *
7 * FILE: vnrou.c *
8 * *
9 * TYPE: continuous multivariate random variate *
10 * METHOD: naive ratio-of-uniforms method *
11 * *
12 * DESCRIPTION: *
13 * Given PDF and (optionally) a bounding rectangle for the acceptance *
14 * region. *
15 * Produce a value x consistent with its density *
16 * The bounding rectangle is computed numerically if it is not given. *
17 * *
18 * REQUIRED: *
19 * pointer to the density function *
20 * OPTIONAL: *
21 * mode of the density *
22 * bounding rectangle of acceptance region *
23 * *
24 *****************************************************************************
25 * *
26 * Copyright (c) 2000-2006 Wolfgang Hoermann and Josef Leydold *
27 * Department of Statistics and Mathematics, WU Wien, Austria *
28 * *
29 * This program is free software; you can redistribute it and/or modify *
30 * it under the terms of the GNU General Public License as published by *
31 * the Free Software Foundation; either version 2 of the License, or *
32 * (at your option) any later version. *
33 * *
34 * This program is distributed in the hope that it will be useful, *
35 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
36 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
37 * GNU General Public License for more details. *
38 * *
39 * You should have received a copy of the GNU General Public License *
40 * along with this program; if not, write to the *
41 * Free Software Foundation, Inc., *
42 * 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA *
43 * *
44 *****************************************************************************
45 * *
46 * REFERENCES: *
47 * [1] Wakefield J.C., Gelfand A.E., Smith A.F.M. *
48 * Efficient generation of random variates via the ratio-of-uniforms *
49 * method. *
50 * Statistics and Computing (1991) 1, pp (129-133) *
51 * *
52 * [2] Hoermann, W., Leydold J., and Derflinger, G. (2004): *
53 * Automatic non-uniform random variate generation, Springer, Berlin. *
54 * Section 2.4, Algorithm 2.9 (RoU), p.35 *
55 * *
56 *****************************************************************************/
57
58 /*---------------------------------------------------------------------------*/
59
60 #include <unur_source.h>
61 #include <distr/distr.h>
62 #include <distr/distr_source.h>
63 #include <distr/cvec.h>
64 #include <utils/fmax_source.h>
65 #include <utils/hooke_source.h>
66 #include <utils/matrix_source.h>
67 #include <utils/unur_fp_source.h>
68 #include <utils/mrou_rectangle_struct.h>
69 #include <utils/mrou_rectangle_source.h>
70 #include <urng/urng.h>
71 #include "unur_methods_source.h"
72 #include "x_gen_source.h"
73 #include "vnrou.h"
74 #include "vnrou_struct.h"
75
76 #ifdef UNUR_ENABLE_INFO
77 # include <tests/unuran_tests.h>
78 #endif
79
80 /*---------------------------------------------------------------------------*/
81 /* Variants: */
82
83 #define VNROU_VARFLAG_VERIFY 0x002u /* run verify mode */
84
85 /*---------------------------------------------------------------------------*/
86 /* Debugging flags */
87 /* bit 01 ... pameters and structure of generator (do not use here) */
88 /* bits 02-12 ... setup */
89 /* bits 13-24 ... adaptive steps */
90 /* bits 25-32 ... trace sampling */
91
92 #define VNROU_DEBUG_REINIT 0x00000010u /* print parameters after reinit */
93
94 /*---------------------------------------------------------------------------*/
95 /* Flags for logging set calls */
96
97 #define VNROU_SET_U 0x001u /* set u values of bounding rectangle */
98 #define VNROU_SET_V 0x002u /* set v values of bounding rectangle */
99 #define VNROU_SET_R 0x008u /* set r-parameter */
100
101 /*---------------------------------------------------------------------------*/
102
103 #define GENTYPE "VNROU" /* type of generator */
104
105 /*---------------------------------------------------------------------------*/
106
107 static struct unur_gen *_unur_vnrou_init( struct unur_par *par );
108 /*---------------------------------------------------------------------------*/
109 /* Initialize new generator. */
110 /*---------------------------------------------------------------------------*/
111
112 static int _unur_vnrou_reinit( struct unur_gen *gen );
113 /*---------------------------------------------------------------------------*/
114 /* Reinitialize generator. */
115 /*---------------------------------------------------------------------------*/
116
117 static struct unur_gen *_unur_vnrou_create( struct unur_par *par );
118 /*---------------------------------------------------------------------------*/
119 /* create new (almost empty) generator object. */
120 /*---------------------------------------------------------------------------*/
121
122 static struct unur_gen *_unur_vnrou_clone( const struct unur_gen *gen );
123 /*---------------------------------------------------------------------------*/
124 /* copy (clone) generator object. */
125 /*---------------------------------------------------------------------------*/
126
127 static void _unur_vnrou_free( struct unur_gen *gen);
128 /*---------------------------------------------------------------------------*/
129 /* destroy generator object. */
130 /*---------------------------------------------------------------------------*/
131
132 static int _unur_vnrou_sample_cvec( struct unur_gen *gen, double *vec );
133 static int _unur_vnrou_sample_check( struct unur_gen *gen, double *vec );
134 /*---------------------------------------------------------------------------*/
135 /* sample from generator. */
136 /*---------------------------------------------------------------------------*/
137
138 static int _unur_vnrou_rectangle( struct unur_gen *gen );
139 /*---------------------------------------------------------------------------*/
140 /* compute (minimal) bounding rectangle. */
141 /*---------------------------------------------------------------------------*/
142
143 #ifdef UNUR_ENABLE_LOGGING
144 /*---------------------------------------------------------------------------*/
145 /* the following functions print debugging information on output stream, */
146 /* i.e., into the LOG file if not specified otherwise. */
147 /*---------------------------------------------------------------------------*/
148 static void _unur_vnrou_debug_init( const struct unur_gen *gen );
149
150 /*---------------------------------------------------------------------------*/
151 /* print after generator has been initialized has completed. */
152 /*---------------------------------------------------------------------------*/
153 #endif
154
155 #ifdef UNUR_ENABLE_INFO
156 static void _unur_vnrou_info( struct unur_gen *gen, int help );
157 /*---------------------------------------------------------------------------*/
158 /* create info string. */
159 /*---------------------------------------------------------------------------*/
160 #endif
161
162 /*---------------------------------------------------------------------------*/
163 /* abbreviations */
164
165 #define DISTR_IN distr->data.cvec /* data for distribution object */
166
167 #define PAR ((struct unur_vnrou_par*)par->datap) /* data for parameter object */
168 #define GEN ((struct unur_vnrou_gen*)gen->datap) /* data for generator object */
169 #define DISTR gen->distr->data.cvec /* data for distribution in generator object */
170 #define SAMPLE gen->sample.cvec /* pointer to sampling routine */
171 #define PDF(x) _unur_cvec_PDF((x),(gen->distr)) /* call to PDF */
172
173 /*---------------------------------------------------------------------------*/
174
175 #define _unur_vnrou_getSAMPLE(gen) \
176 ( ((gen)->variant & VNROU_VARFLAG_VERIFY) \
177 ? _unur_vnrou_sample_check : _unur_vnrou_sample_cvec )
178
179 /*---------------------------------------------------------------------------*/
180
181 /*****************************************************************************/
182 /** Public: User Interface (API) **/
183 /*****************************************************************************/
184
185 struct unur_par *
unur_vnrou_new(const struct unur_distr * distr)186 unur_vnrou_new( const struct unur_distr *distr )
187 /*----------------------------------------------------------------------*/
188 /* get default parameters */
189 /* */
190 /* parameters: */
191 /* distr ... pointer to distribution object */
192 /* */
193 /* return: */
194 /* default parameters (pointer to structure) */
195 /* */
196 /* error: */
197 /* return NULL */
198 /*----------------------------------------------------------------------*/
199 {
200 struct unur_par *par;
201
202 /* check arguments */
203 _unur_check_NULL( GENTYPE,distr,NULL );
204
205 /* check distribution */
206 if (distr->type != UNUR_DISTR_CVEC) {
207 _unur_error(GENTYPE,UNUR_ERR_DISTR_INVALID,""); return NULL; }
208 COOKIE_CHECK(distr,CK_DISTR_CVEC,NULL);
209
210 if (DISTR_IN.pdf == NULL) {
211 _unur_error(GENTYPE,UNUR_ERR_DISTR_REQUIRED,"PDF");
212 return NULL;
213 }
214
215 /* allocate structure */
216 par = _unur_par_new( sizeof(struct unur_vnrou_par) );
217 COOKIE_SET(par,CK_VNROU_PAR);
218
219 /* copy input */
220 par->distr = distr; /* pointer to distribution object */
221
222 /* set default values */
223 PAR->r = 1.; /* r-parameter of the generalized method */
224 PAR->vmax = 0.; /* v-boundary of bounding rectangle (unknown) */
225 PAR->umin = NULL; /* u-boundary of bounding rectangle (unknown) */
226 PAR->umax = NULL; /* u-boundary of bounding rectangle (unknown) */
227 par->method = UNUR_METH_VNROU; /* method and default variant */
228 par->variant = 0u; /* default variant */
229 par->set = 0u; /* inidicate default parameters */
230 par->urng = unur_get_default_urng(); /* use default urng */
231 par->urng_aux = NULL; /* no auxilliary URNG required */
232 par->debug = _unur_default_debugflag; /* set default debugging flags */
233
234 /* routine for starting generator */
235 par->init = _unur_vnrou_init;
236
237 return par;
238
239 } /* end of unur_vnrou_new() */
240
241 /*****************************************************************************/
242
243 int
unur_vnrou_set_u(struct unur_par * par,double * umin,double * umax)244 unur_vnrou_set_u( struct unur_par *par, double *umin, double *umax )
245 /*----------------------------------------------------------------------*/
246 /* Sets left and right u-boundary of bounding rectangle. */
247 /* */
248 /* parameters: */
249 /* par ... pointer to parameter for building generator object */
250 /* umin ... left boundary of rectangle */
251 /* umax ... right boundary of rectangle */
252 /* */
253 /* return: */
254 /* UNUR_SUCCESS ... on success */
255 /* error code ... on error */
256 /*----------------------------------------------------------------------*/
257 {
258 int d; /* index used in dimension loops (0 <= d < dim) */
259
260 /* check arguments */
261 _unur_check_NULL( GENTYPE, par, UNUR_ERR_NULL );
262 _unur_check_par_object( par, VNROU );
263 _unur_check_NULL( GENTYPE, umin, UNUR_ERR_NULL );
264 _unur_check_NULL( GENTYPE, umax, UNUR_ERR_NULL );
265
266 /* check new parameter for generator */
267 for (d=0; d<par->distr->dim; d++) {
268 if (!_unur_FP_greater(umax[d],umin[d])) {
269 _unur_warning(GENTYPE,UNUR_ERR_PAR_SET,"umax <= umin");
270 return UNUR_ERR_PAR_SET;
271 }
272 }
273
274 /* set values */
275 PAR->umin = umin;
276 PAR->umax = umax;
277
278 /* changelog */
279 par->set |= VNROU_SET_U;
280
281 return UNUR_SUCCESS;
282
283 } /* end of unur_vnrou_set_u() */
284
285 /*---------------------------------------------------------------------------*/
286
287 int
unur_vnrou_chg_u(struct unur_gen * gen,double * umin,double * umax)288 unur_vnrou_chg_u( struct unur_gen *gen, double *umin, double *umax )
289 /*----------------------------------------------------------------------*/
290 /* Sets left and right u-boundary of bounding rectangle. */
291 /* */
292 /* parameters: */
293 /* gen ... pointer to generator object */
294 /* umin ... left boundary of rectangle */
295 /* umax ... right boundary of rectangle */
296 /* */
297 /* return: */
298 /* UNUR_SUCCESS ... on success */
299 /* error code ... on error */
300 /*----------------------------------------------------------------------*/
301 {
302 int d; /* index used in dimension loops (0 <= d < dim) */
303
304 /* check arguments */
305 _unur_check_NULL( GENTYPE, gen, UNUR_ERR_NULL );
306 _unur_check_gen_object( gen, VNROU, UNUR_ERR_GEN_INVALID );
307 _unur_check_NULL( GENTYPE, umin, UNUR_ERR_NULL );
308 _unur_check_NULL( GENTYPE, umax, UNUR_ERR_NULL );
309
310 /* check new parameter for generator */
311 for (d=0; d<GEN->dim; d++) {
312 if (!_unur_FP_greater(umax[d],umin[d])) {
313 _unur_warning(GENTYPE,UNUR_ERR_PAR_SET,"umax <= umin");
314 return UNUR_ERR_PAR_SET;
315 }
316 }
317
318 /* set values */
319 memcpy(GEN->umin, umin, GEN->dim * sizeof(double));
320 memcpy(GEN->umax, umax, GEN->dim * sizeof(double));
321
322 /* changelog */
323 gen->set |= VNROU_SET_U;
324
325 return UNUR_SUCCESS;
326
327 } /* end of unur_vnrou_chg_u() */
328
329 /*---------------------------------------------------------------------------*/
330
331 int
unur_vnrou_set_v(struct unur_par * par,double vmax)332 unur_vnrou_set_v( struct unur_par *par, double vmax )
333 /*----------------------------------------------------------------------*/
334 /* Sets upper v-boundary of bounding rectangle. */
335 /* */
336 /* parameters: */
337 /* par ... pointer to parameter for building generator object */
338 /* vmax ... upper boundary of rectangle */
339 /* */
340 /* return: */
341 /* UNUR_SUCCESS ... on success */
342 /* error code ... on error */
343 /*----------------------------------------------------------------------*/
344 {
345 /* check arguments */
346 _unur_check_NULL( GENTYPE, par, UNUR_ERR_NULL );
347 _unur_check_par_object( par, VNROU );
348
349 /* check new parameter for generator */
350 if (vmax <= 0.) {
351 _unur_warning(GENTYPE,UNUR_ERR_PAR_SET,"vmax <= 0");
352 return UNUR_ERR_PAR_SET;
353 }
354
355 /* store values */
356 PAR->vmax = vmax;
357
358 /* changelog */
359 par->set |= VNROU_SET_V;
360
361 return UNUR_SUCCESS;
362
363 } /* end of unur_vnrou_set_v() */
364
365 /*---------------------------------------------------------------------------*/
366
367 int
unur_vnrou_chg_v(struct unur_gen * gen,double vmax)368 unur_vnrou_chg_v( struct unur_gen *gen, double vmax )
369 /*----------------------------------------------------------------------*/
370 /* Sets upper v-boundary of bounding rectangle. */
371 /* */
372 /* parameters: */
373 /* gen ... pointer to generator object */
374 /* vmax ... upper boundary of rectangle */
375 /* */
376 /* return: */
377 /* UNUR_SUCCESS ... on success */
378 /* error code ... on error */
379 /*----------------------------------------------------------------------*/
380 {
381 /* check arguments */
382 _unur_check_NULL( GENTYPE, gen, UNUR_ERR_NULL );
383 _unur_check_gen_object( gen, VNROU, UNUR_ERR_GEN_INVALID );
384
385 /* check new parameter for generator */
386 if (vmax <= 0.) {
387 _unur_warning(GENTYPE,UNUR_ERR_PAR_SET,"vmax <= 0");
388 return UNUR_ERR_PAR_SET;
389 }
390
391 /* store values */
392 GEN->vmax = vmax;
393
394 /* changelog */
395 gen->set |= VNROU_SET_V;
396
397 return UNUR_SUCCESS;
398
399 } /* end of unur_vnrou_chg_v() */
400
401 /*---------------------------------------------------------------------------*/
402
403 int
unur_vnrou_set_r(struct unur_par * par,double r)404 unur_vnrou_set_r( struct unur_par *par, double r )
405 /*----------------------------------------------------------------------*/
406 /* Set the r-parameter for the generalized ratio-of-uniforms method. */
407 /* */
408 /* parameters: */
409 /* par ... pointer to parameter for building generator object */
410 /* r ... r-parameter */
411 /* */
412 /* return: */
413 /* UNUR_SUCCESS ... on success */
414 /* error code ... on error */
415 /*----------------------------------------------------------------------*/
416 {
417 /* check arguments */
418 _unur_check_NULL( GENTYPE, par, UNUR_ERR_NULL );
419 _unur_check_par_object( par, VNROU );
420
421 /* check new parameter for generator */
422 if (r <= 0.) {
423 _unur_warning(GENTYPE,UNUR_ERR_PAR_SET,"r<=0");
424 return UNUR_ERR_PAR_SET;
425 }
426
427 /* store data */
428 PAR->r = r;
429
430 /* changelog */
431 par->set |= VNROU_SET_R;
432
433 /* o.k. */
434 return UNUR_SUCCESS;
435
436 } /* end of unur_vnrou_set_r() */
437
438 /*---------------------------------------------------------------------------*/
439
440 int
unur_vnrou_set_verify(struct unur_par * par,int verify)441 unur_vnrou_set_verify( struct unur_par *par, int verify )
442 /*----------------------------------------------------------------------*/
443 /* turn verifying of algorithm while sampling on/off */
444 /* */
445 /* parameters: */
446 /* par ... pointer to parameter for building generator object */
447 /* verify ... 0 = no verifying, !0 = verifying */
448 /* */
449 /* return: */
450 /* UNUR_SUCCESS ... on success */
451 /* error code ... on error */
452 /* */
453 /* comment: */
454 /* no verifying is the default */
455 /*----------------------------------------------------------------------*/
456 {
457 /* check arguments */
458 _unur_check_NULL( GENTYPE, par, UNUR_ERR_NULL );
459 _unur_check_par_object( par, VNROU );
460
461 /* we use a bit in variant */
462 par->variant = (verify) ? (par->variant | VNROU_VARFLAG_VERIFY) : (par->variant & (~VNROU_VARFLAG_VERIFY));
463
464 /* o.k. */
465 return UNUR_SUCCESS;
466
467 } /* end of unur_vnrou_set_verify() */
468
469 /*---------------------------------------------------------------------------*/
470
471 int
unur_vnrou_chg_verify(struct unur_gen * gen,int verify)472 unur_vnrou_chg_verify( struct unur_gen *gen, int verify )
473 /*----------------------------------------------------------------------*/
474 /* turn verifying of algorithm while sampling on/off */
475 /* */
476 /* parameters: */
477 /* gen ... pointer to generator object */
478 /* verify ... 0 = no verifying, !0 = verifying */
479 /* */
480 /* return: */
481 /* UNUR_SUCCESS ... on success */
482 /* error code ... on error */
483 /* */
484 /* comment: */
485 /* no verifying is the default */
486 /*----------------------------------------------------------------------*/
487 {
488 /* check input */
489 _unur_check_NULL( GENTYPE, gen, UNUR_ERR_NULL );
490 _unur_check_gen_object( gen, VNROU, UNUR_ERR_GEN_INVALID );
491
492 /* we must not change this switch when sampling has been disabled by
493 using a pointer to the error producing routine */
494 if (SAMPLE == _unur_sample_cvec_error)
495 return UNUR_FAILURE;
496
497 if (verify)
498 /* turn verify bounding rectangle on */
499 gen->variant |= VNROU_VARFLAG_VERIFY;
500 else
501 /* turn verify bounding rectangle off */
502 gen->variant &= ~VNROU_VARFLAG_VERIFY;
503
504 SAMPLE = _unur_vnrou_getSAMPLE(gen);
505
506 /* o.k. */
507 return UNUR_SUCCESS;
508
509 } /* end of unur_vnrou_chg_verify() */
510
511 /*---------------------------------------------------------------------------*/
512
513 double
unur_vnrou_get_volumehat(const struct unur_gen * gen)514 unur_vnrou_get_volumehat( const struct unur_gen *gen )
515 /*----------------------------------------------------------------------*/
516 /* Get the volume of below the hat. */
517 /* For normalized densities, i.e. when the volume below PDF is 1, */
518 /* this value equals the rejection constant for the vnrou method. */
519 /* */
520 /* parameters: */
521 /* gen ... pointer to generator object */
522 /* */
523 /*----------------------------------------------------------------------*/
524 {
525 double vol;
526 int d;
527
528 /* check arguments */
529 _unur_check_NULL( GENTYPE, gen, INFINITY );
530 _unur_check_gen_object( gen, VNROU, INFINITY );
531
532 /* compute volume of bounding rectangle */
533 vol = GEN->vmax;
534 for (d=0; d<GEN->dim; d++) {
535 vol *= (GEN->umax[d]-GEN->umin[d]);
536 }
537 /* compute volume of corresponding hat function */
538 vol *= (GEN->r*GEN->dim+1);
539
540 /* return result */
541 return vol;
542 } /* end of unur_vnrou_get_volumehat() */
543
544
545 /*****************************************************************************/
546 /** Private **/
547 /*****************************************************************************/
548
549 struct unur_gen *
_unur_vnrou_init(struct unur_par * par)550 _unur_vnrou_init( struct unur_par *par )
551 /*----------------------------------------------------------------------*/
552 /* initialize new generator */
553 /* */
554 /* parameters: */
555 /* params ... pointer to paramters for building generator object */
556 /* */
557 /* return: */
558 /* pointer to generator object */
559 /* */
560 /* error: */
561 /* return NULL */
562 /*----------------------------------------------------------------------*/
563 {
564 struct unur_gen *gen;
565
566 /* check arguments */
567 CHECK_NULL(par,NULL);
568
569 /* check input */
570 if ( par->method != UNUR_METH_VNROU ) {
571 _unur_error(GENTYPE,UNUR_ERR_PAR_INVALID,"");
572 return NULL; }
573 COOKIE_CHECK(par,CK_VNROU_PAR,NULL);
574
575 /* create a new empty generator object */
576 gen = _unur_vnrou_create(par);
577 _unur_par_free(par);
578 if (!gen) return NULL;
579
580 /* compute bounding rectangle */
581 if (_unur_vnrou_rectangle(gen)!=UNUR_SUCCESS) {
582 _unur_vnrou_free(gen); return NULL;
583 }
584
585 #ifdef UNUR_ENABLE_LOGGING
586 /* write info into LOG file */
587 if (gen->debug) _unur_vnrou_debug_init(gen);
588 #endif
589
590 return gen;
591
592 } /* end of _unur_vnrou_init() */
593
594 /*---------------------------------------------------------------------------*/
595
596 int
_unur_vnrou_reinit(struct unur_gen * gen)597 _unur_vnrou_reinit( struct unur_gen *gen )
598 /*----------------------------------------------------------------------*/
599 /* re-initialize (existing) generator. */
600 /* */
601 /* parameters: */
602 /* gen ... pointer to generator object */
603 /* */
604 /* return: */
605 /* UNUR_SUCCESS ... on success */
606 /* error code ... on error */
607 /*----------------------------------------------------------------------*/
608 {
609 int rcode;
610
611 /* compute bounding rectangle */
612 if ( (rcode = _unur_vnrou_rectangle(gen))!=UNUR_SUCCESS) {
613 return rcode;
614 }
615
616 /* (re)set sampling routine */
617 SAMPLE = _unur_vnrou_getSAMPLE(gen);
618
619 #ifdef UNUR_ENABLE_LOGGING
620 /* write info into LOG file */
621 if (gen->debug & VNROU_DEBUG_REINIT) _unur_vnrou_debug_init(gen);
622 #endif
623
624 return UNUR_SUCCESS;
625 } /* end of _unur_vnrou_reinit() */
626
627 /*---------------------------------------------------------------------------*/
628
629 struct unur_gen *
_unur_vnrou_create(struct unur_par * par)630 _unur_vnrou_create( struct unur_par *par )
631 /*----------------------------------------------------------------------*/
632 /* allocate memory for generator */
633 /* */
634 /* parameters: */
635 /* par ... pointer to parameter for building generator object */
636 /* */
637 /* return: */
638 /* pointer to (empty) generator object with default settings */
639 /* */
640 /* error: */
641 /* return NULL */
642 /*----------------------------------------------------------------------*/
643 {
644 struct unur_gen *gen;
645
646 /* check arguments */
647 CHECK_NULL(par,NULL); COOKIE_CHECK(par,CK_VNROU_PAR,NULL);
648
649 /* create new generic generator object */
650 gen = _unur_generic_create( par, sizeof(struct unur_vnrou_gen) );
651
652 /* magic cookies */
653 COOKIE_SET(gen,CK_VNROU_GEN);
654
655 /* set generator identifier */
656 gen->genid = _unur_set_genid(GENTYPE);
657
658 /* routines for sampling and destroying generator */
659 SAMPLE = _unur_vnrou_getSAMPLE(gen);
660 gen->destroy = _unur_vnrou_free;
661 gen->clone = _unur_vnrou_clone;
662 gen->reinit = _unur_vnrou_reinit;
663
664 /* copy parameters into generator object */
665 GEN->dim = gen->distr->dim; /* dimension */
666 GEN->r = PAR->r; /* r-parameter of the vnrou method */
667 GEN->vmax = PAR->vmax; /* upper v-boundary of bounding rectangle */
668
669 /* allocate memory for u-boundary arrays */
670 GEN->umin = _unur_xmalloc( GEN->dim * sizeof(double)); /* bounding rectangle */
671 GEN->umax = _unur_xmalloc( GEN->dim * sizeof(double)); /* bounding rectangle */
672
673 if (PAR->umin != NULL) memcpy(GEN->umin, PAR->umin, GEN->dim * sizeof(double));
674 if (PAR->umax != NULL) memcpy(GEN->umax, PAR->umax, GEN->dim * sizeof(double));
675
676 /* get center of the distribution */
677 GEN->center = unur_distr_cvec_get_center(gen->distr);
678
679 #ifdef UNUR_ENABLE_INFO
680 /* set function for creating info string */
681 gen->info = _unur_vnrou_info;
682 #endif
683
684 /* return pointer to (almost empty) generator object */
685 return gen;
686
687 } /* end of _unur_vnrou_create() */
688
689 /*---------------------------------------------------------------------------*/
690
691 struct unur_gen *
_unur_vnrou_clone(const struct unur_gen * gen)692 _unur_vnrou_clone( const struct unur_gen *gen )
693 /*----------------------------------------------------------------------*/
694 /* copy (clone) generator object */
695 /* */
696 /* parameters: */
697 /* gen ... pointer to generator object */
698 /* */
699 /* return: */
700 /* pointer to clone of generator object */
701 /* */
702 /* error: */
703 /* return NULL */
704 /*----------------------------------------------------------------------*/
705 {
706 #define CLONE ((struct unur_vnrou_gen*)clone->datap)
707
708 struct unur_gen *clone;
709
710 /* check arguments */
711 CHECK_NULL(gen,NULL); COOKIE_CHECK(gen,CK_VNROU_GEN,NULL);
712
713 /* create generic clone */
714 clone = _unur_generic_clone( gen, GENTYPE );
715
716 /* allocate memory for u-arrays */
717 CLONE->umin = _unur_xmalloc( GEN->dim * sizeof(double));
718 CLONE->umax = _unur_xmalloc( GEN->dim * sizeof(double));
719
720 /* copy parameters into clone object */
721 memcpy(CLONE->umin, GEN->umin, GEN->dim * sizeof(double));
722 memcpy(CLONE->umax, GEN->umax, GEN->dim * sizeof(double));
723
724 /* copy data */
725 CLONE->center = unur_distr_cvec_get_center(clone->distr);
726
727 return clone;
728
729 #undef CLONE
730 } /* end of _unur_vnrou_clone() */
731
732 /*****************************************************************************/
733
734 int
_unur_vnrou_sample_cvec(struct unur_gen * gen,double * vec)735 _unur_vnrou_sample_cvec( struct unur_gen *gen, double *vec )
736 /*----------------------------------------------------------------------*/
737 /* sample from generator */
738 /* */
739 /* parameters: */
740 /* gen ... pointer to generator object */
741 /* vec ... random vector (result) */
742 /*----------------------------------------------------------------------*/
743 {
744 double U, V;
745 int d, dim; /* index used in dimension loops (0 <= d < dim) */
746
747 /* check arguments */
748 CHECK_NULL(gen,UNUR_ERR_NULL);
749 COOKIE_CHECK(gen,CK_VNROU_GEN,UNUR_ERR_COOKIE);
750
751 dim = GEN->dim;
752
753 while (1) {
754
755 /* generate point uniformly on rectangle */
756 while ( _unur_iszero(V = _unur_call_urng(gen->urng)) );
757 V *= GEN->vmax;
758 for (d=0; d<dim; d++) {
759 U = GEN->umin[d] + _unur_call_urng(gen->urng) * (GEN->umax[d] - GEN->umin[d]);
760 vec[d] = U/pow(V,GEN->r) + GEN->center[d];
761 }
762
763 /* X[] inside domain ? */
764
765 /* accept or reject */
766 if (V <= pow(PDF(vec),1./(GEN->r * dim + 1.)))
767 return UNUR_SUCCESS;
768 }
769
770 } /* end of _unur_vnrou_sample() */
771
772 /*---------------------------------------------------------------------------*/
773
774 int
_unur_vnrou_sample_check(struct unur_gen * gen,double * vec)775 _unur_vnrou_sample_check( struct unur_gen *gen, double *vec )
776 /*----------------------------------------------------------------------*/
777 /* sample from generator and verify that method can be used */
778 /* */
779 /* parameters: */
780 /* gen ... pointer to generator object */
781 /* vec ... random sample vector (return) */
782 /*----------------------------------------------------------------------*/
783 {
784 double U, V;
785 int d, dim; /* index used in dimension loops (0 <= d < dim) */
786 int hat_error;
787
788 double fx,sfx,xfx;
789
790 /* check arguments */
791 CHECK_NULL(gen,UNUR_ERR_NULL);
792 COOKIE_CHECK(gen,CK_VNROU_GEN,UNUR_ERR_COOKIE);
793
794 dim = GEN->dim;
795
796 while (1) {
797 /* generate point uniformly on rectangle */
798 while ( _unur_iszero(V = _unur_call_urng(gen->urng)) );
799 V *= GEN->vmax;
800 for (d=0; d<dim; d++) {
801 U = GEN->umin[d] + _unur_call_urng(gen->urng) * (GEN->umax[d] - GEN->umin[d]);
802 vec[d] = U/pow(V,GEN->r) + GEN->center[d];
803 }
804
805 /* X[] inside domain ? */
806
807 /* evaluate PDF */
808 fx = PDF(vec);
809
810 /* a point on the boundary of the region of acceptance
811 has the coordinates ( (vec[]-center[]) * (fx)^(r/r*dim+1)), fx^(1/r*dim+1) ). */
812 sfx = pow( fx, 1./(GEN->r * dim+1.) );
813 /* check hat */
814 hat_error=0;
815 if ( sfx > (1.+DBL_EPSILON) * GEN->vmax ) hat_error++;
816
817 sfx = pow( fx, GEN->r/(GEN->r * dim + 1.) );
818 for (d=0; d<dim; d++) {
819 xfx = (vec[d]-GEN->center[d]) * sfx;
820 if ( (xfx < (1.+UNUR_EPSILON) * GEN->umin[d])
821 || (xfx > (1.+UNUR_EPSILON) * GEN->umax[d]))
822 hat_error++;
823 }
824
825 if (hat_error>0) _unur_error(gen->genid,UNUR_ERR_GEN_CONDITION,"PDF(x) > hat(x)");
826
827 /* accept or reject */
828 if (V <= pow(PDF(vec),1./( GEN->r * dim + 1.)))
829 return UNUR_SUCCESS;
830 }
831
832 } /* end of _unur_vnrou_sample_check() */
833
834 /*****************************************************************************/
835
836 void
_unur_vnrou_free(struct unur_gen * gen)837 _unur_vnrou_free( struct unur_gen *gen )
838 /*----------------------------------------------------------------------*/
839 /* deallocate generator object */
840 /* */
841 /* parameters: */
842 /* gen ... pointer to generator object */
843 /*----------------------------------------------------------------------*/
844 {
845 /* check arguments */
846 if( !gen ) /* nothing to do */
847 return;
848
849 /* check input */
850 if ( gen->method != UNUR_METH_VNROU ) {
851 _unur_warning(gen->genid,UNUR_ERR_GEN_INVALID,"");
852 return; }
853 COOKIE_CHECK(gen,CK_VNROU_GEN,RETURN_VOID);
854
855 /* we cannot use this generator object any more */
856 SAMPLE = NULL; /* make sure to show up a programming error */
857
858 /* free memory */
859 if (GEN->umin) free(GEN->umin);
860 if (GEN->umax) free(GEN->umax);
861 _unur_generic_free(gen);
862
863 } /* end of _unur_vnrou_free() */
864
865 /*****************************************************************************/
866 /** Auxilliary Routines **/
867 /*****************************************************************************/
868
869 int
_unur_vnrou_rectangle(struct unur_gen * gen)870 _unur_vnrou_rectangle( struct unur_gen *gen )
871 /*----------------------------------------------------------------------*/
872 /* compute universal bounding rectangle */
873 /* */
874 /* parameters: */
875 /* gen ... pointer to generator object */
876 /* */
877 /* return: */
878 /* UNUR_SUCCESS ... on success */
879 /* error code ... on error */
880 /*----------------------------------------------------------------------*/
881 {
882
883 int d; /* index used in dimension loops (0 <= d < dim) */
884 struct MROU_RECTANGLE *rr;
885 int rectangle_compute;
886
887 /* check arguments */
888 CHECK_NULL( gen, UNUR_ERR_NULL );
889 COOKIE_CHECK( gen,CK_VNROU_GEN, UNUR_ERR_COOKIE );
890
891 /* Boundary rectangle is already set */
892 if ((gen->set & VNROU_SET_U) && (gen->set & VNROU_SET_V)) {
893 return UNUR_SUCCESS;
894 }
895
896 /* Allocating and filling mrou_rectangle struct */
897 rr = _unur_mrou_rectangle_new();
898
899 rr->distr = gen->distr;
900 rr->dim = GEN->dim;
901 rr->umin = GEN->umin;
902 rr->umax = GEN->umax;
903 rr->r = GEN->r;
904 rr->center = GEN->center;
905 rr->genid = gen->genid;
906
907 /* calculate bounding rectangle */
908 rectangle_compute = _unur_mrou_rectangle_compute(rr);
909
910 if (!(gen->set & VNROU_SET_V)) {
911 /* user has not provided any upper bound for v */
912 GEN->vmax = rr->vmax;
913 }
914
915 if (!(gen->set & VNROU_SET_U)) {
916 /* user has not provided any bounds for u */
917 for (d=0; d<GEN->dim; d++) {
918 GEN->umin[d] = rr->umin[d];
919 GEN->umax[d] = rr->umax[d];
920 }
921 }
922
923 free(rr);
924
925 if (rectangle_compute != UNUR_SUCCESS)
926 return UNUR_ERR_INF;
927
928 /* o.k. */
929 return UNUR_SUCCESS;
930
931 } /* end of _unur_vnrou_rectangle() */
932
933 /*****************************************************************************/
934 /** Debugging utilities **/
935 /*****************************************************************************/
936
937 /*---------------------------------------------------------------------------*/
938 #ifdef UNUR_ENABLE_LOGGING
939 /*---------------------------------------------------------------------------*/
940
941 void
_unur_vnrou_debug_init(const struct unur_gen * gen)942 _unur_vnrou_debug_init( const struct unur_gen *gen )
943 /*----------------------------------------------------------------------*/
944 /* write info about generator into LOG file */
945 /* */
946 /* parameters: */
947 /* gen ... pointer to generator object */
948 /*----------------------------------------------------------------------*/
949 {
950 FILE *LOG;
951 int d, dim; /* index used in dimension loops (0 <= d < dim) */
952 double vol;
953
954 /* check arguments */
955 CHECK_NULL(gen,RETURN_VOID); COOKIE_CHECK(gen,CK_VNROU_GEN,RETURN_VOID);
956
957 LOG = unur_get_stream();
958 dim = GEN->dim;
959
960 fprintf(LOG,"%s:\n",gen->genid);
961 fprintf(LOG,"%s: type = continuous multivariate random variates\n",gen->genid);
962 fprintf(LOG,"%s: method = vnrou (naive ratio-of-uniforms)\n",gen->genid);
963 fprintf(LOG,"%s:\n",gen->genid);
964
965 _unur_distr_cvec_debug( gen->distr, gen->genid );
966
967 fprintf(LOG,"%s: sampling routine = _unur_vnrou_sample",gen->genid);
968 if (gen->variant & VNROU_VARFLAG_VERIFY) fprintf(LOG,"_check");
969 fprintf(LOG,"()\n%s:\n",gen->genid);
970
971 /* parameters */
972 fprintf(LOG,"%s: r-parameter = %g",gen->genid, GEN->r);
973 _unur_print_if_default(gen,VNROU_SET_R);
974 fprintf(LOG,"\n%s:\n",gen->genid);
975
976 /* print center */
977 _unur_matrix_print_vector( GEN->dim, GEN->center, "center =", LOG, gen->genid, "\t ");
978
979 /* print bounding rectangle */
980 fprintf(LOG,"%s: Rectangle:",gen->genid);
981 if (!((gen->set & VNROU_SET_U) && (gen->set & VNROU_SET_V)))
982 fprintf(LOG,"\t[computed]");
983 else
984 fprintf(LOG,"\t[input]");
985 fprintf(LOG,"\n");
986
987 vol = GEN->vmax;
988 fprintf(LOG,"%s:\tvmax = %g\n",gen->genid, GEN->vmax);
989 for (d=0; d<dim; d++) {
990 vol *= (GEN->umax[d]-GEN->umin[d]);
991 fprintf(LOG,"%s:\tumin[%d],umax[%d] = (%g,%g)\n",gen->genid,
992 d, d, GEN->umin[d], GEN->umax[d]);
993 }
994 fprintf(LOG,"%s:\n",gen->genid);
995 fprintf(LOG,"%s:\tvolume = %g\t(hat = %g)\n",gen->genid, vol, vol*(GEN->r*GEN->dim+1));
996 fprintf(LOG,"%s:\n",gen->genid);
997
998 } /* end of _unur_vnrou_debug_init() */
999
1000 /*---------------------------------------------------------------------------*/
1001 #endif /* end UNUR_ENABLE_LOGGING */
1002 /*---------------------------------------------------------------------------*/
1003
1004
1005 /*---------------------------------------------------------------------------*/
1006 #ifdef UNUR_ENABLE_INFO
1007 /*---------------------------------------------------------------------------*/
1008
1009 void
_unur_vnrou_info(struct unur_gen * gen,int help)1010 _unur_vnrou_info( struct unur_gen *gen, int help )
1011 /*----------------------------------------------------------------------*/
1012 /* create character string that contains information about the */
1013 /* given generator object. */
1014 /* */
1015 /* parameters: */
1016 /* gen ... pointer to generator object */
1017 /* help ... whether to print additional comments */
1018 /*----------------------------------------------------------------------*/
1019 {
1020 struct unur_string *info = gen->infostr;
1021 struct unur_distr *distr = gen->distr;
1022 int samplesize = 10000;
1023 int i;
1024 double hvol;
1025
1026 /* generator ID */
1027 _unur_string_append(info,"generator ID: %s\n\n", gen->genid);
1028
1029 /* distribution */
1030 _unur_string_append(info,"distribution:\n");
1031 _unur_distr_info_typename(gen);
1032 _unur_string_append(info," dimension = %d\n",GEN->dim);
1033 _unur_string_append(info," functions = PDF\n");
1034 _unur_distr_cvec_info_domain(gen);
1035
1036 if ( distr->set & UNUR_DISTR_SET_MODE ) {
1037 _unur_string_append(info," mode = ");
1038 _unur_distr_info_vector( gen, DISTR.mode, GEN->dim);
1039 }
1040 _unur_string_append(info,"\n");
1041
1042 _unur_string_append(info," center = ");
1043 _unur_distr_info_vector( gen, GEN->center, GEN->dim);
1044 if ( !(distr->set & UNUR_DISTR_SET_CENTER) ) {
1045 if ( distr->set & UNUR_DISTR_SET_MODE )
1046 _unur_string_append(info," [= mode]");
1047 else
1048 _unur_string_append(info," [default]");
1049 }
1050 _unur_string_append(info,"\n\n");
1051
1052 /* if (help) { */
1053 /* _unur_string_append(info,"\n"); */
1054 /* } */
1055
1056 /* method */
1057 _unur_string_append(info,"method: VNROU (Naive Ratio-Of-Uniforms)\n");
1058 _unur_string_append(info," r = %g\n", GEN->r);
1059 _unur_string_append(info,"\n");
1060
1061 /* performance */
1062 _unur_string_append(info,"performance characteristics:\n");
1063
1064 _unur_string_append(info," bounding rectangle = ");
1065 for (i=0; i<GEN->dim; i++)
1066 _unur_string_append(info,"%s(%g,%g)", i?"x":"", GEN->umin[i], GEN->umax[i]);
1067 _unur_string_append(info," x (0,%g)\n", GEN->vmax);
1068
1069 hvol = GEN->vmax;
1070 for (i=0; i<GEN->dim; i++)
1071 hvol *= GEN->umax[i] - GEN->umin[i];
1072 _unur_string_append(info," volume(hat) = %g\n", hvol);
1073
1074 _unur_string_append(info," rejection constant ");
1075 if ((distr->set & UNUR_DISTR_SET_PDFVOLUME) && _unur_isone(GEN->r))
1076 _unur_string_append(info,"= %g\n", (GEN->dim + 1.) * hvol / DISTR.volume);
1077 else
1078 _unur_string_append(info,"= %.2f [approx.]\n",
1079 unur_test_count_urn(gen,samplesize,0,NULL)/((1.+GEN->dim)*samplesize));
1080 _unur_string_append(info,"\n");
1081
1082 /* parameters */
1083 if (help) {
1084 _unur_string_append(info,"parameters:\n");
1085
1086 _unur_string_append(info," r = %g %s\n", GEN->r,
1087 (gen->set & VNROU_SET_R) ? "" : "[default]");
1088
1089 _unur_string_append(info," v = %g %s\n", GEN->vmax,
1090 (gen->set & VNROU_SET_V) ? "" : "[numeric.]");
1091
1092 _unur_string_append(info," u = ");
1093 _unur_distr_info_vector( gen, GEN->umin, GEN->dim);
1094 _unur_string_append(info," -- ");
1095 _unur_distr_info_vector( gen, GEN->umax, GEN->dim);
1096 _unur_string_append(info,"%s\n",(gen->set & VNROU_SET_U) ? "" : " [numeric.]");
1097
1098 if (gen->variant & VNROU_VARFLAG_VERIFY)
1099 _unur_string_append(info," verify = on\n");
1100
1101 _unur_string_append(info,"\n");
1102 }
1103
1104 /* Hints */
1105 if (help) {
1106 if ( !(gen->set & VNROU_SET_V) )
1107 _unur_string_append(info,"[ Hint: %s ]\n",
1108 "You can set \"v\" to avoid numerical estimate." );
1109 if ( !(gen->set & VNROU_SET_U) )
1110 _unur_string_append(info,"[ Hint: %s ]\n",
1111 "You can set \"u\" to avoid slow (and inexact) numerical estimates." );
1112 _unur_string_append(info,"\n");
1113 }
1114
1115 } /* end of _unur_vnrou_info() */
1116
1117 /*---------------------------------------------------------------------------*/
1118 #endif /* end UNUR_ENABLE_INFO */
1119 /*---------------------------------------------------------------------------*/
1120