pp_bless('PDL::GSL::RNG'); # make the functions generated go into our namespace, and
# not PDL's namespace
pp_addpm({At=>Top},<<'EOD');
=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());
$a=zeroes(5,5,5)
$rng->get_uniform($a); # inplace
$b=$rng->get_uniform(3,4,5); # creates new pdl
=head1 FUNCTIONS
=head2 new()
=for ref
The new method initializes a new instance of the RNG.
The avaible RNGs are:
slatec, cmrg, gfsr4, minstd, mrg, mt19937, r250, ran0, ran1, ran2,
ran3, rand48, rand, random8_bsd, random8_glibc2, random8_libc5,
random128_bsd, random128_glibc2, random128_libc5, random256_bsd,
random256_glibc2, random256_libc5, random32_bsd, random32_glibc2,
random32_libc5, random64_bsd, random64_glibc2, random64_libc5,
random_bsd, random_glibc2, random_libc5, randu, ranf, ranlux389,
ranlux, ranmar, taus, transputer, tt800, uni32, uni, vax, zuf,
default. The last one (default) uses the enviroment variable
GSL_RNG_TYPE. 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);
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 get_uniform()
=for ref
This function creates a piddle with given dimensions or accept an
existing piddle and fills it. get_uniform() returns values 0<=x<1,
Usage:
=for usage
$piddle = $rng->get_uniform($list_of_integers)
$rng->get_uniform($piddle);
Example:
=for example
$a = zeroes 5,6; $max=100;
$o = $rng->get_uniform(10,10); $rng->get_uniform($a);
=head2 get_uniform_pos()
=for ref
This function creates a piddle with given dimensions or accept an
existing piddle and fills it. get_uniform_pos() returns values 0<x<1,
Usage:
=for usage
$piddle = $rng->get_uniform_pos($list_of_integers)
$rng->get_uniform_pos($piddle);
Example:
=for example
$a = zeroes 5,6;
$o = $rng->get_uniform_pos(10,10); $rng->get_uniform_pos($a);
=head2 get()
=for ref
This function creates a piddle with given dimensions or accept an
existing piddle and fills it. get() returns integer values
beetween a minimum and a maximum specific to evry RNG.
Usage:
=for usage
$piddle = $rng->get($list_of_integers)
$rng->get($piddle);
Example:
=for example
$a = zeroes 5,6;
$o = $rng->get(10,10); $rng->get($a);
=head2 get_int()
=for ref
This function creates a piddle with given dimensions or accept an
existing piddle and fills it. get_int() returns integer values
beetween 0 and $max.
Usage:
=for usage
$piddle = $rng->get($max, $list_of_integers)
$rng->get($max, $piddle);
Example:
=for example
$a = zeroes 5,6; $max=100;
$o = $rng->get(10,10); $rng->get($a);
=head2 ran_gaussian()
=for ref
These functions return random deviates from given distribution.
The general form is
ran_[distrib](args)
where distrib can be any of the ones shown below.
They accept the parameters of the distribution and a specification of
where to put output. This spec can be in form of list of integers that
specify the dimensions of the ouput piddle or an existing piddle that
will be filled with values inplace.
Usage:
=for usage
# gaussian dist
$piddle = $rng->ran_gaussian($sigma,[list of integers]);
$rng->ran_gaussian($sigma,$piddle);
# gaussian tail
$piddle = $rng->ran_ugaussian_tail($tail,[list of integers]);
$rng->ran_ugaussian_tail($tail,$piddle);
# exponential dist
$piddle = $rng->ran_exponential($mu,[list of integers]);
$rng->ran_exponential($mu,$piddle);
# laplacian dist
$piddle = $rng->ran_laplace($mu,[list of integers]);
$rng->ran_laplace($mu,$piddle);
$piddle = $rng->ran_exppow($mu,$a,[list of integers]);
$rng->ran_exppow($mu,$a,$piddle);
$piddle = $rng->ran_cauchy($mu,[list of integers]);
$rng->ran_cauchy($mu,$piddle);
$piddle = $rng->ran_rayleigh($sigma,[list of integers]);
$rng->ran_rayleigh($sigma,$piddle);
$piddle = $rng->ran_rayleigh_tail($a,$sigma,[list of integers]);
$rng->ran_rayleigh_tail($a,$sigma,$piddle);
$piddle = $rng->ran_levy($mu,$a,[list of integers]);
$rng->ran_levy($mu,$a,$piddle);
$piddle = $rng->ran_gamma($a,$b,[list of integers]);
$rng->ran_gamma($a,$b,$piddle);
$piddle = $rng->ran_flat($a,$b,[list of integers]);
$rng->ran_flat($a,$b,$piddle);
$piddle = $rng->ran_lognormal($zeta, $sigma,[list of integers]);
$rng->ran_lognormal($zeta, $sigma,$piddle);
$piddle = $rng->ran_chisq($nu,[list of integers]);
$rng->ran_chisq($nu,$piddle);
$piddle = $rng->ran_fdist($nu1, $nu2,[list of integers]);
$rng->ran_fdist($nu1, $nu2,$piddle);
$piddle = $rng->ran_tdist($nu,[list of integers]);
$rng->ran_tdist($nu,$piddle);
$piddle = $rng->ran_beta($a,$b,[list of integers]);
$rng->ran_beta($a,$b,$piddle);
$piddle = $rng->ran_logistic($m,[list of integers]u)
$rng->ran_logistic($m,$piddleu)
$piddle = $rng->ran_pareto($a,$b,[list of integers]);
$rng->ran_pareto($a,$b,$piddle);
$piddle = $rng->ran_weibull($mu,$a,[list of integers]);
$rng->ran_weibull($mu,$a,$piddle);
$piddle = $rng->ran_gumbel1($a,$b,[list of integers]);
$rng->ran_gumbel1($a,$b,$piddle);
$piddle = $rng->ran_gumbel2($a,$b,[list of integers]);
$rng->ran_gumbel2($a,$b,$piddle);
$piddle = $rng->ran_poisson($mu,[list of integers]);
$rng->ran_poisson($mu,$piddle);
$piddle = $rng->ran_bernoulli($p,[list of integers]);
$rng->ran_bernoulli($p,$piddle);
$piddle = $rng->ran_binomial($p,$n,[list of integers]);
$rng->ran_binomial($p,$n,$piddle);
$piddle = $rng->ran_negative_binomial($p,$n,[list of integers]);
$rng->ran_negative_binomial($p,$n,$piddle);
$piddle = $rng->ran_pascal($p,$n,[list of integers]);
$rng->ran_pascal($p,$n,$piddle);
$piddle = $rng->ran_geometric($p,[list of integers]);
$rng->ran_geometric($p,$piddle);
$piddle = $rng->ran_hypergeometric($n1, $n2, $t,[list of integers]);
$rng->ran_hypergeometric($n1, $n2, $t,$piddle);
$piddle = $rng->ran_logarithmic($p,[list of integers]);
$rng->ran_logarithmic($p,$piddle);
Example:
=for example
$o = $rng->ran_gaussian($sigma,10,10);
$rng->ran_gaussian($sigma,$a);
=head2 ran_gaussian_var()
=for ref
This method is similar to L<ran_[distrib]()|ran_gaussian()> except that it
takes the
parameters of the distribution as a piddle and returns a piddle of equal
dimensions. Of course, you can use the same set of distributions as in
the previous method (see also the L<ran_gaussian|ran_gaussian()> entry above).
Usage:
=for usage
$piddle = $rng->ran_[distribution]($distr_parameters_list,$piddle_dim_list);
$rng->ran_[distribution]($distr_parameters_list,$piddle);
Example:
=for example
$sigma_pdl = rvals zeroes 11,11;
$o = $rng->ran_gaussian_var($sigma_pdl);
=head2 ran_additive_gaussian()
=for ref
Add Gaussian noise of given sigma to a piddle.
Usage:
=for usage
$rng->ran_additive_gaussian($sigma,$piddle);
Example:
=for example
$rng->ran_additive_gaussian(1,$image);
=head2 ran_additive_poisson()
=for ref
Add Poisson noise of given sigma to a piddle.
Usage:
=for usage
$rng->ran_additive_poisson($mu,$piddle);
Example:
=for example
$rng->ran_additive_poisson(1,$image);
=head2 ran_feed_poisson()
=for ref
This method simulates shot noise, taking the values of piddle as
values for mu to be fed in the poissonian RNG.
Usage:
=for usage
$rng->ran_feed_poisson($piddle);
Example:
=for example
$rng->ran_feed_poisson($image);
=head2 ran_bivariate_gaussian()
=for ref
Generates $n bivariate gaussian random deviates.
Usage:
=for usage
$piddle = $rng->ran_bivariate_gaussian($sigma_x,$sigma_y,$rho,$n);
Example:
=for example
$o = $rng->ran_bivariate_gaussian(1,2,0.5,1000);
=head2 ran_dir()
=for ref
Returns C<$n> random vectors in C<$ndim> dimensions.
Usage:
=for usage
$piddle = $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
ran_discrete(). You specify the probability of the integer number
that are returned by ran_discrete().
Usage:
=for usage
$discrete_dist_handle = $rng->ran_discrete_preproc($double_piddle_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);
=head2 ran_discrete()
=for ref
Is used to get the desired samples once a proper handle has been
enstablished (see ran_discrete_preproc()).
Usage:
=for usage
$piddle = $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);
# =head1 Shuffling and choosing.
=head2 ran_shuffle()
=for ref
Shuffles values in piddle
Usage:
=for usage
$rng->ran_shuffle($piddle);
=head2 ran_shuffle_vec()
=for ref
Shuffles values in piddle
Usage:
=for usage
$rng->ran_shuffle_vec(@vec);
=head2 ran_choose_vec()
=for ref
Chooses values from $inpiddle to $outpiddle.
Usage:
=for usage
$rng->ran_choose($inpiddle,$outpiddle);
=head2 ran_choose_vec()
=for ref
Chooses $n values from @vec.
Usage:
=for usage
@choosen = $rng->ran_choose_vec($n,@vec);
# =head1 Caotic maps.
=head2 ran_ver()
=for ref
Returns a piddle with $n values generated by the Verhulst map from $x0 and
paramater $r.
Usage:
=for usage
$rng->ran_ver($x0, $r, $n);
=head2 ran_caos()
=for ref
Returns values from Verhuls map with $r=4.0 and randomly choosen
$x0. The values are scaled by $m.
Usage:
=for usage
$rng->ran_caos($m,$n);
=head1 BUGS
Feedback is welcome. Log bugs in the PDL bug database (the
database is always linked from http://pdl.perl.org).
=head1 SEE ALSO
L<PDL>
The GSL documentation is online at
http://sources.redhat.com/gsl/ref/gsl-ref_toc.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;
$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\') {
gsl_' . $fname . '_meat($var[0],' . $par . '$$obj);
return $var[0];
}
else {
my $p;
$p = zeroes @var;
gsl_' . $fname . '_meat($p,' . $par . '$$obj);
return $p;
}
}
'
}
pp_addpm(<<'EOPM');
use strict;
# PDL::GSL::RNG::nullcreate just creates a null PDL. Used
# for the GSL functions that create PDLs
sub nullcreate{
my ($type,$arg) = @_;
PDL->nullcreate($arg);
}
EOPM
pp_addpm(make_get_sub('get_uniform',''));
pp_addpm(make_get_sub('get_uniform_pos',''));
pp_addpm(make_get_sub('get',''));
pp_addpm(make_get_sub('get_int','$n,'));
pp_addhdr('
#include <string.h>
#include "gsl/gsl_rng.h"
#include "gsl/gsl_randist.h"
');
sub pp_defnd { # hide the docs
my ($name, %hash) = @_;
pp_def($name,%hash,Doc=>undef);
}
pp_defnd(
'gsl_get_uniform_meat',
Pars => '[o]a()',
GenericTypes => [F,D],
OtherPars => 'int rng',
Code => '
$a() = gsl_rng_uniform((gsl_rng *) $COMP(rng));');
pp_defnd(
'gsl_get_uniform_pos_meat',
Pars => '[o]a()',
GenericTypes => [F,D],
OtherPars => 'int rng',
Code => '
$a() = gsl_rng_uniform_pos((gsl_rng *) $COMP(rng));');
pp_defnd(
'gsl_get_meat',
Pars => '[o]a()',
OtherPars => 'int rng',
Code => '
$a() = gsl_rng_get((gsl_rng *) $COMP(rng));');
pp_defnd(
'gsl_get_int_meat',
Pars => '[o]a()',
OtherPars => 'int n; int rng',
Code => '
$a() = gsl_rng_uniform_int((gsl_rng *) $COMP(rng),$COMP(n));');
# randist stuff
sub add_randist {
my ($name,$npar) = @_;
my ($pars1,$fcall1,$arglist);
if ($npar==1) {
$pars1='double a; int rng';
$fcall1='$COMP(a)';
$arglist='$a,';
$pars2='a()';
$fcall2='$a()';
}
if ($npar==2) {
$pars1='double a; double b; int rng';
$fcall1='$COMP(a),$COMP(b)';
$arglist='$a,$b,';
$pars2='a();b()';
$fcall2='$a(),$b()';
}
if ($npar==3) {
$pars1='double a; double b; double c; int rng';
$fcall1='$COMP(a),$COMP(b),$COMP(c)';
$arglist='$a,$b,$c,';
$pars2='a();b();c()';
$fcall2='$a(),$b(),$c()';
}
pp_defnd(
'ran_' . $name . '_meat',
Pars => '[o]x()',
OtherPars => $pars1,
Code =>'
$x() = gsl_ran_' . $name . '((gsl_rng *) $COMP(rng),' . $fcall1 . ');');
pp_addpm('
sub ran_' . $name . ' {
my ($obj,' . $arglist . '@var) = @_;
if (ref($var[0]) eq \'PDL\') {
ran_' . $name . '_meat($var[0],' . $arglist . '$$obj);
return $var[0];
}
else {
my $p;
$p = zeroes @var;
ran_' . $name . '_meat($p,' . $arglist . '$$obj);
return $p;
}
}
');
pp_defnd(
'ran_' . $name . '_var_meat',
Pars => $pars2 . ';[o]x()',
OtherPars => 'int rng',
Code =>'
$x() = gsl_ran_' . $name . '((gsl_rng *) $COMP(rng),' . $fcall2 . ');');
pp_addpm('
sub ran_' . $name . '_var {
my ($obj,@var) = @_;
if (scalar(@var) != ' . $npar . ') {barf("Bad number of parameters!");}
return ran_' . $name . '_var_meat(@var,$$obj);
}
');
# pp_defnd(
# 'ran_' . $name . '_add_meat',
# Pars => '[o]x()',
# OtherPars => $pars1,
# Code =>'
#$x() += gsl_ran_' . $name . '((gsl_rng *) $COMP(rng),' . $fcall1 . ');');
# pp_addpm('
#sub ran_' . $name . '_add {
#my ($obj,' . $arglist . '@var) = @_;
#if (ref($var[0]) eq \'PDL\') {
# PDL::ran_' . $name . '_add_meat($var[0],' . $arglist . '$$obj);
# return $var[0];
#}
#else {
# barf("In add mode you must specify a piddle!");
#}
#}
#');
# if ($npar==1) {
# pp_defnd(
# 'ran_' . $name . '_feed_meat',
# Pars => '[o]x()',
# OtherPars => 'int rng',
# Code =>'
#$x() = gsl_ran_' . $name . '((gsl_rng *) $COMP(rng), $x());');
# pp_addpm('
#sub ran_' . $name . '_feed {
#my ($obj, @var) = @_;
#if (ref($var[0]) eq \'PDL\') {
# PDL::ran_' . $name . '_feed_meat($var[0], $$obj);
# return $var[0];
#}
#else {
# barf("In feed mode you must specify a piddle!");
#}
#}
#');
# }
}
add_randist('gaussian',1);
add_randist('ugaussian_tail',1);
add_randist('exponential',1);
add_randist('laplace',1);
add_randist('exppow',2);
add_randist('cauchy',1);
add_randist('rayleigh',1);
add_randist('rayleigh_tail',2);
add_randist('levy',2);
add_randist('gamma',2);
add_randist('flat',2);
add_randist('lognormal',2);
add_randist('chisq',1);
add_randist('fdist',2);
add_randist('tdist',1);
add_randist('beta',2);
add_randist('logistic',1);
add_randist('pareto',2);
add_randist('weibull',2);
add_randist('gumbel1',2);
add_randist('gumbel2',2);
add_randist('poisson',1);
add_randist('bernoulli',1);
add_randist('binomial',2);
add_randist('negative_binomial',2);
add_randist('pascal',2);
add_randist('geometric',1);
add_randist('hypergeometric',3);
add_randist('logarithmic',1);
# specific rnadist
pp_defnd(
'ran_additive_gaussian_meat',
Pars => ';[o]x()',
OtherPars => 'double sigma; int rng',
Code =>'$x() += gsl_ran_gaussian((gsl_rng *) $COMP(rng), $COMP(sigma));');
pp_addpm('
sub ran_additive_gaussian {
my ($obj,$sigma,@var) = @_;
if (ref($var[0]) eq \'PDL\') {
ran_additive_gaussian_meat($var[0],$sigma,$$obj);
return $var[0];
}
else {
barf("In additive gaussian mode you must specify a piddle!");
}
}
');
pp_defnd(
'ran_additive_poisson_meat',
Pars => ';[o]x()',
OtherPars => 'double sigma; int rng',
Code =>'$x() += gsl_ran_poisson((gsl_rng *) $COMP(rng), $COMP(sigma));');
pp_addpm('
sub ran_additive_poisson {
my ($obj,$sigma,@var) = @_;
if (ref($var[0]) eq \'PDL\') {
ran_additive_poisson_meat($var[0],$sigma,$$obj);
return $var[0];
}
else {
barf("In additive poisson mode you must specify a piddle!");
}
}
');
pp_defnd(
'ran_feed_poisson_meat',
Pars => ';[o]x()',
OtherPars => 'int rng',
Code =>'$x() = gsl_ran_poisson((gsl_rng *) $COMP(rng), $x());');
pp_addpm('
sub ran_feed_poisson {
my ($obj,@var) = @_;
if (ref($var[0]) eq \'PDL\') {
ran_feed_poisson_meat($var[0],$$obj);
return $var[0];
}
else {
barf("In poisson mode you must specify a piddle!");
}
}
');
pp_defnd(
'ran_bivariate_gaussian_meat',
Pars => ';[o]x(n)',
OtherPars => 'double sigma_x; double sigma_y; double rho; int rng',
Code =>'
double xx,yy;
gsl_ran_bivariate_gaussian((gsl_rng *) $COMP(rng), $COMP(sigma_x), $COMP(sigma_y),$COMP(rho), &xx, &yy);
$x(n=>0)=xx;
$x(n=>1)=yy;
');
pp_addpm('
sub ran_bivariate_gaussian {
my ($obj,$sigma_x,$sigma_y,$rho,$n) = @_;
if ($n>0) {
my $p = zeroes(2,$n);
ran_bivariate_gaussian_meat($p,$sigma_x,$sigma_y,$rho,$$obj);
return $p;
}
else {
barf("Not enough parameters for gaussian bivariate!");
}
}
');
pp_defnd(
'ran_dir_2d_meat',
Pars => ';[o]x(n)',
OtherPars => 'int rng',
Code =>'
double xx,yy;
gsl_ran_dir_2d((gsl_rng *) $COMP(rng), &xx, &yy);
$x(n=>0)=xx;
$x(n=>1)=yy;
');
pp_defnd(
'ran_dir_3d_meat',
Pars => ';[o]x(n)',
OtherPars => 'int rng',
Code =>'
double xx,yy,zz;
gsl_ran_dir_3d((gsl_rng *) $COMP(rng), &xx, &yy, &zz);
$x(n=>0)=xx;
$x(n=>1)=yy;
$x(n=>2)=zz;
');
$MAX_DIMENSIONS = 100;
pp_defnd(
'ran_dir_nd_meat',
Pars => ';[o]x(n)',
OtherPars => 'int ns => n; int rng',
Code =>'
double xxx[' . $MAX_DIMENSIONS .'];
gsl_ran_dir_nd((gsl_rng *) $COMP(rng), $COMP(ns), xxx);
loop (n) %{ $x() = xxx[n]; %}');
pp_addpm('
sub ran_dir {
my ($obj,$ndim,$n) = @_;
if ($n>0) {
my $p = zeroes($ndim,$n);
if ($ndim==2) { ran_dir_2d_meat($p,$$obj); }
elsif ($ndim==3) { ran_dir_3d_meat($p,$$obj); }
elsif ($ndim>=4 && $ndim<=' . $MAX_DIMENSIONS . ') { ran_dir_nd_meat($p,$ndim,$$obj); }
else { barf("Bad number of dimensions!"); }
return $p;
}
else {
barf("Not enough parameters for random vectors!");
}
}
');
pp_defnd(
'ran_discrete_meat',
Pars => ';[o]x()',
OtherPars => 'int rng_discrete; int rng',
Code =>'
$x()=gsl_ran_discrete((gsl_rng *) $COMP(rng), (gsl_ran_discrete_t *) $COMP(rng_discrete)); ');
pp_addpm('
sub ran_discrete {
my ($obj, $rdt, @var) = @_;
if (ref($var[0]) eq \'PDL\') {
ran_discrete_meat($var[0], $$rdt, $$obj);
return $var[0];
}
else {
my $p;
$p = zeroes @var;
ran_discrete_meat($p, $$rdt, $$obj);
return $p;
}
}
');
pp_addpm('
sub ran_shuffle_vec {
my ($obj,@in) = @_;
my (@out,$i,$p);
$p = long [0..$#in];
$obj->ran_shuffle($p);
for($i=0;$i<scalar(@in);$i++) {
$out[$p->at($i)]=$in[$i];
}
return @out;
}
');
pp_addpm('
sub ran_choose_vec {
my ($obj,$nout,@in) = @_;
my (@out,$i,$pin,$pout);
$pin = long [0..$#in];
$pout = long [0..($nout-1)];
$obj->ran_choose($pin,$pout);
for($i=0;$i<$nout;$i++) {
$out[$i]=$in[$pout->at($i)];
}
return @out;
}
');
pp_defnd(
'ran_ver_meat',
Pars => ';[o]x(n)',
OtherPars => 'double x0; double r;int ns => n; int rng',
Code =>'
double xx=$COMP(x0);
loop (n) %{ $x() = xx; xx = $COMP(r)*(1-xx)*xx; %}');
pp_defnd(
'ran_caos_meat',
Pars => ';[o]x(n)',
OtherPars => 'double m; int ns => n; int rng',
Code =>'
double xx=gsl_ran_gaussian((gsl_rng *) $COMP(rng),0.1)+0.5;
loop (n) %{ $x() = (xx-0.5)*$COMP(m); xx = 4.0*(1-xx)*xx; %}');
pp_addpm('
sub ran_ver {
my ($obj,$x0,$r,$n) = @_;
if ($n>0) {
my $p = zeroes($n);
ran_ver_meat($p,$x0,$r,$n,$$obj);
return $p;
}
else {
barf("Not enough parameters for ran_ver!");
}
}
');
pp_addpm('
sub ran_caos {
my ($obj,$m,$n) = @_;
if ($n>0) {
my $p = zeroes($n);
ran_caos_meat($p,$m,$n,$$obj);
return $p;
}
else {
barf("Not enough parameters for ran_caos!");
}
}
');
# 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(slatec)
DEF_RNG(cmrg)
DEF_RNG(gfsr4)
DEF_RNG(minstd)
DEF_RNG(mrg)
DEF_RNG(mt19937)
DEF_RNG(r250)
DEF_RNG(ran0)
DEF_RNG(ran1)
DEF_RNG(ran2)
DEF_RNG(ran3)
DEF_RNG(rand48)
DEF_RNG(rand)
DEF_RNG(random8_bsd)
DEF_RNG(random8_glibc2)
DEF_RNG(random8_libc5)
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(random_bsd)
DEF_RNG(random_glibc2)
DEF_RNG(random_libc5)
DEF_RNG(randu)
DEF_RNG(ranf)
DEF_RNG(ranlux389)
DEF_RNG(ranlux)
DEF_RNG(ranmar)
DEF_RNG(taus)
DEF_RNG(transputer)
DEF_RNG(tt800)
DEF_RNG(uni32)
DEF_RNG(uni)
DEF_RNG(vax)
DEF_RNG(zuf)
DEF_RNG(default)
if (rng==NULL) {
barf("Unknown RNG, plese use one of the following: %s", rngs);
}
else
RETVAL = rng;
OUTPUT:
RETVAL
void
set_seed(rng, seed)
gsl_rng * rng
int seed
CODE:
gsl_rng_set(rng,seed);
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 = (gsl_rng *) SvIV((SV*)SvRV(sv));
/* fprintf(stderr,"Freeing %d\n",rng); */
gsl_rng_free((gsl_rng *) rng);
gsl_ran_discrete_t *
ran_discrete_preproc(rng, p)
gsl_rng * rng
pdl * p
CODE:
int n;
if (p->ndims!=1 || p->datatype!=PDL_D) {
barf("Bad input to ran_discrete_preproc!");
}
n = p->dims[0];
PDL->make_physical(p);
RETVAL = gsl_ran_discrete_preproc(n,(double *) p->data);
OUTPUT:
RETVAL
void
ran_shuffle(rng, in)
gsl_rng * rng
pdl * in
CODE:
int size, n;
n = in->nvals;
PDL->make_physical(in);
switch(in->datatype) {
case PDL_B: size=sizeof(PDL_Byte); break;
case PDL_S: size=sizeof(PDL_Short); break;
case PDL_US: size=sizeof(PDL_Ushort); break;
case PDL_L: size=sizeof(PDL_Long); break;
case PDL_F: size=sizeof(PDL_Float); break;
case PDL_D: size=sizeof(PDL_Double); break;
}
gsl_ran_shuffle(rng,(double *) in->data,n,size);
void
ran_choose(rng, in, out)
gsl_rng * rng
pdl * in
pdl * out
CODE:
int size, n,m;
n = in->nvals;
m = out->nvals;
if (in->datatype != out->datatype) barf("Data Types must match for ran_chooser");
PDL->make_physical(in);
PDL->make_physical(out);
switch(in->datatype) {
case PDL_B: size=sizeof(PDL_Byte); break;
case PDL_S: size=sizeof(PDL_Short); break;
case PDL_US: size=sizeof(PDL_Ushort); break;
case PDL_L: size=sizeof(PDL_Long); break;
case PDL_F: size=sizeof(PDL_Float); break;
case PDL_D: size=sizeof(PDL_Double); break;
}
gsl_ran_choose(rng,(double *) out->data, m, (double *) in->data,n,size);
');
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_done();