#!/usr/bin/perl -w
##################
#
# Trying out Math::Matrix
#
# I see that
#   * the pieces of the matrix can be retrieved by $m->[row][column]
#   * the $m->print function is more awkward than a $m->asString would be
#   * I can find an inverse matrix by solving [matrix][inverse]=[identity]
#
# Jim on April 7
#
#################
use strict;
use Math::Matrix;

# === playing around ==================================


my $m = new Math::Matrix ([ rand, rand, rand ],
			  [ rand, rand, rand ],
			  [ rand, rand, rand ]);

print "A 3x3 random matrix is\n"; 
$m->print;

print "which has determinant " . $m->determinant . ".\n";

print "and transpose \n";
$m->transpose->print;

my $m_inv = inverse($m);
print "Its inverse is \n";
$m_inv->print;

print "which gives the product of 'm . m_inv' as\n";
print $m->multiply($m_inv);

print "and the product 'm_inv . m' as\n";
print $m_inv->multiply($m);

print "The acccuracy of the 'm * m_inv' result is\n";
my $nearly_identity = $m->multiply($m_inv);
for my $row (0..2) {
  for my $column(0..2) {
    print " (m * m_inv) [row][column] = " . 
      $nearly_identity->[$row][$column] . "\n";
  }
}

print "\nSo there.\n\n";

# === subroutines =====================================

# inverse($matrix) returns the inverse of a square $matrix.
# If $matrix isn't a Math::Matrix object, or if its determinant is 0, 
# or if it isn't square, then 0 is returned to indicate an error.
sub inverse {
  my $matrix = shift;
  return 0 unless ref($matrix) eq 'Math::Matrix';
  my ($N_rows,$N_columns) = $matrix->size;
  return 0 unless $N_rows == $N_columns;
  return 0 unless $matrix->determinant != 0;
  # Find the inverse by solving [ matrix ] [ inverse ] = [ identity ]
  # which Math::Matrix does by taking [ matrix identity ] as a single matrix.
  return $matrix->concat(identity($N_rows))->solve;
}

# Return true if $n is a positive integer;
sub is_integer {
  my $n = shift;
  if ( int($n)==$n and $n>0){
    return 1;
  }
  else {
    return 0;
  }
}

# identity($n) returns an $n x $n Math::Matrix identity matrix.
# If $n isn't a positive integer it returns 0.
sub identity {
  my $n = shift;    # The size
  return 0 unless is_integer($n);
  my @i;
  push @i, [(0)x$n] for 1..$n;
  $i[$_][$_]=1 for 0..($n-1);
  return new Math::Matrix @i;
}


# === testing ===========================
use Test::More qw(no_plan);
print "--- tests -------------\n";
ok( ( is_integer(2) 
      and not is_integer(2.5) 
      and not is_integer(-1)  
    )                               , "is_integer");

ok( ( identity(0)==0
      and identity(3)->[0][0]==1
      and identity(3)->[0][1]==0
    )                               , "identity");

my $two_by_three     = new Math::Matrix( [1,2,3], [10,20,40] );
my $zero_determinant = new Math::Matrix( [1,2], [2,4] );
my $small            = new Math::Matrix( [1,0], [2,4] );
my $small_inverse    = inverse($small);
#$small->print;
#print "   inverse is \n";
#$small_inverse->print;
ok( ( inverse($two_by_three)==0
      and inverse('a')==0
      and inverse($zero_determinant)==0
      and $small->multiply(inverse($small))->equal(identity(2))
    ),                              , "inverse");
