view release on metacpan or  search on metacpan

lib/Bio/Metabolic/MatrixOps.pm  view on Meta::CPAN

Returns the kernel of matrix $m, i.e. the matrix with linearly independent
column vectors $c satisfying the equation $m x $c = 0.

=cut

sub PDL::Matrix::kernel {
    my $matrix = shift->copy;

    my ( $rem, $perm, $rank ) = $matrix->row_echelon_int();

    my ( $cols, $rows ) = $rem->dims();

    my $vec;
    my @veclist = ();

    my ( $cnt, $r );
    for ( $cnt = $rank ; $cnt < $cols ; $cnt++ ) {

        #    print "Solution number ".eval($cnt-$rank+1).":\n";
        $vec = zeroes($cols);
        set $vec, $cnt, 1;

        for ( $r = $rank - 1 ; $r >= 0 ; $r-- ) {

            #      print "Calculating row $r:\n";
            my ($row) = $rem->slice("($r),:");
            my $sum = inner( $row, $vec );

            # ensure only integers remain
            my ( $gcd, $redsum, $redpivot ) =
              euclid( $sum, $rem->at( $r, $r ) );
            $vec *= $redpivot;
            set $vec, $r, -( $sum * $redpivot / $rem->at( $r, $r ) );

            #      print "$vec\n";
        }

        push( @veclist, $vec );
    }

    #  print "kernel: Rank=$rank,Columns=$cols,vecs:@veclist\n";
    #  my $retmat = cat(@veclist)->xchg(0,1);
    #  $retmat->permrows($perm);
    #  return $retmat;
    return null unless @veclist;
    my $retmat = bless cat(@veclist)->xchg( 0, 1 ), ref($matrix);
    return $retmat->permrows($perm);
}

=head2 method invert # probably obsolete!!!! Check with PDL::Matrix / PDL::MatrixOps

$inv = $m->invert;
Returns the inverse of $m, undef if $m is not invertible.

=cut

sub PDL::Matrix::invert {
    my $matrix = shift->copy;

    #  print "ref(matrix) is ".ref($matrix)."\n";
    #  my ($pkg,$file,$line)=caller();

    my @dims = $matrix->dims;

    if ( $dims[0] != $dims[1] ) {
        croak("Inverse for non-square matrix not defined!\n");
        return undef;
    }

    my $n = $dims[0];    # n x n matrix

    # $inverse = 1
    #  print "matrix $matrix";
    #  print "invert called from package $pkg, file $file, line $line (n=$n)\n";
    #  my $inverse = ref($matrix)->new($n,$n);
    my $inverse = mzeroes( $n, $n );

    #  (my $tmp = $inverse->diagonal(0,1)) .= 1;   # some strange error!!!!
    #  for (my $del=0;$del<$n;$del++) {set($inverse,$del,$del,1)};
    $inverse->diagonal( 0, 1 ) .= 1;

    for ( my $colnr = 0 ; $colnr < $n ; $colnr++ ) {

        # find pivot element
        my $colnz = which( $matrix->slice(":,($colnr)")->sever != 0 );
        my $ppiv  = which( $colnz >= $colnr );
        if ( $ppiv->nelem == 0 ) {
            print "Matrix not invertible\n";
            return undef;
        }
        my $pivotrow = $colnz->at( $ppiv->at(0) );

        if ( $pivotrow != $colnr ) {
            $matrix->xrow( $colnr,  $pivotrow );
            $inverse->xrow( $colnr, $pivotrow );
        }

        my $akk = $matrix->at( $colnr, $colnr );    # this is a_{kk}
        for ( my $rownr = 0 ; $rownr < $n ; $rownr++ ) {
            if ( $rownr != $colnr ) {

                #	my $q = - $matrix->at($colnr,$rownr) / $akk;
                my $q = -$matrix->at( $rownr, $colnr ) / $akk;
                $matrix->addrows( $rownr,  $colnr, $q );
                $inverse->addrows( $rownr, $colnr, $q );
            }
        }

        #    print "Matrix now (column $colnr) : $matrix\n";
        #    print "Inverse now (column $colnr) : $inverse\n";

    }

    #  my $diag = $matrix->diagonal(0,1);
    #  if (!(which($diag==0)->isempty)) {
    #    print "Matrix non inversible!\n";
    #    return undef;
    #  }

    for ( my $d = 0 ; $d < $n ; $d++ ) {
        if ( $matrix->at( $d, $d ) == 0 ) {