#
# file:    Matrix.pm
# purpose: a basic Matrix class
# created: pasha jul 13 2014
#
# matrices are in the following format:
# reference to the two-dimensional array, in the row-major format; counting is 
# started from 1 (the zero element is undefined);
# so, the matrix element m_{ij} is accessed as $matrix->[$i][$j]
#
# methods:
# $self->print()         prints the matrix
# $self->print_elems()
#

package Matrix;

use strict;
use warnings;
use pasha::dieonwarn;
use overload '==' => \&_equal;
use overload '*'  => \&_mult;
#use overload '*'  => \&_mult_scalar;
use overload '+'  => \&_sum;
use overload '-'  => \&_diff;

use constant EPSILON => 1e-12;

# possible args: 'size' => 
sub new ($@)
  {
    my $self = bless ({}, shift);
    my %args = @_;
    $self->set_size ($args{'size'}) if (defined ($args{'size'}));
    return ($self);
  }

sub set_size ($$)
  {
    my ($self, $size) = @_;
    $self->{'_size'} = $size;
  }


# set an element
# args: 1 - i
#       2 - j
#       3 - element value
sub set ($$$$)
  {
    my ($self, $i, $j, $val) = @_;
    $self->{'_matrix'}->[$i][$j] = $val;
  }


# get an element
# args: 1 - i
#       2 - j
sub get ($$$)
  {
    return ($_[0]->{'_matrix'}->[$_[1]][$_[2]]);
  }


sub print ($)
  {
    my $self = shift;
    my $j;
    for (my $i=1; $i <= $self->{'_size'}; $i++)
      {
		for ($j=1; $j <= $self->{'_size'}; $j++)
		  {
			printf ('%.11g', sprintf ("%.11f", $self->{'_matrix'}->[$i][$j]));
			print (' ');
		  }
		print ("\n");
      }

	###### tmp: print something ##########
	#print ('max non1: ' . $self->max_non1() . "\n");
	###### end of tmp ##################
  }





################### invariants of a matrix #########################
sub eigenvals ($)
  {
	my $self = shift;
	if (! defined $self->{'_eigenvals'})
	  {
		# unfold the matrix in one-dimensional array 
		# (still in the row-major format)
		my @flat = (); 
		my @row;
		for my $j (1 .. $self->{'_size'})
		  {
			@row = @{$self->{'_matrix'}->[$j]};
			shift (@row);
			push (@flat, @row);
		  }
		my $eigen = Math::MatrixDecomposition::Eigen->new;
		$eigen->decompose (\@flat);
		$self->{'_eigenvals'} = [];
		for ($eigen->values())
		  {
			push (@{$self->{'_eigenvals'}}, _round_small ($_));
		  }
	  }
	return ($self->{'_eigenvals'});
  }


sub max_eigenval ($)
  {
	my $self = shift;
	if (! defined ($self->{'_max_eigenval'}))
	  {
		for (@{$self->eigenvals()})
		  {
			if (defined ($self->{'_max_eigenval'}))
			  {
				$self->{'_max_eigenval'} = $_ if ($self->{'_max_eigenval'} < $_);
			  }
			else
			  {
				$self->{'_max_eigenval'} = $_;
			  }
		  }
	  }
	return ($self->{'_max_eigenval'});
  }


# sum of squares of eigenvalues
sub sum_sq ($)
  {
	my $self = shift;
	if (! defined ($self->{'_sum_sq'}))
	  {
		$self->{'_sum_sq'} = 0;
		for (@{$self->eigenvals()})
		  {
			$self->{'_sum_sq'} += $_ * $_;
		  }
	  }
	return ($self->{'_sum_sq'});
  }


# l1-norm (sum of absolute values of all elements)
sub l1_norm ($) 
  {
	my $self = shift;
	if (! defined ($self->{'_l1_norm'}))
	  {
		$self->{'_l1_norm'} = 0;
		my $j;
		for my $i (1 .. $self->{'_size'})
		  {
			for $j (1 .. $self->{'_size'})
			  {
				$self->{'_l1_norm'} += abs ($self->{'_matrix'}->[$i][$j]);
			  }
		  } # end of loop thru the matrix coordinates
	  }
	return ($self->{'_l1_norm'});
  }


