/* ===========================================================================
 *
 *                            PUBLIC DOMAIN NOTICE
 *               National Center for Biotechnology Information
 *
 *  This software/database is a "United States Government Work" under the
 *  terms of the United States Copyright Act.  It was written as part of
 *  the author's official duties as a United States Government employee and
 *  thus cannot be copyrighted.  This software/database is freely available
 *  to the public for use. The National Library of Medicine and the U.S.
 *  Government have not placed any restriction on its use or reproduction.
 *
 *  Although all reasonable efforts have been taken to ensure the accuracy
 *  and reliability of the software and data, the NLM and the U.S.
 *  Government do not and cannot warrant the performance or results that
 *  may be obtained by using this software or data. The NLM and the U.S.
 *  Government disclaim all warranties, express or implied, including
 *  warranties of performance, merchantability or fitness for any particular
 *  purpose.
 *
 *  Please cite the author in any work or product based on this material.
 *
 * ===========================================================================
 *
 * Author:  Christiam Camacho
 *
 */

/** @file setup_factory.cpp
 * NOTE: This file contains work in progress and the APIs are likely to change,
 * please do not rely on them until this notice is removed.
 */

#include <ncbi_pch.hpp>
#include <algo/blast/api/setup_factory.hpp>
#include <algo/blast/api/uniform_search.hpp>    // for CSearchDatabase
#include <algo/blast/api/blast_options.hpp>
#include <algo/blast/api/seqsrc_seqdb.hpp>      // for SeqDbBlastSeqSrcInit
#include <algo/blast/api/blast_mtlock.hpp>      // for Blast_DiagnosticsInitMT
#include <algo/blast/api/blast_dbindex.hpp>

#include "blast_aux_priv.hpp"
#include "blast_memento_priv.hpp"
#include "blast_setup.hpp"

// SeqAlignVector building
#include "blast_seqalign.hpp"
#include <objects/seqloc/Seq_interval.hpp>
#include <objects/seqalign/seqalign__.hpp>
#include <serial/iterator.hpp>

// CORE BLAST includes
#include <algo/blast/core/blast_setup.h>
#include <algo/blast/core/blast_hspstream.h>
#include <algo/blast/core/hspfilter_collector.h>
#include <algo/blast/core/hspfilter_besthit.h>
#include <algo/blast/core/hspfilter_culling.h>
#include <algo/blast/core/hspfilter_mapper.h>

#include <sstream>
#include "../core/jumper.h"

/** @addtogroup AlgoBlast
 *
 * @{
 */

BEGIN_NCBI_SCOPE
USING_SCOPE(objects);
BEGIN_SCOPE(blast)

CRef<CBlastRPSInfo>
CSetupFactory::CreateRpsStructures(const string& rps_dbname,
                                   CRef<CBlastOptions> options)
{
    CBlastRPSInfo::EOpenFlags mode = 
             (options->GetCompositionBasedStats() == eNoCompositionBasedStats) 
             ? CBlastRPSInfo::fRpsBlast : CBlastRPSInfo::fRpsBlastWithCBS;
    CRef<CBlastRPSInfo> retval(new CBlastRPSInfo(rps_dbname, mode));
    options->SetMatrixName(retval->GetMatrixName());
    options->SetGapOpeningCost(retval->GetGapOpeningCost());
    options->SetGapExtensionCost(retval->GetGapExtensionCost());
    return retval;
}

/** 
 * @brief Auxiliary function to extract the Seq-ids from the ILocalQueryData
 * and bundle them in a Packed-seqint
 * 
 * @param query_data input query data [in]
 * 
 * @return Packed-seqint with query Seq-locs
 */
static
CRef<CPacked_seqint> s_LocalQueryData2Packed_seqint(ILocalQueryData& query_data)
{
    const int kNumQueries = static_cast<int>(query_data.GetNumQueries());
    if (kNumQueries == 0) {
        return CRef<CPacked_seqint>();
    }

    CRef<CPacked_seqint> retval(new CPacked_seqint);
    for (int i = 0; i < kNumQueries; i++) {
        CConstRef<CSeq_id> id(query_data.GetSeq_loc(i)->GetId());
        if (query_data.GetSeq_loc(i)->IsInt()) {
            retval->AddInterval(query_data.GetSeq_loc(i)->GetInt());
        } else if (id.NotEmpty()) {
            TSeqPos len = 0;
            try { len = static_cast<TSeqPos>(query_data.GetSeqLength(i)); }
            catch (...) { 
                /* exception means that it's an invalid seqid, so we do
                 * nothing, the error message should be captured elsewhere */
            }
            retval->AddInterval(*id, 0, len);
        }
    }

    return retval;
}

