use strict;
use warnings;
use PDL::Types qw(types ppdefs_all);
pp_bless('PDL::GSL::RNG'); # make the functions generated go into our namespace, and
# not PDL's namespace
pp_addpm({At=>'Top'},<<'EOD');
use strict;
use warnings;
=head1 NAME
PDL::GSL::RNG - PDL interface to RNG and randist routines in GSL
=head1 DESCRIPTION
This is an interface to the rng and randist packages present
in the GNU Scientific Library.
=head1 SYNOPSIS
use PDL;
use PDL::GSL::RNG;
$rng = PDL::GSL::RNG->new('taus');
$rng->set_seed(time());
$x=zeroes(5,5,5)
$rng->get_uniform($x); # inplace
$y=$rng->get_uniform(3,4,5); # creates new pdl
=head1 NOMENCLATURE
Throughout this documentation we strive to use the same variables that
are present in the original GSL documentation (see L<See
Also|"SEE-ALSO">). Oftentimes those variables are called C<a> and
C<b>. Since good Perl coding practices discourage the use of Perl
variables C<$a> and C<$b>, here we refer to Parameters C<a> and C<b>
as C<$pa> and C<$pb>, respectively, and Limits (of domain or
integration) as C<$la> and C<$lb>.
=cut
EOD
pp_addpm({At=>'Middle'},<<'EOD');
=head2 new
=for ref
The new method initializes a new instance of the RNG.
The available RNGs are:
=over
=item coveyou
=item cmrg
=item fishman18
=item fishman20
=item fishman2x
=item gfsr4
=item knuthran
=item knuthran2
=item knuthran2002
=item lecuyer21
=item minstd
=item mrg
=item mt19937
=item mt19937_1999
=item mt19937_1998
=item r250
=item ran0
=item ran1
=item ran2
=item ran3
=item rand
=item rand48
=item random128_bsd
=item random128_glibc2
=item random128_libc5
=item random256_bsd
=item random256_glibc2
=item random256_libc5
=item random32_bsd
=item random32_glibc2
=item random32_libc5
=item random64_bsd
=item random64_glibc2
=item random64_libc5
=item random8_bsd
=item random8_glibc2
=item random8_libc5
=item random_bsd
=item random_glibc2
=item random_libc5
=item randu
=item ranf
=item ranlux
=item ranlux389
=item ranlxd1
=item ranlxd2
=item ranlxs0
=item ranlxs1
=item ranlxs2
=item ranmar
=item slatec
=item taus
=item taus2
=item taus113
=item transputer
=item tt800
=item uni
=item uni32
=item vax
=item waterman14
=item zuf
=item default
=back
The last one (default) uses the environment variable GSL_RNG_TYPE.
Note that only a few of these rngs are recommended for general
use. Please check the GSL documentation for more information.
=for usage
Usage:
$blessed_ref = PDL::GSL::RNG->new($RNG_name);
Example:
=for example
$rng = PDL::GSL::RNG->new('taus');
=head2 set_seed
=for ref
Sets the RNG seed.
Usage:
=for usage
$rng->set_seed($integer);
# or
$rng = PDL::GSL::RNG->new('taus')->set_seed($integer);
Example:
=for example
$rng->set_seed(666);
=head2 min
=for ref
Return the minimum value generable by this RNG.
Usage:
=for usage
$integer = $rng->min();
Example:
=for example
$min = $rng->min(); $max = $rng->max();
=head2 max
=for ref
Return the maximum value generable by the RNG.
Usage:
=for usage
$integer = $rng->max();
Example:
=for example
$min = $rng->min(); $max = $rng->max();
=head2 name
=for ref
Returns the name of the RNG.
Usage:
=for usage
$string = $rng->name();
Example:
=for example
$name = $rng->name();
=head2 ran_shuffle
=for ref
Shuffles values in ndarray, treating it as flat.
Usage:
=for usage
$rng->ran_shuffle($ndarray);
=head2 ran_shuffle_vec
=for ref
Returns values in Perl list, shuffled.
Usage:
=for usage
@shuffled = $rng->ran_shuffle_vec(@vec);
=head2 ran_choose
=for ref
Chooses values from C<$inndarray> to C<$outndarray>, treating both as flat.
Usage:
=for usage
$rng->ran_choose($inndarray,$outndarray);
=head2 ran_choose_vec
=for ref
Chooses C<$n> values from C<@vec>.
Usage:
=for usage
@chosen = $rng->ran_choose_vec($n,@vec);
=head2 ran_dir
=for ref
Returns C<$n> random vectors in C<$ndim> dimensions.
Usage:
=for usage
$ndarray = $rng->ran_dir($ndim,$n);
Example:
=for example
$o = $rng->ran_dir($ndim,$n);
=head2 ran_discrete_preproc
=for ref
This method returns a handle that must be used when calling
L</ran_discrete>. You specify the probability of the integer number
that are returned by L</ran_discrete>.
Usage:
=for usage
$discrete_dist_handle = $rng->ran_discrete_preproc($double_ndarray_prob);
Example:
=for example
$prob = pdl [0.1,0.3,0.6];
$ddh = $rng->ran_discrete_preproc($prob);
$o = $rng->ran_discrete($discrete_dist_handle,100);
=cut
EOD
pp_addpm({At=>'Bot'},<<'EOD');
=head1 BUGS
Feedback is welcome. Log bugs in the PDL bug database (the
database is always linked from L<http://pdl.perl.org/>).
=head1 SEE ALSO
L<PDL>
The GSL documentation for random number distributions is online at
L<https://www.gnu.org/software/gsl/doc/html/randist.html>
=head1 AUTHOR
This file copyright (C) 1999 Christian Pellegrin <chri@infis.univ.trieste.it>
Docs mangled by C. Soeller. All rights reserved. There
is no warranty. You are allowed to redistribute this software /
documentation under certain conditions. For details, see the file
COPYING in the PDL distribution. If this file is separated from the
PDL distribution, the copyright notice should be included in the file.
The GSL RNG and randist modules were written by James Theiler.
=cut
EOD
# PP interface to RNG
##############################
#
# make_get_sub generates a wrapper PDL subroutine that handles the
# fill-a-PDL and create-a-PDL cases for each of the GSL functions.
# --CED
#
sub make_get_sub {
my ($fname,$par) =@_;
my $s = '
sub ' . $fname . ' {
my ($obj,' . $par . '@var) = @_;';
if ($par ne '') {
my $ss=$par;
$ss =~ s/,//;
$s .= 'if (!(' . $ss . '>0)) {barf("first parameter must be an int >0")};';
}
$s .= 'if (ref($var[0]) eq \'PDL\') {
_' . $fname . '_int($var[0],' . $par . '$obj);
return $var[0];
}
else {
my $p;
$p = zeroes @var;
_' . $fname . '_int($p,' . $par . '$obj);
return $p;
}
}
'
}
pp_addhdr('
#include <string.h>
#include "gsl/gsl_rng.h"
#include "gsl/gsl_randist.h"
');
use Config;
my %p_type = (
'double' => 'double',
'unsigned int' => $Config{intsize} == 4 ? 'uint' :
$Config{intsize} == 8 ? 'ulonglong' : die("Unknown intsize, not [48]"),
'size_t' => 'SIZE',
);
use Text::ParseWords qw(shellwords);
chomp(my $header = `gsl-config --cflags`);
$header =~ s#\\#/#g; # win32
($header) = map {s/^-I//;$_} grep /^-I/, shellwords $header;
$header .= '/gsl/gsl_randist.h';
my $h = do { open my $fh, "<", $header or die "$header: $!"; local $/ = undef; <$fh> };
my @defs;
DEF: for my $def ($h =~ m/(double\s+gsl_ran_\S*_pdf\s*\(.+?\))/sg) {
local $_ = $def;
s/\n\s+/ /xmsg;
s/const //g;
if (m/^(\w+)\s+(\w+)\s*\(\s*(.+)\s*\)/s) {
my ($out_type, $function, $pars) = ($1, $2, $3);
my @pars = split /,/, $pars;
for (@pars) {
die "Couldn't parse '$_' in '$def'" if !m/^(.+?)(\w+)\s*(\[\s*\])?\s*$/;
my ($type, $par, $vec) = ($1, $2, $3);
next DEF if $type =~ /gsl_ran_discrete_t/;
s/^ | $//g for ($type, $par);
$par = lc $par;
# numbers interfere with $ISGOOD in BadCode
$par =~ s/1/a/g;
$par =~ s/2/b/g;
die("unknown type '$type' in '$def'") if !$p_type{$type};
$_ = [$p_type{$type}, $par, $vec ? 1 : ()];
}
push @defs, [ $out_type, $function, \@pars ];
}
}
@defs = sort { $a->[1] cmp $b->[1] } @defs; # sort by function name
my @export_names;
for my $f (@defs) {
my ($out_type, $function, $pars) = @$f;
my ($p, $code) = gen_def($out_type, $function, @$pars);
(my $pp_name = $function) =~ s/^gsl_//;
push @export_names, $pp_name;
pp_def($pp_name,
HandleBad => 1,
GenericTypes => ['D'],
Pars => $p,
Code => $code,
Doc => '',
);
}
sub gen_def {
my ($out_type, $function, @pars) = @_;
my $SIZEVAR;
my @formal_pars = map {
if ($_->[0] eq 'SIZE') {
die "$function got second sizevar (@$_) but already got '$SIZEVAR'"
if $SIZEVAR;
$SIZEVAR = $_->[1];
()
} else {
my $prefix = $_->[0] eq 'double' ? '' : "$_->[0] ";
my $suffix = !$_->[2] ? '()' :
$SIZEVAR ? "($SIZEVAR)" :
die "$function got vector par (@$_) but no sizevar";
$prefix.$_->[1].$suffix
}
} @pars;
my $pars = join '; ', @formal_pars, '[o]out()';
my @call_pars = map $_->[0] eq 'SIZE' ? "\$SIZE($_->[1])" :
$_->[2] ? "\$P($_->[1])" : "\$$_->[1]()",
@pars;
my $code = join "\n", "PDL_IF_BAD(char anybad=0;,)",
(map $_->[2] ?
"PDL_IF_BAD(loop ($SIZEVAR) %{ if (\$ISBAD($_->[1]())) { anybad=1; break; } %},)" :
"PDL_IF_BAD(if ( \$ISBAD($_->[1]())) anybad=1;,)",
grep $_->[0] ne 'SIZE', @pars),
"PDL_IF_BAD(if (anybad) { \$SETBAD(out()); continue; },)",
"\$out() = $function(@{[join ', ', @call_pars]});\n";
return ($pars, $code);
}
sub pp_defnd { # hide the docs
my ($name, %hash) = @_;
pp_def($name,%hash,Doc=>undef);
}
pp_def('ran_shuffle_1d',
Pars => '[io]a(n)',
GenericTypes => [ppdefs_all()],
OtherPars => 'gsl_rng *rng',
Code => '
gsl_ran_shuffle($COMP(rng), $P(a), $SIZE(n), sizeof($GENERIC()));
',
PMCode => <<'EOF',
sub ran_shuffle_1d { _ran_shuffle_1d_int(@_[1,0]) }
EOF
Doc => <<'EOF',
=for ref
Takes n-dimensional ndarray, and shuffles it along its zero-th dimension.
Usage:
=for usage
$vec2d = sequence(10,10);
$rng->ran_shuffle_1d($vec2d);
EOF
);
pp_def('get_uniform',
Pars => '[o]a()',
GenericTypes => ['F','D'],
OtherPars => 'gsl_rng *rng',
Code => '
$a() = gsl_rng_uniform($COMP(rng));',
PMCode => make_get_sub('get_uniform',''),
Doc => <<'EOF',
=for ref
This function creates an ndarray with given dimensions or accepts an
existing ndarray and fills it. get_uniform() returns values 0<=x<1,
Usage:
=for usage
$ndarray = $rng->get_uniform($list_of_integers)
$rng->get_uniform($ndarray);
Example:
=for example
$x = zeroes 5,6; $max=100;
$o = $rng->get_uniform(10,10); $rng->get_uniform($x);
EOF
);
pp_def('get_uniform_pos',
Pars => '[o]a()',
GenericTypes => ['F','D'],
OtherPars => 'gsl_rng *rng',
Code => '
$a() = gsl_rng_uniform_pos($COMP(rng));',
PMCode => make_get_sub('get_uniform_pos',''),
Doc => <<'EOF',
=for ref
This function creates an ndarray with given dimensions or accepts an
existing ndarray and fills it. get_uniform_pos() returns values 0<x<1,
Usage:
=for usage
$ndarray = $rng->get_uniform_pos($list_of_integers)
$rng->get_uniform_pos($ndarray);
Example:
=for example
$x = zeroes 5,6;
$o = $rng->get_uniform_pos(10,10); $rng->get_uniform_pos($x);
EOF
);
pp_def('get',
Pars => '[o]a()',
GenericTypes => ['F','D'],
OtherPars => 'gsl_rng *rng',
Code => '
$a() = gsl_rng_get($COMP(rng));',
PMCode => make_get_sub('get',''),
Doc => <<'EOF',
=for ref
This function creates an ndarray with given dimensions or accepts an
existing ndarray and fills it. get() returns integer values
between a minimum and a maximum specific to every RNG.
Usage:
=for usage
$ndarray = $rng->get($list_of_integers)
$rng->get($ndarray);
Example:
=for example
$x = zeroes 5,6;
$o = $rng->get(10,10); $rng->get($x);
EOF
);
pp_def('get_int',
Pars => '[o]a()',
GenericTypes => ['F','D'],
OtherPars => 'IV n; gsl_rng *rng',
Code => '
$a() = gsl_rng_uniform_int($COMP(rng),$COMP(n));',
PMCode => make_get_sub('get_int','$n,'),
Doc => <<'EOF',
=for ref
This function creates an ndarray with given dimensions or accepts an
existing ndarray and fills it. get_int() returns integer values
between 0 and $max.
Usage:
=for usage
$ndarray = $rng->get($max, $list_of_integers)
$rng->get($max, $ndarray);
Example:
=for example
$x = zeroes 5,6; $max=100;
$o = $rng->get(10,10); $rng->get($x);
EOF
);
# randist stuff
sub add_randist {
my ($name,$doc,$params) = @_;
my $npar = @$params;
my $pars1=join '; ', (map "double $_", @$params), 'gsl_rng *rng';
my $fcall1=join ',', map "\$COMP($_)", @$params;
my $arglist=join '', map "\$$_,", @$params;
my $pars2=join ';', map "$_()", @$params;
my $fcall2=join ',', map "\$$_()", @$params;
my $arglist2=join ',', map "\$${_}_ndarray", @$params;
pp_def(
'ran_' . $name,
Pars => '[o]output()',
GenericTypes => ['F','D'],
OtherPars => $pars1,
Code =>'
$output() = gsl_ran_' . $name . '($COMP(rng),' . $fcall1 . ');',
PMCode =>'
sub ran_' . $name . ' {
my ($obj,' . $arglist . '@var) = @_;
if (ref($var[0]) eq \'PDL\') {
_ran_' . $name . '_int($var[0],' . $arglist . '$obj);
return $var[0];
}
else {
my $p;
$p = zeroes @var;
_ran_' . $name . '_int($p,' . $arglist . '$obj);
return $p;
}
}
',
Doc => <<EOF,
=for ref
Fills output ndarray with random $doc
Usage:
=for usage
\$ndarray = \$rng->ran_$name(@{[join ', ', map qq{\$$_}, @$params]},[list of integers = output ndarray dims]);
\$rng->ran_$name(@{[join ', ', map qq{\$$_}, @$params]}, \$output_ndarray);
Example:
=for example
\$o = \$rng->ran_$name(@{[join ', ', map qq{\$$_}, @$params]},10,10);
\$rng->ran_$name(@{[join ', ', map qq{\$$_}, @$params]},\$o);
EOF
);
pp_def(
'ran_' . $name . '_var',
Pars => $pars2 . ';[o]output()',
GenericTypes => ['F','D'],
OtherPars => 'gsl_rng *rng',
Code =>'
$output() = gsl_ran_' . $name . '($COMP(rng),' . $fcall2 . ');',
PMCode =>'
sub ran_' . $name . '_var {
my ($obj,@var) = @_;
if (scalar(@var) != ' . $npar . ') {barf("Bad number of parameters!");}
_ran_' . $name . '_var_int(@var,my $x=PDL->null,$obj);
return $x;
}
',
Doc => <<EOF,
=for ref
Similar to L</ran_$name> except that it takes the distribution
parameters as an ndarray and returns an ndarray of equal dimensions.
Usage:
=for usage
\$ndarray = \$rng->ran_${name}_var($arglist2);
EOF
);
}
add_randist('gaussian','values from Gaussian distribution with mean zero and standard deviation C<$sigma>.',[qw(sigma)]);
add_randist('ugaussian_tail', 'variates from the upper tail of a Gaussian distribution with C<standard deviation = 1> (AKA unit Gaussian distribution).',[qw(tail)]);
add_randist('exponential', 'variates from the exponential distribution with mean C<$mu>.',[qw(mu)]);
add_randist('laplace', 'variates from the Laplace distribution with width C<$pa>.',[qw(pa)]);
add_randist('exppow', 'variates from the exponential power distribution with scale parameter C<$pa> and exponent C<$pb>.',[qw(pa pb)]);
add_randist('cauchy', 'variates from the Cauchy distribution with scale parameter C<$pa>.',[qw(pa)]);
add_randist('rayleigh', 'variates from the Rayleigh distribution with scale parameter C<$sigma>.',[qw(sigma)]);
add_randist('rayleigh_tail', 'variates from the tail of the Rayleigh distribution
with scale parameter C<$sigma> and a lower limit of C<$la>.',[qw(x sigma)]);
add_randist('levy', 'variates from the Levy symmetric stable distribution with scale C<$c> and exponent C<$alpha>.',[qw(mu x)]);
add_randist('gamma', 'variates from the gamma distribution.',[qw(pa pb)]);
add_randist('flat', 'variates from the flat (uniform) distribution from C<$la> to C<$lb>.',[qw(la lb)]);
add_randist('lognormal', 'variates from the lognormal distribution with parameters C<$mu> (location) and C<$sigma> (scale).',[qw(mu sigma)]);
add_randist('chisq', 'variates from the chi-squared distribution with C<$nu> degrees of freedom.',[qw(nu)]);
add_randist('fdist', 'variates from the F-distribution with degrees of freedom C<$nu1> and C<$nu2>.',[qw(nu1 nu2)]);
add_randist('tdist', q{variates from the t-distribution (AKA Student's
t-distribution) with C<$nu> degrees of freedom.},[qw(nu)]);
add_randist('beta', 'variates from the beta distribution with parameters C<$pa> and C<$pb>.',[qw(pa pb)]);
add_randist('logistic', 'random variates from the logistic distribution.',[qw(m)]);
add_randist('pareto', 'variates from the Pareto distribution of order C<$pa> and scale C<$lb>.',[qw(pa lb)]);
add_randist('weibull', 'variates from the Weibull distribution with scale C<$pa> and exponent C<$pb>. (Some literature uses C<lambda> for C<$pa> and C<k> for C<$pb>.)',[qw(pa pb)]);
add_randist('gumbel1', 'variates from the Type-1 Gumbel distribution.',[qw(pa pb)]);
add_randist('gumbel2', 'variates from the Type-2 Gumbel distribution.',[qw(pa pb)]);
add_randist('poisson','integer values from the Poisson distribution with mean C<$mu>.',[qw(mu)]);
add_randist('bernoulli', 'values 0 or 1, the result of a Bernoulli trial with probability C<$p>.',[qw(p)]);
add_randist('binomial', 'integer values from the binomial distribution, the number of successes in C<$n> independent trials with probability C<$p>.',[qw(p n)]);
add_randist('negative_binomial', 'integer values from the negative binomial
distribution, the number of failures occurring before C<$n> successes in
independent trials with probability C<$p> of success. Note that C<$n> is
not required to be an integer.',[qw(p n)]);
add_randist('pascal', 'integer values from the Pascal distribution.
The Pascal distribution is simply a negative binomial distribution
(see L</ran_negative_binomial>) with an integer value of C<$n>.',[qw(p n)]);
add_randist('geometric', 'integer values from the geometric distribution, the number of independent trials with probability C<$p> until the first success.',[qw(p)]);
add_randist('hypergeometric', 'integer values from the hypergeometric distribution.
If a population contains C<$n1> elements of type 1 and C<$n2> elements of
type 2 then the hypergeometric distribution gives the probability of obtaining
C<$x> elements of type 1 in C<$t> samples from the population without replacement.',[qw(n1 n2 t)]);
add_randist('logarithmic', 'integer values from the logarithmic distribution.',[qw(p)]);
# specific randist
pp_def('ran_additive_gaussian',
Pars => '[o]x()',
GenericTypes => ['F','D'],
OtherPars => 'double sigma; gsl_rng *rng',
Code =>'$x() += gsl_ran_gaussian($COMP(rng), $COMP(sigma));',
PMCode =>'
sub ran_additive_gaussian {
my ($obj,$sigma,$var) = @_;
barf("In additive gaussian mode you must specify an ndarray!")
if ref($var) ne \'PDL\';
_ran_additive_gaussian_int($var,$sigma,$obj);
return $var;
}
',
Doc => <<'EOF',
=for ref
Add Gaussian noise of given sigma to an ndarray.
Usage:
=for usage
$rng->ran_additive_gaussian($sigma,$ndarray);
Example:
=for example
$rng->ran_additive_gaussian(1,$image);
EOF
);
pp_def('ran_additive_poisson',
Pars => '[o]x()',
GenericTypes => ['F','D'],
OtherPars => 'double sigma; gsl_rng *rng',
Code =>'$x() += gsl_ran_poisson($COMP(rng), $COMP(sigma));',
PMCode =>'
sub ran_additive_poisson {
my ($obj,$sigma,$var) = @_;
barf("In additive poisson mode you must specify an ndarray!")
if ref($var) ne \'PDL\';
_ran_additive_poisson_int($var,$sigma,$obj);
return $var;
}
',
Doc => <<'EOF',
=for ref
Add Poisson noise of given C<$mu> to a C<$ndarray>.
Usage:
=for usage
$rng->ran_additive_poisson($mu,$ndarray);
Example:
=for example
$rng->ran_additive_poisson(1,$image);
EOF
);
pp_def('ran_feed_poisson',
Pars => '[o]x()',
GenericTypes => ['F','D'],
OtherPars => 'gsl_rng *rng',
Code =>'$x() = gsl_ran_poisson($COMP(rng), $x());',
PMCode =>'
sub ran_feed_poisson {
my ($obj,$var) = @_;
barf("In poisson mode you must specify an ndarray!")
if ref($var) ne \'PDL\';
_ran_feed_poisson_int($var,$obj);
return $var;
}
',
Doc => <<'EOF',
=for ref
This method simulates shot noise, taking the values of ndarray as
values for C<$mu> to be fed in the poissonian RNG.
Usage:
=for usage
$rng->ran_feed_poisson($ndarray);
Example:
=for example
$rng->ran_feed_poisson($image);
EOF
);
pp_def('ran_bivariate_gaussian',
Pars => '[o]x(n)',
GenericTypes => ['F','D'],
OtherPars => 'double sigma_x; double sigma_y; double rho; gsl_rng *rng',
Code => <<'EOF',
double xx,yy;
gsl_ran_bivariate_gaussian($COMP(rng), $COMP(sigma_x), $COMP(sigma_y),$COMP(rho), &xx, &yy);
$x(n=>0)=xx;
$x(n=>1)=yy;
EOF
PMCode => <<'EOF',
sub ran_bivariate_gaussian {
my ($obj,$sigma_x,$sigma_y,$rho,$n) = @_;
barf("Not enough parameters for gaussian bivariate!") if $n<=0;
my $p = zeroes(2,$n);
_ran_bivariate_gaussian_int($p,$sigma_x,$sigma_y,$rho,$obj);
return $p;
}
EOF
Doc => <<'EOF',
=for ref
Generates C<$n> bivariate gaussian random deviates.
Usage:
=for usage
$ndarray = $rng->ran_bivariate_gaussian($sigma_x,$sigma_y,$rho,$n);
Example:
=for example
$o = $rng->ran_bivariate_gaussian(1,2,0.5,1000);
EOF
);
pp_defnd('ran_dir_2d',
Pars => '[o]x(n)',
GenericTypes => ['F','D'],
OtherPars => 'gsl_rng *rng',
Code =>'
double xx,yy;
gsl_ran_dir_2d($COMP(rng), &xx, &yy);
$x(n=>0)=xx;
$x(n=>1)=yy;
');
pp_defnd('ran_dir_3d',
Pars => '[o]x(n)',
GenericTypes => ['F','D'],
OtherPars => 'gsl_rng *rng',
Code =>'
double xx,yy,zz;
gsl_ran_dir_3d($COMP(rng), &xx, &yy, &zz);
$x(n=>0)=xx;
$x(n=>1)=yy;
$x(n=>2)=zz;
');
my $MAX_DIMENSIONS = 100;
pp_defnd('ran_dir_nd',
Pars => '[o]x(n)',
GenericTypes => ['F','D'],
OtherPars => 'IV ns => n; gsl_rng *rng',
Code =>'
double xxx[' . $MAX_DIMENSIONS .'];
gsl_ran_dir_nd($COMP(rng), $COMP(ns), xxx);
loop (n) %{ $x() = xxx[n]; %}');
pp_addpm('
sub ran_dir {
my ($obj,$ndim,$n) = @_;
barf("Not enough parameters for random vectors!") if $n<=0;
my $p = zeroes($ndim,$n);
if ($ndim==2) { ran_dir_2d($p,$obj); }
elsif ($ndim==3) { ran_dir_3d($p,$obj); }
elsif ($ndim>=4 && $ndim<=' . $MAX_DIMENSIONS . ') { ran_dir_nd($p,$ndim,$obj); }
else { barf("Bad number of dimensions!"); }
return $p;
}
');
pp_def('ran_discrete',
Pars => '[o]x()',
GenericTypes => ['F','D'],
OtherPars => 'gsl_ran_discrete_t *rng_discrete; gsl_rng *rng',
Code =>'
$x()=gsl_ran_discrete($COMP(rng), $COMP(rng_discrete)); ',
PMCode =>'
sub ran_discrete {
my ($obj, $rdt, @var) = @_;
if (ref($var[0]) eq \'PDL\') {
_ran_discrete_int($var[0], $rdt, $obj);
return $var[0];
}
else {
my $p;
$p = zeroes @var;
_ran_discrete_int($p, $rdt, $obj);
return $p;
}
}
',
Doc => <<'EOF',
=for ref
Is used to get the desired samples once a proper handle has been
enstablished (see ran_discrete_preproc()).
Usage:
=for usage
$ndarray = $rng->ran_discrete($discrete_dist_handle,$num);
Example:
=for example
$prob = pdl [0.1,0.3,0.6];
$ddh = $rng->ran_discrete_preproc($prob);
$o = $rng->ran_discrete($discrete_dist_handle,100);
EOF
);
pp_addpm('
sub ran_shuffle_vec {
my ($obj,@in) = @_;
$obj->ran_shuffle(my $p = PDL->sequence(PDL::indx(), 0+@in));
@in[$p->list];
}
');
pp_addpm('
sub ran_choose_vec {
my ($obj,$nout,@in) = @_;
$obj->ran_choose(PDL->sequence(PDL::indx(), 0+@in),my $pout = PDL->zeroes(PDL::indx(), $nout));
@in[$pout->list];
}
');
pp_def('ran_ver',
Pars => '[o]x(n)',
GenericTypes => ['F','D'],
OtherPars => 'double x0; double r;IV ns => n; gsl_rng *rng',
Code =>'
double xx=$COMP(x0);
loop (n) %{ $x() = xx; xx = $COMP(r)*(1-xx)*xx; %}',
PMCode =>'
sub ran_ver {
my ($obj,$x0,$r,$n) = @_;
barf("Not enough parameters for ran_ver!") if $n<=0;
my $p = zeroes($n);
_ran_ver_int($p,$x0,$r,$n,$obj);
return $p;
}
',
Doc => <<'EOF',
=for ref
Returns an ndarray with C<$n> values generated by the Verhulst map from C<$x0> and
parameter C<$r>.
Usage:
=for usage
$rng->ran_ver($x0, $r, $n);
EOF
);
pp_def('ran_caos',
Pars => '[o]x(n)',
GenericTypes => ['F','D'],
OtherPars => 'double m; IV ns => n; gsl_rng *rng',
Code =>'
double xx=gsl_ran_gaussian($COMP(rng),0.1)+0.5;
loop (n) %{ $x() = (xx-0.5)*$COMP(m); xx = 4.0*(1-xx)*xx; %}',
PMCode =>'
sub ran_caos {
my ($obj,$m,$n) = @_;
barf("Not enough parameters for ran_caos!") if $n<=0;
my $p = zeroes($n);
_ran_caos_int($p,$m,$n,$obj);
return $p;
}
',
Doc => <<'EOF',
=for ref
Returns values from Verhuls map with C<$r=4.0> and randomly chosen
C<$x0>. The values are scaled by C<$m>.
Usage:
=for usage
$rng->ran_caos($m,$n);
EOF
);
# XS function for the RNG object
pp_addxs('','
MODULE = PDL::GSL::RNG PACKAGE = PDL::GSL::RNG
#define DEF_RNG(X) if (!strcmp(TYPE,#X)) rng=gsl_rng_alloc( gsl_rng_ ## X ); strcat(rngs,#X ", ");
gsl_rng *
new (CLASS,TYPE)
char *CLASS
char *TYPE
CODE:
gsl_rng * rng = NULL;
char rngs[5000];
strcpy(rngs,"");
DEF_RNG(borosh13)
DEF_RNG(coveyou)
DEF_RNG(cmrg)
DEF_RNG(fishman18)
DEF_RNG(fishman20)
DEF_RNG(fishman2x)
DEF_RNG(gfsr4)
DEF_RNG(knuthran)
DEF_RNG(knuthran2)
DEF_RNG(knuthran2002)
DEF_RNG(lecuyer21)
DEF_RNG(minstd)
DEF_RNG(mrg)
DEF_RNG(mt19937)
DEF_RNG(mt19937_1999)
DEF_RNG(mt19937_1998)
DEF_RNG(r250)
DEF_RNG(ran0)
DEF_RNG(ran1)
DEF_RNG(ran2)
DEF_RNG(ran3)
DEF_RNG(rand)
DEF_RNG(rand48)
DEF_RNG(random128_bsd)
DEF_RNG(random128_glibc2)
DEF_RNG(random128_libc5)
DEF_RNG(random256_bsd)
DEF_RNG(random256_glibc2)
DEF_RNG(random256_libc5)
DEF_RNG(random32_bsd)
DEF_RNG(random32_glibc2)
DEF_RNG(random32_libc5)
DEF_RNG(random64_bsd)
DEF_RNG(random64_glibc2)
DEF_RNG(random64_libc5)
DEF_RNG(random8_bsd)
DEF_RNG(random8_glibc2)
DEF_RNG(random8_libc5)
DEF_RNG(random_bsd)
DEF_RNG(random_glibc2)
DEF_RNG(random_libc5)
DEF_RNG(randu)
DEF_RNG(ranf)
DEF_RNG(ranlux)
DEF_RNG(ranlux389)
DEF_RNG(ranlxd1)
DEF_RNG(ranlxd2)
DEF_RNG(ranlxs0)
DEF_RNG(ranlxs1)
DEF_RNG(ranlxs2)
DEF_RNG(ranmar)
DEF_RNG(slatec)
DEF_RNG(taus)
DEF_RNG(taus2)
DEF_RNG(taus113)
DEF_RNG(transputer)
DEF_RNG(tt800)
DEF_RNG(uni)
DEF_RNG(uni32)
DEF_RNG(vax)
DEF_RNG(waterman14)
DEF_RNG(zuf)
DEF_RNG(default)
if (rng==NULL) {
barf("Unknown RNG, please use one of the following: %s", rngs);
}
else
RETVAL = rng;
OUTPUT:
RETVAL
void
set_seed(rng, seed)
gsl_rng * rng
int seed
PPCODE:
gsl_rng_set(rng,seed);
XPUSHs(ST(0)); /* return self */
unsigned int
min(rng)
gsl_rng * rng
CODE:
RETVAL = gsl_rng_min(rng);
OUTPUT:
RETVAL
unsigned int
max(rng)
gsl_rng * rng
CODE:
RETVAL = gsl_rng_max(rng);
OUTPUT:
RETVAL
char*
name(rng)
gsl_rng * rng
CODE:
RETVAL = (char *) gsl_rng_name(rng);
OUTPUT:
RETVAL
void
DESTROY(sv)
SV * sv
CODE:
gsl_rng *rng = INT2PTR(gsl_rng *, SvIV((SV*)SvRV(sv)));
gsl_rng_free(rng);
gsl_ran_discrete_t *
ran_discrete_preproc(rng, p)
gsl_rng * rng
pdl * p
CODE:
IV n;
if (p->ndims!=1 || p->datatype!=PDL_D) {
barf("Bad input to ran_discrete_preproc!");
}
PDL->barf_if_error(PDL->make_physical(p));
n = p->dims[0];
RETVAL = gsl_ran_discrete_preproc(n,(double *) p->data);
OUTPUT:
RETVAL
void
ran_shuffle(rng, in)
gsl_rng * rng
pdl * in
CODE:
IV size, n;
PDL->barf_if_error(PDL->make_physical(in));
n = in->nvals;
size = PDL->howbig(in->datatype);
gsl_ran_shuffle(rng,in->data,n,size);
PDL->changed(in, PDL_PARENTDATACHANGED, 0);
void
ran_choose(rng, in, out)
gsl_rng * rng
pdl * in
pdl * out
CODE:
IV size, n,m;
if (in->datatype != out->datatype) barf("Data Types must match for ran_chooser");
PDL->barf_if_error(PDL->make_physical(in));
PDL->barf_if_error(PDL->make_physical(out));
n = in->nvals;
m = out->nvals;
size = PDL->howbig(in->datatype);
gsl_ran_choose(rng,out->data, m, in->data,n,size);
PDL->changed(out, PDL_PARENTDATACHANGED, 0);
');
pp_core_importList(' qw/ zeroes long barf /'); # import just a named list to our namespace, so we don't get warning
# messages like 'warning 'min' redefined at line ...'
pp_export_nothing; # set to not export anything. (This is a OO package, it doesn't need to export any methods.)
pp_add_exported(@export_names); # ... except functions suitable for that
pp_add_boot('gsl_set_error_handler_off();
');
pp_done();