openModeller
Version 1.4.0
|
00001 00027 #include <openmodeller/pre/PreChiSquare.hh> 00028 #include <openmodeller/Sampler.hh> 00029 #include <openmodeller/Log.hh> 00030 #include <openmodeller/Exceptions.hh> 00031 00032 #include <stdio.h> 00033 #include <math.h> 00034 00035 PreChiSquare::PreChiSquare() //constructor 00036 { 00037 } 00038 00039 PreChiSquare::~PreChiSquare() //destructor 00040 { 00041 } 00042 00043 bool 00044 PreChiSquare::checkParameters( const PreParameters& parameters ) const 00045 { 00046 SamplerPtr samplerPtr; 00047 if( ! parameters.retrieve( "Sampler", samplerPtr ) ) 00048 { 00049 Log::instance()->error( "Missing parameter: Sampler. \n" ); 00050 return false; 00051 } 00052 00053 if ( ! samplerPtr->getEnvironment() ) 00054 { 00055 std::string msg = "Sampler has no environment.\n"; 00056 00057 Log::instance()->error( msg.c_str() ); 00058 00059 throw InvalidParameterException( msg ); 00060 } 00061 return true; 00062 } 00063 00064 bool 00065 PreChiSquare::runImplementation() 00066 { 00067 size_t layer1, layer2; 00068 00069 init(); 00070 00071 for( layer1 = 0; layer1 < num_layers; layer1++ ) 00072 { 00073 statistic1.push_back(0); 00074 statistic2.push_back(0); 00075 } 00076 00077 for( layer1 = 0; layer1 < num_layers; layer1++ ) 00078 for(layer2 = layer1+1; layer2 < num_layers; layer2++) 00079 { 00080 setMeasured(layer1, layer2); 00081 setExpected(); 00082 setChicell(); 00083 setStatistic(layer1, layer2); 00084 } 00085 00086 SamplerPtr samplerPtr; 00087 params_.retrieve( "Sampler", samplerPtr ); 00088 00089 size_t aux; 00090 for( unsigned int ind = 0 ; ind < num_layers ; ind++ ) 00091 { 00092 string layer_id = samplerPtr->getEnvironment()->getLayerPath(ind); 00093 00094 PreParameters result; 00095 00096 aux = statistic1[ind] + statistic2[ind]; 00097 00098 result.store( "number of correlated layers", (int)aux ); 00099 00100 result_by_layer_[layer_id] = result; 00101 } 00102 00103 return true; 00104 } 00105 00106 void 00107 PreChiSquare::getAcceptedParameters( stringMap& info) 00108 { 00109 info["Sampler"] = "samplerPtr"; 00110 } 00111 00112 void 00113 PreChiSquare::getLayersetResultSpec ( stringMap& info) 00114 { 00115 00116 } 00117 00118 void 00119 PreChiSquare::getLayerResultSpec ( stringMap& info) 00120 { 00121 info["number of correlated layers"] = "int"; 00122 } 00123 00124 bool 00125 PreChiSquare::init() 00126 { 00127 SamplerPtr samplerPtr; 00128 params_.retrieve( "Sampler", samplerPtr ); 00129 00130 setNlayers(samplerPtr); 00131 00132 my_presences = samplerPtr->getPresences(); 00133 00134 setNpoints(); 00135 setNclass(); 00136 setMinimum(); 00137 setDelta(); 00138 00139 if (nclass > classLimit) 00140 { 00141 Log::instance()->error( "ChiSquare: measured, expected, chicell: number of class > %d\n", classLimit ); 00142 return false; 00143 } 00144 00145 return true; 00146 } 00147 00148 size_t 00149 PreChiSquare::getNpoints() 00150 { 00151 return num_points; 00152 } 00153 00154 void 00155 PreChiSquare::setNpoints() 00156 { 00157 num_points = my_presences->numOccurrences(); 00158 } 00159 00160 size_t 00161 PreChiSquare::getNlayers() 00162 { 00163 return num_layers; 00164 } 00165 00166 void 00167 PreChiSquare::setNlayers(SamplerPtr samplerPtr) 00168 { 00169 num_layers = samplerPtr->numIndependent(); 00170 00171 if ( num_layers < 2 ) 00172 { 00173 std::string msg = "chisquare needs at least 2 layers.\n"; 00174 00175 Log::instance()->error( msg.c_str() ); 00176 00177 throw InvalidParameterException( msg ); 00178 } 00179 } 00180 00181 size_t 00182 PreChiSquare::getNclass() 00183 { 00184 return nclass; 00185 } 00186 00187 void 00188 PreChiSquare::setNclass() 00189 { 00190 nclass = (size_t)floor(sqrt((double)(num_points))); 00191 } 00192 00193 void 00194 PreChiSquare::setMinimum() 00195 { 00196 OccurrencesImpl::iterator it = my_presences->begin(); 00197 OccurrencesImpl::iterator last = my_presences->end(); 00198 00199 Sample sample = (*it)->environment(); 00200 minimum = sample; 00201 ++it; 00202 while ( it != last ) 00203 { 00204 Sample const& sample = (*it)->environment(); 00205 minimum &= sample; 00206 ++it; 00207 } 00208 } 00209 00210 Sample 00211 PreChiSquare::getMinimum() 00212 { 00213 return minimum; 00214 } 00215 00216 void 00217 PreChiSquare::setDelta() 00218 { 00219 OccurrencesImpl::iterator it = my_presences->begin(); 00220 OccurrencesImpl::iterator last = my_presences->end(); 00221 00222 Sample sample = (*it)->environment(); 00223 delta = sample; 00224 ++it; 00225 while ( it != last ) 00226 { 00227 Sample const& sample = (*it)->environment(); 00228 delta |= sample; 00229 ++it; 00230 } 00231 00232 delta -= minimum; 00233 delta /= (Scalar)nclass; 00234 } 00235 00236 Sample 00237 PreChiSquare::getDelta() 00238 { 00239 return delta; 00240 } 00241 00242 void 00243 PreChiSquare::setMeasured(size_t layer1, size_t layer2) 00244 { 00245 OccurrencesImpl::iterator it = my_presences->begin(); 00246 OccurrencesImpl::iterator last = my_presences->end(); 00247 00248 size_t row, col; 00249 00250 for (size_t i = 0; i < classLimit; i++) 00251 for (size_t j = 0; j < classLimit; j++) 00252 measured[i][j] = 0.0; 00253 00254 while ( it != last ) 00255 { 00256 Sample const& sample = (*it)->environment(); 00257 00258 row = (size_t)floor( ( sample[layer1] - minimum[layer1] ) / delta[layer1] ); 00259 if (row == nclass) 00260 row--; 00261 col = (size_t)floor( ( sample[layer2] - minimum[layer2] ) / delta[layer2] ); 00262 if (col == nclass) 00263 col--; 00264 measured[row][col] += 1.0; 00265 00266 ++it; 00267 } 00268 } 00269 void 00270 PreChiSquare::setExpected() 00271 { 00272 size_t i, j; 00273 Scalar col_sum[classLimit], row_sum[classLimit], sum; 00274 00275 for (i = 0; i < nclass; i++) 00276 { 00277 col_sum[i] = measured[i][0]; 00278 for ( j = 1; j < nclass; j++) 00279 col_sum[i] += measured[i][j]; 00280 } 00281 00282 for (j = 0; j < nclass; j++) 00283 { 00284 row_sum[j] = measured[0][j]; 00285 for ( i = 1; i < nclass; i++) 00286 row_sum[j] += measured[i][j]; 00287 } 00288 00289 sum = col_sum[0]; 00290 for (i = 1; i < nclass; i++) 00291 sum += col_sum[i]; 00292 00293 for (i = 0; i < nclass; i++) 00294 for (j = 0; j < nclass; j++) 00295 expected[i][j] = col_sum[i] * row_sum[j] / sum; 00296 } 00297 void 00298 PreChiSquare::setChicell() 00299 { 00300 size_t i, j; 00301 Scalar aux; 00302 00303 for (i = 0; i < nclass; i++) 00304 for (j = 0; j < nclass; j++) 00305 if (expected[i][j] == 0.0) 00306 chicell[i][j] = 0.0; 00307 else 00308 { 00309 aux = expected[i][j] - measured[i][j]; 00310 chicell[i][j] = (aux * aux) / expected[i][j]; 00311 } 00312 } 00313 00314 void 00315 PreChiSquare::setStatistic(size_t layer1, size_t layer2) 00316 { 00317 size_t i, j; 00318 Scalar chi=0.0; 00319 Scalar delimita[classLimit] = {0, 3.8415, 5.9915, 7.8147, 9.4877, 11.0705, 12.5916, 14.0671, 15.5073, 16.9190, 18.3070, 19.6751, 21.0261, 22.3620, 23.6848, 24.9958 }; 00320 00321 for (i = 0; i < nclass; i++) 00322 for (j = 0; j < nclass; j++) 00323 chi += chicell[i][j]; 00324 00325 if (chi < delimita[nclass - 1]) 00326 { 00327 statistic1[layer1] += 1; 00328 statistic2[layer2] += 1; 00329 } 00330 }