GNU Octave  3.8.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
oct-rand.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 2003-2013 John W. Eaton
4 
5 This file is part of Octave.
6 
7 Octave is free software; you can redistribute it and/or modify it
8 under the terms of the GNU General Public License as published by the
9 Free Software Foundation; either version 3 of the License, or (at your
10 option) any later version.
11 
12 Octave is distributed in the hope that it will be useful, but WITHOUT
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
15 for more details.
16 
17 You should have received a copy of the GNU General Public License
18 along with Octave; see the file COPYING. If not, see
19 <http://www.gnu.org/licenses/>.
20 
21 */
22 
23 #ifdef HAVE_CONFIG_H
24 #include <config.h>
25 #endif
26 
27 #include <map>
28 #include <vector>
29 
30 #include <stdint.h>
31 
32 #include "data-conv.h"
33 #include "f77-fcn.h"
34 #include "lo-error.h"
35 #include "lo-ieee.h"
36 #include "lo-mappers.h"
37 #include "mach-info.h"
38 #include "oct-locbuf.h"
39 #include "oct-rand.h"
40 #include "oct-time.h"
41 #include "randgamma.h"
42 #include "randmtzig.h"
43 #include "randpoisson.h"
44 #include "singleton-cleanup.h"
45 
46 extern "C"
47 {
48  F77_RET_T
49  F77_FUNC (dgennor, DGENNOR) (const double&, const double&, double&);
50 
51  F77_RET_T
52  F77_FUNC (dgenunf, DGENUNF) (const double&, const double&, double&);
53 
54  F77_RET_T
55  F77_FUNC (dgenexp, DGENEXP) (const double&, double&);
56 
57  F77_RET_T
58  F77_FUNC (dignpoi, DIGNPOI) (const double&, double&);
59 
60  F77_RET_T
61  F77_FUNC (dgengam, DGENGAM) (const double&, const double&, double&);
62 
63  F77_RET_T
64  F77_FUNC (setall, SETALL) (const int32_t&, const int32_t&);
65 
66  F77_RET_T
67  F77_FUNC (getsd, GETSD) (int32_t&, int32_t&);
68 
69  F77_RET_T
70  F77_FUNC (setsd, SETSD) (const int32_t&, const int32_t&);
71 
72  F77_RET_T
73  F77_FUNC (setcgn, SETCGN) (const int32_t&);
74 }
75 
77 
79  : current_distribution (uniform_dist), use_old_generators (false),
80  rand_states ()
81 {
83 
85 }
86 
87 bool
89 {
90  bool retval = true;
91 
92  if (! instance)
93  {
94  instance = new octave_rand ();
95 
96  if (instance)
98  }
99 
100  if (! instance)
101  {
102  (*current_liboctave_error_handler)
103  ("unable to create octave_rand object!");
104 
105  retval = false;
106  }
107 
108  return retval;
109 }
110 
111 double
113 {
114  union d2i { double d; int32_t i[2]; };
115  union d2i u;
116 
118 
119  switch (ff)
120  {
122  F77_FUNC (getsd, GETSD) (u.i[1], u.i[0]);
123  break;
124  default:
125  F77_FUNC (getsd, GETSD) (u.i[0], u.i[1]);
126  break;
127  }
128 
129  return u.d;
130 }
131 
132 static int32_t
133 force_to_fit_range (int32_t i, int32_t lo, int32_t hi)
134 {
135  assert (hi > lo && lo >= 0 && hi > lo);
136 
137  i = i > 0 ? i : -i;
138 
139  if (i < lo)
140  i = lo;
141  else if (i > hi)
142  i = i % hi;
143 
144  return i;
145 }
146 
147 void
149 {
150  use_old_generators = true;
151 
152  int i0, i1;
153  union d2i { double d; int32_t i[2]; };
154  union d2i u;
155  u.d = s;
156 
158 
159  switch (ff)
160  {
162  i1 = force_to_fit_range (u.i[0], 1, 2147483563);
163  i0 = force_to_fit_range (u.i[1], 1, 2147483399);
164  break;
165  default:
166  i0 = force_to_fit_range (u.i[0], 1, 2147483563);
167  i1 = force_to_fit_range (u.i[1], 1, 2147483399);
168  break;
169  }
170 
171  F77_FUNC (setsd, SETSD) (i0, i1);
172 }
173 
174 void
176 {
177  use_old_generators = true;
179 }
180 
182 octave_rand::do_state (const std::string& d)
183 {
184  return rand_states[d.empty () ? current_distribution : get_dist_id (d)];
185 }
186 
187 void
188 octave_rand::do_state (const ColumnVector& s, const std::string& d)
189 {
190  use_old_generators = false;
191 
192  int old_dist = current_distribution;
193 
194  int new_dist = d.empty () ? current_distribution : get_dist_id (d);
195 
196  ColumnVector saved_state;
197 
198  if (old_dist != new_dist)
199  saved_state = get_internal_state ();
200 
201  set_internal_state (s);
202 
203  rand_states[new_dist] = get_internal_state ();
204 
205  if (old_dist != new_dist)
206  rand_states[old_dist] = saved_state;
207 }
208 
209 void
210 octave_rand::do_reset (const std::string& d)
211 {
212  use_old_generators = false;
213 
214  int old_dist = current_distribution;
215 
216  int new_dist = d.empty () ? current_distribution : get_dist_id (d);
217 
218  ColumnVector saved_state;
219 
220  if (old_dist != new_dist)
221  saved_state = get_internal_state ();
222 
224  rand_states[new_dist] = get_internal_state ();
225 
226  if (old_dist != new_dist)
227  rand_states[old_dist] = saved_state;
228 }
229 
230 std::string
232 {
233  std::string retval;
234 
235  switch (current_distribution)
236  {
237  case uniform_dist:
238  retval = "uniform";
239  break;
240 
241  case normal_dist:
242  retval = "normal";
243  break;
244 
245  case expon_dist:
246  retval = "exponential";
247  break;
248 
249  case poisson_dist:
250  retval = "poisson";
251  break;
252 
253  case gamma_dist:
254  retval = "gamma";
255  break;
256 
257  default:
258  (*current_liboctave_error_handler)
259  ("rand: invalid distribution ID = %d", current_distribution);
260  break;
261  }
262 
263  return retval;
264 }
265 
266 void
267 octave_rand::do_distribution (const std::string& d)
268 {
269  int id = get_dist_id (d);
270 
271  switch (id)
272  {
273  case uniform_dist:
275  break;
276 
277  case normal_dist:
279  break;
280 
281  case expon_dist:
283  break;
284 
285  case poisson_dist:
287  break;
288 
289  case gamma_dist:
291  break;
292 
293  default:
294  (*current_liboctave_error_handler)
295  ("rand: invalid distribution ID = %d", id);
296  break;
297  }
298 }
299 
300 void
302 {
304 
305  F77_FUNC (setcgn, SETCGN) (uniform_dist);
306 }
307 
308 void
310 {
312 
313  F77_FUNC (setcgn, SETCGN) (normal_dist);
314 }
315 
316 void
318 {
320 
321  F77_FUNC (setcgn, SETCGN) (expon_dist);
322 }
323 
324 void
326 {
328 
329  F77_FUNC (setcgn, SETCGN) (poisson_dist);
330 }
331 
332 void
334 {
336 
337  F77_FUNC (setcgn, SETCGN) (gamma_dist);
338 }
339 
340 
341 double
343 {
344  double retval = 0.0;
345 
346  if (use_old_generators)
347  {
348  switch (current_distribution)
349  {
350  case uniform_dist:
351  F77_FUNC (dgenunf, DGENUNF) (0.0, 1.0, retval);
352  break;
353 
354  case normal_dist:
355  F77_FUNC (dgennor, DGENNOR) (0.0, 1.0, retval);
356  break;
357 
358  case expon_dist:
359  F77_FUNC (dgenexp, DGENEXP) (1.0, retval);
360  break;
361 
362  case poisson_dist:
363  if (a < 0.0 || ! xfinite (a))
364  retval = octave_NaN;
365  else
366  {
367  // workaround bug in ignpoi, by calling with different Mu
368  F77_FUNC (dignpoi, DIGNPOI) (a + 1, retval);
369  F77_FUNC (dignpoi, DIGNPOI) (a, retval);
370  }
371  break;
372 
373  case gamma_dist:
374  if (a <= 0.0 || ! xfinite (a))
375  retval = octave_NaN;
376  else
377  F77_FUNC (dgengam, DGENGAM) (1.0, a, retval);
378  break;
379 
380  default:
381  (*current_liboctave_error_handler)
382  ("rand: invalid distribution ID = %d", current_distribution);
383  break;
384  }
385  }
386  else
387  {
388  switch (current_distribution)
389  {
390  case uniform_dist:
391  retval = oct_randu ();
392  break;
393 
394  case normal_dist:
395  retval = oct_randn ();
396  break;
397 
398  case expon_dist:
399  retval = oct_rande ();
400  break;
401 
402  case poisson_dist:
403  retval = oct_randp (a);
404  break;
405 
406  case gamma_dist:
407  retval = oct_randg (a);
408  break;
409 
410  default:
411  (*current_liboctave_error_handler)
412  ("rand: invalid distribution ID = %d", current_distribution);
413  break;
414  }
415 
416  save_state ();
417  }
418 
419  return retval;
420 }
421 
422 float
424 {
425  float retval = 0.0;
426 
427  if (use_old_generators)
428  {
429  double da = a;
430  double dretval = 0.0;
431  switch (current_distribution)
432  {
433  case uniform_dist:
434  F77_FUNC (dgenunf, DGENUNF) (0.0, 1.0, dretval);
435  break;
436 
437  case normal_dist:
438  F77_FUNC (dgennor, DGENNOR) (0.0, 1.0, dretval);
439  break;
440 
441  case expon_dist:
442  F77_FUNC (dgenexp, DGENEXP) (1.0, dretval);
443  break;
444 
445  case poisson_dist:
446  if (da < 0.0 || ! xfinite (a))
447  dretval = octave_NaN;
448  else
449  {
450  // workaround bug in ignpoi, by calling with different Mu
451  F77_FUNC (dignpoi, DIGNPOI) (da + 1, dretval);
452  F77_FUNC (dignpoi, DIGNPOI) (da, dretval);
453  }
454  break;
455 
456  case gamma_dist:
457  if (da <= 0.0 || ! xfinite (a))
458  retval = octave_NaN;
459  else
460  F77_FUNC (dgengam, DGENGAM) (1.0, da, dretval);
461  break;
462 
463  default:
464  (*current_liboctave_error_handler)
465  ("rand: invalid distribution ID = %d", current_distribution);
466  break;
467  }
468  retval = dretval;
469  }
470  else
471  {
472  switch (current_distribution)
473  {
474  case uniform_dist:
475  retval = oct_float_randu ();
476  break;
477 
478  case normal_dist:
479  retval = oct_float_randn ();
480  break;
481 
482  case expon_dist:
483  retval = oct_float_rande ();
484  break;
485 
486  case poisson_dist:
487  // Keep poisson distribution in double precision for accuracy
488  retval = oct_randp (a);
489  break;
490 
491  case gamma_dist:
492  retval = oct_float_randg (a);
493  break;
494 
495  default:
496  (*current_liboctave_error_handler)
497  ("rand: invalid distribution ID = %d", current_distribution);
498  break;
499  }
500 
501  save_state ();
502  }
503 
504  return retval;
505 }
506 
509 {
510  Array<double> retval;
511 
512  if (n > 0)
513  {
514  retval.clear (n, 1);
515 
516  fill (retval.capacity (), retval.fortran_vec (), a);
517  }
518  else if (n < 0)
519  (*current_liboctave_error_handler) ("rand: invalid negative argument");
520 
521  return retval;
522 }
523 
526 {
527  Array<float> retval;
528 
529  if (n > 0)
530  {
531  retval.clear (n, 1);
532 
533  fill (retval.capacity (), retval.fortran_vec (), a);
534  }
535  else if (n < 0)
536  (*current_liboctave_error_handler) ("rand: invalid negative argument");
537 
538  return retval;
539 }
540 
541 NDArray
542 octave_rand::do_nd_array (const dim_vector& dims, double a)
543 {
544  NDArray retval;
545 
546  if (! dims.all_zero ())
547  {
548  retval.clear (dims);
549 
550  fill (retval.capacity (), retval.fortran_vec (), a);
551  }
552 
553  return retval;
554 }
555 
558 {
559  FloatNDArray retval;
560 
561  if (! dims.all_zero ())
562  {
563  retval.clear (dims);
564 
565  fill (retval.capacity (), retval.fortran_vec (), a);
566  }
567 
568  return retval;
569 }
570 
571 // Make the random number generator give us a different sequence every
572 // time we start octave unless we specifically set the seed. The
573 // technique used below will cycle monthly, but it it does seem to
574 // work ok to give fairly different seeds each time Octave starts.
575 
576 void
578 {
579  octave_localtime tm;
580  int stored_distribution = current_distribution;
581  F77_FUNC (setcgn, SETCGN) (uniform_dist);
582 
583  int hour = tm.hour () + 1;
584  int minute = tm.min () + 1;
585  int second = tm.sec () + 1;
586 
587  int32_t s0 = tm.mday () * hour * minute * second;
588  int32_t s1 = hour * minute * second;
589 
590  s0 = force_to_fit_range (s0, 1, 2147483563);
591  s1 = force_to_fit_range (s1, 1, 2147483399);
592 
593  F77_FUNC (setall, SETALL) (s0, s1);
594  F77_FUNC (setcgn, SETCGN) (stored_distribution);
595 }
596 
597 void
599 {
601 
603 
605 
607  s = get_internal_state ();
609 
611  s = get_internal_state ();
612  rand_states[expon_dist] = s;
613 
615  s = get_internal_state ();
617 
619  s = get_internal_state ();
620  rand_states[gamma_dist] = s;
621 }
622 
625 {
626  ColumnVector s (MT_N + 1);
627 
628  OCTAVE_LOCAL_BUFFER (uint32_t, tmp, MT_N + 1);
629 
630  oct_get_state (tmp);
631 
632  for (octave_idx_type i = 0; i <= MT_N; i++)
633  s.elem (i) = static_cast<double> (tmp[i]);
634 
635  return s;
636 }
637 
638 void
640 {
642 }
643 
644 int
645 octave_rand::get_dist_id (const std::string& d)
646 {
647  int retval = unknown_dist;
648 
649  if (d == "uniform" || d == "rand")
650  retval = uniform_dist;
651  else if (d == "normal" || d == "randn")
652  retval = normal_dist;
653  else if (d == "exponential" || d == "rande")
654  retval = expon_dist;
655  else if (d == "poisson" || d == "randp")
656  retval = poisson_dist;
657  else if (d == "gamma" || d == "randg")
658  retval = gamma_dist;
659  else
661  ("rand: invalid distribution '%s'", d.c_str ());
662 
663  return retval;
664 }
665 
666 // Guarantee reproducible conversion of negative initialization values to
667 // random number algorithm. Note that Matlab employs slightly different rules.
668 // 1) Seed saturates at 2^32-1 for any value larger than that.
669 // 2) NaN, Inf are translated to 2^32-1.
670 // 3) -Inf is translated to 0.
671 static uint32_t
672 double2uint32 (double d)
673 {
674  uint32_t u;
675  static const double TWOUP32 = std::numeric_limits<uint32_t>::max() + 1.0;
676 
677  if (! xfinite (d))
678  u = 0;
679  else
680  {
681  d = fmod (d, TWOUP32);
682  if (d < 0)
683  d += TWOUP32;
684  u = static_cast<uint32_t> (d);
685  }
686 
687  return u;
688 }
689 
690 void
692 {
693  octave_idx_type len = s.length ();
694  octave_idx_type n = len < MT_N + 1 ? len : MT_N + 1;
695 
696  OCTAVE_LOCAL_BUFFER (uint32_t, tmp, MT_N + 1);
697 
698  for (octave_idx_type i = 0; i < n; i++)
699  tmp[i] = double2uint32 (s.elem (i));
700 
701  if (len == MT_N + 1 && tmp[MT_N] <= MT_N && tmp[MT_N] > 0)
702  oct_set_state (tmp);
703  else
704  oct_init_by_array (tmp, len);
705 }
706 
707 void
709 {
710  if (dist != current_distribution)
711  {
712  current_distribution = dist;
713 
715  }
716 }
717 
718 #define MAKE_RAND(len) \
719  do \
720  { \
721  double val; \
722  for (volatile octave_idx_type i = 0; i < len; i++) \
723  { \
724  octave_quit (); \
725  RAND_FUNC (val); \
726  v[i] = val; \
727  } \
728  } \
729  while (0)
730 
731 void
732 octave_rand::fill (octave_idx_type len, double *v, double a)
733 {
734  if (len < 1)
735  return;
736 
737  switch (current_distribution)
738  {
739  case uniform_dist:
740  if (use_old_generators)
741  {
742 #define RAND_FUNC(x) F77_FUNC (dgenunf, DGENUNF) (0.0, 1.0, x)
743  MAKE_RAND (len);
744 #undef RAND_FUNC
745  }
746  else
747  oct_fill_randu (len, v);
748  break;
749 
750  case normal_dist:
751  if (use_old_generators)
752  {
753 #define RAND_FUNC(x) F77_FUNC (dgennor, DGENNOR) (0.0, 1.0, x)
754  MAKE_RAND (len);
755 #undef RAND_FUNC
756  }
757  else
758  oct_fill_randn (len, v);
759  break;
760 
761  case expon_dist:
762  if (use_old_generators)
763  {
764 #define RAND_FUNC(x) F77_FUNC (dgenexp, DGENEXP) (1.0, x)
765  MAKE_RAND (len);
766 #undef RAND_FUNC
767  }
768  else
769  oct_fill_rande (len, v);
770  break;
771 
772  case poisson_dist:
773  if (use_old_generators)
774  {
775  if (a < 0.0 || ! xfinite (a))
776 #define RAND_FUNC(x) x = octave_NaN;
777  MAKE_RAND (len);
778 #undef RAND_FUNC
779  else
780  {
781  // workaround bug in ignpoi, by calling with different Mu
782  double tmp;
783  F77_FUNC (dignpoi, DIGNPOI) (a + 1, tmp);
784 #define RAND_FUNC(x) F77_FUNC (dignpoi, DIGNPOI) (a, x)
785  MAKE_RAND (len);
786 #undef RAND_FUNC
787  }
788  }
789  else
790  oct_fill_randp (a, len, v);
791  break;
792 
793  case gamma_dist:
794  if (use_old_generators)
795  {
796  if (a <= 0.0 || ! xfinite (a))
797 #define RAND_FUNC(x) x = octave_NaN;
798  MAKE_RAND (len);
799 #undef RAND_FUNC
800  else
801 #define RAND_FUNC(x) F77_FUNC (dgengam, DGENGAM) (1.0, a, x)
802  MAKE_RAND (len);
803 #undef RAND_FUNC
804  }
805  else
806  oct_fill_randg (a, len, v);
807  break;
808 
809  default:
810  (*current_liboctave_error_handler)
811  ("rand: invalid distribution ID = %d", current_distribution);
812  break;
813  }
814 
815  save_state ();
816 
817  return;
818 }
819 
820 void
821 octave_rand::fill (octave_idx_type len, float *v, float a)
822 {
823  if (len < 1)
824  return;
825 
826  switch (current_distribution)
827  {
828  case uniform_dist:
829  if (use_old_generators)
830  {
831 #define RAND_FUNC(x) F77_FUNC (dgenunf, DGENUNF) (0.0, 1.0, x)
832  MAKE_RAND (len);
833 #undef RAND_FUNC
834  }
835  else
836  oct_fill_float_randu (len, v);
837  break;
838 
839  case normal_dist:
840  if (use_old_generators)
841  {
842 #define RAND_FUNC(x) F77_FUNC (dgennor, DGENNOR) (0.0, 1.0, x)
843  MAKE_RAND (len);
844 #undef RAND_FUNC
845  }
846  else
847  oct_fill_float_randn (len, v);
848  break;
849 
850  case expon_dist:
851  if (use_old_generators)
852  {
853 #define RAND_FUNC(x) F77_FUNC (dgenexp, DGENEXP) (1.0, x)
854  MAKE_RAND (len);
855 #undef RAND_FUNC
856  }
857  else
858  oct_fill_float_rande (len, v);
859  break;
860 
861  case poisson_dist:
862  if (use_old_generators)
863  {
864  double da = a;
865  if (da < 0.0 || ! xfinite (a))
866 #define RAND_FUNC(x) x = octave_NaN;
867  MAKE_RAND (len);
868 #undef RAND_FUNC
869  else
870  {
871  // workaround bug in ignpoi, by calling with different Mu
872  double tmp;
873  F77_FUNC (dignpoi, DIGNPOI) (da + 1, tmp);
874 #define RAND_FUNC(x) F77_FUNC (dignpoi, DIGNPOI) (da, x)
875  MAKE_RAND (len);
876 #undef RAND_FUNC
877  }
878  }
879  else
880  oct_fill_float_randp (a, len, v);
881  break;
882 
883  case gamma_dist:
884  if (use_old_generators)
885  {
886  double da = a;
887  if (da <= 0.0 || ! xfinite (a))
888 #define RAND_FUNC(x) x = octave_NaN;
889  MAKE_RAND (len);
890 #undef RAND_FUNC
891  else
892 #define RAND_FUNC(x) F77_FUNC (dgengam, DGENGAM) (1.0, da, x)
893  MAKE_RAND (len);
894 #undef RAND_FUNC
895  }
896  else
897  oct_fill_float_randg (a, len, v);
898  break;
899 
900  default:
901  (*current_liboctave_error_handler)
902  ("rand: invalid distribution ID = %d", current_distribution);
903  break;
904  }
905 
906  save_state ();
907 
908  return;
909 }