/***************************************************************************
 *   Copyright (C) 2005 by Andreas Pokorny                                 *
 *   andreas.pokorny@biozentrum.uni-wuerzburg.de                           *
 *                                                                         *
 *   This file is part of profdist and cbcanalyzer                         *
 *                                                                         *
 *   Both profdist and cbcanalyzer are free software; you can redistribute *
 *   it and/or modify it under the terms of the GNU General Public License *
 *   as published by the Free Software Foundation; either version 2 of the *
 *   License, or (at your option) any later version.                       *
 *                                                                         *
 *   Profdist and cbcanalyzer are distributed in the hope that it will be  *
 *   useful, but WITHOUT ANY WARRANTY; without even the implied warranty   *
 *   of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the      *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU General Public License     *
 *   along with this program; if not, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/

#include "correction.h"
#include "distance.h"
#include "fixed_matrix.h"
namespace profdist {
//--------------------- Distance computation -----------------------------------------------
void compute_distance( AlignCode const& source, profdist::distance_matrix & matrix, profdist::rate_matrix const& Q, profdist::CorrectionModel model, profdist::ExpMType expmType, profdist::FZero fMax, ProgressSink & prg )
{
  std::size_t num_seq = source.get_num_sequences()
    , max_steps = ( num_seq - 1 ) * ( num_seq - 2 ) / 2
    , num_steps = 0;
  double value = 0.0;
  bool cont = true;
  matrix.resize( num_seq, num_seq, 0.0 );

  for ( std::size_t i = 0; i < num_seq - 1; ++i )
  {
    profdist::count_matrix const& A = source.get_matrix(i);
    // compute the distance between Sequence 1 and i + 2 with SubstMatrix N_1_i+2
    value = correction( A , Q, model, expmType, fMax );

    // set the values in the Distancematrix
    matrix( 0, i + 1 ) = value;
    matrix( i + 1, 0 ) = value;

    for( std::size_t j = i + 1; j < num_seq - 1; ++j )
    {
      // Initialialize SubstMatrix N_i+2_j+2 with the diagonale from N_1_i+2
      profdist::count_matrix B( 0U );
      for( int k = 0; k < 4; ++k )
        B[k][k] = A[k][k];

      // i+1 and j+1 since we start at the second  sequence
      AlignCode::const_diff_iterator it_b = source.begin_difference( j )
        , it_a = source.begin_difference( i )
        , a_end = source.end_difference( i )
        , b_end = source.end_difference( j )
        ;

      while ( it_a != a_end && it_b != b_end  )
      {
        if ( it_a->first == it_b->first ) // both are different to first sequence
        {
          ++B[ it_a->second ][ it_b->second ];
          ++it_a;
          ++it_b;
        }
        else if( it_a->first < it_b->first ) // Unterschied festgestellt und wird in Unterschiedsliste aufgenommen
        {
          ++B[ it_a->second][ source.get_reference_element( it_a->first ) ];
          ++it_a;
        }
        else
        {
          ++B[ it_b->second][ source.get_reference_element( it_b->first ) ];
          --B[ source.get_reference_element( it_b->first )][ source.get_reference_element( it_b->first ) ];
          ++it_b;
        }
      } //end of while

      while ( it_a != a_end  )
      {
        ++B[ it_a->second][ source.get_reference_element( it_a->first ) ];
        ++it_a;
      }

      while ( it_b != b_end )
      {
        ++B[ it_b->second][ source.get_reference_element( it_b->first ) ];
        --B[ source.get_reference_element( it_b->first )][ source.get_reference_element( it_b->first ) ];
        ++it_b;
      }

      value = correction( B, Q, model, expmType, fMax );

      matrix( i + 1, j + 1 ) = value;
      matrix( j + 1, i + 1 ) = value;

      if( ++num_steps == max_steps )
        cont = prg.update( num_steps, "Distance matrix calculated!" );
      else
        cont = prg.update( num_steps );

      if ( !cont )
        return;
    } // end for j
  } //end for  i

  return;
}

void compute_distance( AlignCode const& source, profdist::distance_matrix & matrix, profdist::rate_matrix const& Q, profdist::CorrectionModel model, profdist::ExpMType expmType, profdist::FZero fMax )
{
  std::size_t num_seq = source.get_num_sequences()
    , max_steps = ( num_seq - 1 ) * ( num_seq - 2 ) / 2
    , num_steps = 0;
  double value = 0.0;
  bool cont = true;
  matrix.resize( num_seq, num_seq, 0.0 );

  for ( std::size_t i = 0; i < num_seq - 1; ++i )
  {
    profdist::count_matrix const& A = source.get_matrix(i);
    // compute the distance between Sequence 1 and i + 2 with SubstMatrix N_1_i+2
    value = correction( A , Q, model, expmType, fMax );

    // set the values in the Distancematrix
    matrix( 0, i + 1 ) = matrix( i + 1, 0 ) = value;

    for( std::size_t j = i + 1; j < num_seq - 1; ++j )
    {
      // Initialialize SubstMatrix N_i+2_j+2 with the diagonale from N_1_i+2
      profdist::count_matrix B( 0U );
      for( int k = 0; k < 4; ++k )
        B[k][k]= A[k][k];

      // i+1 and j+1 since we start at the second  sequence
      AlignCode::const_diff_iterator it_b = source.begin_difference( j )
        , it_a = source.begin_difference( i )
        , a_end = source.end_difference( i )
        , b_end = source.end_difference( j )
        ;

      while ( it_a != a_end && it_b != b_end  )
      {
        if ( it_a->first == it_b->first ) // both are different to first sequence
        {
          ++B[ it_a->second][ it_b->second ];
          ++it_a;
          ++it_b;
        }
        else if( it_a->first < it_b->first ) // Unterschied festgestellt und wird in Unterschiedsliste aufgenommen
        {
          ++B[ it_a->second][ source.get_reference_element( it_a->first ) ];
          ++it_a;
        }
        else
        {
          ++B[ it_b->second][ source.get_reference_element( it_b->first ) ];
          --B[ source.get_reference_element( it_b->first )][ source.get_reference_element( it_b->first ) ];
          ++it_b;
        }
      } //end of while

      while ( it_a != a_end  )
      {
        ++B[ it_a->second][ source.get_reference_element( it_a->first ) ];
        ++it_a;
      }

      while ( it_b != b_end )
      {
        ++B[ it_b->second][ source.get_reference_element( it_b->first ) ];
        --B[ source.get_reference_element( it_b->first )][ source.get_reference_element( it_b->first ) ];
        ++it_b;
      }

      value = correction( B, Q, model, expmType, fMax );

      matrix( i + 1,  j + 1 ) = value;
      matrix( j + 1,  i + 1 ) = value;

    } // end for j
  } //end for  i

  return;
}


void compute_distance( Profile const& source, profdist::distance_matrix & matrix, profdist::rate_matrix const& Q, profdist::CorrectionModel model, profdist::ExpMType expmType, profdist::FZero fMax )
{
  std::size_t num_profiles = source.get_num_profiles();
  matrix.resize( num_profiles, num_profiles, 0.0 );
  std::size_t num_sites = source.get_num_sites(); 
  
  for( std::size_t i = 0; i != num_profiles - 1; ++i )
    for( std::size_t j = i+1; j != num_profiles; ++j ) {
      
      Profile::ConstSingleProfile i_prof( source[i] )
        , j_prof( source[j] );
      
      
      profdist::small_float_count_matrix cm( 0.0 );
      for( std::size_t site_index = 0; site_index != num_sites; ++site_index )
        for( std::size_t A = 0 ; A < 4; ++ A) {
          double val = i_prof.get_value( site_index, A );
          for( std::size_t B = 0 ; B < 4; ++ B) 
            cm[A][B] = cm[A][B] + val* j_prof.get_value( site_index, B );
        }

      matrix(i,j) = matrix(j,i) = correction( cm , Q, model, expmType, fMax );
    }
}

std::istream& read_rate_matrix( std::istream& in, profdist::rate_matrix & matrix )
{
  in >> std::ws;
  std::size_t num;
  if(in >> num && num == 4 ) 
  {
    std::size_t i = 0; 
    in >>  matrix[0][0] >> matrix[0][1] >> matrix[0][2] >> matrix[0][3]
      >> matrix[1][0] >> matrix[1][1] >> matrix[1][2] >> matrix[1][3]
      >> matrix[2][0] >> matrix[2][1] >> matrix[2][2] >> matrix[2][3]
      >> matrix[3][0] >> matrix[3][1] >> matrix[3][2] >> matrix[3][3];
  }
  else 
    in.fail();
  return in;
}

}


