PLearn 0.1
ProbSparseMatrix.cc
Go to the documentation of this file.
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 :
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines