1 // Copyright 2008 Michiaki Hamada
2 // Modified 2015 Yutaro Konta
3 
4 // This class calculates the scale factor (lambda), and the letter
5 // probabilities, that are implicit in a scoring matrix.  The
6 // calculation might fail for various reasons, putting it into a
7 // "bad/undefined" state.
8 
9 // If the score matrix is symmetric, then the two sets of letter
10 // probabilities should be identical.  With this code, they might
11 // differ minutely from exact identity.
12 
13 #ifndef LAMBDA_CALCULATOR_HH
14 #define LAMBDA_CALCULATOR_HH
15 
16 #include <vector>
17 
18 namespace cbrc{
19 
20 typedef const int *const_int_ptr;
21 
22 class LambdaCalculator{
23  public:
LambdaCalculator()24   LambdaCalculator() { setBad(); }
25 
26   void calculate(const const_int_ptr *matrix, int alphSize);
27 
28   // Put us in the bad/undefined state.
29   void setBad();
30 
31   // Are we in the bad/undefined state?
isBad() const32   bool isBad() const { return (lambda_ < 0); }
33 
34   // The scale factor.  In the bad/undefined state, it is negative.
lambda() const35   double lambda() const { return lambda_; }
36 
37   // The probabilities of letters corresponding to matrix rows (1st index).
38   // In the bad/undefined state, it is a null pointer.
letterProbs1() const39   const double *letterProbs1() const {return isBad() ? 0 : &letterProbs1_[0];}
40 
41   // The probabilities of letters corresponding to matrix columns (2nd index).
42   // In the bad/undefined state, it is a null pointer.
letterProbs2() const43   const double *letterProbs2() const {return isBad() ? 0 : &letterProbs2_[0];}
44 
45  private:
46   double lambda_;
47   std::vector<double> letterProbs1_;
48   std::vector<double> letterProbs2_;
49 
50   bool find_ub(double **matrix, int alpha_size, double *ub);
51   bool binary_search(double** matrix, int alpha_size, double lb, double ub, std::vector<double>& letprob1, std::vector<double>& letprob2, double* lambda, int maxiter);
52   double calculate_lambda(double** matrix, int alpha_size, std::vector<double>& letprob1, std::vector<double>& letprob2, int maxiter, int max_boundary_search_iter, double lb_ratio);
53   bool check_lambda(double** matrix, double lambda, int alpha_size, std::vector<double>& letprob1, std::vector<double>& letprob2);
54 };
55 
56 }  // end namespace
57 
58 #endif
59