lib/Math/GSL.pm

   1 package Math::GSL;
   2 use base qw(Exporter);
   3 use base qw(DynaLoader);
   4 use strict;
   5 use warnings;
   6 use Math::GSL::Machine qw/:all/;
   7 use Math::GSL::Const qw/:all/;
   8 use Math::GSL::Errno qw/:all/;
   9 use Math::GSL::Vector qw/:file/;
  10 use Carp qw/croak/;
  11 use Config;
  12 use Test::More;
  13 our @EXPORT = qw();
  14 our @EXPORT_OK = qw( ok_similar ok_status  is_similar
  15                      is_similar_relative verify verify_results
  16                      $GSL_MODE_DEFAULT $GSL_PREC_DOUBLE
  17                      $GSL_PREC_SINGLE $GSL_PREC_APPROX
  18                      is_windows gsl_inf gsl_fopen gsl_fclose gsl_nan
  19                    );
  20
  21 our %EXPORT_TAGS = (
  22                      all => \@EXPORT_OK,
  23                    );
  24
  25 our ($GSL_PREC_DOUBLE, $GSL_PREC_SINGLE, $GSL_PREC_APPROX ) = 0..2;
  26 our $GSL_MODE_DEFAULT = $GSL_PREC_DOUBLE;
  27 our $VERSION = '0.09_01';
  28
  29 =head1 NAME
  30
  31 Math::GSL - Perl interface to the  GNU Scientific Library (GSL)
  32
  33 =head1 VERSION
  34
  35 Version 0.09_01
  36
  37 =cut
  38
  39 =head1 SYNOPSIS
  40
  41     use Math::GSL::Matrix qw/:all/;
  42     my $matrix = Math::GSL::Matrix->new(5,5);   # 5x5 zero matrix
  43     $matrix->set_col(0, [1..5])
  44            ->set_row(2, [5..9]);
  45     my @matrix = $matrix->as_list;              # matrix as Perl list
  46     my $gsl_matrix = $matrix->raw;              # underlying GSL object
  47
  48 Each GSL subsystem has it's own module. For example, the random number generator
  49 subsystem is Math::GSL::RNG. Many subsystems have a more Perlish and
  50 object-oriented frontend which can be used, as the above example shows. The raw
  51 GSL object is useful for using the low-level GSL functions, which in the case of
  52 the Matrix subsytem, would be of the form gsl_matrix_* . Each module has further
  53 documentation about the low-level C functions as well as using the more
  54 intuitive (but slightly slower) object-oriented interface.
  55
  56 =head1 SUBSYSTEMS
  57
  58     Math::GSL::BLAS             - Linear Algebra Functions
  59     Math::GSL::BSpline          - BSplines
  60     Math::GSL::CBLAS            - Linear Algebra Functions
  61     Math::GSL::CDF              - Cumulative Distribution Functions
  62     Math::GSL::Chebyshev        - Chebyshev Polynomials
  63     Math::GSL::Combination      - Combinatoric Functions
  64     Math::GSL::Complex          - Complex Numbers
  65     Math::GSL::Const            - Various Constants
  66     Math::GSL::DHT              - Discrete Hilbert Transform
  67     Math::GSL::Deriv            - Numerical Derivative
  68     Math::GSL::Eigen            - Eigenvalues and Eigenvectors
  69     Math::GSL::Errno            - Error Handling
  70     Math::GSL::FFT              - Fast Fourier Transform
  71     Math::GSL::Fit              - Curve Fitting
  72     Math::GSL::Heapsort         - Sorting Heaps
  73     Math::GSL::Histogram        - Histograms
  74     Math::GSL::Histogram2D      - 2D Histograms
  75     Math::GSL::Integration      - Numerical Integration
  76     Math::GSL::Interp           - Interpolation
  77     Math::GSL::Linalg           - Linear Algebra
  78     Math::GSL::Machine          - Machine Specific Information
  79     Math::GSL::Matrix           - NxM Matrices
  80     Math::GSL::Min              - Minimization
  81     Math::GSL::Monte            - Monte Carlo Integrations
  82     Math::GSL::Multifit         - Multivariable Fitting
  83     Math::GSL::Multimin         - Multivariable Minimization
  84     Math::GSL::Multiroots       - Muiltvariable Root Finding
  85     Math::GSL::NTuple           - N Tuples
  86     Math::GSL::ODEIV            - Ordinary Differential Equation Solvers (Initial Value Problems)
  87     Math::GSL::Permutation      - Permutations
  88     Math::GSL::Poly             - Polynmials
  89     Math::GSL::PowInt           - Integer Power Functions
  90     Math::GSL::QRNG             - Quasi-Random Number Generators
  91     Math::GSL::RNG              - Random Number Generators
  92     Math::GSL::Randist          - Random Number Distributions
  93     Math::GSL::Roots            - Root Finding Algorithms
  94     Math::GSL::SF               - Special Functions
  95     Math::GSL::Siman            - Simulated Annealing
  96     Math::GSL::Sort             - Sorting
  97     Math::GSL::Spline           - Splines
  98     Math::GSL::Statistics       - Statistics Functions
  99     Math::GSL::Sum              - Summation
 100     Math::GSL::Sys              - Sytem utility functions
 101     Math::GSL::Vector           - N-dimensional Vectors
 102     Math::GSL::Wavelet          - Basic Wavelets
 103     Math::GSL::Wavelet2D        - 2D Wavelets
 104
 105
 106 =head1 AUTHORS
 107
 108 Jonathan Leto, C<< <jonathan@leto.net> >> and Thierry Moisan C<< <thierry.moisan@gmail.com> >>
 109
 110 =head1 BUGS
 111
 112 This software is still in active development, we appreciate your detailed bug reports and
 113 documentation patches. Please report any bugs or feature requests to the authors directly.
 114
 115
 116 =head1 SUPPORT
 117
 118 You can find documentation for this module with the perldoc command.
 119
 120     perldoc Math::GSL
 121
 122 or online at L<http://leto.net/code/Math-GSL/>
 123
 124 =over 4
 125
 126 =item * AnnoCPAN: Annotated CPAN documentation
 127
 128 L<http://annocpan.org/dist/Math::GSL>
 129
 130 =item * CPAN Ratings
 131
 132 L<http://cpanratings.perl.org/d/Math::GSL>
 133
 134 =item * Search CPAN
 135
 136 L<http://search.cpan.org/dist/Math::GSL>
 137
 138 =back
 139
 140 =head1 DEVELOPMENT
 141
 142 If you would like the help develop Math::GSL, email the authors
 143 and do
 144
 145     git clone http://leto.net/code/Math-GSL.git
 146     cd Math-GSL
 147     git checkout -b bleed   # create new local branch
 148     git pull origin bleed   # pull in remote bleed
 149
 150 to get the latest source, which is a two-headed beast with branches master and
 151 bleed. The master branch is our stable branch, which is periodically sync-ed
 152 with bleed. To view the latest source code online, go to
 153 L<http://leto.net/gitweb/>.  The latest version of Git can be found at
 154 L<http://git.or.cz>.
 155
 156 =head1 ACKNOWLEDGEMENTS
 157
 158 Thanks to PDX.pm, The Perl Foundation and everyone at Google who makes
 159 the Summer of Code happen each year. You rock.
 160
 161 =head1 DEDICATION
 162
 163 This Perl module is dedicated in memory of Nick Ing-Simmons.
 164
 165 =head1 COPYRIGHT & LICENSE
 166
 167 Copyright 2008 Jonathan Leto, Thierry Moisan all rights reserved.
 168
 169 This program is free software; you can redistribute it and/or modify it
 170 under the same terms as Perl itself.
 171
 172 =cut
 173
 174 sub new
 175 {
 176     my ($self,$args) = @_;
 177     my $class = ref $self || $self || 'Math::GSL';
 178     my $this = { };
 179     bless $this, $class;
 180 }
 181
 182 sub subsystems
 183 {
 184     # causing havoc to cpantesters
 185     my @disable = qw(NTuple);
 186     return sort qw/
 187         Diff         Machine      Statistics
 188         Eigen        Matrix       Poly
 189         BSpline      Errno        PowInt
 190         CBLAS        FFT          Min           IEEEUtils
 191         CDF          Fit          QRNG
 192         Chebyshev    Monte        RNG           Vector
 193         Heapsort     Multifit     Randist       Roots
 194         Combination  Histogram    Multimin      Wavelet
 195         Complex      Histogram2D  Multiroots    Wavelet2D
 196         Const        Siman        Sum           Sys
 197         Integration  Sort
 198         DHT          Interp       ODEIV         SF
 199         Deriv        Linalg       Permutation   Spline
 200     /;
 201 }
 202
 203
 204 sub is_similar {
 205     my ($x,$y, $eps, $similarity_function) = @_;
 206     $eps ||= 1e-8;
 207     if (ref $x eq 'ARRAY' && ref $y eq 'ARRAY') {
 208         if ( $#$x != $#$y ){
 209             warn "is_similar(): argument of different lengths, $#$x != $#$y !!!";
 210             return 0;
 211         } else {
 212             map {
 213                     my $delta = abs($x->[$_] - $y->[$_]);
 214                     if($delta > $eps){
 215                         warn "\n\tElements start differing at index $_, delta = $delta\n";
 216                         warn qq{\t\t\$x->[$_] = } . $x->[$_] . "\n";
 217                         warn qq{\t\t\$y->[$_] = } . $y->[$_] . "\n";
 218                         return 0;
 219                     }
 220                 } (0..$#$x);
 221             return 1;
 222         }
 223     } else {
 224         if( ref $similarity_function eq 'CODE') {
 225                $similarity_function->($x,$y) <= $eps ? return 1 : return 0;
 226         } else {
 227             abs($x-$y) <= $eps ? return 1 : return 0;
 228         }
 229     }
 230 }
 231 sub ok_status {
 232     my ($got, $expected) = @_;
 233     $expected ||= $GSL_SUCCESS;
 234     ok( $got == $expected, gsl_strerror(int($got)) );
 235 }
 236 sub ok_similar {
 237     my ($x,$y, $msg, $eps) = @_;
 238     ok(is_similar($x,$y,$eps), $msg);
 239 }
 240
 241 sub is_similar_relative {
 242     my ($x,$y, $eps) = @_;
 243     return is_similar($x,$y,$eps, sub { abs( ($_[0] - $_[1])/abs($_[1]) ) } );
 244 }
 245
 246 # this is a huge hack
 247 sub verify_results
 248 {
 249     my ($results,$class) = @_;
 250     # GSL uses a factor of 100
 251     my $factor   = 20;
 252     my $eps      = 2048*$Math::GSL::Machine::GSL_DBL_EPSILON; # TOL3
 253     my ($x,$res);
 254
 255     croak "Usage: verify_results(%results, \$class)" unless $class;
 256     while (my($code,$expected)=each %$results){
 257         my $r        = Math::GSL::SF::gsl_sf_result_struct->new;
 258         my $status   = eval qq{${class}::$code};
 259
 260         ok(0, qq{'$code' died} ) if !defined $status;
 261
 262         if ( defined $r && $code =~ /_e\(.*\$r/) {
 263             $x   = $r->{val};
 264             $eps = $factor*$r->{err};
 265             $res = abs($x-$expected);
 266
 267             if ($ENV{DEBUG} ){
 268                 _dump_result($r);
 269                 print "got $code = $x\n";
 270                 printf "expected   : %.18g\n", $expected ;
 271                 printf "difference : %.18g\n", $res;
 272                 printf "unexpected error of %.18g\n", $res-$eps if ($res-$eps>0);
 273             }
 274             if ($x =~ /nan|inf/i) {
 275                     ok( $expected eq $x, "'$expected'?='$x'" );
 276             } else {
 277                 ok( $res <= $eps, "$code ?= $x,\nres= +-$res, eps=$eps" );
 278             }
 279         }
 280     }
 281 }
 282 sub verify
 283 {
 284     my ($results,$class) = @_;
 285     croak "Usage: verify(%results, \$class)" unless $class;
 286     while (my($code,$result)=each %$results){
 287         my $x = eval qq{${class}::$code};
 288         ok(0, $@) if $@;
 289
 290         my ($expected,$eps);
 291         if (ref $result){
 292             ($expected,$eps)=@$result;
 293         } else {
 294             ($expected,$eps)=($result,1e-8);
 295         }
 296         my $res = abs($x - $expected);
 297         if ($x =~ /nan|inf/i ){
 298             ok( $expected eq $x, "'$expected' ?='$x'" );
 299         } else {
 300             ok( $res <= $eps, "$code ?= $x,\nres= +-$res, eps=$eps" );
 301         }
 302     }
 303 }
 304
 305 sub gsl_inf
 306 {
 307     is_windows() ?  unpack 'd', scalar reverse pack 'H*', '7FF0000000000000'
 308                  : q{inf}
 309 }
 310 sub gsl_nan
 311 {
 312     is_windows() ? unpack 'd', scalar reverse pack 'H*', '7FF8000000000000'
 313                  : q{nan}
 314              }
 315 sub gsl_fopen
 316 {
 317     my ($file, $mode) = @_;
 318     $mode .= '+b' if (is_windows() and $mode !~ /\+b/);
 319     return Math::GSL::Vector::fopen($file, $mode);
 320 }
 321
 322 sub gsl_fclose
 323 {
 324     my $file = shift;
 325     return Math::GSL::Vector::fclose($file);
 326 }
 327
 328
 329 sub _dump_result($)
 330 {
 331     my $r=shift;
 332     printf "result->err: %.18g\n", $r->{err};
 333     printf "result->val: %.18g\n", $r->{val};
 334 }
 335 sub is_windows() { $^O =~ /MSWin32/i }
 336
 337 42;