1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
|
#ifndef INCLUDED_NJN_LOCALMAXSTAT
#define INCLUDED_NJN_LOCALMAXSTAT
/* $Id: $
* ===========================================================================
*
* PUBLIC DOMAIN NOTICE
* National Center for Biotechnology Information
*
* This software/database is a "United States Government Work" under the
* terms of the United States Copyright Act. It was written as part of
* the author's offical duties as a United States Government employee and
* thus cannot be copyrighted. This software/database is freely available
* to the public for use. The National Library of Medicine and the U.S.
* Government have not placed any restriction on its use or reproduction.
*
* Although all reasonable efforts have been taken to ensure the accuracy
* and reliability of the software and data, the NLM and the U.S.
* Government do not and cannot warrant the performance or results that
* may be obtained by using this software or data. The NLM and the U.S.
* Government disclaim all warranties, express or implied, including
* warranties of performance, merchantability or fitness for any particular
* purpose.
*
* Please cite the author in any work or product based on this material.
*
* ===========================================================================*/
/*****************************************************************************
File name: njn_localmaxstat.hpp
Author: John Spouge
Contents: Random walk parameters
******************************************************************************/
#include <math.h>
#include <vector>
#include <assert.h>
#include <ostream>
namespace Njn {
class LocalMaxStat {
// calculates the statistical parameters for the local maximum in a random walk
//
// The scores are uniqued and
// with the correspondence to probabilities maintained, placed in ascending order.
public:
// The following subroutines control the time for the dynamic programming computation.
// The default time_ = 0.0 permits the computation to run forever.
static void setTime (double time_ = 0.0) {assert (time_ >= 0.0); s_time = time_;} // set time for the dynamic programming computation
static double getTime () {return s_time;} // get time for the dynamic programming computation
// For an object o,
// if the computation is terminated before it finishes,
// o.getTerminated () == true.
inline LocalMaxStat (
size_t dimension_ = 0, // #(distinct values)
const long int *score_ = 0, // scores in increasing order
const double *prob_ = 0) // probability of corresponding value
: d_dimension (0), d_score_p (0), d_prob_p (0),
d_lambda (0.0), d_k (0.0), d_c (0.0), d_thetaMin (0.0), d_rMin (0.0),
d_delta (0), d_thetaMinusDelta (0.0),
d_mu (0.0), d_sigma (0.0), d_muAssoc (0.0), d_sigmaAssoc (0.0),
d_meanWDLE (0.0), d_terminated (false)
{
copy (dimension_, score_, prob_);
}
inline LocalMaxStat (const LocalMaxStat &localMaxStat_) // random walk parameters
: d_dimension (0), d_score_p (0), d_prob_p (0),
d_lambda (0.0), d_k (0.0), d_c (0.0), d_thetaMin (0.0), d_rMin (0.0),
d_delta (0), d_thetaMinusDelta (0.0),
d_mu (0.0), d_sigma (0.0), d_muAssoc (0.0), d_sigmaAssoc (0.0),
d_meanWDLE (0.0), d_terminated (false)
{
copy (localMaxStat_);
}
inline ~LocalMaxStat () {free2 ();}
inline operator bool () // ? is the object ready for computation ?
const {
return d_dimension != 0;
}
inline LocalMaxStat &operator= (const LocalMaxStat &localMaxStat_) // random walk parameters
{
if (this != &localMaxStat_) copy (localMaxStat_);
return *this;
}
void copy (
size_t dimension_, // #(distinct values) of scores & probabilities (which are paired)
const long int *score_, // scores in increasing order
const double *prob_); // probabilities
inline void copy (const LocalMaxStat &localMaxStat_)
{
copy (localMaxStat_.getDimension (), localMaxStat_.getScore (), localMaxStat_.getProb (),
localMaxStat_.getLambda (), localMaxStat_.getK (), localMaxStat_.getC (),
localMaxStat_.getThetaMin (), localMaxStat_.getRMin (),
localMaxStat_.getDelta (), localMaxStat_.getThetaMinusDelta (),
localMaxStat_.getMu (), localMaxStat_.getSigma (), localMaxStat_.getMuAssoc (), localMaxStat_.getSigmaAssoc (),
localMaxStat_.getMeanWDLE (), localMaxStat_.getTerminated ());
}
void copy (
size_t dimension_, // #(distinct values) of scores & probabilities (which are paired)
const long int *score_, // scores in increasing order
const double *prob_, // probabilities
double lambda_, // lambda for associated random walk
double k_, // k for random walk : exponential prefactor
double c_, // c for random walk : exponential prefactor (global alignment)
double thetaMin_, // theta for minimum expectation (exp (theta * score))
double rMin_, // minimum expectation (exp (theta * score))
long int delta_, // span
double thetaMinusDelta_, // renewal span parameter
double mu_, // step mean for random walk
double sigma_, // step standard deviation for random walk
double muAssoc_, // step mean for associated random walk (relative entropy)
double sigmaAssoc_, // step standard deviation for associated random walk
double meanDLE_, // expected renewal length
bool terminated_ = false); // ? Was the dynamic programming computation terminated prematurely ?
inline std::ostream &out (std::ostream &ostr_) const {return ostr_;} // output
double getR (double theta_) const; // r (theta_) : dominant eigenvalue for theta_
double getA () const {return getMuAssoc () == 0 ? HUGE_VAL : 1.0 / getMuAssoc ();} // expected [length / y] for achieving y
double getAlpha () const {return getSigmaAssoc () * getSigmaAssoc () * getA () * getA () * getA ();} // var [length] / y for achieving y
inline size_t getDimension () const {return d_dimension;} // #(distinct values) of scores & probabilities (which are paired)
inline const long int *getScore () const {return d_score_p;} // scores in increasing order
inline const double *getProb () const {return d_prob_p;} // probabilities
inline double getLambda () const {return d_lambda;} // lambda for associated random walk
inline double getK () const {return d_k;} // k for random walk : exponential prefactor
inline double getC () const {return d_c;} // c for random walk : exponential prefactor (global alignment)
inline double getThetaMin () const {return d_thetaMin;} // theta for minimum expectation (exp (theta * score))
inline double getRMin () const {return d_rMin;} // minimum expectation (exp (theta * score))
inline long int getDelta () const {return d_delta;} // span
inline double getThetaMinusDelta () const {return d_thetaMinusDelta;} // renewal span parameter
inline double getMu () const {return d_mu;} // step mean for random walk
inline double getSigma () const {return d_sigma;} // step standard deviation for random walk
inline double getMuAssoc () const {return d_muAssoc;} // step mean for associated random walk (relative entropy)
inline double getSigmaAssoc () const {return d_sigmaAssoc;} // step standard deviation for associated random walk
inline double getMeanWDLE () const {return d_meanWDLE;} // expected renewal length for weak ladder epochs
inline bool getTerminated () const {return d_terminated;} // ? Was the dynamic programming computation terminated prematurely ?
private:
// random walk distribution
size_t d_dimension; // #(distinct values) of scores & probabilities (which are paired)
long int *d_score_p; // scores in increasing order
double *d_prob_p; // probabilities
// Karlin-Altschul parameters
double d_lambda; // lambda for associated random walk
double d_k; // k for random walk : exponential prefactor
double d_c; // c for random walk : exponential prefactor (global alignment)
double d_thetaMin; // theta for minimum expectation (exp (theta * score))
double d_rMin; // minimum expectation (exp (theta * score))
long int d_delta; // span
double d_thetaMinusDelta; // renewal span parameter
double d_mu; // step mean for random walk
double d_sigma; // step standard deviation for random walk
double d_muAssoc; // step mean for associated random walk (relative entropy)
double d_sigmaAssoc; // step standard deviation for associated random walk
double d_meanWDLE; // expected renewal length for weak ladder epochs
bool d_terminated; // ? Was the dynamic programming computation terminated prematurely ?
void init (size_t dimension_);
void free2 ();
void clear ();
void dynProgCalc ();
// k for random walk : exponential prefactor
// expected renewal length for weak ladder epochs
static double s_time;
};
}
#endif
|
