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