sub l2_norm ($)
  {
	my $self = shift;
	if (! defined ($self->{'_l2_norm'}))
	  {
		$self->{'_l2_norm'} = 0;
		my $j;
		for my $i (1 .. $self->{'_size'})
		  {
			for $j (1 .. $self->{'_size'})
			  {
				$self->{'_l2_norm'} += 
				  $self->{'_matrix'}->[$i][$j] * $self->{'_matrix'}->[$i][$j];
			  }
		  } # end of loop thru the matrix coordinates
		$self->{'_l2_norm'} = sqrt ($self->{'_l2_norm'});
	  }
	return ($self->{'_l2_norm'});
  }


# entropy
sub entropy ($)
  {
	my $self = shift;
	if (! defined ($self->{'_entropy'}))
	  {
		$self->{'_entropy'} = 0; 
		for (@{$self->eigenvals()})
		  {
			$self->{'_entropy'} -= $_ * log($_) if ($_ > 0);
		  }
	  }
	return ($self->{'_entropy'});
  }
################### end of invariants of a matrix ##################



################### operations on matrices #########################
# return: 1 if equal, 0 if not, undef if cannot compare 
sub _equal ($$)
  {
    my ($left, $right) = @_;
    return (undef) if ($left->{'_size'} != $right->{'_size'});
	my $j;
    for my $i (1 .. $left->{'_size'})
	  {
		for $j (1 .. $left->{'_size'})
		  {
			return (0) if 
			  (abs ($left->{'_matrix'}->[$i][$j] - $right->{'_matrix'}->[$i][$j]) > 
			   EPSILON);
		  }
      } # end of loop thru the matrix coordinates
    return (1);
  }


# sum of two matrices
sub _sum ($$)
  {
	my ($A, $B) = @_;
	if ($A->{'_size'} != $B->{'_size'})
	  {
		return (undef);
	  }
	my $S = new Matrix ('size' => $A->{'_size'});
	my $j;
	for my $i (1 .. $A->{'_size'})
	  {
		for $j (1 .. $A->{'_size'})
		  {
			$S->{'_matrix'}->[$i][$j] = 
			  $A->{'_matrix'}->[$i][$j] + $B->{'_matrix'}->[$i][$j];
		  }
	  }
	return ($S);
  }


# multiply matrix on a scalar
sub _mult_scalar ($$)
  {
	my ($l, $A) = @_;
	my $P = new Matrix ('size' => $A->{'_size'});
	my $j;
	for my $i (1 .. $A->{'_size'})
	  {
		for $j (1 .. $A->{'_size'})
		  {
			$P->{'_matrix'}->[$i][$j] = $l * $A->{'_matrix'}->[$i][$j];
		  }
	  }
	return ($P);
  }


# difference of two matrices 
sub _diff ($$)
  {
	my ($A, $B) = @_;
	return (_sum ($A, _mult_scalar (-1, $B)));
  }


# product of two matrices
sub _mult ($$)
  {
	my ($A, $B) = @_;
	if ($A->{'_size'} != $B->{'_size'})
	  {
		print (STDERR "different sizes of matrix factors\n");
		return (undef);
	  }
	my $P = new Matrix ('size' => $A->{'_size'});
	my ($j, $k);
	for my $i (1 .. $A->{'_size'})
	  {
		for $j (1 .. $A->{'_size'})
		  {
			$P->{'_matrix'}->[$i][$j] = 0;
			for $k (1 .. $A->{'_size'})
			  {
				$P->{'_matrix'}->[$i][$j] += 
				  $A->{'_matrix'}->[$i][$k] * $B->{'_matrix'}->[$k][$j];
			  }
		  }
	  }
	return ($P);
  }
################### end of operations on matrices ##################


sub _round_small ($)
  {
    return ((($_[0] < EPSILON) && ($_[0] > - EPSILON)) ? 0 : $_[0]);
  }

1;

__END__