BlastScoreBlk*
CSetupFactory::CreateScoreBlock(const CBlastOptionsMemento* opts_memento,
                                CRef<ILocalQueryData> query_data,
                                BlastSeqLoc** lookup_segments,
                                TSearchMessages& search_messages,
                                TSeqLocInfoVector* masked_query_regions,
                                const CBlastRPSInfo* rps_info)
{
    _ASSERT(opts_memento);

    double rps_scale_factor(1.0);

    if (rps_info) {
        _ASSERT(Blast_ProgramIsRpsBlast(opts_memento->m_ProgramType));
        rps_scale_factor = rps_info->GetScalingFactor();
    }

    CBlast_Message blast_msg;
    CBlastMaskLoc core_masked_query_regions;

    BlastQueryInfo* query_info = query_data->GetQueryInfo();
    BLAST_SequenceBlk* queries = query_data->GetSequenceBlk();

    // reset composition based statistics for the identity scoring matrix
    // (we do not have all necessary data to compute CBS)
    if (opts_memento->m_ScoringOpts->matrix &&
        strcasecmp("identity", opts_memento->m_ScoringOpts->matrix) == 0 &&
        opts_memento->m_ExtnOpts->compositionBasedStats > 0) {

        opts_memento->m_ExtnOpts->compositionBasedStats = 0;
        ostringstream os;
        os << "Composition-based statistics cannot be used with the IDENTITY "
           << "matrix, resetting the composition-based statistics option to 0";
        search_messages.AddMessageAllQueries(eBlastSevWarning, 0, os.str());
    }

    BlastScoreBlk* retval(0);
    Int2 status = BLAST_MainSetUp(opts_memento->m_ProgramType,
                                  opts_memento->m_QueryOpts,
                                  opts_memento->m_ScoringOpts,
                                  queries,
                                  query_info,
                                  rps_scale_factor,
                                  lookup_segments,
                                  &core_masked_query_regions,
                                  &retval,
                                  &blast_msg,
                                  &BlastFindMatrixPath);

    Blast_Message2TSearchMessages(blast_msg.Get(), query_info, search_messages);
    if (status != 0 && 
        (!(blast_msg.Get()) || (blast_msg.Get() && blast_msg.Get()->severity == eBlastSevError)))
    {
        	string msg;
        	if (search_messages.HasMessages()) {
            		msg = search_messages.ToString();
        	} else {
            		msg = "BLAST_MainSetUp failed (" + NStr::IntToString(status) + 
            			" error code)";
                }
                // Clean up leaks before we throw
                retval = BlastScoreBlkFree(retval);
                *lookup_segments =  BlastSeqLocFree(*lookup_segments);
                NCBI_THROW(CBlastException, eCoreBlastError, msg);
    }

    if (masked_query_regions) {
        CRef<CPacked_seqint> query_locations = 
            s_LocalQueryData2Packed_seqint(*query_data);
        Blast_GetSeqLocInfoVector(opts_memento->m_ProgramType,
                                  *query_locations,
                                  core_masked_query_regions,
                                  *masked_query_regions);
    }

    return retval;
}

LookupTableWrap*
CSetupFactory::CreateLookupTable(CRef<ILocalQueryData> query_data,
                                 const CBlastOptionsMemento* opts_memento,
                                 BlastScoreBlk* score_blk,
                                 CRef< CBlastSeqLocWrap > lookup_segments_wrap,
                                 const CBlastRPSInfo* rps_info,
                                 BlastSeqSrc* seqsrc,
                                 size_t num_threads)
{
    BLAST_SequenceBlk* queries = query_data->GetSequenceBlk();
    CBlast_Message blast_msg;
    LookupTableWrap* retval(0);

    BlastSeqLoc * lookup_segments = lookup_segments_wrap->getLocs();

    Int2 status = LookupTableWrapInit_MT(queries,
                                         opts_memento->m_LutOpts,
                                         opts_memento->m_QueryOpts,
                                         lookup_segments,
                                         score_blk,
                                         &retval,
                                         rps_info ? (*rps_info)() : 0,
                                         &blast_msg,
                                         seqsrc,
                                         static_cast<Uint4>(num_threads));
    if (status != 0) {
         TSearchMessages search_messages;
         Blast_Message2TSearchMessages(blast_msg.Get(), 
                                           query_data->GetQueryInfo(), 
                                           search_messages);
         string msg;
         if (search_messages.HasMessages()) {
              msg = search_messages.ToString();
         } else {
              msg = "LookupTableWrapInit failed (" + 
                   NStr::IntToString(status) + " error code)";
         }
         NCBI_THROW(CBlastException, eCoreBlastError, msg);
    }

    // For PHI BLAST, save information about pattern occurrences in query in
    // the BlastQueryInfo structure
    if (Blast_ProgramIsPhiBlast(opts_memento->m_ProgramType)) {
        SPHIPatternSearchBlk* phi_lookup_table
            = (SPHIPatternSearchBlk*) retval->lut;
        status = Blast_SetPHIPatternInfo(opts_memento->m_ProgramType,
                                phi_lookup_table,
                                queries,
                                lookup_segments,
                                query_data->GetQueryInfo(), 
                                &blast_msg);
        if (status != 0) {  
             TSearchMessages search_messages;
             Blast_Message2TSearchMessages(blast_msg.Get(), 
                                           query_data->GetQueryInfo(), 
                                           search_messages);
             string msg;
             if (search_messages.HasMessages()) {
                 msg = search_messages.ToString();
             } else {
                 msg = "Blast_SetPHIPatternInfo failed (" + 
                     NStr::IntToString(status) + " error code)";
             }
             NCBI_THROW(CBlastException, eCoreBlastError, msg);
        }
    }

    if (seqsrc) {
        GetDbIndexSetQueryInfoFn()( retval, lookup_segments_wrap);
    }

    return retval;
}

