#!/usr/local/bin/perl -w
##---------------------------------------------------------------------------##
##  File:
##      @(#) Matrix.pm
##  Author:
##      Robert M. Hubley   rhubley@systemsbiology.org
##  Description:
##      An for holding a generic biological sequence similarity
##      search result.
##
#******************************************************************************
#* Copyright (C) Institute for Systems Biology 2003-2004 Developed by
#* Arian Smit and Robert Hubley.
#*
#* This work is licensed under the Open Source License v2.1.  To view a copy
#* of this license, visit http://www.opensource.org/licenses/osl-2.1.php or
#* see the license.txt file contained in this distribution.
#*
#******************************************************************************
# Implementation Details:
#
# bless( {
#        }
#     }, 'Matrix' );
#
###############################################################################
# ChangeLog
#
#     $Log$
#
###############################################################################
# To Do:
#

=head1 NAME

Matrix

=head1 SYNOPSIS

use Matrix

Usage: 

    my $matrix = Matrix->new();

  or 

    my $matrix = Matrix->new( fileName => foo.matrix );

=head1 DESCRIPTION

A class for storing a RepeatMasker matrix

=head1 SEE ALSO

=over 4

=back

=head1 COPYRIGHT

Copyright 2011 Institute for Systems Biology

=head1 AUTHOR

Robert Hubley <rhubley@systemsbiology.org>

=head1 INSTANCE METHODS

=cut 

package Matrix;
use strict;
use Data::Dumper;
use Carp;
use vars qw(@ISA @EXPORT @EXPORT_OK %EXPORT_TAGS $VERSION);

require Exporter;

# Matrix Type
use constant crossmatchMatrix => 1;
use constant blastMatrix      => 2;

@ISA = qw(Exporter);

@EXPORT = qw();

@EXPORT_OK = qw();

%EXPORT_TAGS = ( all => [ @EXPORT_OK ] );

my $CLASS = "Matrix";

##-------------------------------------------------------------------------##
## Constructor:
##-------------------------------------------------------------------------##
sub new {
  my $class          = shift;
  my %nameValuePairs = @_;

  # Create ourself as a hash
  my $this = {};

  # Bless this hash in the name of the father, the son...
  bless $this, $class;

  # Allow import of values
  if ( %nameValuePairs ) {

    # Special parameters first
    if ( defined $nameValuePairs{'fileName'} ) {
      $this->parseFromFile( $nameValuePairs{'fileName'},
                            $nameValuePairs{'matrixType'} );
    }
    else {

      # Default calls to set routines
      while ( my ( $name, $value ) = each( %nameValuePairs ) ) {
        my $method = "set" . _ucFirst( $name );
        unless ( $this->can( $method ) ) {
          croak( "Matrix::add: Instance variable $name doesn't exist." . "" );
        }
        $this->$method( $value );
      }
    }
  }

  return $this;
}

##-------------------------------------------------------------------------##

=head2 clone()

  Use: my $newObj = $obj->clone();

  Clone a Matrix *duplicating* all the values of the old
  object in the new one.

=cut

##-------------------------------------------------------------------------##
sub clone {
  my $this = shift;

  my %newHash = %{$this};
  my $newObj  = \%newHash;

  bless $newObj, ref( $this );

  return $newObj;
}

##-------------------------------------------------------------------------##
## Get and Set Methods
##-------------------------------------------------------------------------##

##-------------------------------------------------------------------------##

=head2 get_setMatrixName()

  Use: my $value    = getMatrixName( );
  Use: my $oldValue = setMatrixName( $value );

  Get/Set the name of the matrix.

=cut

##-------------------------------------------------------------------------##
sub getMatrixName {
  my $obj = shift;

  my $value = $obj->{'matrixName'};

  return $value;
}

sub setMatrixName {
  my $obj      = shift;
  my $value    = shift;
  my $oldValue = undef;

  $oldValue = $obj->{'matrixName'};
  $obj->{'matrixName'} = $value;

  return $oldValue;
}

##-------------------------------------------------------------------------##

=head2 get_setAlphabetRef()

  Use: my $alphabetArrayRef = getAlphabetRef();
  Use: my $oldAlphabetArrayRef = setAlphabetRef( $alphabetArrayRef );

  Get/Set the alphabet of the matrix.

