2 use Scalar::Util 'blessed';
5 use Math::GSL::Errno qw/$GSL_SUCCESS/;
6 use Math::GSL::BLAS qw/gsl_blas_ddot/;
7 use Math::GSL::Test qw/is_similar/;
10 '*' => \&_multiplication,
19 @EXPORT_all = qw/fopen fclose
20 gsl_vector_alloc gsl_vector_calloc gsl_vector_alloc_from_block gsl_vector_alloc_from_vector
21 gsl_vector_free gsl_vector_view_array gsl_vector_const_view_array gsl_vector_view_array_with_stride
22 gsl_vector_const_view_array_with_stride gsl_vector_subvector gsl_vector_subvector_wi gsl_vector_subvector_with_stride
23 gsl_vector_const_subvec gsl_vector_const_subvec gsl_vector_get gsl_vector_set
24 gsl_vector_ptr gsl_vector_const_ptr gsl_vector_set_zero gsl_vector_set_all
25 gsl_vector_set_basis gsl_vector_fread gsl_vector_fwrite gsl_vector_fscanf
26 gsl_vector_fprintf gsl_vector_memcpy gsl_vector_reverse gsl_vector_swap
27 gsl_vector_swap_elements gsl_vector_max gsl_vector_min gsl_vector_minmax
28 gsl_vector_max_index gsl_vector_min_index gsl_vector_minmax_index
29 gsl_vector_add gsl_vector_sub gsl_vector_mul gsl_vector_div
30 gsl_vector_scale gsl_vector_add_constant gsl_vector_isnull
31 gsl_vector_ispos gsl_vector_isneg gsl_vector_isnonneg
32 gsl_vector_float_alloc gsl_vector_float_calloc gsl_vector_float_alloc_from_block
33 gsl_vector_float_alloc_from_vector gsl_vector_float_free gsl_vector_float_view_array
34 gsl_vector_float_view_array_with_stride gsl_vector_float_const_view_array gsl_vector_float_const_view_array_with_stride
35 gsl_vector_float_subvector gsl_vector_float_subvector_with_stride gsl_vector_float_const_subvector
36 gsl_vector_float_const_subvector_with_stride gsl_vector_float_get gsl_vector_float_set gsl_vector_float_ptr
37 gsl_vector_float_const_ptr gsl_vector_float_set_zero gsl_vector_float_set_all gsl_vector_float_set_basis
38 gsl_vector_float_fread gsl_vector_float_fwrite gsl_vector_float_fscanf gsl_vector_float_fprintf
39 gsl_vector_float_memcpy gsl_vector_float_reverse gsl_vector_float_swap gsl_vector_float_swap_elements
40 gsl_vector_float_max gsl_vector_float_min gsl_vector_float_minmax gsl_vector_float_max_index gsl_vector_float_min_index
41 gsl_vector_float_minmax_index gsl_vector_float_add gsl_vector_float_sub gsl_vector_float_mul gsl_vector_float_div gsl_vector_float_scale
42 gsl_vector_float_add_constant gsl_vector_float_isnull gsl_vector_float_ispos gsl_vector_float_isneg gsl_vector_float_isnonneg
45 @EXPORT_file =qw/ fopen fclose/;
46 @EXPORT_OK = (@EXPORT_all, @EXPORT_file);
47 %EXPORT_TAGS = ( file => \@EXPORT_file, all => \@EXPORT_all );
51 Math::GSL::Vector - Functions concerning vectors
55 use Math::GSL::Vector qw/:all/;
56 my $vec1 = Math::GSL::Vector->new([1, 7, 94, 15 ]);
58 my $vec3 = Math::GSL::Vector>new(10); # 10 element zero vector
59 my $vec4 = $vec1 + $vec2;
61 # set the element at index 1 to 9
62 # and the element at index 3 to 8
63 $vec3->set([ 1, 3 ], [ 9, 8 ]);
65 my @vec = $vec2->as_list; # return elements as Perl list
67 my $dot_product = $vec1 * $vec2;
68 my $length = $vec2->length;
69 my $first = $vec1->get(0);
74 =head1 Objected Oriented Interface to GSL Math::GSL::Vector
76 =head2 Math::GSL::Vector->new()
78 Creates a new Vector of the given size.
80 my $vector = Math::GSL::Vector->new(3);
82 You can also create and set directly the values of the vector like this :
84 my $vector = Math::GSL::Vector->new([2,4,1]);
89 my ($class, $values) = @_;
90 my $length = $#$values;
93 if ( ref $values eq 'ARRAY' ){
94 die __PACKAGE__.'::new($x) - $x must be a nonempty array reference' if $length == -1;
95 $vector = gsl_vector_alloc($length+1);
96 map { gsl_vector_set($vector, $_, $values->[$_] ) } (0 .. $length);
97 $this->{_length} = $length+1;
98 } elsif ( (int($values) == $values) && ($values > 0)) {
99 $vector = gsl_vector_alloc($values);
100 gsl_vector_set_zero($vector);
101 $this->{_length} = $values;
103 die __PACKAGE__.'::new($x) - $x must be an int or array reference';
105 $this->{_vector} = $vector;
110 Get the underlying GSL vector object created by SWIG, useful for using gsl_vector_* functions which do not have an OO counterpart.
112 my $vector = Math::GSL::Vector->new(3);
113 my $gsl_vector = $vector->raw;
114 my $stuff = gsl_vector_get($gsl_vector, 1);
120 return $self->{_vector};
125 Exchanges the values in the vectors $v with $w by copying.
127 my $v = Math::GSL::Vector->new([1..5]);
128 my $w = Math::GSL::Vector->new([3..7]);
134 my ($self,$other) = @_;
135 croak "Math::GSL::Vector: \$v->swap(\$w) : \$w must be a Math::GSL::Vector"
136 unless ref $other eq 'Math::GSL::Vector';
137 gsl_vector_swap( $self->raw, $other->raw );
143 Reverse the elements in the vector.
151 gsl_vector_reverse($self->raw);
157 Returns the minimum value contained in the vector.
159 my $vector = Math::GSL::Vector->new([2,4,1]);
160 my $minimum = $vector->min;
166 return gsl_vector_min($self->raw);
171 Returns the minimum value contained in the vector.
173 my $vector = Math::GSL::Vector->new([2,4,1]);
174 my $maximum = $vector->max;
180 return gsl_vector_max($self->raw);
185 Returns the number of elements contained in the vector.
187 my $vector = Math::GSL::Vector->new([2,4,1]);
188 my $length = $vector->length;
192 sub length { my $self=shift; $self->{_length} }
196 Returns the p-norm of $v, which defaults to the Euclidean (p=2) norm when no argument is given.
198 my $euclidean_distance = $v->norm;
208 map { $norm += $p == 1 ? abs $_ : $_ ** $p } ($self->as_list);
209 return $norm ** (1 / $p);
214 Divide each element of a vector by it's norm, hence creating a unit vector. Returns the vector for chaining.
215 If you just want the value of the norm without changing the vector, use C<norm()>. The default value for C<$p>
216 is 2, which gives the familiar Euclidean distance norm.
218 my $unit_vector = $vector->normalize(2);
222 my $unit_vector = $vector->normalize;
230 my $norm = $self->norm($p);
231 return $self if ($norm == 0);
233 my $status = gsl_vector_scale( $self->raw, 1/$norm );
234 croak "Math::GSL::Vector - could not scale vectr" unless $status == $GSL_SUCCESS;
240 Gets the content of a Math::GSL::Vector object as a Perl list.
242 my $vector = Math::GSL::Vector->new(3);
244 my @values = $vector->as_list;
250 $self->get( [ 0 .. $self->length - 1 ] );
255 Gets the value of an of a Math::GSL::Vector object.
257 my $vector = Math::GSL::Vector->new(3);
259 my @values = $vector->get(2);
261 You can also enter an array of indices to receive their corresponding values:
263 my $vector = Math::GSL::Vector->new(3);
265 my @values = $vector->get([0,2]);
270 my ($self, $indices) = @_;
271 return map { gsl_vector_get($self->raw, $_ ) } @$indices ;
276 Sets values of an of a Math::GSL::Vector object.
278 my $vector = Math::GSL::Vector->new(3);
279 $vector->set([1,2], [8,23]);
281 This sets the second and third value to 8 and 23.
286 my ($self, $indices, $values) = @_;
287 die (__PACKAGE__.'::set($indices, $values) - $indices and $values must be array references of the same length')
288 unless ( ref $indices eq 'ARRAY' && ref $values eq 'ARRAY' && $#$indices == $#$values );
290 map { gsl_vector_set($self->{_vector}, $indices->[$_], $values->[$_] ) } (0..$#$indices);
297 Returns a copy of the vector, which has the same length and values but resides at a different location in memory.
299 my $vector = Math::GSL::Vector->new([10 .. 20]);
300 my $copy = $vector->copy;
307 my $copy = Math::GSL::Vector->new( $self->length );
308 if ( gsl_vector_memcpy($copy->raw, $self->raw) != $GSL_SUCCESS ) {
309 croak "Math::GSL - error copying memory, aborting";
314 sub _multiplication {
315 my ($left,$right) = @_;
316 my $lcopy = $left->copy;
318 if ( blessed $right && $right->isa('Math::GSL::Vector') ) {
319 return $lcopy->dot_product($right);
321 gsl_vector_scale($lcopy->raw, $right);
327 my ($left, $right, $flip) = @_;
330 my $lcopy = $left->copy;
331 gsl_vector_scale($lcopy->raw, -1 );
332 gsl_vector_add_constant($lcopy->raw, $right);
335 return _addition($left, -1.0*$right);
345 my ($left, $right, $flip) = @_;
347 my $lcopy = $left->copy;
349 if ( blessed $right && $right->isa('Math::GSL::Vector') && blessed $left && $left->isa('Math::GSL::Vector') ) {
350 if ( $left->length == $right->length ) {
351 gsl_vector_add($lcopy->raw, $right->raw);
353 croak "Math::GSL - addition of vectors must be called with two objects vectors and must have the same length";
356 gsl_vector_add_constant($lcopy->raw, $right);
362 my ($left,$right) = @_;
364 if ( blessed $right && $right->isa('Math::GSL::Vector') &&
365 $left->length == $right->length ) {
366 my @l = $left->as_list;
367 my @r = $right->as_list;
368 map { $sum += $l[$_] * $r[$_] } (0..$#l);
371 croak "dot_product() must be called with two vectors";
376 my ($left,$right) = @_;
378 my ($status, $product) = gsl_blas_ddot($left->raw,$right->raw);
379 croak sprintf "Math::GSL::dot_product - %s", gsl_strerror($status) if ($status != $GSL_SUCCESS);
384 my ($left,$right) = @_;
386 return 0 if ($left->length != $right->length);
388 return is_similar( [$left->as_list ], [$right->as_list ]);
392 my ($left, $right) = @_;
393 return !_equal($left,$right);
398 Here is a list of all the functions included in this module :
402 =item C<gsl_vector_alloc($x)> - create a vector of size $x
404 =item C<gsl_vector_calloc($x)> - create a vector of size $x and initializes all the elements of the vector to zero
406 =item C<gsl_vector_alloc_from_block>
408 =item C<gsl_vector_alloc_from_vector>
410 =item C<gsl_vector_free($v)> - free a previously allocated vector $v
412 =item C<gsl_vector_view_array($base, $n)> - This function returns a vector view of an array reference $base. The start of the new vector is given by $base and has $n elements. Mathematically, the i-th element of the new vector v' is given by, v'(i) = $base->[i] where the index i runs from 0 to $n-1. The array containing the elements of v is not owned by the new vector view. When the view goes out of scope the original array will continue to exist. The original memory can only be deallocated by freeing the original pointer base. Of course, the original array should not be deallocated while the view is still in use.
414 =item C<gsl_vector_const_view_array($base, $n)> - This function is equivalent to gsl_vector_view_array but can be used for arrays which are declared const.
416 =item C<gsl_vector_view_array_with_stride($base, $stride, $n)> - This function returns a vector view of an array reference $base with an additional $stride argument. The subvector is formed in the same way as for gsl_vector_view_array but the new vector has $n elements with a step-size of $stride from one element to the next in the original array. Mathematically, the i-th element of the new vector v' is given by, v'(i) = $base->[i*$stride] where the index i runs from 0 to $n-1. Note that the view gives direct access to the underlying elements of the original array. A vector view $view can be passed to any subroutine which takes a vector argument just as a directly allocated vector would be, using $view->{vector}.
418 =item C<gsl_vector_const_view_array_with_stride($base, $stride, $n)> - This function is equivalent to gsl_vector_view_array_with_stride but can be used for arrays which are declared const.
420 =item C<gsl_vector_subvector($v, $offset, $n)> - return a vector_view type which contains a subvector of $v, with a size of $size, starting from the $offset position
422 =item C<gsl_vector_subvector_with_stride($v, $offset, $stride, $size)> - return a vector_view type which contains a subvector of $v, with a size of $size, starting from the $offset position and with a $stride step between each element of $v
424 =item C<gsl_vector_const_subvector>
426 =item C<gsl_vector_get($v, $i)> - return the $i-th element of a vector $v
428 =item C<gsl_vector_set($v, $i, $x)> - return the vector $v with his $i-th element set to $x
430 =item C<gsl_vector_ptr>
432 =item C<gsl_vector_const_ptr>
434 =item C<gsl_vector_set_zero($v)> - set all the elements of $v to 0
436 =item C<gsl_vector_set_all($v, $x)> - set all the elements of $v to $x
438 =item C<gsl_vector_set_basis($v, $i)> - set all the elements of $v to 0 except for the $i-th element which is set to 1 and return 0 if the operation succeded, 1 otherwise.
440 =item C<gsl_vector_fread($file, $v)> - This function reads into the vector $v from the open stream $file opened with gsl_fopen function from the Math::GSL module in binary format. The vector $v must be preallocated with the correct length since the function uses the size of $v to determine how many bytes to read. The return value is 0 for success and 1 if there was a problem reading from the file.
442 =item C<gsl_vector_fwrite($file, $v)> - This function writes the elements of the vector $v to the stream $file opened with gsl_fopen function from the Math::GSL module in binary format. The return value is 0 for success and 1 if there was a problem writing to the file. Since the data is written in the native binary format it may not be portable between different architectures.
444 =item C<gsl_vector_fscanf($file, $v)> This function reads formatted data from the stream $file opened with gsl_fopen function from the Math::GSL module into the vector $v. The vector $v must be preallocated with the correct length since the function uses the size of $v to determine how many numbers to read. The function returns 0 for success and 1 if there was a problem reading from the file.
446 =item C<gsl_vector_fprintf($file, $v, $format)> -This function writes the elements of the vector $v line-by-line to the stream $file opened with gsl_fopen function from the Math::GSL module using the format specifier $format, which should be one of the "%g", "%e" or "%f" formats for floating point numbers and "%d" for integers. The function returns 0 for success and 1 if there was a problem writing to the file.
448 =item C<gsl_vector_memcpy($dest, $src)> - This function copies the elements of the vector $src into the vector $dest and return 0 if the opertaion succeded, 1 otherwise. The two vectors must have the same length.
450 =item C<gsl_vector_reverse($v)> - reverse the order of the elements of the vector $v and return 0 if the opertaion succeded, 1 otherwise
452 =item C<gsl_vector_swap($v, $v2)> - swap the values of the vectors $v and $v2 and return 0 if the opertaion succeded, 1 otherwise
454 =item C<gsl_vector_swap_elements($v, $i, $j)> - permute the elements at position $i and $j in the vector $v and return 0 if the operation succeded, 1 otherwise.
456 =item C<gsl_vector_max($v)> - return the maximum value in the vector $v
458 =item C<gsl_vector_min($v)> - return the minimum value in the vector $v
460 =item C<gsl_vector_minmax($v)> - return two values, the first is the minimum value in the vector $v and the second is the maximum value.
462 =item C<gsl_vector_max_index($v)> - return the position of the maximum value in the vector $v
464 =item C<gsl_vector_min_index($v)> - return the position of the minimum value in the vector $v
466 =item C<gsl_vector_minmax_index> - return two values, the first is the position of the minimum value in the vector $v and the second is the position of the maximum value.
468 =item C<gsl_vector_add($v, $v2)> - add the elements of $v2 to the elements of $v, the two vectors must have the same length and return 0 if the operation succeded, 1 otherwise.
470 =item C<gsl_vector_sub($v, $v2)> - substract the elements of $v2 from the elements of $v, the two vectors must have the same length and return 0 if the operation succeded, 1 otherwise.
472 =item C<gsl_vector_mul($v, $v2)> - multiply the elements of $v by the elements of $v2, the two vectors must have the same length and return 0 if the operation succeded, 1 otherwise.
474 =item C<gsl_vector_div($v, $v2)> - divides the elements of $v by the elements of $v2, the two vectors must have the same length and return 0 if the operation succeded, 1 otherwise.
476 =item C<gsl_vector_scale($v, $x)> - multiplty the elements of the vector $v by a constant $x and return 0 if the operation succeded, 1 otherwise.
478 =item C<gsl_vector_add_constant($v, $x)> - add a constant $x to the elements of the vector $v and return 0 if the operation succeded, 1 otherwise.
480 =item C<gsl_vector_isnull($v)> - verify if all the elements of the vector $v are null, return 0 if it's the case, 1 otherwise.
482 =item C<gsl_vector_ispos($v)> - verify if all the elements of the vector $v are positive, return 0 if it's the case, 1 otherwise.
484 =item C<gsl_vector_isneg($v)> - verify if all the elements of the vector $v are negative, return 0 if it's the case, 1 otherwise.
486 =item C<gsl_vector_isnonneg($v)> - verify if all the elements the vector $v are not negative, return 0 if it's the case, 1 otherwise.
490 Precision on the vector_view type : every modification you'll make on a vector_view will also modify the original vector.
491 For example, the following code will zero the even elements of the vector $v of length $size, while leaving the odd elements untouched :
495 =item C<$v_even= gsl_vector_subvector_with_stride ($v, 0, 2, $size/2);>
497 =item C<gsl_vector_set_zero ($v_even-E<gt>{vector});>
501 For more informations on the functions, we refer you to the GSL offcial documentation:
502 L<http://www.gnu.org/software/gsl/manual/html_node/>
504 Tip : search on google: site:http://www.gnu.org/software/gsl/manual/html_node/ name_of_the_function_you_want
508 Here is an example using both interfaces.
510 use Math::GSL::Vector qw/:all/;
512 print "We'll create this vector : [0,1,4,9,16] \n";
513 my $vector = Math::GSL::Vector->new([0,1,4,9,16]);
514 my ($min, $max) = gsl_vector_minmax_index($vector->raw);
516 print "We then check the index value of the maximum and minimum values of the vector. \n";
517 print "The index of the maximum should be 4 and we received $max \n";
518 print "The index of the minimum should be 0 and we received $min \n";
519 print "We'll then swap the first and the third elements of the vector \n";
521 gsl_vector_swap_elements($vector->raw, 0, 3);
522 my @got = $vector->as_list;
523 print "The vector should now be like this : [9,1,4,0,16] \n";
524 print "and we received : [ @got ]\n";
528 Jonathan Leto <jonathan@leto.net> and Thierry Moisan <thierry.moisan@gmail.com>
530 =head1 COPYRIGHT AND LICENSE
532 Copyright (C) 2008-2009 Jonathan Leto and Thierry Moisan
534 This program is free software; you can redistribute it and/or modify it
535 under the same terms as Perl itself.