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 Carp qw/croak/;
  10 use Config;
  11 use Data::Dumper;
  12 use Test::More;
  13 use Scalar::Util qw/looks_like_number/;
  14 our @EXPORT = qw();
  15 our @EXPORT_OK = qw( ok_similar ok_status  is_similar
  16                      is_similar_relative verify verify_results
  17                      $GSL_MODE_DEFAULT $GSL_PREC_DOUBLE
  18                      $GSL_PREC_SINGLE $GSL_PREC_APPROX
  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.07';
  28
  29 =head1 NAME
  30
  31 Math::GSL - Perl interface to the  GNU Scientific Library (GSL) using SWIG
  32
  33 =head1 VERSION
  34
  35 Version 0.07
  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     return sort qw/
 185         Diff         Machine      Statistics
 186         Block        Eigen        Matrix        Poly
 187         BSpline      Errno        PowInt
 188         CBLAS        FFT          Min           IEEEUtils
 189         CDF          Fit          QRNG
 190         Chebyshev    Monte        RNG           Vector
 191         Heapsort     Multifit     Randist       Roots
 192         Combination  Histogram    Multimin      Wavelet
 193         Complex      Histogram2D  Multiroots    Wavelet2D
 194         Const        Siman        Sum           Sys
 195         Integration  NTuple       Sort
 196         DHT          Interp       ODEIV         SF
 197         Deriv        Linalg       Permutation   Spline
 198     /;
 199 }
 200
 201
 202 sub is_similar {
 203     my ($x,$y, $eps, $similarity_function) = @_;
 204     $eps ||= 1e-8;
 205     if (ref $x eq 'ARRAY' && ref $y eq 'ARRAY') {
 206         if ( $#$x != $#$y ){
 207             warn "is_similar(): argument of different lengths, $#$x != $#$y !!!";
 208             return 0;
 209         } else {
 210             map {
 211                     my $delta = abs($x->[$_] - $y->[$_]);
 212                     if($delta > $eps){
 213                         warn "\n\tElements start differing at index $_, delta = $delta\n";
 214                         warn qq{\t\t\$x->[$_] = } . $x->[$_] . "\n";
 215                         warn qq{\t\t\$y->[$_] = } . $y->[$_] . "\n";
 216                         return 0;
 217                     }
 218                 } (0..$#$x);
 219             return 1;
 220         }
 221     } else {
 222         if( ref $similarity_function eq 'CODE') {
 223                $similarity_function->($x,$y) <= $eps ? return 1 : return 0;
 224         } else {
 225             abs($x-$y) <= $eps ? return 1 : return 0;
 226         }
 227     }
 228 }
 229 sub ok_status {
 230     my ($got, $expected) = @_;
 231     $expected ||= $GSL_SUCCESS;
 232     ok( $got == $expected, gsl_strerror(int($got)) );
 233 }
 234 sub ok_similar {
 235     my ($x,$y, $msg, $eps) = @_;
 236     ok(is_similar($x,$y,$eps), $msg);
 237 }
 238
 239 sub is_similar_relative {
 240     my ($x,$y, $eps) = @_;
 241     return is_similar($x,$y,$eps, sub { abs( ($_[0] - $_[1])/abs($_[1]) ) } );
 242 }
 243
 244 # this is a huge hack
 245 sub verify_results
 246 {
 247     my ($results,$class) = @_;
 248     # GSL uses a factor of 100
 249     my $factor   = 20;
 250     my $eps      = 2048*$Math::GSL::Machine::GSL_DBL_EPSILON; # TOL3
 251     my ($x,$res);
 252
 253     croak "Usage: verify_results(%results, \$class)" unless $class;
 254     while (my($code,$expected)=each %$results){
 255         my $r        = Math::GSL::SF::gsl_sf_result_struct->new;
 256         my $status   = eval qq{${class}::$code};
 257
 258         ok(0, qq{'$code' died} ) if !defined $status;
 259
 260         if ( defined $r && $code =~ /_e\(.*\$r/) {
 261             $x   = $r->{val};
 262             $eps = $factor*$r->{err};
 263             $res = abs($x-$expected);
 264
 265             if ($ENV{DEBUG} ){
 266                 _dump_result($r);
 267                 print "got $code = $x\n";
 268                 printf "expected   : %.18g\n", $expected ;
 269                 printf "difference : %.18g\n", $res;
 270                 printf "unexpected error of %.18g\n", $res-$eps if ($res-$eps>0);
 271             }
 272
 273             if ($x =~ /nan|inf/i) {
 274                     ok( $expected eq $x, "'$expected'?='$x'" );
 275             } else {
 276                 ok( $res <= $eps, "$code ?= $x,\nres= +-$res, eps=$eps" );
 277             }
 278         }
 279     }
 280 }
 281 sub verify
 282 {
 283     my ($results,$class) = @_;
 284     croak "Usage: verify(%results, \$class)" unless $class;
 285     while (my($code,$result)=each %$results){
 286         my $x = eval qq{${class}::$code};
 287         ok(0, $@) if $@;
 288
 289         my ($expected,$eps);
 290         if (ref $result){
 291             ($expected,$eps)=@$result;
 292         } else {
 293             ($expected,$eps)=($result,1e-8);
 294         }
 295         my $res = abs($x - $expected);
 296         if ($x =~ /nan|inf/i ){
 297             ok( $expected eq $x, "'$expected' ?='$x'" );
 298         } else {
 299             ok( $res <= $eps, "$code ?= $x,\nres= +-$res, eps=$eps" );
 300         }
 301     }
 302 }
 303
 304 sub _dump_result($)
 305 {
 306     my $r=shift;
 307     printf "result->err: %.18g\n", $r->{err};
 308     printf "result->val: %.18g\n", $r->{val};
 309 }
 310
 311 42;