File Coverage

blib/lib/Math/Matrix/MaybeGSL.pm

Criterion	Covered	Total	%
statement	115	206	55.8
branch	24	62	38.7
condition	8	19	42.1
subroutine	35	44	79.5
pod	2	2	100.0
total	184	333	55.2

# ABSTRACT: Uniform use of Math::MatrixReal and Math::GSL::Matrix.

use strict;
use warnings;

package Math::Matrix::MaybeGSL;
$Math::Matrix::MaybeGSL::VERSION = '0.008';
use parent 'Exporter';
our @EXPORT = qw{Matrix};

use overload
       '*=' => '_assign_multiply',
        '*' => '_multiply',
       '+=' => '_assign_add',
        '+' => '_add',
 'fallback' =>   undef;

sub _choose_matrix_module {
    return 'Math::GSL::Matrix' if $INC{'Math/GSL/Matrix.pm'};
    return 'Math::MatrixReal'  if $INC{'Math/MatrixReal.pm'};

    my @err;

    return 'Math::GSL::Matrix' if eval {
        require Math::GSL;
        require Math::GSL::Matrix;
        1;
    };
    push @err, "Error loading Math::GSL::Matrix: $@";

    return 'Math::MatrixReal' if eval { require Math::MatrixReal; 1; };
    push @err, "Error loading Math::MatrixReal: $@";

    die join( "\n", "Couldn't load a Matrix module:", @err );
}

sub Matrix { __PACKAGE__ }

sub _call {
    my ($method, $obj, @args) = @_;
    $obj->{matrix}->$method(@args);
}

sub isGSL {
    our $matrix_module;
    return $matrix_module eq "Math::GSL::Matrix";
}

BEGIN {
    our $matrix_module = _choose_matrix_module();
    my %functions
      = (
         'any' => {
                   new => sub {
                       my (undef, $rows, $cols) = @_;
                       return _new( $matrix_module->new($rows, $cols) );
                   },
                   dim     => sub { _call(dim => @_) },
                   each    => sub { _new(_call(each => @_)) },
                   as_list => sub { _call(as_list => @_) },
                   det     => sub { _call(det => @_)} ,
                  },
         'Math::GSL::Matrix' => {
            assign        => sub { _call(set_elem => ($_[0], $_[1]-1, $_[2]-1, $_[3])); },
            element       => sub { _call(get_elem => ($_[0], $_[1]-1, $_[2]-1, $_[3])); },
            new_from_cols => sub { _new(_gsl_new_from_cols($_[1])) },
            new_from_rows => sub { _new(_gsl_new_from_rows($_[1])) },
            vconcat       => sub { _new(_call(vconcat => $_[0], $_[1]{matrix})) },
            hconcat       => sub { _new(_call(hconcat => $_[0], $_[1]{matrix})) },
            write         => sub { _gsl_write(@_) },
            read          => sub { _gsl_read($_[1]) },
            max           => sub {
                if (wantarray) {
                    my ($v, @pos) = _call(max => @_);
                    return ($v, map { $_ + 1 } @pos);
                } else {
                    return scalar(_call(max => @_));
                };
            },
            min           => sub {
                if (wantarray) {
                    my ($v, @pos) = _call(min => @_);
                    return ($v, map { $_ + 1 } @pos);
                } else {
                    return scalar(_call(min => @_));
                };
            },
            row           => sub { _new(_call(row => $_[0], $_[1]-1)) },
            find_zeros    => sub { _gsl_find_zeros(@_) },
            transpose     => sub { _gsl_transpose(@_) },
           },
         'Math::MatrixReal' => {
            assign        => sub { _call(assign        => @_); },
            element       => sub { _call(element       => @_); },
            new_from_cols => sub { _new( $matrix_module->new_from_cols($_[1])) },
            new_from_rows => sub { _new( $matrix_module->new_from_rows($_[1])) },
            vconcat       => sub { _new( ~((~$_[0]{matrix}) . (~$_[1]{matrix})) ) },
            hconcat       => sub { _new(     $_[0]{matrix}  .   $_[1]{matrix}   ) },
            write         => sub { _mreal_write(@_) },
            read          => sub { _mreal_read($_[1]) },
            max           => sub { _mreal_max($_[0]{matrix}) },
            min           => sub { _mreal_min($_[0]{matrix}) },
            row           => sub { _new( $_[0]{matrix}->row($_[1]) ) },
            find_zeros    => sub { _mreal_find_zeros(@_) },
            transpose     => sub { _new( ~$_[0]{matrix} ) },
                               },
	);

    no strict 'refs';

    for my $func (keys %{$functions{$matrix_module}}) {
        # Use Sub::Install later?
        $_ = __PACKAGE__ . "::$func";
        *$_ = $functions{$matrix_module}{$func};
    }
    for my $func (keys %{$functions{any}}) {
        # Use Sub::Install later?
        $_ = __PACKAGE__ . "::$func";
        *$_ = $functions{any}{$func};
    }

}