BlastDiagnostics*
CSetupFactory::CreateDiagnosticsStructure()
{
    return Blast_DiagnosticsInit();
}

BlastDiagnostics*
CSetupFactory::CreateDiagnosticsStructureMT()
{
    return Blast_DiagnosticsInitMT(Blast_CMT_LOCKInit());
}

BlastHSPStream*
CSetupFactory::CreateHspStream(const CBlastOptionsMemento* opts_memento,
                               size_t number_of_queries,
                               BlastHSPWriter* writer)
{
    _ASSERT(opts_memento);
    return BlastHSPStreamNew(opts_memento->m_ProgramType, 
                             opts_memento->m_ExtnOpts, TRUE,
                             static_cast<Int4>(number_of_queries), writer);
}

BlastHSPWriter*
CSetupFactory::CreateHspWriter(const CBlastOptionsMemento* opts_memento,
                               BLAST_SequenceBlk* query,
                               BlastQueryInfo* query_info)
{
    BlastHSPWriterInfo* writer_info = NULL;
    
    const BlastHSPFilteringOptions* filt_opts =
        opts_memento->m_HitSaveOpts->hsp_filt_opt;
    if(filt_opts && (filt_opts->best_hit) && filt_opts->best_hit->score_edge == 0) {
    	//We want to only do score edge in prelimiary stage, skip if score edge in 0
    	filt_opts = NULL;
    }

    if (Blast_ProgramIsMapping(opts_memento->m_ProgramType)) {

        BlastHSPMapperParams* params =
            BlastHSPMapperParamsNew(opts_memento->m_HitSaveOpts,
                                    opts_memento->m_ScoringOpts);

        writer_info = BlastHSPMapperInfoNew(params);
    }
    else if (filt_opts) {
        if (filt_opts->best_hit && (filt_opts->best_hit_stage & ePrelimSearch))
        {
            BlastHSPBestHitParams* params = 
                BlastHSPBestHitParamsNew(opts_memento->m_HitSaveOpts,
                     filt_opts->best_hit,
                     opts_memento->m_ExtnOpts->compositionBasedStats, 
                     opts_memento->m_ScoringOpts->gapped_calculation);
            // Disable overhang in prelimiary stage
            params->overhang = 0;
            writer_info = BlastHSPBestHitInfoNew(params);
        }
        else if (filt_opts->culling_opts && 
                 (filt_opts->culling_stage & ePrelimSearch))
        {
            BlastHSPCullingParams* params = 
                BlastHSPCullingParamsNew(opts_memento->m_HitSaveOpts,
                     filt_opts->culling_opts,
                     opts_memento->m_ExtnOpts->compositionBasedStats,
                     opts_memento->m_ScoringOpts->gapped_calculation);
            if(params->culling_max > 1){
            	params->culling_max += 3;
            }
            writer_info = BlastHSPCullingInfoNew(params);
        }
        else {
        	BlastHSPCollectorParams * params =
        			BlastHSPCollectorParamsNew(opts_memento->m_HitSaveOpts,
                       opts_memento->m_ExtnOpts->compositionBasedStats,
                       opts_memento->m_ScoringOpts->gapped_calculation);
        	writer_info = BlastHSPCollectorInfoNew(params);

        }
    } else {
        /* Use the collector filtering algorithm as the default */
        BlastHSPCollectorParams * params = 
            BlastHSPCollectorParamsNew(opts_memento->m_HitSaveOpts, 
                       opts_memento->m_ExtnOpts->compositionBasedStats, 
                       opts_memento->m_ScoringOpts->gapped_calculation);
        writer_info = BlastHSPCollectorInfoNew(params);
    }
    
    BlastHSPWriter* retval = BlastHSPWriterNew(&writer_info, query_info, query);
    _ASSERT(writer_info == NULL);
    return retval;
}

