/*
 * Copyright (c) Medical Research Council 2002. All rights reserved.
 *
 * Permission to use, copy, modify and distribute this software and its
 * documentation for any purpose is hereby granted without fee, provided that
 * this copyright and notice appears in all copies.
 *
 * This file was written as part of the Staden Package at the MRC Laboratory
 * of Molecular Biology, Hills Road, Cambridge, CB2 2QH, United Kingdom.
 *
 * MRC disclaims all warranties with regard to this software.
 *
 */


#include <cassert>
#include <cstdio>                       // For std::printf(), std::sprintf()
#include <cstring>                      // For std::strcpy(), std::memset()
#include <cstdlib>                      // For std::min()
#include <new>                          // For std::bad_alloc
#include <algorithm>                    // For std::max()
#include <list.hpp>                     // For list template
#include <align.hpp>                    // For Alignment object
#include <trace.hpp>                    // For Trace object
#include <mutscan.hpp>                  // For helpers
#include <mutationtag_utils.hpp>        // For CopyTags(), CompTags(), SortTags(), PruneTags()
#include <mutscan_parameters.hpp>       // For MutScanParameter object
#include <mutscan_preprocess.hpp>       // For MutScanPreprocessor object
#include <mutscan_analysis.hpp>         // For MutScanAnalyser object



// #define VERBOSE_DEBUG



/**
   Initialises an empty mutscan_t structure.
*/
void MutScanInit( mutscan_t* ms )
{
   assert(ms != NULL);
   MutScanParameters Parameter;
   std::memset( ms, 0, sizeof(mutscan_t) );
   for( int n=0; n<MUTSCAN_PARAMETERS; n++ )
       ms->Parameter[n] = Parameter[n].Default();
   ms->Initialised = 1;
}



/**
   Frees up all memory allocated by the trace alignment algorithm
   in the mutscan_t structure.
*/
void MutScanDestroy( mutscan_t* ms )
{
   assert(ms != NULL);
   assert(ms->Initialised);
   try
   {
      // Delete all data
      MutScanDestroyResults( ms );
   }
   catch(...)
   {
      // Shouldn't happen, but we musn't throw exceptions outside dll boundary
      assert(0);
   }
}



/**
   Gets the value of the mutscan parameter 'p'.
*/
double MutScanGetParameter( mutscan_t* ms, mutscan_parameter_t p )
{
   assert(ms != NULL);
   assert(ms->Initialised);
   assert(p<MUTSCAN_PARAMETERS);
   return (p<MUTSCAN_PARAMETERS) ? ms->Parameter[p] : 0.0;
}



/**
   Sets the value of the mutscan parameter 'p'.
*/
void MutScanSetParameter( mutscan_t* ms, mutscan_parameter_t p, double v )
{
   assert(ms != NULL);
   assert(ms->Initialised);
   assert(p<MUTSCAN_PARAMETERS);
   if( p < MUTSCAN_PARAMETERS )
       ms->Parameter[p] = v;
}



/**
   Sets the specified reference trace. The clip points are in base units.
   The caller retains ownership of the Read structure and sequence data.
*/
void MutScanSetReference( mutscan_t* ms, mutlib_strand_t d, Read* r, int ql, int qr )
{
   assert(r != NULL);
   assert(ms != NULL);
   assert(ms->Initialised);
   ms->ReferenceTrace[d].ClipL  = ql;
   ms->ReferenceTrace[d].ClipR  = qr;
   ms->ReferenceTrace[d].Trace  = r;
   ms->ReferenceTrace[d].Strand = d;
   ms->ReferenceTrace[d].New    = 1;
}



/**
   Sets the input traces. The clip points are in base units. The caller retains
   ownership of the Read structure data.
*/
void MutScanSetInput( mutscan_t* ms, mutlib_strand_t d, Read* r, int ql, int qr )
{
   assert(r != NULL);
   assert(ms != NULL);
   assert(ms->Initialised);
   MutScanDestroyResults( ms );
   ms->InputTrace.ClipL  = ql;
   ms->InputTrace.ClipR  = qr;
   ms->InputTrace.Trace  = r;
   ms->InputTrace.Strand = d;
   ms->InputTrace.New    = 1;
}



/**
   Returns the result code generated by the Execute() function. This can be
   used to query the result later without having to save the value returned
   by Execute().
*/
mutlib_result_t MutScanGetResultCode( mutscan_t* ms )
{
   assert(ms != NULL);
   assert(ms->Initialised);
   return ms->ResultCode;
}



/**
   Returns a read-only user friendly error message corresponding to the
   last result code. It may contain more useful information such as a
   filename or the particular intput that was invalid. This can be displayed
   to the user if required.
*/
const char* MutScanGetResultString( mutscan_t* ms )
{
   assert(ms != NULL);
   assert(ms->Initialised);
   return ms->ResultString;
}



