1 /*****************************************************************************
2 * *
3 * UNURAN -- Universal Non-Uniform Random number generator *
4 * *
5 *****************************************************************************
6 * *
7 * FILE: hist.c *
8 * *
9 * TYPE: continuous univariate random variate *
10 * METHOD: generate from histogram *
11 * *
12 * DESCRIPTION: *
13 * Given histogram of observed sample. *
14 * Produce a value x consistent with this histogram (use inversion) *
15 * *
16 * REQUIRED: *
17 * pointer to histogram *
18 * *
19 *****************************************************************************
20 * *
21 * Copyright (c) 2000-2006 Wolfgang Hoermann and Josef Leydold *
22 * Department of Statistics and Mathematics, WU Wien, Austria *
23 * *
24 * This program is free software; you can redistribute it and/or modify *
25 * it under the terms of the GNU General Public License as published by *
26 * the Free Software Foundation; either version 2 of the License, or *
27 * (at your option) any later version. *
28 * *
29 * This program is distributed in the hope that it will be useful, *
30 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
31 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
32 * GNU General Public License for more details. *
33 * *
34 * You should have received a copy of the GNU General Public License *
35 * along with this program; if not, write to the *
36 * Free Software Foundation, Inc., *
37 * 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA *
38 * *
39 *****************************************************************************/
40
41 /*---------------------------------------------------------------------------*/
42
43 #include <unur_source.h>
44 #include <distr/distr.h>
45 #include <distr/distr_source.h>
46 #include <urng/urng.h>
47 #include "unur_methods_source.h"
48 #include "x_gen_source.h"
49 #include "hist.h"
50 #include "hist_struct.h"
51
52 /*---------------------------------------------------------------------------*/
53 /* Variants: none */
54
55 /*---------------------------------------------------------------------------*/
56 /* Debugging flags */
57 /* bit 01 ... pameters and structure of generator (do not use here) */
58 /* bits 02-12 ... setup */
59 /* bits 13-24 ... adaptive steps */
60 /* bits 25-32 ... trace sampling */
61
62 #define HIST_DEBUG_PRINTHIST 0x00000100u
63
64 /*---------------------------------------------------------------------------*/
65 /* Flags for logging set calls */
66
67 /*---------------------------------------------------------------------------*/
68
69 #define GENTYPE "HIST" /* type of generator */
70
71 /*---------------------------------------------------------------------------*/
72
73 static struct unur_gen *_unur_hist_init( struct unur_par *par );
74 /*---------------------------------------------------------------------------*/
75 /* Initialize new generator. */
76 /*---------------------------------------------------------------------------*/
77
78 static struct unur_gen *_unur_hist_create( struct unur_par *par );
79 /*---------------------------------------------------------------------------*/
80 /* create new (almost empty) generator object. */
81 /*---------------------------------------------------------------------------*/
82
83 static struct unur_gen *_unur_hist_clone( const struct unur_gen *gen );
84 /*---------------------------------------------------------------------------*/
85 /* copy (clone) generator object. */
86 /*---------------------------------------------------------------------------*/
87
88 static void _unur_hist_free( struct unur_gen *gen);
89 /*---------------------------------------------------------------------------*/
90 /* destroy generator object. */
91 /*---------------------------------------------------------------------------*/
92
93 static double _unur_hist_sample( struct unur_gen *gen );
94 /*---------------------------------------------------------------------------*/
95 /* sample from generator */
96 /*---------------------------------------------------------------------------*/
97
98 static int _unur_hist_create_tables( struct unur_gen *gen );
99 /*---------------------------------------------------------------------------*/
100 /* create (allocate) tables */
101 /*---------------------------------------------------------------------------*/
102
103 static int _unur_hist_make_guidetable( struct unur_gen *gen );
104 /*---------------------------------------------------------------------------*/
105 /* create table for indexed search */
106 /*---------------------------------------------------------------------------*/
107
108 #ifdef UNUR_ENABLE_LOGGING
109 /*---------------------------------------------------------------------------*/
110 /* the following functions print debugging information on output stream, */
111 /* i.e., into the LOG file if not specified otherwise. */
112 /*---------------------------------------------------------------------------*/
113 static void _unur_hist_debug_init( const struct unur_gen *gen );
114
115 /*---------------------------------------------------------------------------*/
116 /* print after generator has been initialized has completed. */
117 /*---------------------------------------------------------------------------*/
118 #endif
119
120 #ifdef UNUR_ENABLE_INFO
121 static void _unur_hist_info( struct unur_gen *gen, int help );
122 /*---------------------------------------------------------------------------*/
123 /* create info string. */
124 /*---------------------------------------------------------------------------*/
125 #endif
126
127 /*---------------------------------------------------------------------------*/
128 /* abbreviations */
129
130 #define DISTR_IN distr->data.cemp /* data for distribution object */
131
132 #define PAR ((struct unur_hist_par*)par->datap) /* data for parameter object */
133 #define GEN ((struct unur_hist_gen*)gen->datap) /* data for generator object */
134 #define DISTR gen->distr->data.cemp /* data for distribution in generator object */
135
136 #define SAMPLE gen->sample.cont /* pointer to sampling routine */
137
138 /*---------------------------------------------------------------------------*/
139 /* constants */
140
141 /*---------------------------------------------------------------------------*/
142
143 #define _unur_hist_getSAMPLE(gen) (_unur_hist_sample)
144
145 /*---------------------------------------------------------------------------*/
146
147 /*****************************************************************************/
148 /** Public: User Interface (API) **/
149 /*****************************************************************************/
150
151 struct unur_par *
unur_hist_new(const struct unur_distr * distr)152 unur_hist_new( const struct unur_distr *distr )
153 /*----------------------------------------------------------------------*/
154 /* get default parameters */
155 /* */
156 /* parameters: */
157 /* distr ... pointer to distribution object */
158 /* */
159 /* return: */
160 /* default parameters (pointer to structure) */
161 /* */
162 /* error: */
163 /* return NULL */
164 /*----------------------------------------------------------------------*/
165 {
166 struct unur_par *par;
167
168 /* check arguments */
169 _unur_check_NULL( GENTYPE,distr,NULL );
170
171 /* check distribution */
172 if (distr->type != UNUR_DISTR_CEMP) {
173 _unur_error(GENTYPE,UNUR_ERR_DISTR_INVALID,""); return NULL; }
174 COOKIE_CHECK(distr,CK_DISTR_CEMP,NULL);
175
176 if (DISTR_IN.hist_prob == NULL || !(distr->set & UNUR_DISTR_SET_DOMAIN)) {
177 _unur_error(GENTYPE,UNUR_ERR_DISTR_REQUIRED,"histogram"); return NULL; }
178
179 /* allocate structure */
180 par = _unur_par_new( sizeof(struct unur_hist_par) );
181 COOKIE_SET(par,CK_HIST_PAR);
182
183 /* copy input */
184 par->distr = distr; /* pointer to distribution object */
185
186 /* set default values */
187 par->method = UNUR_METH_HIST; /* method */
188 par->variant = 0u; /* default variant */
189
190 par->set = 0u; /* inidicate default parameters */
191 par->urng = unur_get_default_urng(); /* use default urng */
192 par->urng_aux = NULL; /* no auxilliary URNG required */
193
194 par->debug = _unur_default_debugflag; /* set default debugging flags */
195
196 /* routine for starting generator */
197 par->init = _unur_hist_init;
198
199 return par;
200
201 } /* end of unur_hist_new() */
202
203 /*****************************************************************************/
204 /** Private **/
205 /*****************************************************************************/
206
207 struct unur_gen *
_unur_hist_init(struct unur_par * par)208 _unur_hist_init( struct unur_par *par )
209 /*----------------------------------------------------------------------*/
210 /* initialize new generator */
211 /* */
212 /* parameters: */
213 /* par ... pointer to paramters for building generator object */
214 /* */
215 /* return: */
216 /* pointer to generator object */
217 /* */
218 /* error: */
219 /* return NULL */
220 /*----------------------------------------------------------------------*/
221 {
222 struct unur_gen *gen;
223
224 /* check arguments */
225 CHECK_NULL(par,NULL);
226
227 /* check input */
228 if ( par->method != UNUR_METH_HIST ) {
229 _unur_error(GENTYPE,UNUR_ERR_PAR_INVALID,"");
230 return NULL; }
231 COOKIE_CHECK(par,CK_HIST_PAR,NULL);
232
233 /* create a new empty generator object */
234 gen = _unur_hist_create(par);
235 _unur_par_free(par);
236 if (!gen) return NULL;
237
238 /* compute guide table */
239 if ( (_unur_hist_create_tables(gen) != UNUR_SUCCESS) ||
240 (_unur_hist_make_guidetable(gen) != UNUR_SUCCESS) ) {
241 _unur_hist_free(gen); return NULL;
242 }
243
244 #ifdef UNUR_ENABLE_LOGGING
245 /* write info into LOG file */
246 if (gen->debug) _unur_hist_debug_init(gen);
247 #endif
248
249 return gen;
250
251 } /* end of _unur_hist_init() */
252
253 /*---------------------------------------------------------------------------*/
254
255 struct unur_gen *
_unur_hist_create(struct unur_par * par)256 _unur_hist_create( struct unur_par *par )
257 /*----------------------------------------------------------------------*/
258 /* allocate memory for generator */
259 /* */
260 /* parameters: */
261 /* par ... pointer to parameter for building generator object */
262 /* */
263 /* return: */
264 /* pointer to (empty) generator object with default settings */
265 /* */
266 /* error: */
267 /* return NULL */
268 /*----------------------------------------------------------------------*/
269 {
270 struct unur_gen *gen;
271
272 /* check arguments */
273 CHECK_NULL(par,NULL); COOKIE_CHECK(par,CK_HIST_PAR,NULL);
274
275 /* create new generic generator object */
276 gen = _unur_generic_create( par, sizeof(struct unur_hist_gen) );
277
278 /* magic cookies */
279 COOKIE_SET(gen,CK_HIST_GEN);
280
281 /* set generator identifier */
282 gen->genid = _unur_set_genid(GENTYPE);
283
284 /* routines for sampling and destroying generator */
285 SAMPLE = _unur_hist_getSAMPLE(gen);
286 gen->destroy = _unur_hist_free;
287 gen->clone = _unur_hist_clone;
288
289 /* make sure that the domain coincides with bin data */
290 if (DISTR.hist_bins) {
291 DISTR.hmin = DISTR.hist_bins[0];
292 DISTR.hmax = DISTR.hist_bins[DISTR.n_hist];
293 }
294
295 /* copy observed data into generator object */
296 GEN->n_hist = DISTR.n_hist; /* size of histogram */
297 GEN->prob = DISTR.hist_prob; /* probabilities of bins */
298 GEN->hmin = DISTR.hmin; /* lower ... */
299 GEN->hmax = DISTR.hmax; /* ... and upper bound */
300 GEN->hwidth = (DISTR.hmax - DISTR.hmin) / DISTR.n_hist;
301 GEN->bins = (DISTR.hist_bins) ? DISTR.hist_bins : NULL;
302
303 /* set all pointers to NULL */
304 GEN->sum = 0.;
305 GEN->cumpv = NULL;
306 GEN->guide_table = NULL;
307
308 #ifdef UNUR_ENABLE_INFO
309 /* set function for creating info string */
310 gen->info = _unur_hist_info;
311 #endif
312
313 /* return pointer to (almost empty) generator object */
314 return gen;
315
316 } /* end of _unur_hist_create() */
317
318 /*---------------------------------------------------------------------------*/
319
320 struct unur_gen *
_unur_hist_clone(const struct unur_gen * gen)321 _unur_hist_clone( const struct unur_gen *gen )
322 /*----------------------------------------------------------------------*/
323 /* copy (clone) generator object */
324 /* */
325 /* parameters: */
326 /* gen ... pointer to generator object */
327 /* */
328 /* return: */
329 /* pointer to clone of generator object */
330 /* */
331 /* error: */
332 /* return NULL */
333 /*----------------------------------------------------------------------*/
334 {
335 #define CLONE ((struct unur_hist_gen*)clone->datap)
336
337 struct unur_gen *clone;
338
339 /* check arguments */
340 CHECK_NULL(gen,NULL); COOKIE_CHECK(gen,CK_HIST_GEN,NULL);
341
342 /* create generic clone */
343 clone = _unur_generic_clone( gen, GENTYPE );
344
345 /* copy histrogram into generator object */
346 CLONE->prob = clone->distr->data.cemp.hist_prob; /* probabilities of bins */
347 CLONE->bins = clone->distr->data.cemp.hist_bins; /* location of bins */
348
349 /* copy data for distribution */
350 CLONE->cumpv = _unur_xmalloc( GEN->n_hist * sizeof(double) );
351 memcpy( CLONE->cumpv, GEN->cumpv, GEN->n_hist * sizeof(double) );
352 CLONE->guide_table = _unur_xmalloc( GEN->n_hist * sizeof(int) );
353 memcpy( CLONE->guide_table, GEN->guide_table, GEN->n_hist * sizeof(int) );
354
355 return clone;
356
357 #undef CLONE
358 } /* end of _unur_hist_clone() */
359
360 /*---------------------------------------------------------------------------*/
361
362 void
_unur_hist_free(struct unur_gen * gen)363 _unur_hist_free( struct unur_gen *gen )
364 /*----------------------------------------------------------------------*/
365 /* deallocate generator object */
366 /* */
367 /* parameters: */
368 /* gen ... pointer to generator object */
369 /*----------------------------------------------------------------------*/
370 {
371
372 /* check arguments */
373 if( !gen ) /* nothing to do */
374 return;
375
376 /* check input */
377 if ( gen->method != UNUR_METH_HIST ) {
378 _unur_warning(gen->genid,UNUR_ERR_GEN_INVALID,"");
379 return; }
380 COOKIE_CHECK(gen,CK_HIST_GEN,RETURN_VOID);
381
382 /* we cannot use this generator object any more */
383 SAMPLE = NULL; /* make sure to show up a programming error */
384
385 /* free two auxiliary tables */
386 if (GEN->guide_table) free(GEN->guide_table);
387 if (GEN->cumpv) free(GEN->cumpv);
388
389 /* free memory */
390 _unur_generic_free(gen);
391
392 } /* end of _unur_hist_free() */
393
394 /*****************************************************************************/
395
396 double
_unur_hist_sample(struct unur_gen * gen)397 _unur_hist_sample( struct unur_gen *gen )
398 /*----------------------------------------------------------------------*/
399 /* sample from generator */
400 /* */
401 /* parameters: */
402 /* gen ... pointer to generator object */
403 /* */
404 /* return: */
405 /* double (sample from random variate) */
406 /* */
407 /* error: */
408 /* return INFINITY */
409 /*----------------------------------------------------------------------*/
410 {
411 double U;
412 int J;
413
414 /* check arguments */
415 CHECK_NULL(gen,INFINITY); COOKIE_CHECK(gen,CK_HIST_GEN,INFINITY);
416
417 /* sample from U(0,1) */
418 U = _unur_call_urng(gen->urng);
419
420 /* look up in guide table ... */
421 J = GEN->guide_table[(int)(U * GEN->n_hist)];
422 /* ... and search */
423 U *= GEN->sum;
424 while (GEN->cumpv[J] < U) J++;
425
426 /* reuse of uniform random number: U ~ U(0,1) */
427 U = (U - (J ? GEN->cumpv[J-1] : 0.)) / GEN->prob[J];
428
429 /* sample uniformly from bin */
430 if (GEN->bins)
431 return (U * GEN->bins[J+1] + (1.-U) * GEN->bins[J]);
432 else
433 return (GEN->hmin + (U+J)*GEN->hwidth);
434
435 } /* end of _unur_hist_sample() */
436
437
438 /*****************************************************************************/
439 /** Auxilliary Routines **/
440 /*****************************************************************************/
441
442 int
_unur_hist_create_tables(struct unur_gen * gen)443 _unur_hist_create_tables( struct unur_gen *gen )
444 /*----------------------------------------------------------------------*/
445 /* create (allocate) tables */
446 /* */
447 /* parameters: */
448 /* gen ... pointer to generator object */
449 /* */
450 /* return: */
451 /* UNUR_SUCCESS ... on success */
452 /* error code ... on error */
453 /*----------------------------------------------------------------------*/
454 {
455 /* allocation for cummulated probabilities */
456 GEN->cumpv = _unur_xrealloc( GEN->cumpv, GEN->n_hist * sizeof(double) );
457
458 /* allocate memory for the guide table */
459 GEN->guide_table = _unur_xrealloc( GEN->guide_table, GEN->n_hist * sizeof(int) );
460
461 /* o.k. */
462 return UNUR_SUCCESS;
463 } /* end of _unur_hist_create_tables() */
464
465 /*---------------------------------------------------------------------------*/
466
467 int
_unur_hist_make_guidetable(struct unur_gen * gen)468 _unur_hist_make_guidetable( struct unur_gen *gen )
469 /*----------------------------------------------------------------------*/
470 /* create guide table */
471 /* */
472 /* parameters: */
473 /* gen ... pointer to generator object */
474 /* */
475 /* return: */
476 /* UNUR_SUCCESS ... on success */
477 /* error code ... on error */
478 /*----------------------------------------------------------------------*/
479 {
480 double *pv; /* pointer to probability vector */
481 int n_pv; /* length of probability vector */
482 double pvh; /* aux variable for computing cumulated sums */
483 double gstep; /* step size when computing guide table */
484 int i,j;
485
486 /* probability vector */
487 pv = GEN->prob;
488 n_pv = GEN->n_hist;
489
490 /* computation of cumulated probabilities */
491 for( i=0, pvh=0.; i<n_pv; i++ ) {
492 GEN->cumpv[i] = ( pvh += pv[i] );
493 /* ... and check probability vector */
494 if (pv[i] < 0.) {
495 _unur_error(gen->genid,UNUR_ERR_GEN_DATA,"probability < 0");
496 return UNUR_ERR_GEN_DATA;
497 }
498 }
499 GEN->sum = GEN->cumpv[n_pv-1];
500
501 /* computation of guide-table */
502 gstep = GEN->sum / GEN->n_hist;
503 pvh = 0.;
504 for( j=0, i=0; j<GEN->n_hist;j++ ) {
505 while (GEN->cumpv[i] < pvh)
506 i++;
507 if (i >= n_pv) {
508 _unur_warning(gen->genid,UNUR_ERR_ROUNDOFF,"guide table");
509 break;
510 }
511 GEN->guide_table[j] = i;
512 pvh += gstep;
513 }
514
515 /* if there has been an round off error, we have to complete the guide table */
516 for( ; j<GEN->n_hist; j++ )
517 GEN->guide_table[j] = n_pv - 1;
518 /* o.k. */
519 return UNUR_SUCCESS;
520
521 } /* end of _unur_hist_make_urntable() */
522
523
524 /*****************************************************************************/
525 /** Debugging utilities **/
526 /*****************************************************************************/
527
528 /*---------------------------------------------------------------------------*/
529 #ifdef UNUR_ENABLE_LOGGING
530 /*---------------------------------------------------------------------------*/
531
532 void
_unur_hist_debug_init(const struct unur_gen * gen)533 _unur_hist_debug_init( const struct unur_gen *gen )
534 /*----------------------------------------------------------------------*/
535 /* write info about generator into LOG file */
536 /* */
537 /* parameters: */
538 /* gen ... pointer to generator object */
539 /*----------------------------------------------------------------------*/
540 {
541 FILE *LOG;
542
543 /* check arguments */
544 CHECK_NULL(gen,RETURN_VOID); COOKIE_CHECK(gen,CK_HIST_GEN,RETURN_VOID);
545
546 LOG = unur_get_stream();
547
548 fprintf(LOG,"%s:\n",gen->genid);
549 fprintf(LOG,"%s: type = continuous univariate random variates\n",gen->genid);
550 fprintf(LOG,"%s: method = HIST (HISTogram of empirical distribution)\n",gen->genid);
551 fprintf(LOG,"%s:\n",gen->genid);
552
553 _unur_distr_cemp_debug( gen->distr, gen->genid, (gen->debug & HIST_DEBUG_PRINTHIST));
554
555 fprintf(LOG,"%s: sampling routine = _unur_hist_sample()\n",gen->genid);
556 fprintf(LOG,"%s:\n",gen->genid);
557
558 fprintf(LOG,"%s:\n",gen->genid);
559
560 } /* end of _unur_hist_debug_init() */
561
562 /*---------------------------------------------------------------------------*/
563 #endif /* end UNUR_ENABLE_LOGGING */
564 /*---------------------------------------------------------------------------*/
565
566
567 /*---------------------------------------------------------------------------*/
568 #ifdef UNUR_ENABLE_INFO
569 /*---------------------------------------------------------------------------*/
570
571 void
_unur_hist_info(struct unur_gen * gen,int help)572 _unur_hist_info( struct unur_gen *gen, int help )
573 /*----------------------------------------------------------------------*/
574 /* create character string that contains information about the */
575 /* given generator object. */
576 /* */
577 /* parameters: */
578 /* gen ... pointer to generator object */
579 /* help ... whether to print additional comments */
580 /*----------------------------------------------------------------------*/
581 {
582 struct unur_string *info = gen->infostr;
583
584 /* generator ID */
585 _unur_string_append(info,"generator ID: %s\n\n", gen->genid);
586
587 /* distribution */
588 _unur_string_append(info,"distribution:\n");
589 _unur_distr_info_typename(gen);
590 _unur_string_append(info," functions = DATA [histogram of size=%d]\n", DISTR.n_hist);
591 _unur_string_append(info,"\n");
592
593 /* if (help) { */
594 /* _unur_string_append(info,"\n"); */
595 /* } */
596
597 /* method */
598 _unur_string_append(info,"method: HIST (HISTogramm of empirical distribution)\n");
599 _unur_string_append(info,"\n");
600
601 /* performance */
602 /* _unur_string_append(info,"performance characteristics:\n"); */
603 /* _unur_string_append(info,"\n"); */
604
605 /* parameters */
606 if (help) {
607 _unur_string_append(info,"parameters: none\n");
608 _unur_string_append(info,"\n");
609 }
610
611 /* Hints */
612 /* if (help) { */
613 /* _unur_string_append(info,"\n"); */
614 /* } */
615
616 } /* end of _unur_hist_info() */
617
618 /*---------------------------------------------------------------------------*/
619 #endif /* end UNUR_ENABLE_INFO */
620 /*---------------------------------------------------------------------------*/
621