1 /* $Id: $
2 * ===========================================================================
3 *
4 * PUBLIC DOMAIN NOTICE
5 * National Center for Biotechnology Information
6 *
7 * This software/database is a "United States Government Work" under the
8 * terms of the United States Copyright Act. It was written as part of
9 * the author's offical duties as a United States Government employee and
10 * thus cannot be copyrighted. This software/database is freely available
11 * to the public for use. The National Library of Medicine and the U.S.
12 * Government have not placed any restriction on its use or reproduction.
13 *
14 * Although all reasonable efforts have been taken to ensure the accuracy
15 * and reliability of the software and data, the NLM and the U.S.
16 * Government do not and cannot warrant the performance or results that
17 * may be obtained by using this software or data. The NLM and the U.S.
18 * Government disclaim all warranties, express or implied, including
19 * warranties of performance, merchantability or fitness for any particular
20 * purpose.
21 *
22 * Please cite the author in any work or product based on this material.
23 *
24 * ===========================================================================*/
25
26 /*****************************************************************************
27
28 File name: njn_localmaxstat.cpp
29
30 Author: John Spouge
31
32 Contents:
33
34 ******************************************************************************/
35
36 #include "sls_basic.hpp"
37 #include "njn_localmaxstat.hpp"
38 #include "njn_memutil.hpp"
39 #include "njn_dynprogproblim.hpp"
40 #include "njn_function.hpp"
41 #include "njn_integer.hpp"
42 #include "njn_localmaxstatutil.hpp"
43
44 using namespace Njn;
45
46 double LocalMaxStat::s_time = 0.0;
47
init(size_t dimension_)48 void LocalMaxStat::init (size_t dimension_)
49 {
50 if (dimension_ > 0)
51 {
52 d_score_p = new long int [dimension_];
53 d_prob_p = new double [dimension_];
54 }
55
56 d_dimension = dimension_;
57 }
58
free2()59 void LocalMaxStat::free2 ()
60 {
61 if (getDimension () > 0)
62 {
63 delete [] d_score_p; d_score_p = 0;
64 delete [] d_prob_p; d_prob_p = 0;
65 }
66
67 d_dimension = 0;
68 }
69
clear()70 void LocalMaxStat::clear ()
71 {
72 free2 ();
73 init (0);
74
75 d_lambda = 0.0;
76 d_k = 0.0;
77 d_c = 0.0;
78 d_thetaMin = 0.0;
79 d_rMin = 0.0;
80 d_delta = 0;
81 d_thetaMinusDelta = 0.0;
82 d_mu = 0.0;
83 d_sigma = 0.0;
84 d_muAssoc = 0.0;
85 d_sigmaAssoc = 0.0;
86 d_meanWDLE = 0.0;
87 d_terminated = false;
88 }
89
copy(size_t dimension_,const long int * score_,const double * prob_,double lambda_,double k_,double c_,double thetaMin_,double rMin_,long int delta_,double thetaMinusDelta_,double mu_,double sigma_,double muAssoc_,double sigmaAssoc_,double meanLength_,bool terminated_)90 void LocalMaxStat::copy (
91 size_t dimension_, // #(distinct values) of scores & probabilities (which are paired)
92 const long int *score_, // scores
93 const double *prob_, // probabilities
94 double lambda_, // lambda for associated random walk
95 double k_, // k for random walk : exponential prefactor
96 double c_, // c for random walk : exponential prefactor (global alignment)
97 double thetaMin_, // theta for minimum expectation (exp (theta * score))
98 double rMin_, // minimum expectation (exp (theta * score))
99 long int delta_, // span
100 double thetaMinusDelta_, // renewal span parameter
101 double mu_, // n_step mean for random walk
102 double sigma_, // n_step standard deviation for random walk
103 double muAssoc_, // n_step mean for associated random walk (relative entropy)
104 double sigmaAssoc_, // n_step standard deviation for associated random walk
105 double meanLength_, // expected renewal length
106 bool terminated_) // ? Was the dynamic programming computation terminated prematurely ?
107 {
108 free2 ();
109 init (dimension_);
110
111 memcpy (d_score_p, score_, sizeof (long int) * getDimension ());
112 memcpy (d_prob_p, prob_, sizeof (double) * getDimension ());
113
114 d_lambda = lambda_;
115 d_k = k_;
116 d_c = c_;
117 d_thetaMin = thetaMin_;
118 d_rMin = rMin_;
119 d_delta = delta_;
120 d_thetaMinusDelta = thetaMinusDelta_;
121 d_mu = mu_;
122 d_sigma = sigma_;
123 d_muAssoc = muAssoc_;
124 d_sigmaAssoc = sigmaAssoc_;
125 d_meanWDLE = meanLength_;
126 d_terminated = terminated_;
127 }
128
copy(size_t dimension_,const long int * score_,const double * prob_)129 void LocalMaxStat::copy (
130 size_t dimension_, // #(distinct values) of scores & probabilities (which are paired)
131 const long int *score_, // scores in increasing order
132 const double *prob_) // corresponding probabilities
133 {
134 if (dimension_ == 0)
135 {
136 clear ();
137 return;
138 }
139
140 if (! LocalMaxStatUtil::isLogarithmic (dimension_, score_, prob_))
141 {
142 //IoUtil::abort ("LocalMaxStat::copy : ! isLogarithmic");
143 throw Sls::error("Error - you have exceeded the calculation time or memory limit.\nThe error might indicate that the regime is linear or too close to linear to permit efficient computation.\nPossible solutions include changing the randomization seed, or increasing the allowed calculation time and the memory limit.\n",3);
144 }
145
146 size_t i = 0;
147 /*sls deleted size_t j = 0;*/
148 /*sls deleted long int iter = 0;*/
149 /*sls deleted long int value = 0;*/
150
151 free2 ();
152 init (dimension_);
153
154 memcpy (d_score_p, score_, sizeof (long int) * getDimension ());
155 memcpy (d_prob_p, prob_, sizeof (double) * getDimension ());
156
157 d_mu = LocalMaxStatUtil::mu (getDimension (), getScore (), getProb ());
158 d_sigma = 0.0;
159
160 for (i = 0; i < dimension_; i++)
161 {
162 d_sigma += static_cast <double> (score_ [i]) * static_cast <double> (score_ [i]) * prob_ [i];
163 }
164
165 d_sigma -= getMu () * getMu ();
166 d_sigma = Function::psqrt (getSigma ());
167
168 // calculate lambda
169
170 d_lambda = LocalMaxStatUtil::lambda (getDimension (), getScore (), getProb ());
171 d_muAssoc = LocalMaxStatUtil::muAssoc (getDimension (), getScore (), getProb (), getLambda ());
172 d_sigmaAssoc = 0.0;
173
174 for (i = 0; i < getDimension (); i++)
175 {
176 d_sigmaAssoc += static_cast <double> (getScore () [i]) * static_cast <double> (getScore () [i]) *
177 getProb () [i] * exp (getLambda () * static_cast <double> (getScore () [i]));
178 }
179
180 d_sigmaAssoc -= getMuAssoc () * getMuAssoc ();
181 d_sigmaAssoc = Function::psqrt (d_sigmaAssoc);
182
183 d_thetaMin = LocalMaxStatUtil::thetaMin (getDimension (), getScore (), getProb (), getLambda ());
184 d_rMin = LocalMaxStatUtil::rMin (getDimension (), getScore (), getProb (), getLambda (), getThetaMin ());
185
186 d_delta = LocalMaxStatUtil::delta (getDimension (), getScore ());
187 d_thetaMinusDelta = LocalMaxStatUtil::thetaMinusDelta (getLambda (), getDimension (), getScore ());
188
189 dynProgCalc ();
190 }
191
getR(double theta_) const192 double LocalMaxStat::getR (double theta_) const
193 {
194 return LocalMaxStatUtil::r (d_dimension, d_score_p, d_prob_p, theta_);
195 }
196
dynProgCalc()197 void LocalMaxStat::dynProgCalc ()
198 // k for random walk : exponential prefactor
199 // expected renewal length for weak ladder epochs
200 {
201 double eSumAlpha_ = 0.0;
202 double eOneMinusExpSumAlpha_ = 0.0;
203 LocalMaxStatUtil::descendingLadderEpoch (getDimension (), getScore (), getProb (),
204 &eSumAlpha_, &eOneMinusExpSumAlpha_, false,
205 getLambda (), getMu (), getMuAssoc (), getThetaMin (), getRMin (), getTime (), &d_terminated);
206
207 if (getTerminated ()) return;
208
209 // fluctuation sum quantities
210 double ratio = eOneMinusExpSumAlpha_ / eSumAlpha_;
211 d_k = getMu () * getMu () / getThetaMinusDelta () / getMuAssoc () * ratio * ratio;
212 d_meanWDLE = eSumAlpha_ / getMu ();
213 d_c = getK () * getMeanWDLE () / eOneMinusExpSumAlpha_;
214 }
215
216