BlastHSPPipe*
CSetupFactory::CreateHspPipe(const CBlastOptionsMemento* opts_memento,
                             BlastQueryInfo* query_info)
{
    _ASSERT(opts_memento);

    BlastHSPPipe* retval = NULL;
    BlastHSPPipeInfo* pipe_info = NULL;
    
    const BlastHSPFilteringOptions* filt_opts =
        opts_memento->m_HitSaveOpts->hsp_filt_opt;
    if (filt_opts) {
        if (filt_opts->best_hit && 
            (filt_opts->best_hit_stage & eTracebackSearch)) {
            BlastHSPBestHitParams* params = 
                BlastHSPBestHitParamsNew(opts_memento->m_HitSaveOpts,
                     filt_opts->best_hit,
                     opts_memento->m_ExtnOpts->compositionBasedStats,
                     opts_memento->m_ScoringOpts->gapped_calculation);
            BlastHSPPipeInfo_Add(&pipe_info, 
                                 BlastHSPBestHitPipeInfoNew(params));
        } else if (filt_opts->culling_opts &&
                   (filt_opts->culling_stage & eTracebackSearch)) {
            BlastHSPCullingParams* params = 
                BlastHSPCullingParamsNew(opts_memento->m_HitSaveOpts,
                     filt_opts->culling_opts,
                     opts_memento->m_ExtnOpts->compositionBasedStats,
                     opts_memento->m_ScoringOpts->gapped_calculation);
            BlastHSPPipeInfo_Add(&pipe_info, 
                                 BlastHSPCullingPipeInfoNew(params));
        }
    } else {
        ; /* the default is to use no pipes */
    }
    
    retval = BlastHSPPipeNew(&pipe_info, query_info);
    _ASSERT(pipe_info == NULL);
    return retval;
}

BlastSeqSrc*
CSetupFactory::CreateBlastSeqSrc(const CSearchDatabase& db)
{
    return CreateBlastSeqSrc(db.GetSeqDb(), 
                             db.GetFilteringAlgorithm(),
                             db.GetMaskType());
}

BlastSeqSrc*
CSetupFactory::CreateBlastSeqSrc(CSeqDB * db, int filt_algo,
                                 ESubjectMaskingType mask_type)
{
    BlastSeqSrc* retval = SeqDbBlastSeqSrcInit(db, filt_algo, mask_type);
    char* error_str = BlastSeqSrcGetInitError(retval);
    if (error_str) {
        string msg(error_str);
        sfree(error_str);
        retval = BlastSeqSrcFree(retval);
        NCBI_THROW(CBlastException, eSeqSrcInit, msg);
    }
    return retval;
}

void
CSetupFactory::InitializeMegablastDbIndex(CRef<CBlastOptions> options)
{
    _ASSERT(options->GetUseIndex());

    if (options->GetMBIndexLoaded()) {
        return;
    }

    string errstr = "";
    bool partial( false );

    if( options->GetProgramType() != eBlastTypeBlastn
        && options->GetProgramType() != eBlastTypeMapping) {
        errstr = "Database indexing is available for blastn only.";
    }
    else if( options->GetMBTemplateLength() > 0 ) {
        errstr = "Database indexing is not available for discontiguous ";
        errstr += "searches.";
    }
    else {
        errstr = DbIndexInit( 
                options->GetIndexName(), 
                options->GetIsOldStyleMBIndex(), partial, options->GetWordSize() );
    }

    if( errstr != "" ) {
            ERR_POST_EX(1, 1, Warning << errstr << " Database index will not be used." );
            options->SetUseIndex( false );
            return;
    }

    options->SetMBIndexLoaded();
    options->SetLookupTableType( 
            partial ? eMixedMBLookupTable : eIndexedMBLookupTable );
}

SInternalData::SInternalData()
{
    m_Queries = 0;
    m_QueryInfo = 0;
}

SDatabaseScanData::SDatabaseScanData()
   : kNoPhiBlastPattern(-1)
{
   m_NumPatOccurInDB = kNoPhiBlastPattern; 
}

END_SCOPE(blast)
END_NCBI_SCOPE

/* @} */

