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