gatk-3.8/public/VectorPairHMM/src/main/c++/template.h

/*Copyright (c) 2012 The Broad Institute

*Permission is hereby granted, free of charge, to any person
*obtaining a copy of this software and associated documentation
*files (the "Software"), to deal in the Software without
*restriction, including without limitation the rights to use,
*copy, modify, merge, publish, distribute, sublicense, and/or sell
*copies of the Software, and to permit persons to whom the
*Software is furnished to do so, subject to the following
*conditions:

*The above copyright notice and this permission notice shall be
*included in all copies or substantial portions of the Software.

*THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
*EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
*OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
*NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
*HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
*WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
*FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR
*THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/


#ifndef TEMPLATES_H_
#define TEMPLATES_H_

#include "headers.h"

#define MM 0
#define GapM 1
#define MX 2
#define XX 3
#define MY 4
#define YY 5

//#define MROWS  500
//#define MCOLS  1000

#define CAT(X,Y) X####Y
#define CONCAT(X,Y) CAT(X,Y)

#define ALIGNED __attribute__((aligned(32)))

typedef union __attribute__((aligned(32))) {
        ALIGNED __m256 ALIGNED d;
        ALIGNED __m128i ALIGNED s[2];
        ALIGNED float  ALIGNED f[8];
        ALIGNED __m256i ALIGNED i;
} ALIGNED mix_F ALIGNED;

typedef union __attribute__((aligned(32))) {
        ALIGNED __m128 ALIGNED d;
        ALIGNED __m64 ALIGNED s[2];
        ALIGNED float  ALIGNED f[4];
        ALIGNED __m128i ALIGNED i;
} ALIGNED mix_F128 ALIGNED;

typedef union ALIGNED {
  __m128i vec ;
  __m128 vecf ;
  uint32_t masks[4] ;
} MaskVec_F ;

typedef union ALIGNED {
  __m64 vec ;
  __m64 vecf ;
  uint32_t masks[2] ;
} MaskVec_F128 ;

typedef union ALIGNED
{
        ALIGNED __m128i ALIGNED i;
        ALIGNED __m128 ALIGNED f;
} ALIGNED IF_128f ALIGNED;

typedef union ALIGNED
{
        ALIGNED int    ALIGNED i;
        ALIGNED float  ALIGNED f;
} ALIGNED IF_32 ALIGNED;

typedef union __attribute__((aligned(32))) {
        ALIGNED __m256d ALIGNED d;
        ALIGNED __m128i ALIGNED s[2];
        ALIGNED double  ALIGNED f[4];
        ALIGNED __m256i ALIGNED i;
} ALIGNED mix_D ALIGNED;

typedef union __attribute__((aligned(32))) {
        ALIGNED __m128d ALIGNED d;
        ALIGNED __m64 ALIGNED s[2];
        ALIGNED double  ALIGNED f[2];
        ALIGNED __m128i ALIGNED i;
} ALIGNED mix_D128 ALIGNED;

typedef union ALIGNED {
  __m128i vec ;
  __m128d vecf ;
  uint64_t masks[2] ;
} MaskVec_D ;

typedef union ALIGNED {
  __m64 vec ;
  __m64 vecf ;
  uint64_t masks[1] ;
} MaskVec_D128 ;

typedef union ALIGNED
{
        ALIGNED __m128i ALIGNED i;
        ALIGNED __m128d ALIGNED f;
} ALIGNED IF_128d ALIGNED;

typedef union ALIGNED
{
        ALIGNED int64_t ALIGNED i;
        ALIGNED double  ALIGNED f;
} ALIGNED IF_64 ALIGNED;


#define MAX_QUAL 254
#define MAX_JACOBIAN_TOLERANCE 8.0
#define JACOBIAN_LOG_TABLE_STEP 0.0001
#define JACOBIAN_LOG_TABLE_INV_STEP (1.0 / JACOBIAN_LOG_TABLE_STEP)
#define MAXN 70000
#define LOG10_CACHE_SIZE  (4*MAXN)  // we need to be able to go up to 2*(2N) when calculating some of the coefficients
#define JACOBIAN_LOG_TABLE_SIZE ((int) (MAX_JACOBIAN_TOLERANCE / JACOBIAN_LOG_TABLE_STEP) + 1)

template<class NUMBER>
struct ContextBase
{
  public:
    NUMBER ph2pr[128];
    NUMBER INITIAL_CONSTANT;
    NUMBER LOG10_INITIAL_CONSTANT;
    NUMBER RESULT_THRESHOLD; 

    static bool staticMembersInitializedFlag;
    static NUMBER jacobianLogTable[JACOBIAN_LOG_TABLE_SIZE];
    static NUMBER matchToMatchProb[((MAX_QUAL + 1) * (MAX_QUAL + 2)) >> 1];

    static void initializeStaticMembers()
    {
      if(!staticMembersInitializedFlag)
      {
        //Order of calls important - Jacobian first, then MatchToMatch
        initializeJacobianLogTable();
        initializeMatchToMatchProb();
        staticMembersInitializedFlag = true;
      }
    }

    static void deleteStaticMembers()
    {
      if(staticMembersInitializedFlag)
      {
        staticMembersInitializedFlag = false;
      }
    }

    //Called only once during library load - don't bother to optimize with single precision fp
    static void initializeJacobianLogTable()
    {
      for (int k = 0; k < JACOBIAN_LOG_TABLE_SIZE; k++) {
        jacobianLogTable[k] = (NUMBER)(log10(1.0 + pow(10.0, -((double) k) * JACOBIAN_LOG_TABLE_STEP)));
      }
    }

    //Called only once per library load - don't bother optimizing with single fp
    static void initializeMatchToMatchProb()
    {
      double LN10 = log(10);
      double INV_LN10 = 1.0/LN10;
      for (int i = 0, offset = 0; i <= MAX_QUAL; offset += ++i)
        for (int j = 0; j <= i; j++) {
          double log10Sum = approximateLog10SumLog10(-0.1*i, -0.1*j);
          double matchToMatchLog10 =
            log1p(-std::min(1.0,pow(10,log10Sum))) * INV_LN10;
          matchToMatchProb[offset + j] = (NUMBER)(pow(10,matchToMatchLog10));
        }
    }
    //Called during computation - use single precision where possible
    static int fastRound(NUMBER d) {
      return (d > ((NUMBER)0.0)) ? (int) (d + ((NUMBER)0.5)) : (int) (d - ((NUMBER)0.5));
    }
    //Called during computation - use single precision where possible
    static NUMBER approximateLog10SumLog10(NUMBER small, NUMBER big) {
      // make sure small is really the smaller value
      if (small > big) {
        NUMBER t = big;
        big = small;
        small = t;
      }

      if (isinf(small) == -1 || isinf(big) == -1)
        return big;

      NUMBER diff = big - small;
      if (diff >= ((NUMBER)MAX_JACOBIAN_TOLERANCE))
        return big;

      // OK, so |y-x| < tol: we use the following identity then:
      // we need to compute log10(10^x + 10^y)
      // By Jacobian logarithm identity, this is equal to
      // max(x,y) + log10(1+10^-abs(x-y))
      // we compute the second term as a table lookup with integer quantization
      // we have pre-stored correction for 0,0.1,0.2,... 10.0
      int ind = fastRound((NUMBER)(diff * ((NUMBER)JACOBIAN_LOG_TABLE_INV_STEP))); // hard rounding
      return big + jacobianLogTable[ind];
    }
};

template<class NUMBER>
struct Context : public ContextBase<NUMBER>
{};

template<>
struct Context<double> : public ContextBase<double>
{
  Context():ContextBase<double>()
  {
    for (int x = 0; x < 128; x++)
      ph2pr[x] = pow(10.0, -((double)x) / 10.0);

    INITIAL_CONSTANT = ldexp(1.0, 1020.0);
    LOG10_INITIAL_CONSTANT = log10(INITIAL_CONSTANT);
    RESULT_THRESHOLD = 0.0;
  }

  double LOG10(double v){ return log10(v); }
  inline double POW(double b, double e) { return pow(b,e); }

  static double _(double n){ return n; }
  static double _(float n){ return ((double) n); }
};

template<>
struct Context<float> : public ContextBase<float>
{
  Context() : ContextBase<float>()
  {
    for (int x = 0; x < 128; x++)
    {
      ph2pr[x] = powf(10.f, -((float)x) / 10.f);
    }

    INITIAL_CONSTANT = ldexpf(1.f, 120.f);
    LOG10_INITIAL_CONSTANT = log10f(INITIAL_CONSTANT);
    RESULT_THRESHOLD = ldexpf(1.f, -110.f);
  }

  float LOG10(float v){ return log10f(v); }
  inline float POW(float b, float e) { return powf(b,e); }

  static float _(double n){ return ((float) n); }
  static float _(float n){ return n; }
};

#define SET_MATCH_TO_MATCH_PROB(output, insQual, delQual)                       \
{                                                                               \
  int minQual = delQual;                                                        \
  int maxQual = insQual;                                                        \
  if (insQual <= delQual)                                                       \
  {                                                                             \
    minQual = insQual;                                                          \
    maxQual = delQual;                                                          \
  }                                                                             \
  (output) = (MAX_QUAL < maxQual) ?                                             \
  ((NUMBER)1.0) - ctx.POW(((NUMBER)10), ctx.approximateLog10SumLog10(((NUMBER)-0.1)*minQual, ((NUMBER)-0.1)*maxQual))       \
  : ctx.matchToMatchProb[((maxQual * (maxQual + 1)) >> 1) + minQual];           \
}

typedef struct
{
        int rslen, haplen;
	/*int *q, *i, *d, *c;*/
	/*int q[MROWS], i[MROWS], d[MROWS], c[MROWS];*/
	char *q, *i, *d, *c;
        char *hap, *rs;
	int *ihap;
	int *irs;
} testcase;

int normalize(char c);
int read_testcase(testcase *tc, FILE* ifp=0);


#define MIN_ACCEPTED 1e-28f
#define NUM_DISTINCT_CHARS 5
#define AMBIG_CHAR 4

class ConvertChar {

  static uint8_t conversionTable[255] ;

public:

  static void init() {
    assert (NUM_DISTINCT_CHARS == 5) ;
    assert (AMBIG_CHAR == 4) ;

    conversionTable['A'] = 0 ;
    conversionTable['C'] = 1 ;
    conversionTable['T'] = 2 ;
    conversionTable['G'] = 3 ;
    conversionTable['N'] = 4 ;
  }

  static inline uint8_t get(uint8_t input) {
    return conversionTable[input] ;
  }

};


#endif