/**
   Returns the number of tags generated by the mutation scanner and
   available to be read.
*/
int MutScanGetTagCount( mutscan_t* ms )
{
   assert(ms != NULL);
   assert(ms->Initialised);
   return ms->TagCount;
}



/**
   Returns a pointer to the 'nth' tag in the list. If the tag item is empty,
   a null pointer is returned.
*/
mutlib_tag_t* MutScanGetTag( mutscan_t* ms, int n )
{
   assert(ms != NULL);
   assert(ms->Initialised);
   assert(n<ms->TagCount);
   if( n < ms->TagCount )
   {
      assert(ms->Tag != NULL);
      return &ms->Tag[n];
   }
   return 0;
}



/**
   Executes the trace alignment algorithm.
*/
mutlib_result_t MutScanExecute( mutscan_t* ms )
{
    enum { STATE_INITIALISE, STATE_VALIDATE_INPUT, STATE_TRACE_ALIGN,
           STATE_TRACE_DIFFERENCE, STATE_TRACE_PREPROCESS, STATE_MUTATION_ANALYSIS,
           STATE_MUTATION_POSITION, STATE_COVERAGE_TAG, STATE_MUTATION_TAG,
           STATE_EXIT };
    int                 n;
    tracealign_t        ta;
    mutlib_result_t     Result;
    mutlib_strand_t     Strand = MUTLIB_STRAND_FORWARD; // silence warning
    MutScanParameters   Parameter;
    Alignment           Aligner;
    MutScanAnalyser     Analyser;
    int                 SamplesPerBase    = 1;
    MutScanPreprocessor Preprocessor[2];
    Trace               AlignedTrace[2];
    int                 AlignedTraceClipL[2];
    int                 AlignedTraceClipR[2];
    Trace*              DifferenceTrace   = 0;
    bool                Done              = false;
    int                 State             = STATE_INITIALISE;
    const int           BASE_MARGIN_LEFT  = 2;
    const int           BASE_MARGIN_RIGHT = 6;
    assert(ms != NULL);
    assert(ms->Initialised);
    try
    {
        while(!Done)
        {
            switch(State)
            {
                case STATE_INITIALISE:
                    // Destroy old results
                    TraceAlignInit( &ta );
                    MutScanDestroyResults( ms );
                    Strand              = ms->InputTrace.Strand;
                    ms->ResultCode      = MUTLIB_RESULT_SUCCESS;
                    ms->ResultString    = new char [256];
                    ms->ResultString[0] = 0;
                    State = STATE_VALIDATE_INPUT;
                    break;



                case STATE_VALIDATE_INPUT:
                    // Check input values, adjust rhs clip point to catch indels
                    if( ms->InputTrace.New && (ms->InputTrace.ClipR>0) )
                        ms->InputTrace.ClipR += BASE_MARGIN_RIGHT;
                    for( n=0; n<MUTSCAN_PARAMETERS; n++ )
                        Parameter[n].Value( ms->Parameter[n] );
                    if( MutScanValidateInput(ms,Parameter) )
                        State = STATE_EXIT;
                    else
                        State = STATE_TRACE_ALIGN;
                    ms->InputTrace.New = 0;
                    break;



                case STATE_TRACE_ALIGN:
                    // Align the reference and input traces
                    TraceAlignSetReference( &ta, Strand, ms->ReferenceTrace[Strand].Trace, ms->ReferenceTrace[Strand].ClipL, ms->ReferenceTrace[Strand].ClipR );
                    TraceAlignSetInput( &ta, Strand, ms->InputTrace.Trace, ms->InputTrace.ClipL, ms->InputTrace.ClipR );
                    if( TraceAlignExecute(&ta) != MUTLIB_RESULT_SUCCESS )
                    {
                        ms->ResultCode = TraceAlignGetResultCode( &ta );
                        std::strcpy( ms->ResultString, TraceAlignGetResultString(&ta) );
                        State = STATE_EXIT;
                        break;
                    }
                    for( int n=0; n<2; n++ )
                        AlignedTrace[n].Wrap( TraceAlignGetAlignment(&ta, static_cast<mutlib_input_t>(n), &AlignedTraceClipL[n], &AlignedTraceClipR[n]), false );
                    State = STATE_TRACE_PREPROCESS;
                    break;



                case STATE_TRACE_PREPROCESS:
                    // Preprocess the traces. We scan beyond the right hand clip point
                    // to find indels. Where possible, we also start scanning from the
                    // left margin to avoid discontinuties at the edges of the alignment.
                    {State = STATE_TRACE_DIFFERENCE;
                    SamplesPerBase = AlignedTrace[MUTLIB_INPUT_REFERENCE].IntervalMode();
                    int SampleMarginLeft  = BASE_MARGIN_LEFT * SamplesPerBase;
                    for( int n=0; n<2; n++ )
                    {
                        int SampleMarginRight = AlignedTrace[n].Samples() - SamplesPerBase;
                        Result = Preprocessor[n].Execute( ms, AlignedTrace[n], n, SampleMarginLeft, SampleMarginRight );
                        if( Result != MUTLIB_RESULT_SUCCESS )
                        {
                            State = STATE_EXIT;
                            break;
                        }
                    }
                    break;}



                case STATE_TRACE_DIFFERENCE:
                    // Subtract the two aligned traces
                    DifferenceTrace = AlignedTrace[MUTLIB_INPUT].Subtract( AlignedTrace[MUTLIB_INPUT_REFERENCE] );
                    if( DifferenceTrace )
                    {
                        DifferenceTrace->Floor( 35 );
                        DifferenceTrace->FloorHalfwaves();
                        DifferenceTrace->FloorNarrowPeaks( SamplesPerBase/2 );
                        DifferenceTrace->FillGaps();
                        #ifdef VERBOSE_DEBUG
                        char buffer[256];
                        std::strcpy( buffer, AlignedTrace[MUTLIB_INPUT].Name() );
                        std::strcat( buffer, "_diff.ztr" );
                        DifferenceTrace->SaveAs( buffer );
                        #endif
                    }
                    State = STATE_MUTATION_ANALYSIS;
                    break;



                case STATE_MUTATION_ANALYSIS:
                    // Search for and analyse the traces for mutations
                    Result = Analyser.Execute( ms, Preprocessor, AlignedTrace, DifferenceTrace );
                    if( Result != MUTLIB_RESULT_SUCCESS )
                        State = STATE_EXIT;
                    else
                        State = STATE_MUTATION_POSITION;
                    break;



                case STATE_MUTATION_POSITION: {
                    // Map mutation positions back onto original traces
                    MutationTag* pTag = Analyser.MutationTagList.First();
                    while( pTag )
                    {
                        if( !pTag->Marked() )
                        {
                            // Trace function returns zero-based baseno, so we have to add 1
                            int Base = AlignedTrace[MUTLIB_INPUT].BaseNumberFromSamplePosition( pTag->Position(1) );
                            Base += 1;
                            Base += AlignedTraceClipL[MUTLIB_INPUT];
                            pTag->Position( 0, Base );
                        }
                        pTag = Analyser.MutationTagList.Next();
                    }
                    State = STATE_COVERAGE_TAG;
                    break; }


                case STATE_COVERAGE_TAG: {
                    // Create a coverage mutation tag
                    MutationTag* pTag = new MutationTag( "MCOV" );
                    pTag->Strand( Strand );
                    pTag->Position( 0, AlignedTraceClipL[MUTLIB_INPUT]+1 );
                    pTag->Position( 1, AlignedTraceClipR[MUTLIB_INPUT]-1 );
                    Analyser.MutationTagList.Append( pTag );
                    State = STATE_MUTATION_TAG;
                    break; }



                case STATE_MUTATION_TAG:
                    // Create a mutation tags array for output
                    if( Analyser.MutationTagList.Count() > 0 )
                    {
                        // Copy tags from list into a C-style array
                        SimpleArray<mutlib_tag_t> OutTag;
                        OutTag.Create( Analyser.MutationTagList.Count() );
                        CopyTags( OutTag, Analyser.MutationTagList );
                        if( (Strand==MUTLIB_STRAND_REVERSE) && (Parameter[MUTSCAN_PARAMETER_COMPLEMENT_TAGS].Value()>0.0) )
                            CompTags( OutTag );
                        SortTags( OutTag );
                        PruneTags( OutTag );
                        if( OutTag.Length() > (ms->Parameter[MUTSCAN_PARAMETER_ALIGNFAIL_THRESHOLD]+1) )
                        {
                            ms->ResultCode = MUTLIB_RESULT_ALIGNMENT_FAILURE;
                            std::sprintf( ms->ResultString, "Trace alignment failed for %s\n", ms->InputTrace.Trace->trace_name );
                        }
                        ms->Tag      = OutTag.Raw();
                        ms->TagCount = OutTag.Length();
                        OutTag.AutoDestroy( false );
                    }
                    State = STATE_EXIT;
                    break;



                case STATE_EXIT:
                    // Quit
                    Done = true;
                    break;
            }
        }
    }
    catch( std::bad_alloc& )
    {
        // Memory allocation error
        ms->ResultCode = MUTLIB_RESULT_OUT_OF_MEMORY;
        std::strcpy( ms->ResultString, "Not enough memory available to complete the operation.\n" );
    }
    catch(...)
    {
        // Unknown exceptions
        ms->ResultCode = MUTLIB_RESULT_UNEXPECTED_EXCEPTION;
        std::strcpy( ms->ResultString, "An unexpected fatal exception has occurred, please "
                     "report the details to staden-package@mrc-lmb.cam.ac.uk.\n" );
    }



    // Exit
    if( DifferenceTrace )
        delete DifferenceTrace;
    TraceAlignDestroy( &ta );
    return ms->ResultCode;
}