pod/VectorComplex.pod

   1 %perlcode %{
   2 use Scalar::Util 'blessed';
   3 use Data::Dumper;
   4 use Carp qw/croak/;
   5 use Math::GSL::Errno qw/:all/;
   6 use Math::GSL::BLAS qw/gsl_blas_ddot/;
   7 use Math::GSL::Complex qw/:all/;
   8 use Math::GSL::Test qw/is_status_ok/;
   9 use Math::Complex;
  10 use overload
  11     '*'      => \&_multiplication,
  12     '+'      => \&_addition,
  13     '-'      => \&_subtract,
  14     fallback => 1,
  15 ;
  16
  17 @EXPORT_all  = qw/fopen fclose
  18                  gsl_vector_complex_alloc gsl_vector_complex_calloc gsl_vector_complex_alloc_from_block gsl_vector_complex_alloc_from_vector
  19                  gsl_vector_complex_free gsl_vector_complex_view_array gsl_vector_complex_view_array_with_stride gsl_vector_complex_const_view_array
  20                  gsl_vector_complex_const_view_array_with_stride gsl_vector_complex_subvector gsl_vector_complex_subvector_with_stride
  21                  gsl_vector_complex_const_subvector gsl_vector_complex_const_subvector_with_stride gsl_vector_complex_real gsl_vector_complex_imag
  22                  gsl_vector_complex_const_real gsl_vector_complex_const_imag gsl_vector_complex_get gsl_vector_complex_set
  23                  gsl_vector_complex_ptr gsl_vector_complex_const_ptr gsl_vector_complex_set_zero gsl_vector_complex_set_all
  24                  gsl_vector_complex_set_basis gsl_vector_complex_fread gsl_vector_complex_fwrite gsl_vector_complex_fscanf
  25                  gsl_vector_complex_fprintf gsl_vector_complex_memcpy gsl_vector_complex_reverse gsl_vector_complex_swap
  26                  gsl_vector_complex_swap_elements gsl_vector_complex_isnull gsl_vector_complex_ispos gsl_vector_complex_isneg
  27 /;
  28 @EXPORT_file =qw/ fopen fclose/;
  29 @EXPORT_OK = (@EXPORT_all, @EXPORT_file);
  30 %EXPORT_TAGS = ( file => \@EXPORT_file, all => \@EXPORT_all );
  31
  32 =head1 NAME
  33
  34 Math::GSL::VectorComplex - Complex Vectors
  35
  36 =head1 SYNOPSIS
  37
  38     use Math::GSL::VectorComplex qw/:all/;
  39     my $vec1 = Math::GSL::VectorComplex->new([1 + 2*i, 7*i, 5, -3 ]);
  40     my $vec2 = $vec1 * 5;
  41     my $vec3 = Math::GSL::Vector>new(10);   # 10 element zero vector
  42     my $vec4 = $vec1 + $vec2;
  43
  44     # set the element at index 1 to -i
  45     # and the element at index 3 to i
  46     $vec3->set([ 1, -i ], [ 9, i ]);
  47
  48     my @vec = $vec2->as_list;               # return elements as Perl list
  49
  50     my $dot_product = $vec1 * $vec2;
  51     my $length      = $vec2->length;
  52     my $first       = $vec1->get(0);
  53
  54
  55 =cut
  56
  57 =head1 Objected Oriented Interface to GSL Math::GSL::VectorComplex
  58
  59 =head2 new()
  60
  61 Creates a new Vector of the given size.
  62
  63     my $vector = Math::GSL::VectorComplex->new(3);
  64
  65 You can also create and set directly the values of the vector like this :
  66
  67    my $vector = Math::GSL::VectorComplex->new([2,4,1]);
  68
  69 =cut
  70
  71 sub new {
  72     my ($class, $values) = @_;
  73     my $length  = $#$values;
  74     my $this = {};
  75     my $vector;
  76
  77     # we expect $values to have Math::Complex objects
  78     @$values = map { gsl_complex_rect(Re($_), Im($_)) } @$values;
  79
  80     if ( ref $values eq 'ARRAY' ){
  81         die __PACKAGE__.'::new($x) - $x must be a nonempty array reference' if $length == -1;
  82         $vector  = gsl_vector_complex_alloc($length+1);
  83         map { gsl_vector_complex_set($vector, $_, $values->[$_] ) }  (0 .. $length);
  84         $this->{_length} = $length+1;
  85     } elsif ( (int($values) == $values) && ($values > 0)) {
  86         $vector  = gsl_vector_complex_alloc($values);
  87         gsl_vector_complex_set_zero($vector);
  88         $this->{_length} = $values;
  89     } else {
  90         die __PACKAGE__.'::new($x) - $x must be an int or array reference';
  91     }
  92     $this->{_vector} = $vector;
  93     bless $this, $class;
  94 }
  95
  96
  97 =head2 raw()
  98
  99 Get the underlying GSL vector object created by SWIG, useful for using gsl_vector_* functions which do not have an OO counterpart.
 100
 101     my $vector    = Math::GSL::VectorComplex->new(3);
 102     my $gsl_vector = $vector->raw;
 103     my $stuff      = gsl_vector_get($gsl_vector, 1);
 104
 105 =cut
 106
 107 sub raw {
 108     my $self = shift;
 109     return $self->{_vector};
 110 }
 111
 112 =head2 min()
 113
 114 Returns the minimum value contained in the vector.
 115
 116    my $vector = Math::GSL::VectorComplex->new([2,4,1]);
 117    my $minimum = $vector->min;
 118
 119 =cut
 120
 121 sub min {
 122     my $self=shift;
 123     return gsl_vector_min($self->raw);
 124 }
 125
 126 =head2 max()
 127
 128 Returns the minimum value contained in the vector.
 129
 130    my $vector = Math::GSL::VectorComplex->new([2,4,1]);
 131    my $maximum = $vector->max;
 132
 133 =cut
 134
 135 sub max {
 136     my $self=shift;
 137     return gsl_vector_max($self->raw);
 138 }
 139
 140 =head2 length()
 141
 142 Returns the number of elements contained in the vector.
 143
 144    my $vector = Math::GSL::VectorComplex->new([2,4,1]);
 145    my $length = $vector->length;
 146
 147 =cut
 148
 149 sub length { my $self=shift; $self->{_length} }
 150
 151 =head2  as_list()
 152
 153 Gets the content of a Math::GSL::Vector object as a Perl list.
 154
 155     my $vector = Math::GSL::VectorComplex->new(3);
 156     ...
 157     my @values = $vector->as_list;
 158 =cut
 159
 160 sub as_list {
 161     my $self=shift;
 162     # this is wrong
 163     return map { cplxe( gsl_complex_abs($_), gsl_complex_arg($_) ) } $self->get( [ 0 .. $self->length - 1  ] );
 164 }
 165
 166 =head2  get()
 167
 168 Gets the value of an of a Math::GSL::Vector object.
 169
 170     my $vector = Math::GSL::VectorComplex->new(3);
 171     ...
 172     my @values = $vector->get(2);
 173
 174 You can also enter an array of indices to receive their corresponding values:
 175
 176     my $vector = Math::GSL::VectorComplex->new(3);
 177     ...
 178     my @values = $vector->get([0,2]);
 179
 180 =cut
 181
 182 sub get {
 183     my ($self, $indices) = @_;
 184     return  map {  gsl_vector_complex_get($self->raw, $_ ) } @$indices ;
 185 }
 186
 187 =head2 reverse()
 188
 189 Returns the a vector with the elements in reversed order.
 190
 191     use Math::Complex;
 192     my $v1 = Math::GSL::VectorComplex->new([ 1, 2, 3*i]);
 193     my $v2 = $v1->reverse;
 194
 195 =cut
 196
 197 sub reverse {
 198     my $self = shift;
 199     my $copy = $self->copy();
 200     unless ( is_status_ok( gsl_vector_complex_reverse( $copy->raw )) ) {
 201         die( __PACKAGE__.": error reversing vector " . gsl_strerror($status) );
 202     }
 203     return $copy;
 204 }
 205
 206 =head2  set()
 207
 208 Sets values of an of a Math::GSL::Vector object.
 209
 210     my $vector = Math::GSL::VectorComplex->new(3);
 211     $vector->set([1,2], [8,23]);
 212
 213 This sets the second and third value to 8 and 23.
 214
 215 =cut
 216
 217 sub set {
 218     my ($self, $indices, $values) = @_;
 219     die (__PACKAGE__.'::set($indices, $values) - $indices and $values must be array references of the same length')
 220         unless ( ref $indices eq 'ARRAY' && ref $values eq 'ARRAY' &&  $#$indices == $#$values );
 221     eval {
 222         map {  gsl_vector_complex_set($self->{_vector}, $indices->[$_], $values->[$_] ) } (0..$#$indices);
 223     };
 224     # better error handling?
 225     warn $@ if $@;
 226     return;
 227 }
 228
 229 =head2 copy()
 230
 231 Returns a copy of the vector, which has the same length and values but resides at a different location in memory.
 232
 233     my $vector = Math::GSL::VectorComplex->new([10 .. 20]);
 234     my $copy   = $vector->copy;
 235
 236 =cut
 237
 238 sub copy {
 239     my $self = shift;
 240     my $copy = Math::GSL::VectorComplex->new( $self->length );
 241     my $status = gsl_vector_complex_memcpy($copy->raw, $self->raw);
 242     if ( $status != $GSL_SUCCESS ) {
 243         croak "Math::GSL - error copying memory, aborting. $! status=$status";
 244     }
 245     return $copy;
 246 }
 247
 248 =head2 swap()
 249
 250 Exchanges the values in the vectors $v with $w by copying.
 251
 252     my $v = Math::GSL::VectorComplex->new([1..5]);
 253     my $w = Math::GSL::VectorComplex->new([3..7]);
 254     $v->swap( $w );
 255
 256 =cut
 257
 258 sub swap() {
 259     my ($self,$other) = @_;
 260     croak "Math::GSL::VectorComplex : \$v->swap(\$w) - \$w must be a Math::GSL::VectorComplex"
 261         unless ref $other eq 'Math::GSL::VectorComplex';
 262     gsl_vector_complex_swap( $self->raw, $other->raw );
 263     return $self;
 264 }
 265
 266 sub _multiplication {
 267     my ($left,$right) = @_;
 268     my $lcopy = $left->copy;
 269
 270     if ( blessed $right && $right->isa(__PACKAGE__) ) {
 271         return $lcopy->dot_product($right);
 272     } else {
 273         # will be in upcoming gsl 1.12
 274         # gsl_vector_complex_scale($lcopy->raw, $right);
 275     }
 276     return $lcopy;
 277 }
 278
 279 sub _subtract {
 280     my ($left, $right, $flip) = @_;
 281
 282     if ($flip) {
 283         my $lcopy = $left->copy;
 284         # will be in upcoming gsl 1.12
 285         # gsl_vector_complex_scale($lcopy->raw, -1 );
 286         gsl_vector_add_constant($lcopy->raw, $right);
 287         return $lcopy;
 288     } else {
 289         return _addition($left, -1.0*$right);
 290     }
 291 }
 292
 293 sub _addition {
 294     my ($left, $right, $flip) = @_;
 295     my $lcopy = $left->copy;
 296
 297     if ( blessed $right && $right->isa('Math::GSL::Vector') && blessed $left && $left->isa('Math::GSL::Vector') ) {
 298         if ( $left->length == $right->length ) {
 299             gsl_vector_complex_add($lcopy->raw, $right->raw);
 300         } else {
 301             croak "Math::GSL - addition of vectors must be called with two objects vectors and must have the same length";
 302         }
 303     } else {
 304         gsl_vector_complex_add_constant($lcopy->raw, $right);
 305     }
 306     return $lcopy;
 307 }
 308
 309 sub dot_product_pp {
 310     my ($left,$right) = @_;
 311     my $sum=0;
 312     if ( blessed $right && $right->isa('Math::GSL::Vector') && $left->length == $right->length ) {
 313          my @l = $left->as_list;
 314          my @r = $right->as_list;
 315          map { $sum += $l[$_] * $r[$_] } (0..$#l);
 316         return $sum;
 317     } else {
 318         croak "dot_product() must be called with two vectors";
 319     }
 320 }
 321
 322 sub dot_product {
 323     my ($left,$right) = @_;
 324
 325     my ($status, $product) = gsl_blas_ddot($left->raw,$right->raw);
 326     croak sprintf "Math::GSL::dot_product - %s", gsl_strerror($status) if ($status != $GSL_SUCCESS);
 327     return $product;
 328 }
 329
 330 =head1 AUTHORS
 331
 332 Jonathan Leto <jonathan@leto.net> and Thierry Moisan <thierry.moisan@gmail.com>
 333
 334 =head1 COPYRIGHT AND LICENSE
 335
 336 Copyright (C) 2008 Jonathan Leto and Thierry Moisan
 337
 338 This program is free software; you can redistribute it and/or modify it
 339 under the same terms as Perl itself.
 340
 341 =cut
 342
 343 %}