#include "variance/SingletonEvent.hh"
#include "common/SequenceUtils.hh"
#include "variance/SampleStats.hh"

/**
 * Project     : CASAVA
 * Module      : $RCSfile: SingletonEvent.cpp,v $
 * @author     : Tony Cox
 * Copyright   : Copyright (c) Illumina 2008, 2009. All rights reserved.
 *
 ** This software is covered by the "Illumina Genome Analyzer Software
 ** License Agreement" and the "Illumina Source Code License Agreement",
 ** and certain third party copyright/licenses, and any user of this
 ** source file is bound by the terms therein (see accompanying files
 ** Illumina_Genome_Analyzer_Software_License_Agreement.pdf and
 ** Illumina_Source_Code_License_Agreement.pdf and third party
 ** copyright/license notices).
 *
 */

//using namespace ga::common;
using namespace casava::common;

namespace ca {
namespace variance_detection {


  SingletonNotMatched::SingletonNotMatched
    ( const Alignment& alignment,
      const Alignment& alignmentPrevious,
      map <string, SampleStats>& sampleStats ) :
      Singleton( alignment.getData().c_str(),
		 alignment.getQuality().c_str(), alignment ),
      //      alignment_(alignment),
      alignmentPrevious_(alignmentPrevious)
  {
    Match::Strand strand = Match::None;
    std::string sampleName;
    //    Match thisMatch;
    sampleName=SequenceUtils::getSampleName(alignment_);

    standardDeviation_=sampleStats[sampleName].offsetStats_.getStandardDeviation();
    // RJS : removed this assertion in favour of handling standardDeviation_
    //       of 0 where it is used (populateMetric).
    // assert(standardDeviation_!=0);

    strand=alignmentPrevious_.getMatch().getStrand();

    Match thisMatch(alignment_.getMatch());

    if (strand == Match::Forward)
    {
      minPos_=alignmentPrevious_.getMatch().getPosition()
	+sampleStats[sampleName].minOffset_;
      maxPos_=alignmentPrevious_.getMatch().getPosition()
	+sampleStats[sampleName].maxOffset_;
      thisMatch.setStrand(Match::Reverse);
    } // if
    else
    {
      minPos_=alignmentPrevious_.getMatch().getPosition()
	-sampleStats[sampleName].maxOffset_;
      maxPos_=alignmentPrevious_.getMatch().getPosition()
	-sampleStats[sampleName].minOffset_;
      thisMatch.setStrand(Match::Forward);
    } // ~else

    alignment_.setMatch(thisMatch);

    // make slighltly more conservative assumption than read lies in
    // [expected min start of read 1, expected end of read 2]
    // instead of [expected min midpt read 1, expected max endpt read 2]
    maxPos_+=strlen(alignment_.getData().c_str())-1;

    // in absence of other information, guess the position of variant causing
    // the non-match is being in the middle of the reas
    //    minPos_+=strlen(alignment_.getData().c_str())/2;
    //  maxPos_+=strlen(alignment_.getData().c_str())/2;

    // put everything in the forward strand
    if (strand == Match::Forward)
    {
       SequenceUtils::reverseComp(read_);
    } // ~if
  } // c'tor


SingletonNotMatched::SingletonNotMatched
( const Alignment& alignment,
  const int minPos,
  const int maxPos,
  const double standardDeviation ) :
  Singleton( alignment.getData().c_str(),
	     alignment.getQuality().c_str(), alignment ),
  standardDeviation_(standardDeviation)
{
  minPos_=minPos;
  maxPos_=maxPos;
} // ~c'tor




void SingletonNotMatched::populateBounds
( vector<SingletonEvent>& eventStarts,
  vector<SingletonEvent>& eventEnds )
{
    eventStarts.push_back(SingletonEvent(minPos_, this));
    eventEnds.push_back(SingletonEvent(maxPos_, this));
} // ~SingletonNotMatched::populateBounds
  
void SingletonNotMatched::populateMetric( IndelMetricStore& indelMetric )
{
    // If stats are invalid for this read, fall back to uniform distribution
    // over range minPos to maxPos.
    if (standardDeviation_ == 0) {
        const float increment((minPos_ == maxPos_)
                              ? 0
                              : (1.0 / (maxPos_ - minPos_ + 1)));

        for (int i(minPos_); i<=maxPos_; i++) {
            indelMetric[i] += increment;
        }

        return;
    }

    // This computes the PDF of the Gaussian distribution within the interval
    double midPoint=(double)(minPos_+maxPos_)/2.0;
    double c1(1.0/(standardDeviation_*2.506628275)); // 2.506... = root of 2*pi
    double c2(-1.0/(2.0*pow(standardDeviation_,2)));

    double increment;
    for (int i(minPos_);i<=maxPos_;i++)
    {
        increment=midPoint-(double)i;
        increment*=increment;
        increment*=c2;
        increment=exp(increment);
        increment*=c1;
        indelMetric[i]+=increment;
    } // ~for
} // ~SingletonNotMatched::populateMetric( IndelMetricStore& indelMetric )


  SingletonSemiAligned::SingletonSemiAligned
  ( const Alignment& alignment ) :
    Singleton( alignment.getData().c_str(), alignment.getQuality().c_str(),
	       alignment )
    //    alignment_(alignment)
  {

    if (alignment.getMatch().getStrand() == Match::Forward)
    {
      minPos_=alignment.getMatch().getPosition();
      maxPos_=alignment.getMatch().getPosition();
    } // if
    else
    {
      minPos_=alignment.getMatch().getPosition();
      maxPos_=alignment.getMatch().getPosition();
      SequenceUtils::reverseComp(read_);
    } // ~else
    maxPos_+=strlen(alignment.getData().c_str())-1; // add read length

  } // c'tor





  void SingletonSemiAligned::populateBounds
  ( vector<SingletonEvent>& eventStarts,
    vector<SingletonEvent>& eventEnds )
  {
    eventStarts.push_back(SingletonEvent(minPos_, this));
    eventEnds.push_back(SingletonEvent(maxPos_, this));
  } // ~SingletonSemiAligned::populateBounds
  void SingletonSemiAligned::populateMetric( IndelMetricStore& indelMetric )
  {
    float increment;
    // for this we assume the variant has equal chance of sitting anywhere
    // in the read, give each a prob of 1/(read length)
    if (minPos_==maxPos_)
      increment=0;
    else
      increment=1.0/(maxPos_-minPos_+1);
    for (int i(minPos_);i<=maxPos_;i++)
      indelMetric[i]+=increment;

  } // ~SingletonSemiAligned::populateMetric( IndelMetricStore& indelMetric )

}
} // end namespace casava{ namespace { applications