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::Histograma       - 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::Mode             - GSL Precision Modes
  82     Math::GSL::Monte            - Monte Carlo Integrations
  83     Math::GSL::Multifit         - Multivariable Fitting
  84     Math::GSL::Multimin         - Multivariable Minimization
  85     Math::GSL::Multiroots       - Muiltvariable Root Finding
  86     Math::GSL::NTuple           - N Tuples
  87     Math::GSL::ODEIV            - Ordinary Differential Equation Solvers (Initial Value Problems)
  88     Math::GSL::Permutation      - Permutations
  89     Math::GSL::Poly             - Polynmials
  90     Math::GSL::PowInt           - Integer Power Functions
  91     Math::GSL::QRNG             - Quasi-Random Number Generators
  92     Math::GSL::RNG              - Random Number Generators
  93     Math::GSL::Randist          - Random Number Distributions
  94     Math::GSL::Roots            - Root Finding Algorithms
  95     Math::GSL::SF               - Special Functions
  96     Math::GSL::Siman            - Simulated Annealing
  97     Math::GSL::Sort             - Sorting
  98     Math::GSL::Spline           - Splines
  99     Math::GSL::Statistics       - Statistics Functions
 100     Math::GSL::Sum              - Summation
 101     Math::GSL::Sys              - Sytem utility functions
 102     Math::GSL::Vector           - N-dimensional Vectors
 103     Math::GSL::Wavelet          - Basic Wavelets
 104     Math::GSL::Wavelet2D        - 2D Wavelets
 105
 106
 107 =head1 AUTHORS
 108
 109 Jonathan Leto, C<< <jonathan@leto.net> >> and Thierry Moisan C<< <thierry.moisan@gmail.com> >>
 110
 111 =head1 BUGS
 112
 113 This software is still in active development, we appreciate your detailed bug reports.
 114 Please report any bugs or feature requests to the authors directly.
 115
 116
 117 =head1 SUPPORT
 118
 119 You can find documentation for this module with the perldoc command.
 120
 121     perldoc Math::GSL
 122
 123 or online at L<http://leto.net/code/Math-GSL/>
 124
 125 =over 4
 126
 127 =item * AnnoCPAN: Annotated CPAN documentation
 128
 129 L<http://annocpan.org/dist/Math::GSL>
 130
 131 =item * CPAN Ratings
 132
 133 L<http://cpanratings.perl.org/d/Math::GSL>
 134
 135 =item * Search CPAN
 136
 137 L<http://search.cpan.org/dist/Math::GSL>
 138
 139 =back
 140
 141 =head1 DEVELOPMENT
 142
 143 If you would like the help develop Math::GSL, email the authors
 144 and do
 145
 146     git clone http://leto.net/code/Math-GSL.git
 147     cd Math-GSL
 148     git checkout -b bleed   # create new local branch
 149     git pull origin bleed   # pull in remote bleed
 150
 151 to get the latest source, which is a two-headed beast with branches master and
 152 bleed. The master branch is our stable branch, which is periodically sync-ed
 153 with bleed. To view the latest source code online, go to
 154 L<http://leto.net/gitweb/>.  The latest version of Git can be found at
 155 L<http://git.or.cz>.
 156
 157 =head1 ACKNOWLEDGEMENTS
 158
 159
 160 =head1 COPYRIGHT & LICENSE
 161
 162 Copyright 2008 Jonathan Leto, Thierry Moisan all rights reserved.
 163
 164 This program is free software; you can redistribute it and/or modify it
 165 under the same terms as Perl itself.
 166
 167 =cut
 168
 169 sub new
 170 {
 171     my ($self,$args) = @_;
 172     my $class = ref $self || $self || 'Math::GSL';
 173     my $this = { };
 174     bless $this, $class;
 175 }
 176
 177 sub subsystems
 178 {
 179     return sort qw/
 180         Diff         Machine      Statistics
 181         Block        Eigen        Matrix        Poly
 182         BSpline      Errno        PowInt
 183         CBLAS        FFT          Min           IEEEUtils
 184         CDF          Fit          Mode          QRNG
 185         Chebyshev    Monte        RNG           Vector
 186         Heapsort     Multifit     Randist       Roots
 187         Combination  Histogram    Multimin      Wavelet
 188         Complex      Histogram2D  Multiroots    Wavelet2D
 189         Const        Siman        Sum           Sys
 190         Integration  NTuple       Sort
 191         DHT          Interp       ODEIV         SF
 192         Deriv        Linalg       Permutation   Spline
 193     /;
 194 }
 195
 196
 197 sub is_similar {
 198     my ($x,$y, $eps, $similarity_function) = @_;
 199     $eps ||= 1e-8;
 200     if (ref $x eq 'ARRAY' && ref $y eq 'ARRAY') {
 201         if ( $#$x != $#$y ){
 202             warn "is_similar(): argument of different lengths, $#$x != $#$y !!!";
 203             return 0;
 204         } else {
 205             map {
 206                     my $delta = abs($x->[$_] - $y->[$_]);
 207                     if($delta > $eps){
 208                         warn "\n\tElements start differing at index $_, delta = $delta\n";
 209                         warn qq{\t\t\$x->[$_] = } . $x->[$_] . "\n";
 210                         warn qq{\t\t\$y->[$_] = } . $y->[$_] . "\n";
 211                         return 0;
 212                     }
 213                 } (0..$#$x);
 214             return 1;
 215         }
 216     } else {
 217         if( ref $similarity_function eq 'CODE') {
 218                $similarity_function->($x,$y) <= $eps ? return 1 : return 0;
 219         } else {
 220             abs($x-$y) <= $eps ? return 1 : return 0;
 221         }
 222     }
 223 }
 224 sub ok_status {
 225     my ($got, $expected) = @_;
 226     ok( $got == $expected, gsl_strerror($expected) );
 227 }
 228 sub ok_similar {
 229     my ($x,$y, $msg, $eps) = @_;
 230     ok(is_similar($x,$y,$eps), $msg);
 231 }
 232
 233 sub is_similar_relative {
 234     my ($x,$y, $eps) = @_;
 235     return is_similar($x,$y,$eps, sub { abs( ($_[0] - $_[1])/abs($_[1]) ) } );
 236 }
 237
 238 # this is a huge hack
 239 sub verify_results
 240 {
 241     my ($results,$class) = @_;
 242     # GSL uses a factor of 100
 243     my $factor   = 20;
 244     my $eps      = 2048*$Math::GSL::Machine::GSL_DBL_EPSILON; # TOL3
 245     my ($x,$res);
 246
 247     croak "Usage: verify_results(%results, \$class)" unless $class;
 248     while (my($code,$expected)=each %$results){
 249         my $r        = Math::GSL::SF::gsl_sf_result_struct->new;
 250         my $status   = eval qq{${class}::$code};
 251
 252         ok(0, qq{'$code' died} ) if !defined $status;
 253
 254         if ( defined $r && $code =~ /_e\(.*\$r/) {
 255             $x   = $r->{val};
 256             $eps = $factor*$r->{err};
 257             $res = abs($x-$expected);
 258
 259             if ($ENV{DEBUG} ){
 260                 _dump_result($r);
 261                 print "got $code = $x\n";
 262                 printf "expected   : %.18g\n", $expected ;
 263                 printf "difference : %.18g\n", $res;
 264                 printf "unexpected error of %.18g\n", $res-$eps if ($res-$eps>0);
 265             }
 266
 267             if ($x =~ /nan|inf/i) {
 268                     ok( $expected eq $x, "'$expected'?='$x'" );
 269             } else {
 270                 ok( $res <= $eps, "$code ?= $x,\nres= +-$res, eps=$eps" );
 271             }
 272         }
 273     }
 274 }
 275 sub verify
 276 {
 277     my ($results,$class) = @_;
 278     croak "Usage: verify(%results, \$class)" unless $class;
 279     while (my($code,$result)=each %$results){
 280         my $x = eval qq{${class}::$code};
 281         ok(0, $@) if $@;
 282
 283         my ($expected,$eps);
 284         if (ref $result){
 285             ($expected,$eps)=@$result;
 286         } else {
 287             ($expected,$eps)=($result,1e-8);
 288         }
 289         my $res = abs($x - $expected);
 290         if ($x =~ /nan|inf/i ){
 291             ok( $expected eq $x, "'$expected' ?='$x'" );
 292         } else {
 293             ok( $res <= $eps, "$code ?= $x,\nres= +-$res, eps=$eps" );
 294         }
 295     }
 296 }
 297
 298 sub _dump_result($)
 299 {
 300     my $r=shift;
 301     printf "result->err: %.18g\n", $r->{err};
 302     printf "result->val: %.18g\n", $r->{val};
 303 }
 304
 305 42;