/***************************************************************************
 *   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 <algorithm>
#include <iterator>
#include "profile.h"

Profile::Profile( std::size_t sites, std::size_t prof, std::size_t rows) 
  : num_sites( sites ), num_profiles( prof ), num_rows( rows ), profile_array(num_profiles * num_sites* num_rows, 0.0 )
{
}

Profile::Profile( AlignCode const& source, tree_types::profile_map const& tree, std::size_t rows )
  : num_sites(source.get_num_sites()), num_profiles(tree.size()), num_rows( rows ), profile_array( num_profiles*num_sites*num_rows, 0.0 )
{
  inner_create_profiles( source, tree );
}

Profile::Profile( AlignCode const& source, std::list<std::list<std::size_t> > const& profiles, std::size_t rows )
  : num_sites(source.get_num_sites()), num_profiles(profiles.size()), num_rows( rows ), profile_array( num_profiles*num_sites*num_rows, 0.0 )
{
  inner_create_profiles( source, profiles );
}

void Profile::init( std::size_t sites, std::size_t prof, std::size_t rows )
{
  num_sites = sites; 
  num_profiles = prof;
  num_rows = rows;
  profile_array.resize( num_profiles*num_sites*num_rows);
  std::memset( &(profile_array[0]), 0, sizeof(double)*profile_array.size() );
}

void Profile::init( AlignCode const& source, tree_types::profile_map const& tree, std::size_t rows )
{
  init( source.get_num_sites(), tree.size(), rows );
  inner_create_profiles( source, tree );
}
void Profile::init( AlignCode const& source, std::list< std::list<std::size_t> >  const& profiles, std::size_t rows )
{
  init( source.get_num_sites(), profiles.size(), rows );
  inner_create_profiles( source, profiles );
}

void Profile::inner_create_profiles( AlignCode const& source, std::list< std::list<std::size_t> >  const& profiles )
{
  using namespace std;
  size_t prof_index = 0;
  for( list<list<size_t> >::const_iterator it = profiles.begin(), e = profiles.end();  it != e; ++it, ++prof_index )
  {
    SingleProfile profile = get_profile( prof_index );

    size_t profile_size = it->size();
    // first gather information from the reference sequence 
    for( size_t i = 0; i != num_sites; ++i ) 
      profile.get_value(i, source.get_reference_element( i ) ) = profile_size;

    handle_differences( profile, source, it->begin(), it->end() );
    profile.normalize( profile_size );
  }
}
void Profile::process_difference( SingleProfile & profile, AlignCode const& code, std::size_t index )
{
  for( AlignCode::const_diff_iterator diff_it = code.begin_difference( index ), 
      diff_end = code.end_difference( index );
      diff_it != diff_end; ++diff_it ) {
    // .. add the difference to the reference sequence
    ++(profile.get_value( diff_it->first, diff_it->second) );
    // .. and remove the reference entries at that position:
    --(profile.get_value( diff_it->first, code.get_reference_element( diff_it->first ))); 
  }
}

void Profile::inner_create_profiles( AlignCode const& source, tree_types::profile_map const& tree )
{
  for( tree_types::profile_map::const_iterator it = tree.begin(), e = tree.end();  it != e; ++it )
  {
    SingleProfile profile = get_profile( it->first );
    
    if( it->second->is_leaf() ) {  // just a single sequence to process here 

      // first set as if we would process the reference sequence
      for( std::size_t i = 0; i != num_sites; ++i ) 
        profile.get_value(i, source.get_reference_element( i )) = 1;

      // if we in fact process a different then ..
      if( it->second->get_reference_position() != 0 )
      {
        // - 1 because difference vectors starts at zero but stores elements 1 to N-1 
        process_difference( profile, source, it->second->get_reference_position() - 1 );
      }
    }
    else {
      std::size_t set_size = it->second->get_split_set().size();

      // first gather information from the reference sequence 
      for( std::size_t i = 0; i != num_sites; ++i ) 
        profile.get_value(i, source.get_reference_element( i ) ) = set_size;


      handle_differences( profile, source, it->second->get_split_set().begin(), it->second->get_split_set().end() );

      profile.normalize( set_size );
    }
  }
}

void Profile::refine( tree_types::profile_map const& list )
{
  Profile empty( get_num_sites(), list.size(), num_rows ); // creating empty profile 
  empty.swap(*this); // swap to keep current and make "this" empty

  Profile const& backup( empty );
  
  for( tree_types::profile_map::const_iterator it = list.begin(), e = list.end();  it != e; ++it  )
  {
    SingleProfile dest(get_profile( it->first ));
    CountedSets::set_type const& splits = it->second->get_split_set();
    
    if( splits.size() )
    {
      for( Node::set_type::const_iterator set_it = splits.begin(), set_e = splits.end(); set_it != set_e;  ++set_it ) 
        dest += backup.get_profile( *set_it  );
      dest.normalize( splits.size() );
    }
    else // single node => there was no change in that step => just copy the data
    {
      ConstSingleProfile source( backup.get_profile( it->second->get_reference_position() ) );
      dest = source;
    }
  }
}

void Profile::SingleProfile::normalize( std::size_t num_sequences_in_profile )
{
  for( std::size_t i = 0; i != num_sites*num_rows; ++i ) 
    data[i] /= num_sequences_in_profile;
}

Profile::SingleProfile::SingleProfile( value_type* d, std::size_t rows, std::size_t sites )
  : data( d ), num_rows( rows ), num_sites( sites ) { }


Profile::SingleProfile Profile::get_profile(std::size_t profile_index ) 
{
  return SingleProfile( &*(profile_array.begin() + profile_index * num_sites * num_rows), num_rows, num_sites);
}

Profile::SingleProfile Profile::operator[]( std::size_t profile_index ) 
{
  return SingleProfile( &*(profile_array.begin() + profile_index * num_sites * num_rows), num_rows, num_sites );
}

Profile::ConstSingleProfile::ConstSingleProfile( Profile::SingleProfile const& ref )
  : data( ref.data ), num_rows( ref.num_rows ), num_sites( ref.num_sites ) { }
Profile::ConstSingleProfile::ConstSingleProfile( value_type const* d, std::size_t rows, std::size_t sites )
  : data( d ), num_rows( rows ), num_sites( sites ) { }

Profile::ConstSingleProfile Profile::get_profile(std::size_t profile_index )  const
{
  return ConstSingleProfile( &*(profile_array.begin() + profile_index * num_sites * num_rows), num_rows, num_sites);
}

Profile::ConstSingleProfile Profile::operator[]( std::size_t profile_index ) const
{
  return ConstSingleProfile( &*(profile_array.begin() + profile_index * num_sites * num_rows), num_rows, num_sites );
}

Profile::SingleProfile& Profile::SingleProfile::operator=( ConstSingleProfile const& right )
{
  std::memcpy( data, right.data, num_sites*sizeof(Profile::value_type)*num_rows );
  return *this;
}
Profile::SingleProfile& Profile::SingleProfile::operator+=( ConstSingleProfile const& right )
{
  std::size_t i = num_rows*num_sites;
  do {
    --i;
    data[i] += right.data[i];
  }
  while( i );
  return *this;
}

const std::size_t Profile::get_num_sites() const
{
  return num_sites;
}
const std::size_t Profile::get_num_profiles() const
{
  return num_profiles;
}

const std::size_t Profile::get_num_rows() const
{
  return num_rows;
}

/*std::ostream& operator << ( std::ostream & out, Profile const& obj )
{
  // WARNING: This code depends on the Profile memory layout
  out << obj.get_num_sequences() << ' ' << obj.get_num_sites() << '\n';
  
  Profile::const_iterator b = obj.begin();
  for( Profile::const_iterator it = obj.begin(), e = obj.end(); 
      it != e; ++it )
  {
    out << setw(5) << *it;
    if( ( it - b  ) % obj.get_num_rows() == 0 ) 
      out << '\n';
  }
  out << std::endl;
  return out;
}
*/

void Profile::swap( Profile& rhs )
{
  std::swap(num_sites,rhs.num_sites);
  std::swap(num_profiles,rhs.num_profiles);
  std::swap(num_rows,rhs.num_rows );
  std::swap(profile_array,rhs.profile_array);
}

void swap( Profile &lhs, Profile &rhs )
{
  lhs.swap(rhs);
}

std::ostream & operator<<( std::ostream & o, Profile const& obj )
{
  for( std::size_t i = 0; i < obj.get_num_profiles(); ++i ) {
    o << i << ".\n";
    Profile::ConstSingleProfile prof = obj.get_profile( i );
    for( std::size_t j = 0; j < obj.get_num_rows(); ++j ) {
      for( std::size_t k = 0; k < obj.get_num_sites(); ++k ) 
        o << prof.get_value(k,j)<< '|';
      o << std::endl;
    }
  }
  return o;

}


