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/
  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_OK = (@EXPORT_all);
  29 %EXPORT_TAGS = ( all => \@EXPORT_all );
  30
  31 =head1 NAME
  32
  33 Math::GSL::VectorComplex - Complex Vectors
  34
  35 =head1 SYNOPSIS
  36
  37     use Math::GSL::VectorComplex qw/:all/;
  38     my $vec1 = Math::GSL::VectorComplex->new([1 + 2*i, 7*i, 5, -3 ]);
  39     my $vec2 = $vec1 * 5;
  40     my $vec3 = Math::GSL::Vector>new(10);   # 10 element zero vector
  41     my $vec4 = $vec1 + $vec2;
  42
  43     # set the element at index 1 to -i
  44     # and the element at index 3 to i
  45     $vec3->set([ 1, -i ], [ 9, i ]);
  46
  47     my @vec = $vec2->as_list;               # return elements as Perl list
  48
  49     my $dot_product = $vec1 * $vec2;
  50     my $length      = $vec2->length;
  51     my $first       = $vec1->get(0);
  52
  53
  54 =cut
  55
  56 =head1 Objected Oriented Interface to GSL Math::GSL::VectorComplex
  57
  58 =head2 new()
  59
  60 Creates a new Vector of the given size.
  61
  62     my $vector = Math::GSL::VectorComplex->new(3);
  63
  64 You can also create and set directly the values of the vector like this :
  65
  66    my $vector = Math::GSL::VectorComplex->new([2,4,1]);
  67
  68 =cut
  69
  70 sub new {
  71     my ($class, $values) = @_;
  72     my $length  = $#$values;
  73     my $this = {};
  74     my $vector;
  75
  76     # we expect $values to have Math::Complex objects
  77     @$values = map { gsl_complex_rect(Re($_), Im($_)) } @$values;
  78
  79     if ( ref $values eq 'ARRAY' ){
  80         die __PACKAGE__.'::new($x) - $x must be a nonempty array reference' if $length == -1;
  81         $vector  = gsl_vector_complex_alloc($length+1);
  82         map { gsl_vector_complex_set($vector, $_, $values->[$_] ) }  (0 .. $length);
  83         $this->{_length} = $length+1;
  84     } elsif ( (int($values) == $values) && ($values > 0)) {
  85         $vector  = gsl_vector_complex_alloc($values);
  86         gsl_vector_complex_set_zero($vector);
  87         $this->{_length} = $values;
  88     } else {
  89         die __PACKAGE__.'::new($x) - $x must be an int or array reference';
  90     }
  91     $this->{_vector} = $vector;
  92     bless $this, $class;
  93 }
  94
  95
  96 =head2 raw()
  97
  98 Get the underlying GSL vector object created by SWIG, useful for using gsl_vector_* functions which do not have an OO counterpart.
  99
 100     my $vector    = Math::GSL::VectorComplex->new(3);
 101     my $gsl_vector = $vector->raw;
 102     my $stuff      = gsl_vector_get($gsl_vector, 1);
 103
 104 =cut
 105
 106 sub raw {
 107     my $self = shift;
 108     return $self->{_vector};
 109 }
 110
 111 =head2 min()
 112
 113 Returns the minimum value contained in the vector.
 114
 115    my $vector = Math::GSL::VectorComplex->new([2,4,1]);
 116    my $minimum = $vector->min;
 117
 118 =cut
 119
 120 sub min {
 121     my $self=shift;
 122     return gsl_vector_min($self->raw);
 123 }
 124
 125 =head2 max()
 126
 127 Returns the minimum value contained in the vector.
 128
 129    my $vector = Math::GSL::VectorComplex->new([2,4,1]);
 130    my $maximum = $vector->max;
 131
 132 =cut
 133
 134 sub max {
 135     my $self=shift;
 136     return gsl_vector_max($self->raw);
 137 }
 138
 139 =head2 length()
 140
 141 Returns the number of elements contained in the vector.
 142
 143    my $vector = Math::GSL::VectorComplex->new([2,4,1]);
 144    my $length = $vector->length;
 145
 146 =cut
 147
 148 sub length { my $self=shift; $self->{_length} }
 149
 150 =head2  as_list()
 151
 152 Gets the content of a Math::GSL::Vector object as a Perl list.
 153
 154     my $vector = Math::GSL::VectorComplex->new(3);
 155     ...
 156     my @values = $vector->as_list;
 157 =cut
 158
 159 sub as_list {
 160     my $self=shift;
 161     # this is wrong
 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_complex_set($self->{_vector}, $indices->[$_], $values->[$_] ) } (0..$#$indices);
 222     };
 223     # better error handling?
 224     warn $@ if $@;
 225     return;
 226 }
 227
 228 =head2 copy()
 229
 230 Returns a copy of the vector, which has the same length and values but resides at a different location in memory.
 231
 232     my $vector = Math::GSL::VectorComplex->new([10 .. 20]);
 233     my $copy   = $vector->copy;
 234
 235 =cut
 236
 237 sub copy {
 238     my $self = shift;
 239     my $copy = Math::GSL::VectorComplex->new( $self->length );
 240     my $status = gsl_vector_complex_memcpy($copy->raw, $self->raw);
 241     if ( $status != $GSL_SUCCESS ) {
 242         croak "Math::GSL - error copying memory, aborting. $! status=$status";
 243     }
 244     return $copy;
 245 }
 246
 247 =head2 swap()
 248
 249 Exchanges the values in the vectors $v with $w by copying.
 250
 251     my $v = Math::GSL::VectorComplex->new([1..5]);
 252     my $w = Math::GSL::VectorComplex->new([3..7]);
 253     $v->swap( $w );
 254
 255 =cut
 256
 257 sub swap() {
 258     my ($self,$other) = @_;
 259     croak "Math::GSL::VectorComplex : \$v->swap(\$w) - \$w must be a Math::GSL::VectorComplex"
 260         unless ref $other eq 'Math::GSL::VectorComplex';
 261     gsl_vector_complex_swap( $self->raw, $other->raw );
 262     return $self;
 263 }
 264
 265 sub _multiplication {
 266     my ($left,$right) = @_;
 267     my $lcopy = $left->copy;
 268
 269     if ( blessed $right && $right->isa(__PACKAGE__) ) {
 270         return $lcopy->dot_product($right);
 271     } else {
 272         # will be in upcoming gsl 1.12
 273         # gsl_vector_complex_scale($lcopy->raw, $right);
 274     }
 275     return $lcopy;
 276 }
 277
 278 sub _subtract {
 279     my ($left, $right, $flip) = @_;
 280
 281     if ($flip) {
 282         my $lcopy = $left->copy;
 283         # will be in upcoming gsl 1.12
 284         # gsl_vector_complex_scale($lcopy->raw, -1 );
 285         gsl_vector_add_constant($lcopy->raw, $right);
 286         return $lcopy;
 287     } else {
 288         return _addition($left, -1.0*$right);
 289     }
 290 }
 291
 292 sub _addition {
 293     my ($left, $right, $flip) = @_;
 294     my $lcopy = $left->copy;
 295
 296     if ( blessed $right && $right->isa('Math::GSL::Vector') && blessed $left && $left->isa('Math::GSL::Vector') ) {
 297         if ( $left->length == $right->length ) {
 298             gsl_vector_complex_add($lcopy->raw, $right->raw);
 299         } else {
 300             croak "Math::GSL - addition of vectors must be called with two objects vectors and must have the same length";
 301         }
 302     } else {
 303         gsl_vector_complex_add_constant($lcopy->raw, $right);
 304     }
 305     return $lcopy;
 306 }
 307
 308 sub dot_product_pp {
 309     my ($left,$right) = @_;
 310     my $sum=0;
 311     if ( blessed $right && $right->isa('Math::GSL::Vector') && $left->length == $right->length ) {
 312          my @l = $left->as_list;
 313          my @r = $right->as_list;
 314          map { $sum += $l[$_] * $r[$_] } (0..$#l);
 315         return $sum;
 316     } else {
 317         croak "dot_product() must be called with two vectors";
 318     }
 319 }
 320
 321 sub dot_product {
 322     my ($left,$right) = @_;
 323
 324     my ($status, $product) = gsl_blas_ddot($left->raw,$right->raw);
 325     croak sprintf "Math::GSL::dot_product - %s", gsl_strerror($status) if ($status != $GSL_SUCCESS);
 326     return $product;
 327 }
 328
 329 =head1 AUTHORS
 330
 331 Jonathan Leto <jonathan@leto.net> and Thierry Moisan <thierry.moisan@gmail.com>
 332
 333 =head1 COPYRIGHT AND LICENSE
 334
 335 Copyright (C) 2008-2009 Jonathan Leto and Thierry Moisan
 336
 337 This program is free software; you can redistribute it and/or modify it
 338 under the same terms as Perl itself.
 339
 340 =cut
 341
 342 %}