PLearn 0.1
|
00001 // -*- C++ -*- 00002 00003 // PLearn (A C++ Machine Learning Library) 00004 // Copyright (C) 1999-2002 Christian Jauvin 00005 // 00006 00007 // Redistribution and use in source and binary forms, with or without 00008 // modification, are permitted provided that the following conditions are met: 00009 // 00010 // 1. Redistributions of source code must retain the above copyright 00011 // notice, this list of conditions and the following disclaimer. 00012 // 00013 // 2. Redistributions in binary form must reproduce the above copyright 00014 // notice, this list of conditions and the following disclaimer in the 00015 // documentation and/or other materials provided with the distribution. 00016 // 00017 // 3. The name of the authors may not be used to endorse or promote 00018 // products derived from this software without specific prior written 00019 // permission. 00020 // 00021 // THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``AS IS'' AND ANY EXPRESS OR 00022 // IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES 00023 // OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN 00024 // NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, 00025 // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED 00026 // TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR 00027 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF 00028 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 00029 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 00030 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 00031 // 00032 // This file is part of the PLearn library. For more information on the PLearn 00033 // library, go to the PLearn Web site at www.plearn.org 00034 00037 #include "ProbSparseMatrix.h" 00038 00039 namespace PLearn { 00040 00041 ProbSparseMatrix::ProbSparseMatrix(int n_rows, int n_cols, string name, int mode, bool double_access) : DoubleAccessSparseMatrix<real>(n_rows, n_cols, name, mode, double_access) 00042 { 00043 } 00044 00045 void ProbSparseMatrix::incr(int i, int j, real inc, bool warning) 00046 { 00047 if (inc <= 0.0 && warning) 00048 PLWARNING("incrementing value by: %g", inc); 00049 DoubleAccessSparseMatrix<real>::incr(i, j, inc); 00050 } 00051 00052 void ProbSparseMatrix::set(int i, int j, real value, bool warning) 00053 { 00054 if (value <= 0.0 && warning) 00055 PLWARNING("setting value: %g", value); 00056 DoubleAccessSparseMatrix<real>::set(i, j, value); 00057 } 00058 00059 bool ProbSparseMatrix::checkCondProbIntegrity() 00060 { 00061 real sum = 0.0; 00062 if (mode == ROW_WISE) 00063 { 00064 for (int i = 0; i < height; i++) 00065 { 00066 map<int, real>& row_i = rows[i]; 00067 sum = 0.0; 00068 for (map<int, real>::iterator it = row_i.begin(); it != row_i.end(); ++it) 00069 sum += it->second; 00070 if (fabs(sum - 1.0) > 1e-4 && (sum > 0.0)) 00071 return false; 00072 } 00073 return true; 00074 } else 00075 { 00076 for (int j = 0; j < width; j++) 00077 { 00078 map<int, real>& col_j = cols[j]; 00079 sum = 0.0; 00080 for (map<int, real>::iterator it = col_j.begin(); it != col_j.end(); ++it) 00081 sum += it->second; 00082 if (fabs(sum - 1.0) > 1e-4 && (sum > 0.0)) 00083 return false; 00084 } 00085 return true; 00086 } 00087 } 00088 00089 bool ProbSparseMatrix::checkJointProbIntegrity() 00090 { 00091 return (fabs(sumOfElements() - 1.0) > 1e-4); 00092 } 00093 00094 void ProbSparseMatrix::normalizeCond(ProbSparseMatrix& nXY, bool clear_nXY) 00095 { 00096 if (mode == ROW_WISE && (nXY.getMode() == ROW_WISE || nXY.isDoubleAccessible())) 00097 { 00098 clear(); 00099 int nXY_height = nXY.getHeight(); 00100 for (int i = 0; i < nXY_height; i++) 00101 { 00102 real sum_row_i = nXY.sumRow(i); 00103 map<int, real>& row_i = nXY.getRow(i); 00104 for (map<int, real>::iterator it = row_i.begin(); it != row_i.end(); ++it) 00105 { 00106 int j = it->first; 00107 real nij = it->second; 00108 real pij = nij / sum_row_i; 00109 if (pij > 0.0) 00110 set(i, j, pij); 00111 } 00112 } 00113 if (clear_nXY) 00114 nXY.clear(); 00115 } else if (mode == COLUMN_WISE && (nXY.getMode() == COLUMN_WISE || nXY.isDoubleAccessible())) 00116 { 00117 clear(); 00118 int nXY_width = nXY.getWidth(); 00119 for (int j = 0; j < nXY_width; j++) 00120 { 00121 real sum_col_j = nXY.sumCol(j); 00122 map<int, real>& col_j = nXY.getCol(j); 00123 for (map<int, real>::iterator it = col_j.begin(); it != col_j.end(); ++it) 00124 { 00125 int i = it->first; 00126 real nij = it->second; 00127 real pij = nij / sum_col_j; 00128 if (pij > 0.0) 00129 set(i, j, pij); 00130 } 00131 } 00132 if (clear_nXY) 00133 nXY.clear(); 00134 } else 00135 { 00136 PLERROR("pXY and nXY accessibility modes must match"); 00137 } 00138 } 00139 00140 void ProbSparseMatrix::normalizeCond() 00141 { 00142 if (mode == ROW_WISE) 00143 { 00144 for (int i = 0; i < height; i++) 00145 { 00146 real sum_row_i = sumRow(i); 00147 map<int, real>& row_i = rows[i]; 00148 for (map<int, real>::iterator it = row_i.begin(); it != row_i.end(); ++it) 00149 { 00150 int j = it->first; 00151 real nij = it->second; 00152 real pij = nij / sum_row_i; 00153 if (pij > 0.0) 00154 set(i, j, pij); 00155 } 00156 } 00157 } else 00158 { 00159 for (int j = 0; j < width; j++) 00160 { 00161 real sum_col_j = sumCol(j); 00162 map<int, real>& col_j = cols[j]; 00163 for (map<int, real>::iterator it = col_j.begin(); it != col_j.end(); ++it) 00164 { 00165 int i = it->first; 00166 real nij = it->second; 00167 real pij = nij / sum_col_j; 00168 if (pij > 0.0) 00169 set(i, j, pij); 00170 } 00171 } 00172 } 00173 } 00174 00177 void ProbSparseMatrix::normalizeJoint(ProbSparseMatrix& nXY, bool clear_nXY) 00178 { 00179 clear(); 00180 real sum_nXY = nXY.sumOfElements(); 00181 if (nXY.getMode() == ROW_WISE) 00182 { 00183 int nXY_height = nXY.getHeight(); 00184 for (int i = 0; i < nXY_height; i++) 00185 { 00186 map<int, real>& row_i = nXY.getRow(i); 00187 for (map<int, real>::iterator it = row_i.begin(); it != row_i.end(); ++it) 00188 { 00189 int j = it->first; 00190 real nij = it->second; 00191 real pij = nij / sum_nXY; 00192 if (pij > 0.0) 00193 set(i, j, pij); 00194 } 00195 } 00196 } else if (nXY.getMode() == COLUMN_WISE) 00197 { 00198 int nXY_width = nXY.getWidth(); 00199 for (int j = 0; j < nXY_width; j++) 00200 { 00201 map<int, real>& col_j = nXY.getCol(j); 00202 for (map<int, real>::iterator it = col_j.begin(); it != col_j.end(); ++it) 00203 { 00204 int i = it->first; 00205 real nij = it->second; 00206 real pij = nij / sum_nXY; 00207 if (pij > 0.0) 00208 set(i, j, pij); 00209 } 00210 } 00211 } 00212 if (clear_nXY) 00213 nXY.clear(); 00214 } 00215 00216 void ProbSparseMatrix::normalizeJoint() 00217 { 00218 if (mode == ROW_WISE) 00219 { 00220 for (int i = 0; i < height; i++) 00221 { 00222 map<int, real>& row_i = rows[i]; 00223 for (map<int, real>::iterator it = row_i.begin(); it != row_i.end(); ++it) 00224 { 00225 int j = it->first; 00226 real nij = it->second; 00227 real pij = nij / sumOfElements(); 00228 if (pij > 0.0) 00229 set(i, j, pij); 00230 } 00231 } 00232 } else 00233 { 00234 for (int j = 0; j < width; j++) 00235 { 00236 map<int, real>& col_j = cols[j]; 00237 for (map<int, real>::iterator it = col_j.begin(); it != col_j.end(); ++it) 00238 { 00239 int i = it->first; 00240 real nij = it->second; 00241 real pij = nij / sumOfElements(); 00242 if (pij > 0.0) 00243 set(i, j, pij); 00244 } 00245 } 00246 } 00247 } 00248 00250 real ProbSparseMatrix::euclidianDistance( ProbSparseMatrix &p) 00251 { 00252 real distance=0; 00253 real diff; 00254 if(p.getHeight()!=height || p.getWidth() != width) PLERROR("euclidianDistance: matrix dimensions do not match "); 00255 if (mode == ROW_WISE){ 00256 // go thru the first matrix 00257 for (int i = 0; i < height; i++){ 00258 map<int, real>& row_i = rows[i]; 00259 for (map<int, real>::iterator it = row_i.begin(); it != row_i.end(); ++it){ 00260 diff = it->second-p.get(i, it->first); 00261 distance +=sqrt(diff*diff); 00262 } 00263 } 00264 // go thru the second one 00265 for (int i = 0; i < p.getHeight(); i++){ 00266 map<int, real>& row_i = p.getRow(i); 00267 for (map<int, real>::iterator it = row_i.begin(); it != row_i.end(); ++it){ 00268 // if the value exists in the first matrix, it has already been included in the distance 00269 if(exists(i,it->first))continue; 00270 // no value in the first matrix 00271 diff = p.get(i,it->first); 00272 distance +=sqrt(diff*diff); 00273 } 00274 } 00275 }else{ 00276 // go thru the first matrix 00277 for (int j = 0; j < width; j++){ 00278 map<int, real>& col_j = cols[j]; 00279 for (map<int, real>::iterator it = col_j.begin(); it != col_j.end(); ++it){ 00280 diff = it->second-p.get(it->first,j); 00281 distance +=sqrt(diff*diff); 00282 } 00283 } 00284 // go thru the second one 00285 for (int j = 0; j < width; j++){ 00286 map<int, real>& col_j = p.getCol(j); 00287 for (map<int, real>::iterator it = col_j.begin(); it != col_j.end(); ++it){ 00288 // if the value exists in the first matrix, it has already been included in the distance 00289 if(exists(it->first,j))continue; 00290 // no value in the first matrix 00291 diff = p.get(it->first,j); 00292 distance +=sqrt(diff*diff); 00293 } 00294 } 00295 } 00296 return(distance); 00297 } 00298 00299 00300 00301 void ProbSparseMatrix::iterativeProportionalFittingStep( ProbSparseMatrix& p,Vec& lineMarginal, Vec& colMarginal) 00302 // one step of proportional iterative fitting on the matrix with lineMarginal and colMarginal 00303 { 00304 real newVal; 00305 real sum_row_i; 00306 real sum_col_j; 00307 00308 if(p.getHeight()!=lineMarginal.size() || p.getWidth() != colMarginal.size()) PLERROR("iterativeProportionalFittingStep: matrix dimension does not match marginal vectors dimensions"); 00309 if(p.getHeight()!=height || p.getWidth() != width) PLERROR("iterativeProportionalFittingStep: new matrix dimension does not match old matrix dimensions"); 00310 if(p.mode!=mode) PLERROR("iterativeProportionalFittingStep: Matrices access mode must match"); 00311 if (mode == ROW_WISE){ 00312 Vec sum_col(width); 00313 // First pass 00314 for (int i = 0; i < height; i++){ 00315 map<int, real>& row_i = p.getRow(i); 00316 sum_row_i = p.sumRow(i); 00317 for (map<int, real>::iterator it = row_i.begin(); it != row_i.end(); ++it){ 00318 if(sum_row_i==0)PLERROR("iterativeProportionalFittingStep: line %d is empty",i); 00319 newVal= it->second*lineMarginal[i]/sum_row_i; 00320 // store sum of column for next step 00321 sum_col[it->first]+=newVal; 00322 set(i,it->first,newVal); 00323 } 00324 } 00325 // Second Pass 00326 for (int i = 0; i < height; i++){ 00327 // we use the values set in the matrix at the previous stage 00328 map<int, real>& row_i = rows[i]; 00329 for (map<int, real>::iterator it = row_i.begin(); it != row_i.end(); ++it){ 00330 if(sum_col[it->first]==0)PLERROR("iterativeProportionalFittingStep: column %d is empty",i); 00331 newVal= it->second*colMarginal[it->first]/sum_col[it->first]; 00332 set(i,it->first,newVal); 00333 } 00334 } 00335 00336 }else{ 00337 Vec sum_row(height); 00338 for (int j = 0; j < width; j++){ 00339 map<int, real>& col_j = p.getCol(j); 00340 sum_col_j = p.sumCol(j); 00341 if( col_j.begin()!= col_j.end())cout << " " << colMarginal[j]/sum_col_j<< ":"; 00342 for (map<int, real>::iterator it = col_j.begin(); it != col_j.end(); ++it){ 00343 if(sum_col_j==0){ 00344 PLWARNING("iterativeProportionalFittingStep: column %d is empty",j); 00345 continue; 00346 } 00347 newVal= it->second*colMarginal[j]/sum_col_j; 00348 cout << " " <<it->second<<"="<<newVal; 00349 sum_row[it->first]+=newVal; 00350 set(it->first,j,newVal); 00351 } 00352 if( col_j.begin()!= col_j.end())cout << endl; 00353 } 00354 cout << endl; 00355 // Second Pass 00356 for (int j = 0; j < width; j++){ 00357 // we use the values set in the matrix at the previous stage 00358 map<int, real>& col_j = cols[j]; 00359 // cout << " " << lineMarginal[it->first]/sum_row[it->first]<< ":"; 00360 for (map<int, real>::iterator it = col_j.begin(); it != col_j.end(); ++it){ 00361 if(sum_row[it->first]==0){ 00362 PLWARNING("iterativeProportionalFittingStep: line %d is empty",it->first); 00363 continue; 00364 } 00365 newVal= it->second*lineMarginal[it->first]/sum_row[it->first]; 00366 cout << " " <<it->second<<"="<<newVal; 00367 set(it->first,j,newVal); 00368 } 00369 cout << endl; 00370 } 00371 } 00372 } 00373 00374 void ProbSparseMatrix::add( ProbSparseMatrix& p, ProbSparseMatrix& q) 00375 { 00376 real val; 00377 if(p.getHeight()!=q.getHeight() || p.getWidth() != q.getWidth()) PLERROR("euclidianDistance: matrix dimensions do not match "); 00378 if (mode == ROW_WISE){ 00379 // go thru the first matrix 00380 for (int i = 0; i < p.getHeight(); i++){ 00381 map<int, real>& row_i = p.getRow(i); 00382 for (map<int, real>::iterator it = row_i.begin(); it != row_i.end(); ++it){ 00383 val = it->second+q.get(i,it->first); 00384 set(i,it->first,val); 00385 } 00386 } 00387 // go thru the second one 00388 for (int i = 0; i < q.getHeight(); i++){ 00389 map<int, real>& row_i = q.getRow(i); 00390 for (map<int, real>::iterator it = row_i.begin(); it != row_i.end(); ++it){ 00391 // if the value exists in the first matrix, it has already been seen 00392 if(!p.exists(i,it->first)){ 00393 // no value in the first matrix 00394 set(i,it->first,it->second); 00395 } 00396 } 00397 } 00398 }else{ 00399 // go thru the first matrix 00400 for (int j = 0; j < p.getWidth(); j++){ 00401 map<int, real>& col_j = p.getCol(j); 00402 for (map<int, real>::iterator it = col_j.begin(); it != col_j.end(); ++it){ 00403 val = it->second+q.get(it->first,j); 00404 set(it->first,j,val); 00405 } 00406 } 00407 // go thru the second one 00408 for (int j = 0; j < q.getWidth(); j++){ 00409 map<int, real>& col_j = q.getCol(j); 00410 for (map<int, real>::iterator it = col_j.begin(); it != col_j.end(); ++it){ 00411 // if the value exists in the first matrix, it has already been seen 00412 if(!p.exists(it->first,j)){ 00413 // no value in the first matrix 00414 set(it->first,j,it->second); 00415 } 00416 } 00417 } 00418 } 00419 00420 } 00421 00422 } 00423 00424 00425 /* 00426 Local Variables: 00427 mode:c++ 00428 c-basic-offset:4 00429 c-file-style:"stroustrup" 00430 c-file-offsets:((innamespace . 0)(inline-open . 0)) 00431 indent-tabs-mode:nil 00432 fill-column:79 00433 End: 00434 */ 00435 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=79 :