1 /*****************************************************************************
2 * *
3 * UNURAN -- Universal Non-Uniform Random number generator *
4 * *
5 *****************************************************************************
6 * *
7 * FILE: mixt.c *
8 * *
9 * TYPE: (continuous) univariate random variate *
10 * METHOD: mixture of distributions (meta method) *
11 * *
12 * DESCRIPTION: *
13 * Given an array of components and array of probabilities *
14 * produce a value X consistent with the mixture distribution *
15 * *
16 * REQUIRED: *
17 * pointers to generators of components *
18 * array of probabilities *
19 * *
20 *****************************************************************************
21 * *
22 * Copyright (c) 2010 Wolfgang Hoermann and Josef Leydold *
23 * Department of Statistics and Mathematics, WU Wien, Austria *
24 * *
25 * This program is free software; you can redistribute it and/or modify *
26 * it under the terms of the GNU General Public License as published by *
27 * the Free Software Foundation; either version 2 of the License, or *
28 * (at your option) any later version. *
29 * *
30 * This program is distributed in the hope that it will be useful, *
31 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
32 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
33 * GNU General Public License for more details. *
34 * *
35 * You should have received a copy of the GNU General Public License *
36 * along with this program; if not, write to the *
37 * Free Software Foundation, Inc., *
38 * 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA *
39 * *
40 *****************************************************************************/
41
42 /*---------------------------------------------------------------------------*/
43
44 #include <unur_source.h>
45 #include <distr/distr.h>
46 #include <distr/distr_source.h>
47 #include <distr/cont.h>
48 #include <distr/discr.h>
49 #include <urng/urng.h>
50 #include <utils/unur_fp_source.h>
51 #include "unur_methods_source.h"
52 #include "x_gen.h"
53 #include "x_gen_source.h"
54 #include "dgt.h"
55 #include "dgt_struct.h"
56 #include "mixt.h"
57 #include "mixt_struct.h"
58
59 /*---------------------------------------------------------------------------*/
60 /* Variants: */
61
62 #define MIXT_VARFLAG_INVERSION 0x004u /* use inversion method (if possible) */
63
64 /*---------------------------------------------------------------------------*/
65 /* Debugging flags */
66 /* bit 01 ... pameters and structure of generator (do not use here) */
67 /* bits 02-12 ... setup */
68 /* bits 13-24 ... adaptive steps */
69 /* bits 25-32 ... trace sampling */
70
71 /*---------------------------------------------------------------------------*/
72 /* Flags for logging set calls */
73
74 #define MIXT_SET_USEINVERSION 0x001u /* use inverion method */
75
76 /*---------------------------------------------------------------------------*/
77
78 #define GENTYPE "MIXT" /* type of generator */
79
80 /*---------------------------------------------------------------------------*/
81
82 static struct unur_gen *_unur_mixt_init( struct unur_par *par );
83 /*---------------------------------------------------------------------------*/
84 /* Initialize new generator. */
85 /*---------------------------------------------------------------------------*/
86
87 static struct unur_gen *_unur_mixt_create( struct unur_par *par );
88 /*---------------------------------------------------------------------------*/
89 /* create new (almost empty) generator object. */
90 /*---------------------------------------------------------------------------*/
91
92 static int _unur_mixt_check_par( struct unur_gen *gen );
93 /*---------------------------------------------------------------------------*/
94 /* Check parameters of given distribution and method */
95 /*---------------------------------------------------------------------------*/
96
97 static struct unur_gen *_unur_mixt_clone( const struct unur_gen *gen );
98 /*---------------------------------------------------------------------------*/
99 /* copy (clone) generator object. */
100 /*---------------------------------------------------------------------------*/
101
102 static void _unur_mixt_free( struct unur_gen *gen);
103 /*---------------------------------------------------------------------------*/
104 /* destroy generator object. */
105 /*---------------------------------------------------------------------------*/
106
107 static double _unur_mixt_sample( struct unur_gen *gen );
108 /*---------------------------------------------------------------------------*/
109 /* sample from generator */
110 /*---------------------------------------------------------------------------*/
111
112 static double _unur_mixt_sample_inv( struct unur_gen *gen );
113 /*---------------------------------------------------------------------------*/
114 /* sample from generator by inversion */
115 /*---------------------------------------------------------------------------*/
116
117 static struct unur_gen *_unur_mixt_indexgen( const double *prob, int n_prob );
118 /*---------------------------------------------------------------------------*/
119 /* create generator for index. */
120 /*---------------------------------------------------------------------------*/
121
122 static int _unur_mixt_get_boundary( struct unur_gen *gen );
123 /*---------------------------------------------------------------------------*/
124 /* compute boundary of mixture and check for overlapping domains. */
125 /*---------------------------------------------------------------------------*/
126
127 #ifdef UNUR_ENABLE_LOGGING
128 /*---------------------------------------------------------------------------*/
129 /* the following functions print debugging information on output stream, */
130 /* i.e., into the LOG file if not specified otherwise. */
131 /*---------------------------------------------------------------------------*/
132
133 static void _unur_mixt_debug_init( const struct unur_gen *gen );
134 /*---------------------------------------------------------------------------*/
135 /* print after generator has been initialized has completed. */
136 /*---------------------------------------------------------------------------*/
137 #endif
138
139 #ifdef UNUR_ENABLE_INFO
140 static void _unur_mixt_info( struct unur_gen *gen, int help );
141 /*---------------------------------------------------------------------------*/
142 /* create info string. */
143 /*---------------------------------------------------------------------------*/
144 #endif
145
146 /*---------------------------------------------------------------------------*/
147 /* abbreviations */
148
149 #define PAR ((struct unur_mixt_par*)par->datap) /* data for parameter object */
150 #define GEN ((struct unur_mixt_gen*)gen->datap) /* data for generator object */
151 #define DISTR gen->distr->data.cont /* data for distribution in generator object */
152
153 #define SAMPLE gen->sample.cont /* pointer to sampling routine */
154
155 #define BD_LEFT domain[0] /* left boundary of domain of distribution */
156 #define BD_RIGHT domain[1] /* right boundary of domain of distribution */
157
158 /* shortcuts for auxiliary generators */
159 #define INDEX gen_aux
160
161 #define PROB gen_aux->distr->data.discr.pv
162 #define COMP gen_aux_list
163 #define N_COMP n_gen_aux_list
164
165 /*---------------------------------------------------------------------------*/
166
167 #define _unur_mixt_getSAMPLE(gen) \
168 ( ((gen)->variant & MIXT_VARFLAG_INVERSION) \
169 ? _unur_mixt_sample_inv : _unur_mixt_sample )
170
171 /*---------------------------------------------------------------------------*/
172
173 /*****************************************************************************/
174 /** Public: User Interface (API) **/
175 /*****************************************************************************/
176
177 struct unur_par *
unur_mixt_new(int n,const double * prob,struct unur_gen ** comp)178 unur_mixt_new( int n, const double *prob, struct unur_gen **comp )
179 /*----------------------------------------------------------------------*/
180 /* get default parameters */
181 /* */
182 /* parameters: */
183 /* n ... number of components */
184 /* prob ... probabilities for components */
185 /* comp ... array of pointers to components (generators) */
186 /* */
187 /* return: */
188 /* default parameters (pointer to structure) */
189 /* */
190 /* error: */
191 /* return NULL */
192 /*----------------------------------------------------------------------*/
193 {
194 struct unur_par *par;
195
196 /* check arguments */
197 _unur_check_NULL( GENTYPE, prob, NULL );
198 _unur_check_NULL( GENTYPE, comp, NULL );
199 if (n<1) { _unur_error(GENTYPE,UNUR_ERR_DISTR_DOMAIN,"n < 1"); return NULL; }
200 /* checking type of generator objects in 'comp' is delayed to init */
201
202 /* allocate structure */
203 par = _unur_par_new( sizeof(struct unur_mixt_par) );
204 COOKIE_SET(par,CK_MIXT_PAR);
205
206 /* copy input */
207 par->distr = NULL; /* pointer to distribution object */
208
209 /* copy data */
210 PAR->n_comp = n; /* number of components */
211 PAR->prob = prob; /* probabilities for components */
212 PAR->comp = comp; /* array of pointers to components (generators) */
213
214 /* set default values */
215 par->method = UNUR_METH_MIXT; /* method and default variant */
216 par->variant = 0u; /* default variant */
217 par->set = 0u; /* inidicate default parameters */
218 par->urng = unur_get_default_urng(); /* use default urng */
219 par->urng_aux = NULL; /* no auxilliary URNG required */
220
221 par->debug = _unur_default_debugflag; /* set default debugging flags */
222
223 /* routine for starting generator */
224 par->init = _unur_mixt_init;
225
226 return par;
227
228 } /* end of unur_mixt_new() */
229
230 /*****************************************************************************/
231
232 int
unur_mixt_set_useinversion(struct unur_par * par,int useinversion)233 unur_mixt_set_useinversion( struct unur_par *par, int useinversion )
234 /*----------------------------------------------------------------------*/
235 /* set flag for using inversion method (default: off) */
236 /* */
237 /* parameters: */
238 /* par ... pointer to parameter for building generator object */
239 /* useinversion ... !0 = try inversion method, 0 = no inversion */
240 /* */
241 /* return: */
242 /* UNUR_SUCCESS ... on success */
243 /* error code ... on error */
244 /* */
245 /* comment: */
246 /* no squeeze is the default */
247 /*----------------------------------------------------------------------*/
248 {
249 /* check arguments */
250 _unur_check_NULL( GENTYPE, par, UNUR_ERR_NULL );
251 _unur_check_par_object( par, MIXT );
252
253 /* we use a bit in 'variant' */
254 par->variant = (useinversion)
255 ? (par->variant | MIXT_VARFLAG_INVERSION)
256 : (par->variant & (~MIXT_VARFLAG_INVERSION));
257
258 /* changelog */
259 par->set |= MIXT_SET_USEINVERSION;
260
261 /* o.k. */
262 return UNUR_SUCCESS;
263
264 } /* end of unur_mixt_set_useinversion() */
265
266
267 /*****************************************************************************/
268 /** Private **/
269 /*****************************************************************************/
270
271 struct unur_gen *
_unur_mixt_init(struct unur_par * par)272 _unur_mixt_init( struct unur_par *par )
273 /*----------------------------------------------------------------------*/
274 /* initialize new generator */
275 /* */
276 /* parameters: */
277 /* params pointer to paramters for building generator object */
278 /* */
279 /* return: */
280 /* pointer to generator object */
281 /* */
282 /* error: */
283 /* return NULL */
284 /*----------------------------------------------------------------------*/
285 {
286 struct unur_gen *gen;
287 int i;
288
289 /* check arguments */
290 CHECK_NULL(par,NULL);
291
292 /* check input */
293 if ( par->method != UNUR_METH_MIXT ) {
294 _unur_error(GENTYPE,UNUR_ERR_PAR_INVALID,"");
295 return NULL; }
296 COOKIE_CHECK(par,CK_MIXT_PAR,NULL);
297
298 /* create a new empty generator object */
299 gen = _unur_mixt_create(par);
300 if (!gen) { _unur_par_free(par); return NULL; }
301
302 /* probabilities */
303 gen->INDEX = _unur_mixt_indexgen(PAR->prob,PAR->n_comp);
304
305 /* components */
306 gen->N_COMP = PAR->n_comp; /* number of components */
307 gen->COMP = _unur_xmalloc( gen->N_COMP * sizeof(struct unur_gen *));
308 for (i=0; i<gen->N_COMP; i++)
309 gen->COMP[i] = unur_gen_clone(PAR->comp[i]);
310
311 /* free parameters */
312 _unur_par_free(par);
313
314 /* check parameters */
315 if (_unur_mixt_check_par(gen) != UNUR_SUCCESS) {
316 _unur_mixt_free(gen); return NULL;
317 }
318
319 /* compute boundary of mixture and check for overlapping domains */
320 if ( _unur_mixt_get_boundary(gen) != UNUR_SUCCESS ) {
321 _unur_mixt_free(gen); return NULL;
322 }
323
324 /* set name of distribution */
325 unur_distr_set_name(gen->distr, "(mixture)");
326
327 #ifdef UNUR_ENABLE_LOGGING
328 /* write info into LOG file */
329 if (gen->debug) _unur_mixt_debug_init(gen);
330 #endif
331
332 return gen;
333 } /* end of _unur_mixt_init() */
334
335 /*---------------------------------------------------------------------------*/
336
337 struct unur_gen *
_unur_mixt_create(struct unur_par * par)338 _unur_mixt_create( struct unur_par *par )
339 /*----------------------------------------------------------------------*/
340 /* allocate memory for generator */
341 /* */
342 /* parameters: */
343 /* par ... pointer to parameter for building generator object */
344 /* */
345 /* return: */
346 /* pointer to (empty) generator object with default settings */
347 /* */
348 /* error: */
349 /* return NULL */
350 /*----------------------------------------------------------------------*/
351 {
352 struct unur_gen *gen;
353
354 /* check arguments */
355 CHECK_NULL(par,NULL); COOKIE_CHECK(par,CK_MIXT_PAR,NULL);
356
357 /* create new generic generator object */
358 gen = _unur_generic_create( par, sizeof(struct unur_mixt_gen) );
359
360 /* magic cookies */
361 COOKIE_SET(gen,CK_MIXT_GEN);
362
363 /* set generator identifier */
364 gen->genid = _unur_set_genid(GENTYPE);
365
366 /* object 'par' does not contain a distribution object */
367 /* so we create one. */
368 gen->distr = unur_distr_cont_new();
369
370 /* routines for sampling and destroying generator */
371 SAMPLE = _unur_mixt_getSAMPLE(gen);
372 gen->destroy = _unur_mixt_free;
373 gen->clone = _unur_mixt_clone;
374 gen->reinit = NULL; /* reinit not implemented ! */
375
376 /* copy some parameters into generator object */
377 GEN->is_inversion = (gen->variant & MIXT_VARFLAG_INVERSION) ? TRUE : FALSE;
378
379 /* initialize parameters: none */
380
381 #ifdef UNUR_ENABLE_INFO
382 /* set function for creating info string */
383 gen->info = _unur_mixt_info;
384 #endif
385
386 /* return pointer to (almost empty) generator object */
387 return gen;
388
389 } /* end of _unur_mixt_create() */
390
391 /*---------------------------------------------------------------------------*/
392
393 int
_unur_mixt_check_par(struct unur_gen * gen)394 _unur_mixt_check_par( struct unur_gen *gen )
395 /*----------------------------------------------------------------------*/
396 /* check parameters of given distribution and method */
397 /* */
398 /* parameters: */
399 /* gen ... pointer to generator object */
400 /* */
401 /* return: */
402 /* UNUR_SUCCESS ... on success */
403 /* error code ... on error */
404 /*----------------------------------------------------------------------*/
405 {
406 int i;
407 int type;
408
409 /* check probabilities */
410 if (gen->INDEX == NULL) {
411 _unur_error(gen->genid,UNUR_ERR_GEN_DATA,"invalid probabilities");
412 return UNUR_ERR_GEN_DATA;
413 }
414
415 /* check generator objects */
416 for (i=0; i<gen->N_COMP; i++) {
417
418 /* generators must not be NULL */
419 if (gen->COMP[i] == NULL) {
420 _unur_error(gen->genid,UNUR_ERR_NULL,"component is NULL");
421 return UNUR_ERR_NULL;
422 }
423
424 /* all generators must sample from univariate distributions */
425 type = gen->COMP[i]->method & UNUR_MASK_TYPE;
426 if ( type != UNUR_METH_DISCR &&
427 type != UNUR_METH_CONT &&
428 type != UNUR_METH_CEMP ) {
429 _unur_error(gen->genid,UNUR_ERR_GEN_INVALID,"component not univariate");
430 return UNUR_ERR_GEN_INVALID;
431 }
432
433 /* we only can use inversion method if all generators use inversion method */
434 if (GEN->is_inversion && (! _unur_gen_is_inversion (gen->COMP[i]))) {
435 _unur_error(gen->genid,UNUR_ERR_GEN_INVALID,"component does not implement inversion");
436 return UNUR_ERR_GEN_INVALID;
437 }
438 }
439
440 return UNUR_SUCCESS;
441 } /* end of _unur_mixt_check_par() */
442
443 /*---------------------------------------------------------------------------*/
444
445 struct unur_gen *
_unur_mixt_clone(const struct unur_gen * gen)446 _unur_mixt_clone( const struct unur_gen *gen )
447 /*----------------------------------------------------------------------*/
448 /* copy (clone) generator object */
449 /* */
450 /* parameters: */
451 /* gen ... pointer to generator object */
452 /* */
453 /* return: */
454 /* pointer to clone of generator object */
455 /* */
456 /* error: */
457 /* return NULL */
458 /*----------------------------------------------------------------------*/
459 {
460 #define CLONE ((struct unur_mixt_gen*)clone->datap)
461
462 struct unur_gen *clone;
463
464 /* check arguments */
465 CHECK_NULL(gen,NULL); COOKIE_CHECK(gen,CK_MIXT_GEN,NULL);
466
467 /* create generic clone */
468 clone = _unur_generic_clone( gen, GENTYPE );
469
470 return clone;
471
472 #undef CLONE
473 } /* end of _unur_mixt_clone() */
474
475 /*---------------------------------------------------------------------------*/
476
477 void
_unur_mixt_free(struct unur_gen * gen)478 _unur_mixt_free( struct unur_gen *gen )
479 /*----------------------------------------------------------------------*/
480 /* deallocate generator object */
481 /* */
482 /* parameters: */
483 /* gen ... pointer to generator object */
484 /*----------------------------------------------------------------------*/
485 {
486 /* check arguments */
487 if( !gen ) /* nothing to do */
488 return;
489
490 /* check input */
491 if ( gen->method != UNUR_METH_MIXT ) {
492 _unur_warning(gen->genid,UNUR_ERR_GEN_INVALID,"");
493 return; }
494 COOKIE_CHECK(gen,CK_MIXT_GEN,RETURN_VOID);
495
496 /* we cannot use this generator object any more */
497 SAMPLE = NULL; /* make sure to show up a programming error */
498
499 /* free memory */
500 _unur_generic_free(gen);
501
502 } /* end of _unur_mixt_free() */
503
504 /*****************************************************************************/
505
506 double
_unur_mixt_sample(struct unur_gen * gen)507 _unur_mixt_sample( struct unur_gen *gen )
508 /*----------------------------------------------------------------------*/
509 /* sample from generator */
510 /* */
511 /* parameters: */
512 /* gen ... pointer to generator object */
513 /* */
514 /* return: */
515 /* double (sample from random variate) */
516 /* */
517 /* error: */
518 /* return INFINITY */
519 /*----------------------------------------------------------------------*/
520 {
521 struct unur_gen *comp;
522 int J;
523
524 /* check arguments */
525 CHECK_NULL(gen,INFINITY); COOKIE_CHECK(gen,CK_MIXT_GEN,INFINITY);
526
527 /* sample index */
528 J = unur_sample_discr(gen->INDEX);
529
530 /* get component */
531 comp = gen->COMP[J];
532
533 /* sample from selected component */
534 switch(comp->method & UNUR_MASK_TYPE) {
535 case UNUR_METH_DISCR:
536 return ((double) comp->sample.discr(comp));
537 case UNUR_METH_CONT:
538 case UNUR_METH_CEMP:
539 default:
540 return (comp->sample.cont(comp));
541 }
542
543 } /* end of _unur_mixt_sample() */
544
545 /*---------------------------------------------------------------------------*/
546
547 double
_unur_mixt_sample_inv(struct unur_gen * gen)548 _unur_mixt_sample_inv( struct unur_gen *gen )
549 /*----------------------------------------------------------------------*/
550 /* sample from generator by inversion */
551 /* */
552 /* parameters: */
553 /* gen ... pointer to generator object */
554 /* */
555 /* return: */
556 /* double (sample from random variate) */
557 /* */
558 /* error: */
559 /* return INFINITY */
560 /*----------------------------------------------------------------------*/
561 {
562 double U, recycle;
563 int J;
564
565 /* check arguments */
566 CHECK_NULL(gen,INFINITY); COOKIE_CHECK(gen,CK_MIXT_GEN,INFINITY);
567
568 /* sample index */
569 U = _unur_call_urng(gen->urng);
570 J =unur_dgt_eval_invcdf_recycle( gen->INDEX, U, &recycle );
571
572 /* the resolution of recycle is less than that of U. */
573 /* the values 0. and 1. may be result in INFINITY. */
574 /* thus we make a small perturbation in this case. */
575 if (_unur_iszero(recycle)) recycle = DBL_MIN;
576 if (_unur_isone(recycle)) recycle = 1. - DBL_EPSILON;
577
578 /* sample from component */
579 return unur_quantile(gen->COMP[J], recycle);
580
581 } /* end of _unur_mixt_sample_inv() */
582
583 /*---------------------------------------------------------------------------*/
584
585 double
unur_mixt_eval_invcdf(const struct unur_gen * gen,double u)586 unur_mixt_eval_invcdf( const struct unur_gen *gen, double u )
587 /*----------------------------------------------------------------------*/
588 /* evaluate inverse CDF at u. */
589 /* */
590 /* parameters: */
591 /* gen ... pointer to generator object */
592 /* u ... argument for inverse CDF */
593 /* */
594 /* return: */
595 /* double (inverse CDF) */
596 /* */
597 /* error: */
598 /* return INFINITY */
599 /*----------------------------------------------------------------------*/
600 {
601 double recycle;
602 int J;
603
604 /* check arguments */
605 _unur_check_NULL( GENTYPE, gen, INFINITY );
606 if ( ! (gen->method == UNUR_METH_MIXT && GEN->is_inversion) ) {
607 _unur_error(gen->genid,UNUR_ERR_GEN_INVALID,"");
608 return INFINITY;
609 }
610 COOKIE_CHECK(gen,CK_MIXT_GEN,INFINITY);
611
612 if ( ! (u>0. && u<1.)) {
613 if ( ! (u>=0. && u<=1.)) {
614 _unur_warning(gen->genid,UNUR_ERR_DOMAIN,"U not in [0,1]");
615 }
616 if (u<=0.) return DISTR.domain[0];
617 if (u>=1.) return DISTR.domain[1];
618 return u; /* = NaN */
619 }
620
621 /* get index */
622 J =unur_dgt_eval_invcdf_recycle( gen->INDEX, u, &recycle );
623
624 /* the resolution of recycle is less than that of U. */
625 /* the values 0. and 1. may be result in INFINITY. */
626 /* thus we make a small perturbation in this case. */
627 if (_unur_iszero(recycle)) recycle = DBL_MIN;
628 if (_unur_isone(recycle)) recycle = 1. - DBL_EPSILON;
629
630 /* get from component */
631 return unur_quantile(gen->COMP[J], recycle);
632
633 } /* end of unur_mixt_eval_invcdf() */
634
635
636 /*****************************************************************************/
637 /** Auxilliary Routines **/
638 /*****************************************************************************/
639
640 struct unur_gen *
_unur_mixt_indexgen(const double * prob,int n_prob)641 _unur_mixt_indexgen( const double *prob, int n_prob )
642 /*---------------------------------------------------------------------------*/
643 /* create generator for index using method DGT. */
644 /* */
645 /* parameters: */
646 /* prob ... probability vector */
647 /* n_prob ... length of probability vector */
648 /* */
649 /* return: */
650 /* pointer to generator object */
651 /* */
652 /* error: */
653 /* return NULL */
654 /*---------------------------------------------------------------------------*/
655 {
656 struct unur_distr *distr;
657 struct unur_par *par;
658 struct unur_gen *igen;
659
660 /* create generator */
661 distr = unur_distr_discr_new();
662 unur_distr_discr_set_pv(distr, prob, n_prob);
663 par = unur_dgt_new(distr);
664 igen = unur_init(par);
665
666 /* clear working space */
667 unur_distr_free(distr);
668
669 return igen;
670
671 } /* end of _unur_mixt_indexgen() */
672
673 /*---------------------------------------------------------------------------*/
674
675 int
_unur_mixt_get_boundary(struct unur_gen * gen)676 _unur_mixt_get_boundary( struct unur_gen *gen )
677 /*---------------------------------------------------------------------------*/
678 /* compute boundary of mixture and check for overlapping domains. */
679 /* */
680 /* parameters: */
681 /* gen ... pointer to generator object */
682 /* */
683 /* return: */
684 /* UNUR_SUCCESS ... when everything is o.k. */
685 /* UNUR_NULL ... when a component is NULL (which should not happen) */
686 /* UNUR_ERR_GEN_INVALID .. if domains of components overlap or not ordered */
687 /*---------------------------------------------------------------------------*/
688 {
689 int i;
690 int overlap = FALSE;
691 double comp_left, comp_right;
692 double bd_left, bd_right;
693 struct unur_gen *comp;
694
695 /* initialize boundary values for mixture */
696 bd_left = INFINITY;
697 bd_right = -INFINITY;
698
699 for (i=0; i<gen->N_COMP; i++) {
700 comp = gen->COMP[i];
701
702 /* comp has already been check. but we want to be on the safe side. */
703 CHECK_NULL(comp,UNUR_ERR_NULL);
704
705 /* domain of component [i] */
706 switch (comp->method & UNUR_MASK_TYPE) {
707 case UNUR_METH_CONT:
708 comp_left = comp->distr->data.cont.BD_LEFT;
709 comp_right = comp->distr->data.cont.BD_RIGHT;
710 break;
711
712 case UNUR_METH_DISCR:
713 comp_left = (double) (comp->distr->data.discr.BD_LEFT);
714 comp_right = (double) (comp->distr->data.discr.BD_RIGHT);
715 break;
716
717 default:
718 /* cannot estimate boundary */
719 comp_left = -INFINITY;
720 comp_right = INFINITY;
721 }
722
723 /* check for overlapping domains */
724 if ( _unur_FP_less(comp_left,bd_right) )
725 overlap = TRUE;
726
727 /* update domain of mixture */
728 bd_left = _unur_min(bd_left, comp_left);
729 bd_right = _unur_max(bd_right, comp_right);
730 }
731
732 /* overlap or unordered domains? */
733 if (GEN->is_inversion && overlap) {
734 _unur_error(gen->genid,UNUR_ERR_GEN_INVALID,"domains of components overlap or are unsorted");
735 return UNUR_ERR_GEN_INVALID;
736 }
737
738 /* store boundary */
739 unur_distr_cont_set_domain(gen->distr, bd_left, bd_right);
740
741 /* o.k. */
742 return UNUR_SUCCESS;
743 } /* end of _unur_mixt_get_boundary() */
744
745
746 /*****************************************************************************/
747 /** Debugging utilities **/
748 /*****************************************************************************/
749
750 /*---------------------------------------------------------------------------*/
751 #ifdef UNUR_ENABLE_LOGGING
752 /*---------------------------------------------------------------------------*/
753
754 void
_unur_mixt_debug_init(const struct unur_gen * gen)755 _unur_mixt_debug_init( const struct unur_gen *gen )
756 /*----------------------------------------------------------------------*/
757 /* write info about generator into LOG file */
758 /* */
759 /* parameters: */
760 /* gen ... pointer to generator object */
761 /*----------------------------------------------------------------------*/
762 {
763 FILE *LOG;
764 struct unur_gen *comp;
765 int i;
766
767 /* check arguments */
768 CHECK_NULL(gen,RETURN_VOID); COOKIE_CHECK(gen,CK_MIXT_GEN,RETURN_VOID);
769
770 LOG = unur_get_stream();
771
772 fprintf(LOG,"%s:\n",gen->genid);
773 fprintf(LOG,"%s: type = continuous univariate random variates\n",gen->genid);
774 fprintf(LOG,"%s: method = MIXT (MIXTure of distributions -- meta method)\n",gen->genid);
775 fprintf(LOG,"%s:\n",gen->genid);
776
777 _unur_distr_cont_debug( gen->distr, gen->genid );
778
779 fprintf(LOG,"%s: sampling routine = _unur_mixt_sample",gen->genid);
780 if (GEN->is_inversion) fprintf(LOG,"_inv");
781 fprintf(LOG,"()\n%s:\n",gen->genid);
782
783 fprintf(LOG,"%s: use inversion = %s",gen->genid,
784 (GEN->is_inversion) ? "on" : "off");
785 _unur_print_if_default(gen,MIXT_SET_USEINVERSION);
786 fprintf(LOG,"\n%s:\n",gen->genid);
787
788 /* probabilities */
789 fprintf(LOG,"%s: probabilities (%d) = \n",gen->genid, gen->N_COMP);
790 fprintf(LOG,"%s: %g",gen->genid, (gen->PROB)[0]);
791 for (i=1; i<gen->N_COMP; i++)
792 fprintf(LOG,", %g", (gen->PROB)[i]);
793 fprintf(LOG,"\n%s:\n",gen->genid);
794
795 /* components */
796 fprintf(LOG,"%s: components (%d):\n",gen->genid, gen->N_COMP);
797 for (i=0; i<gen->N_COMP; i++) {
798 comp = gen->COMP[i];
799 fprintf(LOG,"%s: [%d]: %s\n",gen->genid, i, comp->genid);
800 fprintf(LOG,"%s:\t type = ",gen->genid);
801 switch (comp->distr->type) {
802 case UNUR_DISTR_CONT:
803 case UNUR_DISTR_CEMP:
804 fprintf(LOG,"continuous\n");
805 break;
806 case UNUR_DISTR_DISCR:
807 fprintf(LOG,"discrete\n");
808 break;
809 default:
810 fprintf(LOG,"[unknown]\n");
811 }
812 fprintf(LOG,"%s:\t name = %s\n",gen->genid, comp->distr->name);
813 }
814 fprintf(LOG,"%s:\n",gen->genid);
815
816 } /* end of _unur_mixt_debug_init() */
817
818 /*---------------------------------------------------------------------------*/
819 #endif /* end UNUR_ENABLE_LOGGING */
820 /*---------------------------------------------------------------------------*/
821
822
823 /*---------------------------------------------------------------------------*/
824 #ifdef UNUR_ENABLE_INFO
825 /*---------------------------------------------------------------------------*/
826
827 void
_unur_mixt_info(struct unur_gen * gen,int help)828 _unur_mixt_info( struct unur_gen *gen, int help )
829 /*----------------------------------------------------------------------*/
830 /* create character string that contains information about the */
831 /* given generator object. */
832 /* */
833 /* parameters: */
834 /* gen ... pointer to generator object */
835 /* help ... whether to print additional comments */
836 /*----------------------------------------------------------------------*/
837 {
838 struct unur_string *info = gen->infostr;
839 struct unur_gen *comp;
840 int i;
841 double sum;
842
843 /* generator ID */
844 _unur_string_append(info,"generator ID: %s\n\n", gen->genid);
845
846 /* distribution */
847 _unur_string_append(info,"distribution:\n");
848 _unur_distr_info_typename(gen);
849 _unur_string_append(info," domain = (%g, %g)\n", DISTR.domain[0],DISTR.domain[1]);
850 _unur_string_append(info," # components = %d\n", gen->N_COMP);
851
852 if (help) {
853 sum = ((struct unur_dgt_gen*)gen->INDEX->datap)->sum;
854 _unur_string_append(info," probabilities = (%g", gen->PROB[0] / sum);
855 for (i=1; i<gen->N_COMP; i++)
856 _unur_string_append(info,", %g", gen->PROB[i] / sum);
857 _unur_string_append(info,")\n");
858
859 _unur_string_append(info," components = \n");
860 for (i=0; i<gen->N_COMP; i++) {
861 comp = gen->COMP[i];
862 _unur_string_append(info,"\t[%d] %s - ",i, comp->genid);
863 switch (comp->distr->type) {
864 case UNUR_DISTR_CONT:
865 case UNUR_DISTR_CEMP:
866 _unur_string_append(info,"continuous");
867 break;
868 case UNUR_DISTR_DISCR:
869 _unur_string_append(info,"discrete");
870 break;
871 default:
872 _unur_string_append(info,"[unknown]");
873 }
874 _unur_string_append(info,": %s\n",comp->distr->name);
875 }
876 }
877 _unur_string_append(info,"\n");
878
879 /* method */
880 _unur_string_append(info,"method: MIXT (MIXTure of distributions -- meta method)\n");
881 _unur_string_append(info," select component = method DGT\n");
882 _unur_string_append(info," inversion method = %s\n",
883 (GEN->is_inversion) ? "TRUE" : "FALSE");
884 _unur_string_append(info,"\n");
885
886 /* performance */
887 _unur_string_append(info,"performance characteristics: depends on components\n");
888 _unur_string_append(info,"\n");
889
890 /* parameters */
891 if (help) {
892 _unur_string_append(info,"parameters:\n");
893 _unur_string_append(info," useinversion = ");
894 if (gen->variant & MIXT_VARFLAG_INVERSION)
895 _unur_string_append(info,"on\n");
896 else
897 _unur_string_append(info,"off [default]\n");
898 }
899
900 } /* end of _unur_mixt_info() */
901
902 /*---------------------------------------------------------------------------*/
903 #endif /* end UNUR_ENABLE_INFO */
904 /*---------------------------------------------------------------------------*/
905