=cut

##-------------------------------------------------------------------------##
sub getAlphabetRef {
  my $obj = shift;

  my $value = $obj->{'alphabetArrayRef'};

  return $value;
}

sub setAlphabetRef {
  my $obj      = shift;
  my $value    = shift;
  my $oldValue = undef;

  $oldValue = $obj->{'alphabetArrayRef'};
  $obj->{'alphabetArrayRef'} = $value;

  return $oldValue;
}

##-------------------------------------------------------------------------##

=head2 getMatrixValuesRef()

  Use: my $alphabetArrayRef = getMatrixValuesRef();

  Get the matrix.

=cut

##-------------------------------------------------------------------------##
sub getMatrixValuesRef {
  my $obj = shift;

  my $value = $obj->{'matrixValues'};

  return $value;
}

##-------------------------------------------------------------------------##

=head2 getMatrixFreqsRef()

  Use: my $alphabetArrayRef = getMatrixFreqsRef();

  Get the matrix.

=cut

##-------------------------------------------------------------------------##
sub getMatrixFreqsRef {
  my $obj = shift;

  my $value = $obj->{'alphabetFreqArrayRef'};

  return $value;
}

##-------------------------------------------------------------------------##

=head2 getLambda()

  Use: my $lambda = getLambda();

  Get the 

=cut

##-------------------------------------------------------------------------##
sub getLambda {
  my $obj = shift;

  my $value = $obj->{'lambda'};

  return $value;
}

##-------------------------------------------------------------------------##
## General Object Methods
##-------------------------------------------------------------------------##

##-------------------------------------------------------------------------##

=head2

  Use: toString();

  Create an string representation of the matrix.

=cut

##-------------------------------------------------------------------------##
sub toString {
  my $obj = shift;
}

##-------------------------------------------------------------------------##

=head2

  Use: parseFromFile( $csvString );

  Populate object with values stored in CSV format.

=cut

