/*****************************************************************************
#   Copyright (C) 1994-2008 by David Gordon.
#   All rights reserved.                           
#                                                                           
#   This software is part of a beta-test version of the Consed/Autofinish
#   package.  It should not be redistributed or
#   used for any commercial purpose, including commercially funded
#   sequencing, without written permission from the author and the
#   University of Washington.
#   
#   This software is provided ``AS IS'' and any express or implied
#   warranties, including, but not limited to, the implied warranties of
#   merchantability and fitness for a particular purpose, are disclaimed.
#   In no event shall the authors or the University of Washington be
#   liable for any direct, indirect, incidental, special, exemplary, or
#   consequential damages (including, but not limited to, procurement of
#   substitute goods or services; loss of use, data, or profits; or
#   business interruption) however caused and on any theory of liability,
#   whether in contract, strict liability, or tort (including negligence
#   or otherwise) arising in any way out of the use of this software, even
#   if advised of the possibility of such damage.
#
#   Building Consed from source is error prone and not simple which is
#   why I provide executables.  Due to time limitations I cannot
#   provide any assistance in building Consed.  Even if you do not
#   modify the source, you may introduce errors due to using a
#   different version of the compiler, a different version of motif,
#   different versions of other libraries than I used, etc.  For this
#   reason, if you discover Consed bugs, I can only offer help with
#   those bugs if you first reproduce those bugs with an executable
#   provided by me--not an executable you have built.
# 
#   Modifying Consed is also difficult.  Although Consed is modular,
#   some modules are used by many other modules.  Thus making a change
#   in one place can have unforeseen effects on many other features.
#   It may takes months for you to notice these other side-effects
#   which may not seen connected at all.  It is not feasable for me to
#   provide help with modifying Consed sources because of the
#   potentially huge amount of time involved.
#
#*****************************************************************************/
#ifndef dGetErrorsFixedAtBase_included
#define dGetErrorsFixedAtBase_included


#include    "consedParameters.h"
#include    <cmath>
using namespace std;


inline double dGetErrorsFixedAtBase( 
                  const double dConsensusErrors,
                  const double dNewReadErrors,
                  const unsigned char ucStrandAndChemistry,
                  const unsigned char ucReadTypesAlreadyCoveringBase ) {


   // This can be the case if we have some doNotFinish tags.  If
   // we don't return here, then dA == 0, and log10( dA ) will be undefined
   if ( dConsensusErrors == 0.0 )
      return( 0.0 );


   double dErrorsFixed;
   double dT;


   if ( ucStrandAndChemistry & ucReadTypesAlreadyCoveringBase )
      dT = 0.8;
   else
      dT = 1.0;

   double dE = MIN( dConsensusErrors, dNewReadErrors );
   double dF = MAX( dConsensusErrors, dNewReadErrors );

   if ( dF == 1.0 ) {
      if ( dConsensusErrors == 1.0 )
         dErrorsFixed = dConsensusErrors - dNewReadErrors;
      else
         dErrorsFixed = 0.0;
   }
   else {

      double dADenominator = ( 1.0 - dE ) * ( 1.0 - dT* dF ) +
         dE* ( dT*dF + 1.0 - dT );
   
      double dA = dE * ( dT * dF + 1.0 - dT ) / dADenominator;
   
      double dBDenominator = ( 1.0 - dE ) * dT * dF + dE * dT * ( 1.0 - dF );
      
      double dB = dE * dT * ( 1.0 - dF ) / dBDenominator;
   
      double dQA = -10.0 * log10( dA );
      double dQB = -10.0 * log10( dB );
   
      double dQResult = ( ( 1.0 - dE ) * ( 1 - dT * dF ) +
                          dE * ( dT * dF + 1.0 - dT ) ) * dQA
         + ( ( 1.0 - dE ) * dT * dF + dE * dT * ( 1.0 - dF ) ) * dQB;
   
      double dNewConsensusErrors = pow( 10.0, -dQResult / 10.0 );

      dErrorsFixed = dConsensusErrors - dNewConsensusErrors;

   }


   // if at this base there is already a read of the same strand and 
   // chemistry, then only half as many errors are fixed

   assert( dErrorsFixed <= dConsensusErrors );
   assert( dErrorsFixed >= 0 );
   
   return( dErrorsFixed );
}


#endif