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 expected $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     return map { cplxe( gsl_complex_abs($_), gsl_complex_arg($_) ) } $self->get( [ 0 .. $self->length - 1  ] );
 163 }
 164
 165 =head2  get()
 166
 167 Gets the value of an of a Math::GSL::Vector object.
 168
 169     my $vector = Math::GSL::VectorComplex->new(3);
 170     ...
 171     my @values = $vector->get(2);
 172
 173 You can also enter an array of indices to receive their corresponding values:
 174
 175     my $vector = Math::GSL::VectorComplex->new(3);
 176     ...
 177     my @values = $vector->get([0,2]);
 178
 179 =cut
 180
 181 sub get {
 182     my ($self, $indices) = @_;
 183     return  map {  gsl_vector_complex_get($self->raw, $_ ) } @$indices ;
 184 }
 185
 186 =head2 reverse()
 187
 188 Returns the a vector with the elements in reversed order.
 189
 190     use Math::Complex;
 191     my $v1 = Math::GSL::VectorComplex->new([ 1, 2, 3*i]);
 192     my $v2 = $v1->reverse;
 193
 194 =cut
 195
 196 sub reverse {
 197     my $self = shift;
 198     my $copy = $self->copy();
 199     unless ( is_status_ok( gsl_vector_complex_reverse( $copy->raw )) ) {
 200         die( __PACKAGE__.": error reversing vector " . gsl_strerror($status) );
 201     }
 202     return $copy;
 203 }
 204
 205 =head2  set()
 206
 207 Sets values of an of a Math::GSL::Vector object.
 208
 209     my $vector = Math::GSL::VectorComplex->new(3);
 210     $vector->set([1,2], [8,23]);
 211
 212 This sets the second and third value to 8 and 23.
 213
 214 =cut
 215
 216 sub set {
 217     my ($self, $indices, $values) = @_;
 218     die (__PACKAGE__.'::set($indices, $values) - $indices and $values must be array references of the same length')
 219         unless ( ref $indices eq 'ARRAY' && ref $values eq 'ARRAY' &&  $#$indices == $#$values );
 220     eval {
 221         map {  gsl_vector_set($self->{_vector}, $indices->[$_], $values->[$_] ) } (0..$#$indices);
 222     };
 223     return;
 224 }
 225
 226 =head2 copy()
 227
 228 Returns a copy of the vector, which has the same length and values but resides at a different location in memory.
 229
 230     my $vector = Math::GSL::VectorComplex->new([10 .. 20]);
 231     my $copy   = $vector->copy;
 232
 233 =cut
 234
 235
 236 sub copy {
 237     my $self = shift;
 238     my $copy = Math::GSL::VectorComplex->new( $self->length );
 239     if ( gsl_vector_complex_memcpy($copy->raw, $self->raw) != $GSL_SUCCESS ) {
 240         croak "Math::GSL - error copying memory, aborting";
 241     }
 242     return $copy;
 243 }
 244
 245 sub _multiplication {
 246     my ($left,$right) = @_;
 247     my $lcopy = $left->copy;
 248
 249     if ( blessed $right && $right->isa(__PACKAGE__) ) {
 250         return $lcopy->dot_product($right);
 251     } else {
 252         # will be in upcoming gsl 1.12
 253         # gsl_vector_complex_scale($lcopy->raw, $right);
 254     }
 255     return $lcopy;
 256 }
 257
 258 sub _subtract {
 259     my ($left, $right, $flip) = @_;
 260
 261     if ($flip) {
 262         my $lcopy = $left->copy;
 263         gsl_vector_scale($lcopy->raw, -1 );
 264         gsl_vector_add_constant($lcopy->raw, $right);
 265         return $lcopy;
 266     } else {
 267         return _addition($left, -1.0*$right);
 268     }
 269 }
 270
 271 sub _addition {
 272     my ($left, $right, $flip) = @_;
 273
 274     my $lcopy = $left->copy;
 275
 276     if ( blessed $right && $right->isa('Math::GSL::Vector') && blessed $left && $left->isa('Math::GSL::Vector') ) {
 277         if ( $left->length == $right->length ) {
 278             gsl_vector_add($lcopy->raw, $right->raw);
 279         } else {
 280             croak "Math::GSL - addition of vectors must be called with two objects vectors and must have the same length";
 281         }
 282     } else {
 283         gsl_vector_add_constant($lcopy->raw, $right);
 284     }
 285     return $lcopy;
 286 }
 287
 288 sub dot_product_pp {
 289     my ($left,$right) = @_;
 290     my $sum=0;
 291     if ( blessed $right && $right->isa('Math::GSL::Vector') &&
 292          $left->length == $right->length ) {
 293          my @l = $left->as_list;
 294          my @r = $right->as_list;
 295          map { $sum += $l[$_] * $r[$_] } (0..$#l);
 296         return $sum;
 297     } else {
 298         croak "dot_product() must be called with two vectors";
 299     }
 300 }
 301
 302 sub dot_product {
 303     my ($left,$right) = @_;
 304
 305     my ($status, $product) = gsl_blas_ddot($left->raw,$right->raw);
 306     croak sprintf "Math::GSL::dot_product - %s", gsl_strerror($status) if ($status != $GSL_SUCCESS);
 307     return $product;
 308 }
 309
 310 =head1 DESCRIPTION
 311
 312 =head1 EXAMPLES
 313
 314 Here is an example using both interfaces.
 315
 316  use Math::GSL::Vector qw/:all/;
 317
 318  print "We'll create this vector : [0,1,4,9,16] \n";
 319  my $vector = Math::GSL::VectorComplex->new([0,1,4,9,16]);
 320  my ($min, $max) = gsl_vector_minmax_index($vector->raw);
 321
 322  print "We then check the index value of the maximum and minimum values of the vector. \n";
 323  print "The index of the maximum should be 4 and we received $max \n";
 324  print "The index of the minimum should be 0 and we received $min \n";
 325  print "We'll then swap the first and the third elements of the vector \n";
 326
 327  gsl_vector_swap_elements($vector->raw, 0, 3);
 328  my @got = $vector->as_list;
 329  print "The vector should now be like this : [9,1,4,0,16] \n";
 330  print "and we received : [ @got ]\n";
 331
 332 =head1 AUTHORS
 333
 334 Jonathan Leto <jonathan@leto.net> and Thierry Moisan <thierry.moisan@gmail.com>
 335
 336 =head1 COPYRIGHT AND LICENSE
 337
 338 Copyright (C) 2008 Jonathan Leto and Thierry Moisan
 339
 340 This program is free software; you can redistribute it and/or modify it
 341 under the same terms as Perl itself.
 342
 343 =cut
 344
 345 1;
 346 %}
 347