%perlcode %{
use Scalar::Util 'blessed';
use Data::Dumper;
use Carp qw/croak/;
use Math::GSL::Errno qw/$GSL_SUCCESS gsl_strerror/;
use Math::GSL::BLAS qw/gsl_blas_ddot/;
use Math::GSL::Test qw/is_similar/;

use overload
    '*'      => \&_multiplication,
    '+'      => \&_addition,
    '-'      => \&_subtract,
    'abs'    => \&_abs,
    '=='     => \&_equal,
    '!='     => \&_not_equal,
    fallback => 1,
;

@EXPORT_all  = qw/fopen fclose
                 gsl_vector_alloc gsl_vector_calloc gsl_vector_alloc_from_block gsl_vector_alloc_from_vector
                 gsl_vector_free gsl_vector_view_array gsl_vector_const_view_array gsl_vector_view_array_with_stride
                 gsl_vector_const_view_array_with_stride gsl_vector_subvector gsl_vector_subvector_wi gsl_vector_subvector_with_stride
                 gsl_vector_const_subvector gsl_vector_const_subvec gsl_vector_get gsl_vector_set
                 gsl_vector_ptr gsl_vector_const_ptr gsl_vector_set_zero gsl_vector_set_all
                 gsl_vector_set_basis gsl_vector_fread gsl_vector_fwrite gsl_vector_fscanf
                 gsl_vector_fprintf gsl_vector_memcpy gsl_vector_reverse gsl_vector_swap
                 gsl_vector_swap_elements gsl_vector_max gsl_vector_min gsl_vector_minmax
                 gsl_vector_max_index gsl_vector_min_index gsl_vector_minmax_index
                 gsl_vector_add gsl_vector_sub gsl_vector_mul gsl_vector_div
                 gsl_vector_scale gsl_vector_add_constant gsl_vector_isnull
                 gsl_vector_ispos gsl_vector_isneg gsl_vector_isnonneg
                 gsl_vector_float_alloc gsl_vector_float_calloc gsl_vector_float_alloc_from_block
                 gsl_vector_float_alloc_from_vector gsl_vector_float_free gsl_vector_float_view_array
                 gsl_vector_float_view_array_with_stride gsl_vector_float_const_view_array gsl_vector_float_const_view_array_with_stride
                 gsl_vector_float_subvector gsl_vector_float_subvector_with_stride gsl_vector_float_const_subvector
                 gsl_vector_float_const_subvector_with_stride gsl_vector_float_get gsl_vector_float_set gsl_vector_float_ptr
                 gsl_vector_float_const_ptr gsl_vector_float_set_zero gsl_vector_float_set_all gsl_vector_float_set_basis
                 gsl_vector_float_fread gsl_vector_float_fwrite gsl_vector_float_fscanf gsl_vector_float_fprintf
                 gsl_vector_float_memcpy gsl_vector_float_reverse gsl_vector_float_swap gsl_vector_float_swap_elements
                 gsl_vector_float_max gsl_vector_float_min gsl_vector_float_minmax gsl_vector_float_max_index gsl_vector_float_min_index
                 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
                 gsl_vector_float_add_constant gsl_vector_float_isnull gsl_vector_float_ispos gsl_vector_float_isneg gsl_vector_float_isnonneg
/;

@EXPORT_file =qw/ fopen fclose/;
@EXPORT_OK = (@EXPORT_all, @EXPORT_file);
%EXPORT_TAGS = ( file => \@EXPORT_file, all => \@EXPORT_all );

=encoding utf8

=head1 NAME

Math::GSL::Vector - Functions concerning vectors

=head1 SYNOPSIS

    use Math::GSL::Vector qw/:all/;
    my $vec1 = Math::GSL::Vector->new([1, 7, 94, 15 ]);
    my $vec2 = $vec1 * 5;
    my $vec3 = Math::GSL::Vector>new(10);   # 10 element zero vector
    my $vec4 = $vec1 + $vec2;

    # set the element at index 1 to 9
    # and the element at index 3 to 8
    $vec3->set([ 1, 3 ], [ 9, 8 ]);

    my @vec = $vec2->as_list;               # return elements as Perl list

    my $dot_product = $vec1 * $vec2;
    my $length      = $vec2->length;
    my $first       = $vec1->get(0);

    # access raw GSL object to call low-level functions
    my $raw         = $vec1->raw;
    my $element     = gsl_vector_get($raw, 2);

=cut

=head1 Objected Oriented Interface to GSL Math::GSL::Vector

=head2 Math::GSL::Vector->new()

Creates a new Vector of the given size.

    my $vector = Math::GSL::Vector->new(3);

You can also create and set directly the values of the vector like this :

   my $vector = Math::GSL::Vector->new([2,4,1]);

=cut

sub new {
    my ($class, $values) = @_;
    my $length  = $#$values;
    my $this = {};
    my $vector;
    if ( ref $values eq 'ARRAY' ){
        die __PACKAGE__.'::new($x) - $x must be a nonempty array reference' if $length == -1;
        $vector  = gsl_vector_alloc($length+1);
        map { gsl_vector_set($vector, $_, $values->[$_] ) }  (0 .. $length);
        $this->{_length} = $length+1;
    } elsif ( (int($values) == $values) && ($values > 0)) {
        $vector  = gsl_vector_alloc($values);
        gsl_vector_set_zero($vector);
        $this->{_length} = $values;
    } else {
        die __PACKAGE__.'::new($x) - $x must be an int or array reference';
    }
    $this->{_vector} = $vector;
    bless $this, $class;
}
=head2 raw()

Get the underlying GSL vector object created by SWIG, useful for using gsl_vector_* functions which do not have an OO counterpart.

    my $vector    = Math::GSL::Vector->new(3);
    my $gsl_vector = $vector->raw;
    my $stuff      = gsl_vector_get($gsl_vector, 1);

=cut

sub raw {
    my $self = shift;
    return $self->{_vector};
}

=head2 swap()

Exchanges the values in the vectors $v with $w by copying.

    my $v = Math::GSL::Vector->new([1..5]);
    my $w = Math::GSL::Vector->new([3..7]);
    $v->swap( $w );

=cut

sub swap() {
    my ($self,$other) = @_;
    croak "Math::GSL::Vector: \$v->swap(\$w) : \$w must be a Math::GSL::Vector"
        unless ref $other eq 'Math::GSL::Vector';
    gsl_vector_swap( $self->raw, $other->raw );
    return $self;
}

=head2 reverse()

Reverse the elements in the vector.

    $v->reverse;

=cut

sub reverse() {
    my $self = shift;
    gsl_vector_reverse($self->raw);
    return $self;
}

=head2 min()

Returns the minimum value contained in the vector.

   my $vector = Math::GSL::Vector->new([2,4,1]);
   my $minimum = $vector->min;

=cut

sub min {
    my $self=shift;
    return gsl_vector_min($self->raw);
}

=head2 max()

Returns the minimum value contained in the vector.

   my $vector = Math::GSL::Vector->new([2,4,1]);
   my $maximum = $vector->max;

=cut

sub max {
    my $self=shift;
    return gsl_vector_max($self->raw);
}

=head2 length()

Returns the number of elements contained in the vector.

   my $vector = Math::GSL::Vector->new([2,4,1]);
   my $length = $vector->length;

=cut

sub length { my $self=shift; $self->{_length} }

=head2 $v->norm($p)

Returns the p-norm of $v, which defaults to the Euclidean (p=2) norm when no argument is given.

    my $euclidean_distance = $v->norm;

=cut

sub norm($;$)
{
    my ($self,$p) = @_;
    my $norm = 0;
    $p ||= 2;

    map { $norm += $p == 1 ? abs $_ : $_ ** $p } ($self->as_list);
    return $norm **  (1 / $p);
}

=head2 normalize($p)

Divide each element of a vector by its norm, hence creating a unit vector.
Returns the vector for chaining.  If you just want the value of the norm
without changing the vector, use C<norm()>. The default value for C<$p> is 2,
which gives the familiar Euclidean distance norm.

    my $unit_vector = $vector->normalize(2);

is the same as

    my $unit_vector = $vector->normalize;

=cut

sub normalize($;$)
{
    my ($self,$p) = @_;
    $p ||= 2;
    my $norm = $self->norm($p);
    return $self if ($norm == 0);

    my $status = gsl_vector_scale( $self->raw, 1/$norm );
    croak "Math::GSL::Vector - could not scale vectr" unless $status == $GSL_SUCCESS;
    return $self;
}

=head2  as_list()

Gets the content of a Math::GSL::Vector object as a Perl list.

    my $vector = Math::GSL::Vector->new(3);
    ...
    my @values = $vector->as_list;

=cut

sub as_list {
    my $self=shift;
    $self->get( [ 0 .. $self->length - 1  ] );
}

=head2  get()

Gets the value of an of a Math::GSL::Vector object.

    my $vector = Math::GSL::Vector->new(3);
    ...
    my @values = $vector->get(2);

You can also enter an array of indices to receive their corresponding values:

    my $vector = Math::GSL::Vector->new(3);
    ...
    my @values = $vector->get([0,2]);

=cut

sub get {
    my ($self, $indices) = @_;

    return gsl_vector_get($self->raw, $indices) unless ref $indices;

    return  map {  gsl_vector_get($self->raw, $_ ) } @$indices ;
}

=head2  set()

Sets values of an of a Math::GSL::Vector object.

    my $vector = Math::GSL::Vector->new(3);
    $vector->set([1,2], [8,23]);

This sets the second and third value to 8 and 23.

=cut

sub set {
    my ($self, $indices, $values) = @_;
    die (__PACKAGE__.'::set($indices, $values) - $indices and $values must be array references of the same length')
        unless ( ref $indices eq 'ARRAY' && ref $values eq 'ARRAY' &&  $#$indices == $#$values );
    eval {
        map {  gsl_vector_set($self->{_vector}, $indices->[$_], $values->[$_] ) } (0..$#$indices);
    };
    return;
}

=head2 copy()

Returns a copy of the vector, which has the same length and values but resides at a different location in memory.

    my $vector = Math::GSL::Vector->new([10 .. 20]);
    my $copy   = $vector->copy;

=cut


sub copy {
    my $self = shift;
    my $copy = Math::GSL::Vector->new( $self->length );
    if ( gsl_vector_memcpy($copy->raw, $self->raw) != $GSL_SUCCESS ) {
        croak "Math::GSL - error copying memory, aborting";
    }
    return $copy;
}

sub _multiplication {
    my ($left,$right) = @_;
    my $lcopy = $left->copy;

    if ( blessed $right && $right->isa('Math::GSL::Vector') ) {
        return $lcopy->dot_product($right);
    } else {
        gsl_vector_scale($lcopy->raw, $right);
    }
    return $lcopy;
}

sub _subtract {
    my ($left, $right, $flip) = @_;

    if ($flip) {
        my $lcopy = $left->copy;
        gsl_vector_scale($lcopy->raw, -1 );
        gsl_vector_add_constant($lcopy->raw, $right);
        return $lcopy;
    } else {
        return _addition($left, -1.0*$right);
    }
}

sub _abs {
    my $self = shift;
    $self->norm;
}

sub _addition {
    my ($left, $right, $flip) = @_;

    my $lcopy = $left->copy;

    if ( blessed $right && $right->isa('Math::GSL::Vector') && blessed $left && $left->isa('Math::GSL::Vector') ) {
        if ( $left->length == $right->length ) {
            gsl_vector_add($lcopy->raw, $right->raw);
        } else {
            croak "Math::GSL - addition of vectors must be called with two objects vectors and must have the same length";
        }
    } else {
        gsl_vector_add_constant($lcopy->raw, $right);
    }
    return $lcopy;
}

sub dot_product_pp {
    my ($left,$right) = @_;
    my $sum=0;
    if ( blessed $right && $right->isa('Math::GSL::Vector') &&
         $left->length == $right->length ) {
         my @l = $left->as_list;
         my @r = $right->as_list;
         map { $sum += $l[$_] * $r[$_] } (0..$#l);
        return $sum;
    } else {
        croak "dot_product() must be called with two vectors";
    }
}

sub dot_product {
    my ($left,$right) = @_;

    my ($status, $product) = gsl_blas_ddot($left->raw,$right->raw);
    croak sprintf "Math::GSL::dot_product - %s", gsl_strerror($status) if ($status != $GSL_SUCCESS);
    return $product;
}

sub _equal {
    my ($left,$right) = @_;

    return 0 if ($left->length != $right->length);

    return is_similar(  [$left->as_list ], [$right->as_list ]);
}

sub _not_equal {
    my ($left, $right) = @_;
    return !_equal($left,$right);
}

=head1 DESCRIPTION

Here is a list of all the functions included in this module :

=over 1

=item C<gsl_vector_alloc($x)>

Create a vector of size $x

=item C<gsl_vector_calloc($x)>

Create a vector of size $x and initializes all the elements of the vector to zero

=item C<gsl_vector_alloc_from_block>

=item C<gsl_vector_alloc_from_vector>

=item C<gsl_vector_free($v)>

Free a previously allocated vector $v

=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.

=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.

=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}.

=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.

=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

=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

=item C<gsl_vector_const_subvector>


=item C<gsl_vector_get($v, $i)>

Return the $i-th element of a raw vector $v, where $v is a Math::GSL::Vector object. For example:

    my $third_element = gsl_vector_get($v->raw, 3);

=item C<gsl_vector_set($v, $i, $x)>

Return the vector $v with his $i-th element set to $x

=item C<gsl_vector_ptr>

=item C<gsl_vector_const_ptr>

=item C<gsl_vector_set_zero($v)>

 set all the elements of $v to 0

=item C<gsl_vector_set_all($v, $x)>

set all the elements of $v to $x

=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 succeeded, 1 otherwise.

=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.

=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.

=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.

=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.

=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 operation succeeded, 1 otherwise. The two vectors must have the
same length.

=item C<gsl_vector_reverse($v)>

reverse the order of the elements of the vector $v and return 0 if the
operation succeeded, 1 otherwise

=item C<gsl_vector_swap($v, $v2)>

swap the values of the vectors $v and $v2 and return 0 if the operation
succeeded, 1 otherwise

=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 succeeded, 1 otherwise.

=item C<gsl_vector_max($v)>

 return the maximum value in the vector $v

=item C<gsl_vector_min($v)>

return the minimum value in the vector $v

=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.

=item C<gsl_vector_max_index($v)>

return the position of the maximum value in the vector $v

=item C<gsl_vector_min_index($v)>

return the position of the minimum value in the vector $v

=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.

=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 succeeded, 1 otherwise.

=item C<gsl_vector_sub($v, $v2)>

subtract the elements of $v2 from the elements of $v, the two vectors must
have the same length and return 0 if the operation succeeded, 1 otherwise.

=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 succeeded, 1 otherwise.

=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 succeeded, 1 otherwise.

=item C<gsl_vector_scale($v, $x)>

multiplty the elements of the vector $v by a constant $x and return 0 if the
operation succeeded, 1 otherwise.

=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 succeeded, 1 otherwise.

=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.

=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.

=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.

=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.

=back

Precision on the vector_view type : every modification you'll make on a
vector_view will also modify the original vector.  For example, the following
code will zero the even elements of the vector $v of length $size, while
leaving the odd elements untouched :

=over 1

=item C<$v_even= gsl_vector_subvector_with_stride ($v, 0, 2, $size/2);>

=item C<gsl_vector_set_zero ($v_even-E<gt>{vector});>

=back

For more information on the functions, we refer you to the GSL official documentation:
L<http://www.gnu.org/software/gsl/manual/html_node/>

=head1 EXAMPLES

Here is an example using both interfaces.

    use Math::GSL::Vector qw/:all/;

    print "We'll create this vector : [0,1,4,9,16] \n";
    my $vector = Math::GSL::Vector->new([0,1,4,9,16]);
    my ($min, $max) = gsl_vector_minmax_index($vector->raw);

    print "We then check the index value of the maximum and minimum values of the vector. \n";
    print "The index of the maximum should be 4 and we received $max \n";
    print "The index of the minimum should be 0 and we received $min \n";
    print "We'll then swap the first and the third elements of the vector \n";

    gsl_vector_swap_elements($vector->raw, 0, 3);
    my @got = $vector->as_list;
    print "The vector should now be like this : [9,1,4,0,16] \n";
    print "and we received : [ @got ]\n";

=head1 AUTHORS

Jonathan "Duke" Leto <jonathan@leto.net> and Thierry Moisan <thierry.moisan@gmail.com>

=head1 COPYRIGHT AND LICENSE

Copyright (C) 2008-2024 Jonathan "Duke" Leto and Thierry Moisan

This program is free software; you can redistribute it and/or modify it
under the same terms as Perl itself.

=cut
%}