1 #ifndef ALGO_BLAST_GUMBEL_PARAMS__INCLUDED_NJN_LOCALMAXSTAT 2 #define ALGO_BLAST_GUMBEL_PARAMS__INCLUDED_NJN_LOCALMAXSTAT 3 4 /* $Id: njn_localmaxstat.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_localmaxstat.hpp 32 33 Author: John Spouge 34 35 Contents: Random walk parameters 36 37 ******************************************************************************/ 38 39 #include <corelib/ncbitype.h> 40 #include <corelib/ncbistl.hpp> 41 #include <assert.h> 42 #include <math.h> 43 #include <ostream> 44 #include <vector> 45 46 BEGIN_NCBI_SCOPE 47 BEGIN_SCOPE(blast) 48 49 BEGIN_SCOPE(Njn) 50 51 52 class LocalMaxStat { 53 54 // calculates the statistical parameters for the local maximum in a random walk 55 // 56 // The scores are uniqued and 57 // with the correspondence to probabilities maintained, placed in ascending order. 58 59 public: 60 61 // The following subroutines control the time for the dynamic programming computation. 62 // The default time_ = 0.0 permits the computation to run forever. setTime(double time_=0.0)63 static void setTime (double time_ = 0.0) {assert (time_ >= 0.0); s_time = time_;} // set time for the dynamic programming computation getTime()64 static double getTime () {return s_time;} // get time for the dynamic programming computation 65 // For an object o, 66 // if the computation is terminated before it finishes, 67 // o.getTerminated () == true. 68 LocalMaxStat(size_t dimension_=0,const Int4 * score_=0,const double * prob_=0)69 inline LocalMaxStat ( 70 size_t dimension_ = 0, // #(distinct values) 71 const Int4 *score_ = 0, // scores in increasing order 72 const double *prob_ = 0) // probability of corresponding value 73 : d_dimension (0), d_score_p (0), d_prob_p (0), 74 d_lambda (0.0), d_k (0.0), d_c (0.0), d_thetaMin (0.0), d_rMin (0.0), 75 d_delta (0), d_thetaMinusDelta (0.0), 76 d_mu (0.0), d_sigma (0.0), d_muAssoc (0.0), d_sigmaAssoc (0.0), 77 d_meanWDLE (0.0), d_terminated (false) 78 { 79 copy (dimension_, score_, prob_); 80 } 81 LocalMaxStat(const LocalMaxStat & localMaxStat_)82 inline LocalMaxStat (const LocalMaxStat &localMaxStat_) // random walk parameters 83 : d_dimension (0), d_score_p (0), d_prob_p (0), 84 d_lambda (0.0), d_k (0.0), d_c (0.0), d_thetaMin (0.0), d_rMin (0.0), 85 d_delta (0), d_thetaMinusDelta (0.0), 86 d_mu (0.0), d_sigma (0.0), d_muAssoc (0.0), d_sigmaAssoc (0.0), 87 d_meanWDLE (0.0), d_terminated (false) 88 { 89 copy (localMaxStat_); 90 } 91 ~LocalMaxStat()92 inline ~LocalMaxStat () {free ();} 93 operator bool() const94 inline operator bool () // ? is the object ready for computation ? 95 const { 96 return d_dimension != 0; 97 } 98 operator =(const LocalMaxStat & localMaxStat_)99 inline LocalMaxStat &operator= (const LocalMaxStat &localMaxStat_) // random walk parameters 100 { 101 if (this != &localMaxStat_) copy (localMaxStat_); 102 return *this; 103 } 104 105 void copy ( 106 size_t dimension_, // #(distinct values) of scores & probabilities (which are paired) 107 const Int4 *score_, // scores in increasing order 108 const double *prob_); // probabilities 109 copy(const LocalMaxStat & localMaxStat_)110 inline void copy (const LocalMaxStat &localMaxStat_) 111 { 112 copy (localMaxStat_.getDimension (), localMaxStat_.getScore (), localMaxStat_.getProb (), 113 localMaxStat_.getLambda (), localMaxStat_.getK (), localMaxStat_.getC (), 114 localMaxStat_.getThetaMin (), localMaxStat_.getRMin (), 115 localMaxStat_.getDelta (), localMaxStat_.getThetaMinusDelta (), 116 localMaxStat_.getMu (), localMaxStat_.getSigma (), localMaxStat_.getMuAssoc (), localMaxStat_.getSigmaAssoc (), 117 localMaxStat_.getMeanWDLE (), localMaxStat_.getTerminated ()); 118 } 119 120 void copy ( 121 size_t dimension_, // #(distinct values) of scores & probabilities (which are paired) 122 const Int4 *score_, // scores in increasing order 123 const double *prob_, // probabilities 124 double lambda_, // lambda for associated random walk 125 double k_, // k for random walk : exponential prefactor 126 double c_, // c for random walk : exponential prefactor (global alignment) 127 double thetaMin_, // theta for minimum expectation (exp (theta * score)) 128 double rMin_, // minimum expectation (exp (theta * score)) 129 Int4 delta_, // span 130 double thetaMinusDelta_, // renewal span parameter 131 double mu_, // step mean for random walk 132 double sigma_, // step standard deviation for random walk 133 double muAssoc_, // step mean for associated random walk (relative entropy) 134 double sigmaAssoc_, // step standard deviation for associated random walk 135 double meanDLE_, // expected renewal length 136 bool terminated_ = false); // ? Was the dynamic programming computation terminated prematurely ? 137 out(std::ostream & ostr_) const138 inline std::ostream &out (std::ostream &ostr_) const {return ostr_;} // output 139 140 double getR (double theta_) const; // r (theta_) : dominant eigenvalue for theta_ 141 getA() const142 double getA () const {return getMuAssoc () == 0 ? HUGE_VAL : 1.0 / getMuAssoc ();} // expected [length / y] for achieving y 143 getAlpha() const144 double getAlpha () const {return getSigmaAssoc () * getSigmaAssoc () * getA () * getA () * getA ();} // var [length] / y for achieving y 145 getDimension() const146 inline size_t getDimension () const {return d_dimension;} // #(distinct values) of scores & probabilities (which are paired) getScore() const147 inline const Int4 *getScore () const {return d_score_p;} // scores in increasing order getProb() const148 inline const double *getProb () const {return d_prob_p;} // probabilities getLambda() const149 inline double getLambda () const {return d_lambda;} // lambda for associated random walk getK() const150 inline double getK () const {return d_k;} // k for random walk : exponential prefactor getC() const151 inline double getC () const {return d_c;} // c for random walk : exponential prefactor (global alignment) getThetaMin() const152 inline double getThetaMin () const {return d_thetaMin;} // theta for minimum expectation (exp (theta * score)) getRMin() const153 inline double getRMin () const {return d_rMin;} // minimum expectation (exp (theta * score)) getDelta() const154 inline Int4 getDelta () const {return d_delta;} // span getThetaMinusDelta() const155 inline double getThetaMinusDelta () const {return d_thetaMinusDelta;} // renewal span parameter getMu() const156 inline double getMu () const {return d_mu;} // step mean for random walk getSigma() const157 inline double getSigma () const {return d_sigma;} // step standard deviation for random walk getMuAssoc() const158 inline double getMuAssoc () const {return d_muAssoc;} // step mean for associated random walk (relative entropy) getSigmaAssoc() const159 inline double getSigmaAssoc () const {return d_sigmaAssoc;} // step standard deviation for associated random walk getMeanWDLE() const160 inline double getMeanWDLE () const {return d_meanWDLE;} // expected renewal length for weak ladder epochs getTerminated() const161 inline bool getTerminated () const {return d_terminated;} // ? Was the dynamic programming computation terminated prematurely ? 162 163 private: 164 165 // random walk distribution 166 size_t d_dimension; // #(distinct values) of scores & probabilities (which are paired) 167 Int4 *d_score_p; // scores in increasing order 168 double *d_prob_p; // probabilities 169 170 // Karlin-Altschul parameters 171 double d_lambda; // lambda for associated random walk 172 double d_k; // k for random walk : exponential prefactor 173 double d_c; // c for random walk : exponential prefactor (global alignment) 174 double d_thetaMin; // theta for minimum expectation (exp (theta * score)) 175 double d_rMin; // minimum expectation (exp (theta * score)) 176 Int4 d_delta; // span 177 double d_thetaMinusDelta; // renewal span parameter 178 double d_mu; // step mean for random walk 179 double d_sigma; // step standard deviation for random walk 180 double d_muAssoc; // step mean for associated random walk (relative entropy) 181 double d_sigmaAssoc; // step standard deviation for associated random walk 182 double d_meanWDLE; // expected renewal length for weak ladder epochs 183 bool d_terminated; // ? Was the dynamic programming computation terminated prematurely ? 184 185 void init (size_t dimension_); 186 void free (); 187 188 void clear (); 189 190 void dynProgCalc (); 191 // k for random walk : exponential prefactor 192 // expected renewal length for weak ladder epochs 193 194 static double s_time; 195 196 }; 197 198 END_SCOPE(Njn) 199 200 END_SCOPE(blast) 201 END_NCBI_SCOPE 202 203 #endif //! ALGO_BLAST_GUMBEL_PARAMS__INCLUDED_NJN_LOCALMAXSTAT 204