1 #ifndef ALGO_BLAST_GUMBEL_PARAMS__INCLUDED_NJN_LOCALMAXSTATUTIL 2 #define ALGO_BLAST_GUMBEL_PARAMS__INCLUDED_NJN_LOCALMAXSTATUTIL 3 4 /* $Id: njn_localmaxstatutil.hpp 183505 2010-02-18 16:10:58Z boratyng $ 5 * =========================================================================== 6 * 7 * PUBLIC DOMAIN NOTICE 8 * National Center for Biotechnology Information 9 * 10 * This software/database is a "United States Government Work" under the 11 * terms of the United States Copyright Act. It was written as part of 12 * the author's offical duties as a United States Government employee and 13 * thus cannot be copyrighted. This software/database is freely available 14 * to the public for use. The National Library of Medicine and the U.S. 15 * Government have not placed any restriction on its use or reproduction. 16 * 17 * Although all reasonable efforts have been taken to ensure the accuracy 18 * and reliability of the software and data, the NLM and the U.S. 19 * Government do not and cannot warrant the performance or results that 20 * may be obtained by using this software or data. The NLM and the U.S. 21 * Government disclaim all warranties, express or implied, including 22 * warranties of performance, merchantability or fitness for any particular 23 * purpose. 24 * 25 * Please cite the author in any work or product based on this material. 26 * 27 * ===========================================================================*/ 28 29 /***************************************************************************** 30 31 File name: njn_localmaxstatutil.hpp 32 33 Author: John Spouge 34 35 Contents: Random walk parameters 36 37 ******************************************************************************/ 38 39 #include <corelib/ncbistl.hpp> 40 #include <corelib/ncbitype.h> 41 #include <corelib/ncbi_limits.h> 42 43 #include "njn_matrix.hpp" 44 45 46 BEGIN_NCBI_SCOPE 47 BEGIN_SCOPE(blast) 48 49 BEGIN_SCOPE(Njn) 50 BEGIN_SCOPE(LocalMaxStatUtil) 51 52 53 const double REL_TOL = 1.0e-6; 54 55 void flatten ( // allocates memory for linear probabilities and scores 56 size_t dimension_, // dimension of equilProb_ 57 const Int4 *const *scoreMatrix_, // packed scoring matrix [0...dimension_)[0...dimension_) 58 const double *const *prob_, // prob_ [0...dimension_)[0...dimension_) : distribution of scores sum to 1.0 59 size_t *dim_, // dimension of p_ 60 Int4 **score_, // score [0...dim_) in increasing order 61 double **p_, // linear p_ [0...dim_) : distribution of scores 62 size_t dimension2_ = 0); // dimension2 of equilProb_ : defaults to dimension_ 63 // asserts (sum (p_) == 1.0); 64 65 double lambda ( 66 size_t dimMatrix_, // dimension of equilProb_ 67 const Int4 *const *scoreMatrix_, // packed scoring matrix [0...dimension_)[0...dimension_) 68 const double *q_); // q_ [0...dimension_) : distribution of independent letters 69 70 double mu ( 71 size_t dimension_, // #(distinct values) 72 const Int4 *score_, // scores in increasing order 73 const double *prob_); // probability of corresponding value 74 75 double lambda ( 76 size_t dimension_, // #(distinct values) 77 const Int4 *score_, // scores in increasing order 78 const double *prob_); // probability of corresponding value 79 // assumes logarithmic regime 80 81 double muAssoc ( 82 size_t dimension_, // #(distinct values) of scores & probabilities (which are paired) 83 const Int4 *score_, // scores in increasing order 84 const double *prob_, // corresponding probabilities 85 double lambda_ = 0.0); // lambda 86 87 double muPowerAssoc ( 88 size_t dimension_, // #(distinct values) of scores & probabilities (which are paired) 89 const Int4 *score_, // scores in increasing order 90 const double *prob_, // corresponding probabilities 91 double lambda_ = 0.0, // lambda 92 Int4 power_ = 1); // power 93 94 double thetaMin ( // minimizing value for r(theta) 95 size_t dimension_, // #(distinct values) 96 const Int4 *score_, // scores in increasing order 97 const double *prob_, // probability of corresponding value 98 double lambda_ = 0.0); // lambda 99 // assumes logarithmic regime 100 101 double rMin ( // minimum value of r(theta) 102 size_t dimension_, // #(distinct values) 103 const Int4 *score_, // scores in increasing order 104 const double *prob_, // probability of corresponding value 105 double lambda_ = 0.0, // lambda 106 double thetaMin_ = 0.0); // argument of rate 107 // assumes logarithmic regime 108 109 double r ( // r(theta) 110 size_t dimension_, // #(distinct values) 111 const Int4 *score_, // scores in increasing order 112 const double *prob_, // probability of corresponding value 113 double theta_); // argument of rate 114 // assumes logarithmic regime 115 116 Int4 delta ( // theta [minus delta] for ungapped sequence comparison 117 size_t dimension_, // #(distinct values) of scores & probabilities (which are paired) 118 const Int4 *score_); // scores 119 120 double thetaMinusDelta ( // theta [minus delta] for ungapped sequence comparison 121 double lambda_, // lambda, the exponential rate for the local maximum 122 size_t dimension_, // #(distinct values) of scores & probabilities (which are paired) 123 const Int4 *score_); // scores 124 125 void descendingLadderEpoch ( 126 size_t dimension_, // #(distinct values) 127 const Int4 *score_, // values 128 const double *prob_, // probability of corresponding value 129 double *eSumAlpha_ = 0, // expectation (sum [alpha]) 130 double *eOneMinusExpSumAlpha_ = 0, // expectation [1.0 - exp (sum [alpha])] 131 bool isStrict_ = false, // ? is this a strict descending ladder epoch 132 double lambda0_ = 0.0, // lambda for flattened distribution (avoid recomputation) 133 double mu0_ = 0.0, // mean of flattened distribution (avoid recomputation) 134 double muAssoc0_ = 0.0, // mean of associated flattened distribution (avoid recomputation) 135 double thetaMin0_ = 0.0, // thetaMin of flattened distribution (avoid recomputation) 136 double rMin0_ = 0.0, // rMin of flattened distribution (avoid recomputation) 137 double time_ = 0.0, // get time for the dynamic programming computation 138 bool *terminated_ = 0); // ? Was the dynamic programming computation terminated prematurely ? 139 140 void descendingLadderEpochRepeat ( 141 size_t dimension_, // #(distinct values) 142 const Int4 *score_, // values 143 const double *prob_, // probability of corresponding value 144 double *eSumAlpha_ = 0, // expectation (sum [alpha]) 145 double *eOneMinusExpSumAlpha_ = 0, // expectation [1.0 - exp (sum [alpha])] 146 bool isStrict_ = false, // ? is this a strict descending ladder epoch 147 double lambda_ = 0.0, // lambda for repeats : default is lambda0_ below 148 size_t endW_ = 0, // maximum w plus 1 149 double *pAlphaW_ = 0, // probability {alpha = w} : pAlphaW_ [0, wEnd) 150 double *eOneMinusExpSumAlphaW_ = 0, // expectation [1.0 - exp (sum [alpha]); alpha = w] : eOneMinusExpSumAlphaW_ [0, wEnd) 151 double lambda0_ = 0.0, // lambda for flattened distribution (avoid recomputation) 152 double mu0_ = 0.0, // mean of flattened distribution (avoid recomputation) 153 double muAssoc0_ = 0.0, // mean of associated flattened distribution (avoid recomputation) 154 double thetaMin0_ = 0.0, // thetaMin of flattened distribution (avoid recomputation) 155 double rMin0_ = 0.0, // rMin of flattened distribution (avoid recomputation) 156 double time_ = 0.0, // get time for the dynamic programming computation 157 bool *terminated_ = 0); // ? Was the dynamic programming computation terminated prematurely ? 158 // assumes logarithmic regime 159 160 bool isProbDist ( 161 size_t dimension_, // #(distinct values) of scores & probabilities (which are paired) 162 const double *prob_); // corresponding probabilities 163 164 bool isScoreIncreasing ( 165 size_t dimension_, // #(distinct values) 166 const Int4 *score_); // scores in increasing order 167 168 bool isLogarithmic ( 169 size_t dimension_, // #(distinct values) 170 const Int4 *score_, // scores in increasing order 171 const double *prob_); // probability of corresponding value 172 173 END_SCOPE(LocalMaxStatUtil) 174 END_SCOPE(Njn) 175 176 END_SCOPE(blast) 177 END_NCBI_SCOPE 178 179 #endif //! ALGO_BLAST_GUMBEL_PARAMS__INCLUDED_NJN_LOCALMAXSTATUTIL 180