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