#
# file:         CorrMatrix.pm
# purpose:      object representing a correlation matrix
# created:      pasha aug 15 2008
# modified:     pasha jul 8,10-11,13 2014
# modification: inherit from Matrix
#
# methods:
# $self->corr()          computes correlation matrix (of the given matrix - i.e correlation 
#                        between rows)
# $self->iter()          computes iterative correlation
# $self->print_elems()
# $self->random()        generate a random matrix
#

package CorrMatrix;

use strict;
use warnings;
use Math::MatrixDecomposition::Eigen;
use pasha::dieonwarn;

#require Exporter;
#our @ISA = qw(Exporter);
use Matrix;
our @ISA = qw(Matrix);
  
  

# generate a random matrix
# args: 1 - matrix size
sub random ($$)
  {
    my ($self, $size) = @_;
    $self->set_size ($size);

    my $i = 1;    # matrix coordinates
    my $j = 2;

	# loop thru the all coordinates in the upper right triangle of the square 
	# matrix: (1,2) (1,3) ... (1,n) (2,3) ... (2,n) ... (n-1,n)
    while (1)  
      {
		# scale a random number from [0,1] to [-1,1]
		$self->{'_matrix'}->[$j][$i] = $self->{'_matrix'}->[$i][$j] 
		  = 2*rand(1) - 1;
	
		# go to the next pair
		$j++;
		if ($j > $size)
		  {
			$i++;
			last if ($i >= $size);
			$j = $i+1;
		  }
      } # end of loop thru the matrix coordinates
  } # end of random()
  

sub set_size ($$)
  {
    my ($self, $size) = @_;
    $self->SUPER::set_size ($size);
    # set diagonal elements
    for (my $i=1; $i <= $size; $i++)
      {
		$self->{'_matrix'}->[$i][$i] = 1.0;
      }
  }

  
# set an element (as in a symmetrtic matrix, i.e. both (i,j) and (j,i))
# args: 1 - i
#       2 - j
#       3 - element value
sub set ($$$$)
  {
    my ($self, $i, $j, $val) = @_;
    $self->{'_matrix'}->[$j][$i] = $self->{'_matrix'}->[$i][$j] = $val;
  }

  
## set diagonal element
#sub set_diag ($$$)
#  {
#    $_[0]->{'_matrix'}->[$_[1]][$_[1]] = $_[2];
#  }



# print just elements in the row, unfolding upper right triangle row-by-row
sub print_elems ($)
  {
    my $self = shift;
    my $i = 1;  # matrix coordinates
    my $j = 2;
    while (1)
      {
		print (STDOUT _round_small ($self->{'_matrix'}->[$i][$j]) . ' ');
		$j++;
		if ($j > $self->{'_size'})
		  {
			$i++;
			last if ($i >= $self->{'_size'});
			$j = $i+1;
		  }
      } # end of loop thru the matrix coordinates
  }



# compute correlation matrix of the given matrix (i.e. correlations between columns)
# return: pair (return code, ref to the correlation matrix)
# return code = 1 if everything is ok, 0 if correlation matrix cannot be computed
# (i.e. when at least one of the rows is constant (1, ..., 1))
sub corr ($)
  {
    my $self = shift;
    my $corr_matrix = new CorrMatrix ('size' => $self->{'_size'});
    my ($j, $mean_i, $mean_j, $ci, $cj);
    for (my $i=1; $i <= $self->{'_size'}; $i++)
      {
	$mean_i = _mean ($self->{'_matrix'}->[$i]);
	$ci = _covar_self ($self->{'_matrix'}->[$i], $mean_i, $self->{'_size'});
	return (0, undef) if ($ci == 0);
	for ($j=$i+1; $j <= $self->{'_size'}; $j++)
	  {
	    $mean_j = _mean ($self->{'_matrix'}->[$j]);
	    $cj = _covar_self ($self->{'_matrix'}->[$j], $mean_j, $self->{'_size'});
	    return (0, undef) if ($cj == 0);
	    $corr_matrix->set 
	      ($i, $j, 
	       _covar ($self->{'_matrix'}->[$i], $self->{'_matrix'}->[$j], $mean_i, $mean_j, 
		       $self->{'_size'}) /  
	       (sqrt ($ci * $cj)));
	  }
      }
    return (1, $corr_matrix);
  } # end of corr()



# iterate correlations the given number of steps
# args: 1 - number of iterations
#       2 - frequency with which to print intermediate matrices; 
#           0 if not print at all
# return: 1 if converges or matrix is degenerate, 0 if not
sub iter ($$$)
  {
    my ($self, $num_iter, $print) = @_;
    my $matrix = $self;
    $self->print() if ($print != 0);
    my ($rc, $matrix_next);
    for (my $i=1; $i <= $num_iter; $i++)
      { 
		($rc, $matrix_next) = $matrix->corr();
		if (($rc != 1) || ($matrix_next == $matrix))
		  {
			if ($print != 0)
			  {
				print ("\nconverges in " . ($i-1) . " steps to:\n");
				$matrix_next->print();
			  }
			return (1);
		  }
		
		###### tmp: compare something #################
		#if (($i > 1) && ($matrix_next->max_eigenval() < $matrix->max_eigenval()))
		#  {
		#	die ("max eigenval decreases!\n");
		#  }

		#print ("commutator:\n");
		#my $comm = $matrix * $matrix_next - $matrix_next * $matrix;
		#$comm->print();
		###### end of tmp #############################

		$matrix = $matrix_next;
		if ($rc == 0)
		  {
			print ("degenerate matrix\n");
			return (1);
		  }
		if (($print != 0) && ($i % $print == 0))
		  {
			print ("\n${i}:\n");
			$matrix->print();
		  }
      }
    if ($print != 0)
      {
		print ("\ndoes not converge after $num_iter steps:\n");
		$matrix->print();
      }
    return (0); 
  } # end of iter()




# mean of vector (starting from 1)
sub _mean ($)
  {
    my $vector = $_[0];
    my $mean = 0;
    my $n = scalar (@{$vector}) - 1;
    for (my $i=1; $i<=$n; $i++)
      {
		$mean += $vector->[$i];
      }
    return ($mean / $n);
  }


# covariance of two vectors (starting from 1)
# args: 0 - ref to the first vector
#       1 - ref to the second vector
#       2 - mean of the first vector
#       3 - mean of the second vector
#       4 - length of vectors
sub _covar ($$$$$)
  {
    my ($vector_a, $vector_b, $mean_a, $mean_b, $n) = @_;
    my $covar = 0;
    for (my $i=1; $i<=$n; $i++)
      {
		$covar += ($vector_a->[$i] - $mean_a) * ($vector_b->[$i] - $mean_b);
      }
    return ($covar);
  }
  

# covariance of vector with itself
# args: 0 - ref to the vector
#       1 - mean of the vector
#       2 - length of the vector
sub _covar_self ($$$)
  {
    return (_covar ($_[0], $_[0], $_[1], $_[1], $_[2]));
  }


1;

__END__