##-------------------------------------------------------------------------##
sub parseFromFile {
  my $this       = shift;
  my $fileName   = shift;
  my $matrixType = shift;

  $matrixType = Matrix::blastMatrix if ( !defined $matrixType );

  my $subroutine = ( caller( 0 ) )[ 0 ] . "::" . ( caller( 0 ) )[ 3 ];

  open MATRIX, "<$fileName"
      or croak( "$subroutine: Could not open $fileName for reading!" );

  my $row = 0;
  while ( <MATRIX> ) {
    $_ = uc( $_ );
    chomp;

    # FREQS A 0.325 C 0.175 G 0.175 T 0.325
    if ( /FREQS\s+(.*)/ ) {
      $this->{'freqHashRef'} = { split( " ", $1 ) };
    }

    #
    #   A  C  G  T
    # A 1 -1 -1 -1
    if ( /^\s*[A-Z]\s+[A-Z]\s+[A-Z]\s+[A-Z]\s+/ ) {
      s/ //g;
      $this->{'alphabetArrayRef'} = [ split( //, $_ ) ];
    }
    elsif ( $this->{'alphabetArrayRef'}
            && /^\s*([^\d]\s+)?[\d-]+\s+[\d-]+\s+[\d-]+\s+[\d-]+/ )
    {
      my @rowValues = split;
      if ( $1 ) {
        shift @rowValues;
      }
      $this->{'matrixValues'}->[ $row++ ] = [ @rowValues ];
    }

  }
  close MATRIX;

  if ( !defined $this->{'freqHashRef'} ) {
    $this->{'freqHashRef'}->{"A"} = 0.25;
    $this->{'freqHashRef'}->{"C"} = 0.25;
    $this->{'freqHashRef'}->{"G"} = 0.25;
    $this->{'freqHashRef'}->{"T"} = 0.25;
  }

  for ( my $i = 0 ; $i <= $#{ $this->{'alphabetArrayRef'} } ; $i++ ) {
    $this->{'alphabetFreqArrayRef'}->[ $i ] =
        $this->{'freqHashRef'}->{ $this->{'alphabetArrayRef'}->[ $i ] } || 0;
  }

  $this->_calculateLambda();

}

##-------------------------------------------------------------------------##
## Use: my transposeMatrix();
##
## Returns
##      Transposes the matrix values ( ie. rows -> cols ) inside
##      the object and recalculates the lambda value.
##-------------------------------------------------------------------------##
sub transposeMatrix {
  my $this    = shift;
  my @tMatrix = ();
  for ( my $i = 0 ; $i <= $#{ $this->{'alphabetArrayRef'} } ; $i++ ) {
    for ( my $j = 0 ; $j <= $#{ $this->{'alphabetArrayRef'} } ; $j++ ) {
      push @{ $tMatrix[ $i ] }, $this->{'matrixValues'}->[ $j ]->[ $i ];
    }
  }
  $this->{'matrixValues'} = \@tMatrix;
  $this->_calculateLambda();
}

##-------------------------------------------------------------------------##
## Use: my _calculateLambda( $matScoresRef, $matFreqsRef );
##
##      $matScoresRef: Score Matrix
##      $matFreqsRef : Matrix alphabet freq vector
##
## Returns
##      $lambda : The lambda parameter derived from the matrix
##                and the matrix alphabet frequencies.  This
##                is derived from Phil Green's swat/cross_match
##                programs.
##-------------------------------------------------------------------------##
sub _calculateLambda {
  my $this = shift;

  my $subroutine = ( caller( 0 ) )[ 0 ] . "::" . ( caller( 0 ) )[ 3 ];

  my $lambda_upper = 0;
  my $lambda_lower = 0;
  my $lambda       = 0.5;

  my $matFreqsRef  = $this->{'alphabetFreqArrayRef'};
  my $matScoresRef = $this->{'matrixValues'};
  my $sum          = 0;

  do {
    $sum = 0; # 9/8/2021 Fixed a major bug in lambda calculation
    # NOTE: This fix only had an impact on a single RM matrix
    #       ( identity.matrix ) as this is the only RM matrix
    #       that exceeded a lambda of 1.0.
    my $check = 0;
    for ( my $i = 0 ; $i <= $#$matFreqsRef ; $i++ ) {
      for ( my $j = 0 ; $j <= $#$matFreqsRef ; $j++ ) {
        if ( $$matFreqsRef[ $i ] && $$matFreqsRef[ $j ] ) {
          $sum += $$matFreqsRef[ $i ] * $$matFreqsRef[ $j ] *
              exp( $lambda * $$matScoresRef[ $i ][ $j ] );
          $check += $$matFreqsRef[ $i ] * $$matFreqsRef[ $j ];
        }
      }
    }
    croak "$subroutine: Sanity check failed lowerLambda freqSum = $check"
        . " but should have been between close to 1 ( .999 - 1.001 )"
        if (    $check > 1.001
             || $check < .999 );

    if ( $sum < 1.0 ) {
      $lambda_lower = $lambda;
      $lambda *= 2.0;
    }
  } while ( $sum < 1.0 );
  $lambda_upper = $lambda;

  while ( $lambda_upper - $lambda_lower > .00001 ) {
    $lambda = ( $lambda_lower + $lambda_upper ) / 2.0;
    my $sum   = 0;
    my $check = 0;
    for ( my $i = 0 ; $i <= $#$matFreqsRef ; $i++ ) {
      for ( my $j = 0 ; $j <= $#$matFreqsRef ; $j++ ) {
        if ( $$matFreqsRef[ $i ] && $$matFreqsRef[ $j ] ) {
          $sum += $$matFreqsRef[ $i ] * $$matFreqsRef[ $j ] *
              exp( $lambda * $$matScoresRef[ $i ][ $j ] );
          $check += $$matFreqsRef[ $i ] * $$matFreqsRef[ $j ];
        }
      }
    }
    croak "$subroutine: Sanity check failed upperLambda freqSum = $check"
        . " but should have been between close to 1 ( .999 - 1.001 )"
        if (    $check > 1.001
             || $check < .999 );

    if ( $sum >= 1.0 ) {
      $lambda_upper = $lambda;
    }
    else {
      $lambda_lower = $lambda;
    }
  }

  $this->{'lambda'} = $lambda;
}

1;