1 #ifndef INCLUDED_SLS_PVALUES 2 #define INCLUDED_SLS_PVALUES 3 4 /* $Id: $ 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: sls_pvalues.hpp 32 33 Author: Sergey Sheetlin 34 35 Contents: P-values calculation routines 36 37 ******************************************************************************/ 38 39 #include "sls_basic.hpp" 40 41 #include <vector> 42 #include <string> 43 #include <math.h> 44 #include <cstdlib> 45 46 #include <cmath> 47 #include <iostream> 48 49 namespace Sls { 50 51 struct ALP_set_of_parameters 52 { ALP_set_of_parametersSls::ALP_set_of_parameters53 ALP_set_of_parameters() 54 { 55 d_params_flag=false; 56 57 b_I=0; 58 b_I_error=0; 59 60 b_J=0; 61 b_J_error=0; 62 63 beta_I=0; 64 beta_I_error=0; 65 66 beta_J=0; 67 beta_J_error=0; 68 69 tau=0; 70 tau_error=0; 71 72 }; 73 74 double lambda; 75 double lambda_error; 76 77 double C; 78 double C_error; 79 80 81 double K; 82 double K_error; 83 84 double a_I; 85 double a_I_error; 86 87 double a_J; 88 double a_J_error; 89 90 double sigma; 91 double sigma_error; 92 93 double alpha_I; 94 double alpha_I_error; 95 96 double alpha_J; 97 double alpha_J_error; 98 99 double a; 100 double a_error; 101 102 double alpha; 103 double alpha_error; 104 105 double gapless_a; 106 double gapless_a_error; 107 108 double gapless_alpha; 109 double gapless_alpha_error; 110 111 long int G; 112 long int G1; 113 long int G2; 114 115 std::vector<double > m_LambdaSbs; 116 std::vector<double > m_KSbs; 117 std::vector<double > m_CSbs; 118 119 std::vector<double > m_SigmaSbs; 120 121 std::vector<double > m_AlphaISbs; 122 std::vector<double > m_AlphaJSbs; 123 124 std::vector<double > m_AISbs; 125 std::vector<double > m_AJSbs; 126 127 double m_CalcTime; 128 129 bool d_params_flag;//if true, then the parameters are defined and P-values can be calculated 130 131 //intercepts 132 double b_I; 133 double b_I_error; 134 135 double b_J; 136 double b_J_error; 137 138 double beta_I; 139 double beta_I_error; 140 141 double beta_J; 142 double beta_J_error; 143 144 double tau; 145 double tau_error; 146 147 std::vector<double > m_BISbs; 148 std::vector<double > m_BJSbs; 149 150 std::vector<double > m_BetaISbs; 151 std::vector<double > m_BetaJSbs; 152 153 std::vector<double > m_TauSbs; 154 155 //tmp values 156 double vi_y_thr; 157 double vj_y_thr; 158 double c_y_thr; 159 }; 160 161 std::ostream &operator<<(std::ostream &s_, 162 const ALP_set_of_parameters &gumbel_params_); 163 164 std::istream &operator>>(std::istream &s_, 165 ALP_set_of_parameters &gumbel_params_); 166 167 class pvalues{ 168 169 public: 170 171 172 pvalues(); 173 174 ~pvalues(); 175 176 177 public: 178 179 static void get_appr_tail_prob_with_cov( 180 const ALP_set_of_parameters &par_, 181 bool blast_, 182 double y_, 183 double m_, 184 double n_, 185 186 double &P_, 187 double &P_error_, 188 189 double &E_, 190 double &E_error_, 191 192 double &area_, 193 bool &area_is_1_flag_); 194 195 static void compute_tmp_values(ALP_set_of_parameters &par_); 196 197 static void get_appr_tail_prob_with_cov_without_errors( 198 const ALP_set_of_parameters &par_, 199 bool blast_, 200 double y_, 201 double m_, 202 double n_, 203 204 double &P_, 205 206 double &E_, 207 208 double &area_, 209 bool &area_is_1_flag_, 210 bool compute_only_area_=false); 211 212 static void get_P_error_using_splitting_method( 213 const ALP_set_of_parameters &par_, 214 bool blast_, 215 double y_, 216 double m_, 217 double n_, 218 219 double &P_, 220 double &P_error_, 221 222 double &E_, 223 double &E_error_, 224 225 bool &area_is_1_flag_); 226 227 228 public: 229 230 static void compute_intercepts( 231 ALP_set_of_parameters &par_); 232 233 void calculate_P_values( 234 long int Score1, 235 long int Score2, 236 double Seq1Len, 237 double Seq2Len, 238 const ALP_set_of_parameters &ParametersSet, 239 std::vector<double> &P_values, 240 std::vector<double> &P_values_errors, 241 std::vector<double> &E_values, 242 std::vector<double> &E_values_errors); 243 244 void calculate_P_values( 245 double Score, 246 double Seq1Len, 247 double Seq2Len, 248 const ALP_set_of_parameters &ParametersSet, 249 double &P_value, 250 double &P_value_error, 251 double &E_value, 252 double &E_value_error, 253 bool read_Sbs_par_flag=true); 254 ran3()255 static inline double ran3()//generates the next random value 256 { 257 double rand_C=(double)((double)rand()/(double)RAND_MAX); 258 return rand_C; 259 }; 260 standard_normal()261 static inline double standard_normal()//generates standard normal random value using the Box�Muller transform 262 { 263 double r1=0; 264 while(r1==0) 265 { 266 r1=ran3(); 267 }; 268 double r2=0; 269 while(r2==0) 270 { 271 r2=ran3(); 272 }; 273 274 double v1=-2*log(r1); 275 if(v1<0) 276 { 277 v1=0; 278 }; 279 return sqrt(v1)*cos(2*pi*r2); 280 }; 281 282 static //returns "true" if the Gumbel parameters are properly defined and "false" otherwise 283 bool assert_Gumbel_parameters( 284 const ALP_set_of_parameters &par_);//a set of Gumbel parameters 285 286 287 288 289 290 public: 291 292 293 bool blast; 294 double eps; 295 double a_normal; 296 double b_normal; 297 long int N_normal; 298 double h_normal; 299 double *p_normal; 300 301 302 }; 303 } 304 305 #endif //! INCLUDED 306 307