sub _mreal_max {
    my $matrix = shift;
    my ($rs, $cs) = $matrix->dim();
    return $matrix->max() if ($rs == 1 || $cs == 1);

    my ($m, $r, $c, $v) = ($matrix->[0], 1, 1, undef);

    for my $i (1..$rs) {
        for my $j (1..$cs) {
            if (!$v || $v < $m->[$i-1][$j-1]) {
                $r = $i;
                $c = $j;
                $v = $m->[$i-1][$j-1];
            }
        }
    }

    return wantarray ? ($v, $r, $c) : $v;
}

sub _mreal_min {
    my $matrix = shift;
    my ($rs, $cs) = $matrix->dim();
    return $matrix->min() if ($rs == 1 || $cs == 1);

    my ($m, $r, $c, $v) = ($matrix->[0], 1, 1, undef);

    for my $i (1..$rs) {
        for my $j (1..$cs) {
            if (!$v || $v > $m->[$i-1][$j-1]) {
                $r = $i;
                $c = $j;
                $v = $m->[$i-1][$j-1];
            }
        }
    }

    return wantarray ? ($v, $r, $c) : $v;
}

sub _gsl_new_from_cols {
    my $cols = shift;

    my $nr_columns = scalar(@$cols);
    my $nr_rows = 0;
    for my $row (@$cols) {
        $nr_rows = scalar(@$row) if @$row > $nr_rows;
    }
    my $m = Math::GSL::Matrix->new($nr_rows, $nr_columns);
    for my $r (0..$nr_rows - 1) {
        for my $c (0..$nr_columns - 1) {
            $m->set_elem($r, $c, $cols->[$c][$r] || 0);
        }
    }
    return $m;
}

sub _gsl_new_from_rows {
    my $rows = shift;

    my $nr_rows = scalar(@$rows);
    my $nr_columns = 0;
    for my $col (@$rows) {
        $nr_columns = scalar(@$col) if @$col > $nr_columns;
    }
    my $m = Math::GSL::Matrix->new($nr_rows, $nr_columns);
    for my $c (0..$nr_columns - 1) {
        for my $r (0..$nr_rows - 1) {
            $m->set_elem($r, $c, $rows->[$r][$c] || 0);
        }
    }
    return $m;
}

sub _new {
    my $mat = shift;
    return bless { matrix => $mat }, __PACKAGE__;
}

sub _assign_multiply {
    my($object,$argument) = @_;

    return( &_multiply($object,$argument,undef) );
}

sub _multiply {
    my ($object, $argument, $flag) = @_;

    $argument = $argument->{matrix} if ref $argument eq __PACKAGE__;
    $object   = $object->{matrix}   if ref $object   eq __PACKAGE__;

    if ((defined $flag) && $flag) {
    	return _new($argument * $object);
    } else {
    	return _new($object * $argument);
    }
}

sub _assign_add {
    my($object,$argument) = @_;

    return( &_add($object,$argument) );
}

sub _add {
    my ($object, $argument) = @_;

    $object   = $object->{matrix}   if ref $object   eq __PACKAGE__;
    $argument = $argument->{matrix} if ref $argument eq __PACKAGE__;

    return _new($object + $argument);
}

