lib/Math/GSL.pm

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