sub _mreal_write {
    my ($m, $filename) = @_;

    my $matrix = $m->{matrix};

    open my $fh, ">", $filename or
      die "Could not create file '$filename': $!";

    # probably faster than creating a full string in memory
    my ($rows, $cols) = $matrix->dim();

    for my $r (0..$rows-1) {
        for my $c (0..$cols-1) {
            print $fh $matrix->[0][$r][$c];
            print $fh "\t" unless $c == $cols-1;
        }
        print $fh "\n";
    }
    close $fh;
}

sub _mreal_read {
    my $filename = shift;

    my $m = [];

    open my $fh, "<", $filename or
      die "could not open file '$filename': $!";

    while (<$fh>) {
        chomp;
        push @$m, [split /\s+/];
    }

    return _new( Math::MatrixReal->new_from_rows($m) );
}

sub _mreal_find_zeros {
    my ($matrix) = @_;
    my ($rs, $cs) = $matrix->dim();

    my @matches;
    my $pos = 0;
    for ($matrix->as_list()) {
        push @matches, [int($pos/$cs)+1, ($pos % $cs)+1] unless $_;
        $pos++;
    }

    return @matches;
}

sub _gsl_read {
    my $filename = shift;

    die "$filename does not exist" unless -f $filename;

    my $fh = Math::GSL::gsl_fopen($filename, "r");
    die "error opening file $filename for reading" unless $fh;

    my $dim = Math::GSL::Matrix::gsl_matrix_alloc(1, 2);
    my $err = Math::GSL::Matrix::gsl_matrix_fread($fh, $dim);
    die "error reading matrix" if $err;

    my $m = Math::GSL::Matrix::gsl_matrix_alloc(
               Math::GSL::Matrix::gsl_matrix_get($dim, 0, 0),
               Math::GSL::Matrix::gsl_matrix_get($dim, 0, 1));
    $err = Math::GSL::Matrix::gsl_matrix_fread($fh, $m);
    die "error reading matrix" if $err;

    Math::GSL::Matrix::gsl_matrix_free($dim);

    Math::GSL::gsl_fclose($fh);
    _new( Math::GSL::Matrix->new($m) );
}

sub _gsl_write {
    my ($self, $filename) = @_;

    my $fh = Math::GSL::gsl_fopen($filename, "w");
    die "error opening file: $filename" unless $fh;

    # create a temporary matrix with the main matrix dimensions
    my $dim = Math::GSL::Matrix::gsl_matrix_alloc(1, 2);
    my ($rows, $cols) = $self->dim;
    Math::GSL::Matrix::gsl_matrix_set($dim, 0, 0, $rows);
    Math::GSL::Matrix::gsl_matrix_set($dim, 0, 1, $cols);

    my $err = Math::GSL::Matrix::gsl_matrix_fwrite($fh, $dim);
    die "error gsl-writting matrix" if $err;

    Math::GSL::Matrix::gsl_matrix_free($dim);

    $err = Math::GSL::Matrix::gsl_matrix_fwrite($fh, $self->{matrix}->raw);
    die "error gsl-writting matrix" if $err;

    Math::GSL::gsl_fclose($fh);

}

sub _gsl_find_zeros {
    my ($matrix) = @_;
    my ($rs, $cs) = $matrix->dim();

    my $raw_matrix = $matrix->{matrix}->raw;
    my @matches;
    for my $i (0..$rs-1) {
        for my $j (0..$cs-1) {
            next if Math::GSL::Matrix::gsl_matrix_get($raw_matrix, $i, $j);
            push @matches, [$i+1, $j+1];
        }
    }
    return @matches;
}

sub _gsl_transpose {
    my ($matrix) = @_;
    my ($rs, $cs) = $matrix->dim();

    my $result = Math::GSL::Matrix::gsl_matrix_alloc($cs, $rs);
    Math::GSL::Matrix::gsl_matrix_transpose_memcpy($result, $matrix->{matrix}->raw());

    return _new(Math::GSL::Matrix->new($result));
}




1;

__END__