PLearn 0.1
vmatmain.cc
Go to the documentation of this file.
00001 // -*- C++ -*-
00002 
00003 // vmatmain.cc
00004 // Copyright (C) 2002 Pascal Vincent, Julien Keable, Xavier Saint-Mleux, Rejean Ducharme
00005 //
00006 // Redistribution and use in source and binary forms, with or without
00007 // modification, are permitted provided that the following conditions are met:
00008 // 
00009 //  1. Redistributions of source code must retain the above copyright
00010 //     notice, this list of conditions and the following disclaimer.
00011 // 
00012 //  2. Redistributions in binary form must reproduce the above copyright
00013 //     notice, this list of conditions and the following disclaimer in the
00014 //     documentation and/or other materials provided with the distribution.
00015 // 
00016 //  3. The name of the authors may not be used to endorse or promote
00017 //     products derived from this software without specific prior written
00018 //     permission.
00019 // 
00020 // THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``AS IS'' AND ANY EXPRESS OR
00021 // IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
00022 // OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN
00023 // NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
00024 // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
00025 // TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00026 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00027 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00028 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00029 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00030 // 
00031 // This file is part of the PLearn library. For more information on the PLearn
00032 // library, go to the PLearn Web site at www.plearn.org
00033 
00034 /* *******************************************************      
00035  * $Id: vmatmain.cc 10283 2009-07-22 15:20:55Z nouiz $
00036  ******************************************************* */
00037 
00038 #include "vmatmain.h"
00039 #include <commands/PLearnCommands/PLearnCommandRegistry.h>
00040 #include <plearn/base/general.h>
00041 #include <plearn/base/stringutils.h>
00042 #include <plearn/base/lexical_cast.h>
00043 #include <plearn/math/StatsCollector.h>
00044 #include <plearn/math/stats_utils.h>
00045 #include <plearn/vmat/ConcatColumnsVMatrix.h>
00046 #include <plearn/vmat/SelectColumnsVMatrix.h>
00047 #include <plearn/vmat/SubVMatrix.h>
00048 #include <plearn/vmat/VMatLanguage.h>
00049 #include <plearn/vmat/VVMatrix.h>
00050 #include <plearn/vmat/VMat.h>
00051 #include <plearn/vmat/SelectRowsFileIndexVMatrix.h>
00052 #include <plearn/math/TMat_maths.h>
00053 #include <plearn/db/getDataSet.h>
00054 #include <plearn/display/Gnuplot.h>
00055 #include <plearn/io/openFile.h>
00056 #include <plearn/io/load_and_save.h>
00057 #include <plearn/base/HelpSystem.h>
00058 #include <plearn/base/PDate.h>
00059 #include <algorithm>                         // for max
00060 #include <iostream>
00061 #include <iomanip>                           // for setw and such
00062 
00063 namespace PLearn {
00064 using namespace std;
00065 
00077 static void save_vmat_as_csv(VMat source, ostream& destination,
00078                              bool skip_missings, int precision = 12,
00079                              string delimiter = ",",
00080                              bool verbose = true,
00081                              bool convert_date = false)
00082 {
00083     char buffer[1000];
00084   
00085     // First, output the fieldnames in quoted CSV format.  Don't forget
00086     // to quote the quotes
00087     TVec<string> fields = source->fieldNames();
00088     for (int i=0, n=fields.size() ; i<n ; ++i) {
00089         string curfield = fields[i];
00090         search_replace(curfield, "\"", "\\\"");
00091         destination << '"' << curfield << '"';
00092         if (i < n-1)
00093             destination << delimiter;
00094     }
00095     destination << "\n";
00096 
00097     PP<ProgressBar> pb;
00098     if (verbose)
00099         pb = new ProgressBar(cout, "Saving to CSV", source.length());
00100 
00101     // Next, output each line.  Perform missing-value checks if required.
00102     for (int i=0, n=source.length() ; i<n ; ++i) {
00103         if (pb)
00104             pb->update(i+1);
00105         Vec currow = source(i);
00106         if (! skip_missings || ! currow.hasMissing()) {
00107             for (int j=0, m=currow.size() ; j<m ; ++j) {
00108                 string strval="";
00109                 if (convert_date && j==0)
00110                     // Date conversion: add 19000000 to convert from CYYMMDD to
00111                     // YYYYMMDD, and always output without trailing . if not
00112                     // necessary
00113                     sprintf(buffer, "%8f", currow[j] + 19000000.0);
00114                 else if((strval=source->getValString(j,currow[j]))!=""){
00115                     search_replace(strval, "\"", "\\\"");
00116                     strval = "\"" + strval + "\"";
00117                     if(strval.length()>1000-1)
00118                         PLERROR("a value is too big!");
00119                     strncpy(buffer,strval.c_str(),1000);
00120                 }else{
00121                     // Normal processing
00122                     sprintf(buffer, "%#.*f", precision, currow[j]);
00123 
00124                     // strip all trailing zeros and final period
00125                     // there is always a period since sprintf includes # modifier
00126                     char* period = buffer;
00127                     while (*period && *period != '.')
00128                         period++;
00129                     for (char* last = period + strlen(period) - 1 ;
00130                          last >= period && (*last == '0' || *last == '.') ; --last) {
00131                         bool should_break = *last == '.';
00132                         *last = '\0';
00133                         if (should_break)
00134                             break;
00135                     }
00136                 }
00137                 destination << buffer;
00138                 if (j < m-1)
00139                     destination << delimiter;
00140             }
00141             destination << "\n";
00142         }
00143     }
00144 }
00145 
00146   
00156 static void save_vmat_as_arff(VMat source, ostream& destination, TVec<string>& date_columns,
00157                              bool skip_missings, int precision = 12, bool verbose = true)
00158 {
00159     PP<ProgressBar> pb;
00160     if (verbose)
00161         pb = new ProgressBar(cout, "Saving to ARFF", source.length());
00162 
00163     // First, write the ARFF specific file header
00164     destination << "@relation arff-database\n";
00165     TVec<string> fields = source->fieldNames();
00166     int nb_fields = fields.size();
00167     for (int i=0; i<nb_fields; ++i)
00168     {
00169         string curfield = fields[i];
00170         destination << "@attribute " << curfield << " ";
00171         map<string,real> string_map = source->getStringToRealMapping(i);
00172         if (date_columns.contains(curfield))  // e.g. @attribute start_date date "yyyyMMdd"
00173             destination << "date \"yyyyMMdd\"\n";
00174         else if (string_map.empty())  // e.g. @attribute Age numeric
00175             destination << "numeric\n";
00176         else  // e.g. @attribute Country {"Canada", "China", "Columbia"}
00177         {
00178             int nb_keys = string_map.size();
00179             int key_i = 0;
00180             destination << "{";
00181             map<string,real>::iterator it;
00182             for (it = string_map.begin(); it != string_map.end(); ++it)
00183             {
00184                 string key = it->first;
00185                 search_replace(key, "\"", "\\\"");
00186                 destination << '"' << key << '"';
00187                 if (key_i < nb_keys-1)
00188                     destination << ", ";
00189                 key_i++;
00190             }
00191             destination << "}\n";
00192         }
00193     }
00194     destination << "@data\n";
00195 
00196     // Next, output each line.  Perform missing-value checks if required.
00197     const string delimiter = ",";
00198     const int buffer_size = 10000;
00199     char buffer[buffer_size];
00200     for (int i=0, n=source.length(); i<n; ++i) {
00201         if (pb)
00202             pb->update(i+1);
00203 
00204         // Skip missing values?
00205         Vec currow = source(i);
00206         if (skip_missings && currow.hasMissing())
00207             continue;
00208 
00209         for (int j=0, m=currow.size(); j<m; ++j)
00210         {
00211             string strval = "";
00212             // Date field
00213             if (date_columns.contains(fields[j]))
00214             {
00215                 // Date conversion: add 19000000 to convert from CYYMMDD to
00216                 // YYYYMMDD, and always output without trailing . if not
00217                 // necessary
00218                 real curfield = currow[j];
00219                 if (is_missing(curfield)  ||  curfield==0)  // missing value
00220                 {
00221                     strval = "?";
00222                     strncpy(buffer, strval.c_str(), buffer_size);
00223                 }
00224                 else
00225                 {
00226                     if (curfield < 10000000)  // CYYMMDD format
00227                         curfield += 19000000.0;
00228                     sprintf(buffer, "%8d", int(curfield));
00229                 }
00230             }
00231             // String mapped field
00232             else if ((strval=source->getValString(j,currow[j])) != "")
00233             {
00234                 search_replace(strval, "\"", "\\\"");
00235                 strval = "\"" + strval + "\"";
00236                 if (strval.length()>buffer_size-1)
00237                     PLERROR("a value is too big!");
00238                 strncpy(buffer, strval.c_str(), buffer_size);
00239             }
00240             // Numeric field
00241             else
00242             {
00243                 real curfield = currow[j];
00244                 if (is_missing(curfield))  // missing value
00245                 {
00246                     strval = "?";
00247                     strncpy(buffer, strval.c_str(), buffer_size);
00248                 }
00249                 else
00250                 {
00251                     // Normal processing
00252                     sprintf(buffer, "%#.*f", precision, currow[j]);
00253 
00254                     // strip all trailing zeros and final period
00255                     // there is always a period since sprintf includes # modifier
00256                     char* period = buffer;
00257                     while (*period && *period != '.')
00258                         period++;
00259                     for (char* last = period + strlen(period) - 1; last >= period && (*last == '0' || *last == '.'); --last)
00260                     {
00261                         bool should_break = *last == '.';
00262                         *last = '\0';
00263                         if (should_break)
00264                             break;
00265                     }
00266                 }
00267             }
00268             destination << buffer;
00269             if (j < m-1)
00270                 destination << delimiter;
00271         }
00272         destination << "\n";
00273     }
00274 }
00275 
00276   
00279 int print_diff(ostream& out, VMat m1, VMat m2, double tolerance, int verbose)
00280 {
00281     int ndiff = 0;
00282     if(m1.length()!=m2.length() || m1.width()!=m2.width())
00283     {
00284         out << "Size of the two matrices differ: " 
00285             << m1.length() << " x " << m1.width() << "  vs.  "
00286             << m2.length() << " x " << m2.width() << endl;
00287         return -1;
00288     }
00289     if(m1->getFieldInfos()!=m2->getFieldInfos()){
00290         Array<VMField> a1=m1->getFieldInfos();
00291         Array<VMField> a2=m2->getFieldInfos();
00292         if(verbose)
00293             pout << "Field infos differ:";
00294         //compare fieldnames
00295         for(int i=0;i<m1.width();i++){
00296             if(a1[i].name!=a2[i].name){
00297                 ++ndiff;
00298                 if(verbose)
00299                     pout << " " << a1[i].name << "!=" << a2[i].name;
00300             }
00301         }
00302         //compare fieldtype
00303         for(int i=0;i<m1.width();i++){
00304             if(a1[i].fieldtype!=a2[i].fieldtype){
00305                 ++ndiff;
00306                 if(verbose)
00307                     pout << " " << a1[i].fieldtype << "!=" << a2[i].fieldtype;
00308             }
00309         }
00310         if(verbose)
00311             pout <<endl;
00312         
00313     }
00314     int l = m1.length();
00315     int w = m1.width();
00316     Vec v1(w);
00317     Vec v2(w);
00318     for(int i=0; i<l; i++)
00319     {
00320         m1->getRow(i,v1);
00321         m2->getRow(i,v2);
00322         for(int j=0; j<w; j++)
00323         {
00324             if (!is_equal(v1[j], v2[j], 1.0, real(tolerance), real(tolerance)))
00325             {
00326                 if (verbose)
00327                     out << "Elements at " << i << ',' << j << " differ by "
00328                         << v1[j] - v2[j] << endl;
00329                 ++ndiff;
00330             } else if (m1->getValString(j, v1[j]) != m2->getValString(j, v2[j])) {
00331                 if (verbose)
00332                     out << "Elements at " << i << ',' << j << " differ: "
00333                         << "'" << m1->getValString(j, v1[j]) << "' != "
00334                         << "'" << m2->getValString(j, v2[j]) << "'" << endl;
00335                 ++ndiff;
00336             }
00337         }
00338     }
00339     if (!verbose) out << ndiff <<endl;
00340     return ndiff;
00341 }
00342 
00343 void interactiveDisplayCDF(const Array<VMat>& vmats)
00344 {
00345     int k = vmats.size();
00346     int w = vmats[0]->width();
00347 
00348     Array<string> name(k);
00349     pout << ">>>> Dimensions of vmats: \n";
00350     for(int i=0; i<k; i++)
00351     {
00352         name[i] = vmats[i]->getMetaDataDir();
00353         name[i] = name[i].erase(name[i].size()-10);
00354         pout << name[i] << ": \t " << vmats[i]->length() << " x " << vmats[i]->width() << endl;
00355     }
00356 
00357     vmats[0]->printFields(pout);
00358 
00359     Gnuplot gp;
00360 
00361     for(;;)
00362     {
00363         // TVec<RealMapping> ranges = vm->getRanges();
00364 
00365         vector<string> command;
00366         int varnum = -1;
00367         real low = -FLT_MAX; // means autorange
00368         real high = FLT_MAX; // means autorange
00369         do
00370         {
00371             pout << "Field (0.." << w-1 << ") [low high] or exit? " << flush;
00372             command = split(pgetline(cin));
00373             if(command.size()==0)
00374                 vmats[0]->printFields(pout);
00375             else if(command.size()==1&&command[0]=="exit")
00376                 exit(0);
00377             else
00378             {
00379                 varnum = vmats[0]->getFieldIndex(command[0],false);
00380                 if(varnum == -1)
00381                     pout<<"Bad column name or number("<<command[0]<<")"<<endl;
00382                 if(varnum<0 || varnum>=w)
00383                     vmats[0]->printFields(pout);
00384                 else if(command.size()==3)
00385                 {
00386                     low = toreal(command[1]);
00387                     high = toreal(command[2]);
00388                 }
00389             }
00390         } while(varnum<0 || varnum>=w);
00391 
00392 
00393         pout << "\n\n*************************************" << endl;
00394         pout << "** #" << varnum << ": " << vmats[0]->fieldName(varnum) << " **" << endl;
00395         pout << "*************************************" << endl;
00396 
00397         Array<Mat> m(k);
00398 
00399         for(int i=0; i<k; i++)
00400         {
00401             TVec<StatsCollector> stats = vmats[i]->getStats();        
00402             StatsCollector& st = stats[varnum];
00403             m[i] = st.cdf(true);
00404             pout << "[ " << name[i]  << " ]" << endl;
00405             pout << st << endl;
00406         }
00407         // pout << "RANGES: " << endl;
00408         // pout << ranges[varnum];
00409 
00410         if(is_equal(low,-FLT_MAX))
00411             gp << "set xrange [*:*]" << endl;      
00412         else
00413             //The -0.01 and 0.01 is to clearly see the end value.
00414             gp << "set xrange [" << low-0.01 << ":" << high+0.01 << "]" << endl;
00415 
00416         if(k>=4)
00417             gp.plot(m[0],"title '"+name[0]+"'", m[1], "title '" + name[1]+"'", m[2], "title '" + name[2]+"'", m[3], "title '"+name[3]+"'");    
00418         else if(k>=3)
00419             gp.plot(m[0],"title '"+name[0]+"'", m[1], "title '"+name[1]+"'", m[2], "title '"+name[2]+"'");
00420         else if(k>=2)
00421             gp.plot(m[0],"title '"+name[0]+"'", m[1], "title '"+name[1]+"'");
00422         else
00423             gp.plot(m[0],"title '"+name[0]+"'");
00424     }
00425 }
00426 
00427 void displayBasicStats(VMat vm)
00428 {
00429     int nfields = vm.width();
00430     TVec<StatsCollector> stats = vm->getStats(true);        
00431 
00432     // find longest field name
00433     size_t fieldlen = 0;
00434     for (int k=0; k<nfields; ++k)
00435         fieldlen = max(fieldlen, vm->fieldName(k).size());
00436     fieldlen++;
00437   
00438     cout << std::left << setw(6)  << "# "
00439          << setw(int(fieldlen)) << " fieldname " << std::right
00440          << setw(15) << " mean "
00441          << setw(15) << " stddev "
00442          << setw(15) << " min "
00443          << setw(15) << " max "
00444          << setw(15) << " count "
00445          << setw(15) << " nmissing "
00446          << setw(15) << " stderr" << endl; 
00447     for(int k=0; k<nfields; k++)
00448     {
00449         cout << std::left << setw(6)  << k << " " 
00450              << setw(int(fieldlen)) << vm->fieldName(k) << " " << std::right
00451              << setw(15) << stats[k].mean() << " " 
00452              << setw(15) << stats[k].stddev() << " "
00453              << setw(15) << stats[k].min() << " " 
00454              << setw(15) << stats[k].max() << " " 
00455              << setw(15) << stats[k].n() << " " 
00456              << setw(15) << stats[k].nmissing() << " " 
00457              << setw(15) << stats[k].stderror() << " " 
00458              << endl;
00459     }
00460 }
00461 
00462 
00463 void printDistanceStatistics(VMat vm, int inputsize)
00464 {
00465     int l = vm.length();
00466     int w = vm.width();
00467     Vec x1(w);
00468     Vec x2(w);
00469     StatsCollector collector(2);  
00470     ProgressBar pb(cerr, "Computing distance statistics", l-1);
00471     for(int i=0; i<l-1; i++)
00472     {
00473         vm->getRow(i,x1);
00474         vm->getRow(i+1,x2);
00475         real d = L2distance(x1.subVec(0,inputsize),x2.subVec(0,inputsize));
00476         collector.update(d);
00477         pb(i);
00478     }
00479 
00480     pout << "Euclidean distance statistics: " << endl;
00481     pout << collector << endl;
00482 }
00483 
00484 /*
00485   void printConditionalStats(VMat vm, int condfield)
00486   {
00487   cout << "*** Ranges ***" << endl;
00488   TVec<RealMapping> ranges = vm->getRanges();
00489   PP<ConditionalStatsCollector> st = vm->getConditionalStats(condfield);
00490   int w = vm->width();
00491   for(int i=0; i<w; i++)
00492   {
00493   cout << "Field #" << i << ": " << vm->fieldName(i) << endl;
00494   cout << "Ranges: " << ranges[i] << endl;
00495   }
00496   cout << "\n\n------------------------------------------------------------" << endl;
00497   cout << "** Raw counts conditioned on field #" << condfield << " (" << vm->fieldName(condfield) << ") **\n" << endl;
00498   for(int k=0; k<w; k++)
00499   {
00500   cout << "#" << k << " " << vm->fieldName(condfield) << endl;
00501   cout << st->counts[k] << endl;
00502   }
00503   
00504   cout << "\n\n------------------------------------------------------------" << endl;
00505   cout << "** Joint probabilities (percentage) **\n" << endl;
00506   for(int k=0; k<w; k++)
00507   {
00508   TMat<int>& C = st->counts[k];
00509   Mat m(C.length(), C.width());
00510   m << C;
00511   m /= sum(m);
00512   m *= real(100);
00513   cout << "#" << k << " " << vm->fieldName(condfield) << endl;
00514   cout << m << endl;
00515   }
00516 
00517   cout << "\n\n------------------------------------------------------------" << endl;
00518   cout << "** Conditional probabilities conditioned on << " << vm->fieldName(condfield) << "  **\n" << endl;
00519   for(int k=0; k<w; k++)
00520   {
00521   TMat<int>& C = st->counts[k];
00522   Mat m(C.length(), C.width());      
00523   m << C;
00524   normalizeRows(m);
00525   m *= real(100);
00526   cout << "#" << k << " " << vm->fieldName(condfield) << endl;
00527   cout << m << endl;
00528   }
00529 
00530   cout << "\n\n------------------------------------------------------------" << endl;
00531   cout << "** Conditional probabilities conditioned on the other variables **\n" << endl;
00532   for(int k=0; k<w; k++)
00533   {
00534   TMat<int>& C = st->counts[k];
00535   Mat m(C.length(), C.width());      
00536   m << C;
00537   normalizeColumns(m);
00538   m *= real(100);
00539   cout << "#" << k << " " << vm->fieldName(condfield) << endl;
00540   cout << m << endl;
00541   }
00542 
00543 
00544   }
00545 */
00546 
00547 /*
00548   int findNextIndexOfValue(VMat m, int col, real value, int startrow=0)
00549   {
00550   if(m->hasMetaDataDir())
00551   {
00552   string fpath = apppend_slash(m->getMetaDataDir())+"CachedColumns/"+tostring(col);
00553   if(!isfile(filepath))
00554         
00555       
00556   }
00557   }
00558 */
00559 
00560 
00561 void plotVMats(char* defs[], int ndefs)
00562 {
00563     /* defs[] is of format:
00564        { "<dataset0>", "<col0>[:<row0>:<nrows0>]", ..., "<datasetN>", "<colN>[:<rowN>:<nrowsN>]" }
00565     */
00566     int nseries= ndefs/2;
00567     TmpFilenames tmpfnames(nseries, "/tmp/", "_vmat_plot_");
00568     Array<VMat> vmats(nseries);
00569     Array<Vec> series(nseries);
00570     string gp_command= "plot ";
00571     for(int i= 0; i < nseries; ++i)
00572     {
00573         vmats[i]= getDataSet(string(defs[2*i]));
00574 
00575         vector<string> spec= PLearn::split(defs[2*i+1], ":");
00576       
00577         series[i].resize(vmats[i].length());
00578         vmats[i]->getColumn(toint(spec[0]),series[i]);
00579 
00580         if(spec.size() == 3)
00581         {
00582             int row= toint(spec[1]), nrows= toint(spec[2]);
00583             if(row+nrows > series[i].length())
00584                 nrows= series[i].length()-row;
00585             series[i]= series[i].subVec(row, nrows);
00586         }
00587         else if(spec.size() != 1)
00588             PLERROR("in plotVMats: invalid spec for vmat %s: '%s'; sould be '<col>[:<row>:<nrows>]'.",
00589                     defs[2*i], defs[2*i+1]);
00590 
00591         saveGnuplot(tmpfnames[i].c_str(), series[i]);
00592         chmod(tmpfnames[i].c_str(),0777);      
00593         gp_command+= " '" + tmpfnames[i] + "' title '" + defs[2*i] + ' ' + defs[2*i+1] + "' " + tostring(i+1)  +", ";
00594     }
00595     gp_command.resize(gp_command.length()-2);
00596 
00597     Gnuplot gp;
00598     gp << gp_command << endl;
00599   
00600     pout << "Press any key to close GNUplot window and exit." << endl;
00601     cin.get();
00602 }
00603 
00604 
00605 VMat getVMat(const PPath& source, const PPath& indexf)
00606 {
00607     VMat vm= getDataSet(source);
00608     if(indexf != "")
00609         vm= new SelectRowsFileIndexVMatrix(vm, indexf);
00610     return vm;
00611 }
00612 
00613 
00614 int vmatmain(int argc, char** argv)
00615 {
00616   
00617     if(argc<3)
00618     {
00619         // Use the VMatCommand help instead of repeating the same help message twice...
00620         // help message in file commands/PLearnCommands/VMatCommand.cc
00621         pout << HelpSystem::helpOnCommand("vmat") << flush;
00622         exit(0);
00623     }
00624 
00625     PPath indexf= "";
00626     if(string(argv[1])=="-i")
00627     {
00628         indexf= argv[2];
00629         argv+= 2;//skip -i and indexfile name
00630     }
00631 
00632     string command = argv[1];
00633 
00634     if(command=="cdf")
00635     {      
00636         Array<VMat> vmats;
00637         for(int i=2; i<argc; i++)
00638         {
00639             string dbname = argv[i];
00640             VMat vm = getVMat(dbname, indexf);
00641             vmats.append(vm);
00642         }
00643         interactiveDisplayCDF(vmats);
00644     }
00645     /*
00646       else if(command=="cond")
00647       {
00648       string dbname = argv[2];
00649       VMat vm = getDataSet(dbname);
00650       cout << "** Using dataset: " << dbname << " **" << endl;
00651       cout << "Metadata for this dataset in: " << vm->getMetaDataDir() << endl;
00652       int condfield = atoi(argv[3]);
00653       printConditionalStats(vm, condfield);    
00654       }
00655     */
00656     else if(command=="convert")
00657     {
00658         if(argc<4)
00659             PLERROR("Usage: vmat convert <source> <destination> "
00660                     "[--mat_to_mem] [--cols=col1,col2,col3,...] [--save_vmat] [--skip-missings] [--precision=N] [--delimiter=CHAR] [--force_float] [--auto_float]");
00661 
00662         PPath source = argv[2];
00663         PPath destination = argv[3];
00664         bool mat_to_mem = false;
00665         if(source==destination)
00666             PLERROR("You are overwriting the source. This is not allowed!");
00698         TVec<string> columns;
00699         TVec<string> date_columns;
00700         bool skip_missings = false;
00701         int precision = 12;
00702         string delimiter = ",";
00703         bool convert_date = false;
00704         bool save_vmat = false;
00705         bool update = false;
00706         bool force_float = false;
00707         bool auto_float = false;
00708 
00709         string ext = extract_extension(destination);
00710 
00711         for (int i=4 ; i < argc && argv[i] ; ++i) {
00712             string curopt = removeblanks(argv[i]);
00713             if (curopt == "")
00714                 continue;
00715             if (curopt.substr(0,7) == "--cols=") {
00716                 string columns_str = curopt.substr(7);
00717                 columns = split(columns_str, ',');
00718             }
00719             else if (curopt.substr(0,12) == "--date-cols=") {
00720                 string columns_str = curopt.substr(12);
00721                 date_columns = split(columns_str, ',');
00722             }
00723             else if (curopt == "--skip-missings")
00724                 skip_missings = true;
00725             else if (curopt.substr(0,12) == "--precision=") {
00726                 precision = toint(curopt.substr(12));
00727             }
00728             else if (curopt.substr(0,12) == "--delimiter=") {
00729                 if(ext!=".csv")
00730                     PLERROR("Vmat convert: the --delimiter option is supported only with .csv destination file. You have a '%s' extension.",ext.c_str());
00731                 delimiter = curopt.substr(12);
00732             }
00733             else if (curopt == "--convert-date")
00734                 convert_date = true;
00735             else if (curopt =="--mat_to_mem")
00736                 mat_to_mem = true;
00737             else if (curopt == "--save_vmat")
00738                 save_vmat = true;
00739             else if (curopt == "--update")
00740                 update = true;
00741             else if (curopt == "--force_float"){
00742                 PLCHECK(ext==".pmat");
00743                 force_float = true;
00744             }else if (curopt == "--auto_float"){
00745                 PLCHECK(ext==".pmat");
00746                 auto_float = true;
00747             }else
00748                 PLWARNING("VMat convert: unrecognized option '%s'; ignoring it...",
00749                           curopt.c_str());
00750         }
00751         VMat vm = getVMat(source, indexf);
00752 
00753         // If columns specified, select them.  Note: SelectColumnsVMatrix is very
00754         // powerful and allows ranges, etc.
00755         if (columns.size() > 0)
00756             vm = new SelectColumnsVMatrix(vm, columns);
00757 
00758         if (ext != ".csv" && skip_missings)
00759             PLWARNING("Option '--skip-missings' not supported for extension '%s'; ignoring it...",
00760                       ext.c_str());
00761         if(mat_to_mem)
00762             vm.precompute();
00763         if(update && vm->isUpToDate(destination))
00764             pout << "The file is up to date. We don't regenerate it."<<endl;
00765         else if(ext==".amat")
00766             // Save strings as strings so they are not lost.
00767             vm->saveAMAT(destination, true, false, true);
00768         else if(ext==".cmat")
00769             vm->saveCMAT(destination);
00770         else if(ext==".pmat")
00771             vm->savePMAT(destination, force_float, auto_float);
00772         else if(ext==".dmat")
00773             vm->saveDMAT(destination);
00774         else if(ext == ".csv")
00775         {
00776             if (destination == "-.csv")
00777                 save_vmat_as_csv(vm, cout, skip_missings, precision, delimiter, true /*verbose*/, convert_date);
00778             else
00779             {
00780                 ofstream out(destination.c_str());
00781                 save_vmat_as_csv(vm, out, skip_missings, precision, delimiter, true /*verbose*/, convert_date);
00782             }
00783         }
00784         else if(ext == ".arff")
00785         {
00786             ofstream out(destination.c_str());
00787             save_vmat_as_arff(vm, out, date_columns, skip_missings, precision);
00788         }
00789         else if(ext == ".vmat")
00790             PLearn::save(destination,vm);
00791         else
00792         {
00793             cerr << "ERROR: can only convert to .amat .pmat .dmat, .vmat or .csv" << endl
00794                  << "Please specify a destination name with a valid extension " << endl;
00795         }
00796         if(save_vmat && extract_extension(source)==".vmat")
00797             PLearn::save(destination+".metadata/orig.vmat",vm);
00798         else if(save_vmat)
00799             PLWARNING("We haven't saved the original file as it is not a vmat");
00800     }
00801     else if(command=="info")
00802     {
00803         for(int i=2;i<argc;i++){
00804             string dbname = argv[i];
00805             VMat vm = getVMat(dbname, indexf);
00806             if(argc>3)
00807                 pout<<dbname<<endl;
00808             pout<<vm.length()<<" x "<<vm.width()<<endl;
00809             pout << "inputsize: " << vm->inputsize() << endl;
00810             pout << "targetsize: " << vm->targetsize() << endl;
00811             pout << "weightsize: " << vm->weightsize() << endl;
00812             pout << "extrasize: " << vm->extrasize() << endl;
00813             VVMatrix * vvm = dynamic_cast<VVMatrix*>((VMatrix*)vm);
00814             if(vvm!=NULL)
00815             {
00816                 pout<< "Last modification (including dependencies of .vmat): "
00817                     << int32_t(vvm->getMtime()) << endl;
00818                 bool ispre=vvm->isPrecomputedAndUpToDate();
00819                 pout<<"precomputed && uptodate : ";
00820                 if(ispre)
00821                 {
00822                     pout <<"yes : " << vvm->getPrecomputedDataName()<<endl;
00823                     pout<< "timestamp of precom. data : "
00824                         <<int32_t(getDataSetDate(vvm->getPrecomputedDataName()))
00825                         << endl;
00826                 }
00827                 else pout <<"no"<<endl;
00828             }
00829         }
00830     }
00831     else if(command=="fields")
00832     {
00833         bool add_info = true;
00834         bool transpose = false;
00835         if (argc >= 4) {
00836             add_info = !(string(argv[3]) == "name_only");
00837         }
00838         if (argc >= 5) {
00839             transpose = (string(argv[4]) == "transpose");
00840         }
00841         string dbname = argv[2];
00842         VMat vm = getVMat(dbname, indexf);
00843         if (add_info) {
00844             pout<<"FieldNames: ";
00845             if (!transpose) {
00846                 pout << endl;
00847             }
00848         }
00849         for(int i=0;i<vm.width();i++) {
00850             if (add_info) {
00851                 pout << i << ": ";
00852             }
00853             pout << vm->fieldName(i);
00854             if (transpose) {
00855                 pout << " ";
00856             } else {
00857                 pout << endl;
00858             }
00859         }
00860         if (transpose) {
00861             // It misses a carriage return after everything is displayed.
00862             pout << endl;
00863         }
00864     }
00865     else if(command=="fieldinfo")
00866     {
00867         if (argc < 4)
00868             PLERROR("The 'fieldinfo' subcommand requires more parameters, please check the help");
00869         string dbname = argv[2];
00870         string fieldname_or_num = argv[3];
00871 
00872         bool print_binning = false;
00873         if (argc == 5) {
00874             if (argv[4] == string("--bin"))
00875                 print_binning = true;
00876             else
00877                 PLERROR("vmat fieldinfo: unrecognized final argument; can be '--bin' "
00878                         "to print the binning");
00879         }
00880         
00881         VMat vm = getVMat(dbname, indexf);
00882         vm->printFieldInfo(pout, fieldname_or_num, print_binning);
00883     }
00884     else if(command=="stats")
00885     {
00886         string dbname = argv[2];
00887         VMat vm = getVMat(dbname, indexf);
00888         displayBasicStats(vm);
00889     }
00890     else if(command=="gendef")
00891     {
00892         string dbname = argv[2];
00893         TVec<int> bins(argc-3);
00894         for(int i=3;i<argc;i++)
00895             bins[i-3]=toint(argv[i]);
00896       
00897         VMat vm = getVMat(dbname, indexf);
00898         TVec<StatsCollector> sc = vm->getStats();
00899         // write stats file in metadatadir
00900         string name = vm->getMetaDataDir()+"/stats.def";
00901         ofstream sfile(name.c_str());
00902         for(int i=0;i<sc.size();i++)
00903         {
00904             sfile<<"DEFINE @"<<vm->fieldName(i)<<".mean "<<tostring(sc[i].mean())<<endl;
00905             sfile<<"DEFINE @"<<vm->fieldName(i)<<".stddev "<<tostring(sc[i].stddev())<<endl;
00906             sfile<<"DEFINE @"<<vm->fieldName(i)<<".stderr "<<tostring(sc[i].stderror())<<endl;
00907             sfile<<"DEFINE @"<<vm->fieldName(i)<<".min "<<tostring(sc[i].min())<<endl;
00908             sfile<<"DEFINE @"<<vm->fieldName(i)<<".max "<<tostring(sc[i].max())<<endl;
00909             sfile<<"DEFINE @"<<vm->fieldName(i)<<".normalized @"<<vm->fieldName(i)<<" @"<<vm->fieldName(i)<<".mean - @"<<
00910                 vm->fieldName(i)<<".stddev /"<<endl;
00911             sfile<<"DEFINE @"<<vm->fieldName(i)<<".sum "<<tostring(sc[i].sum())<<endl;
00912             sfile<<"DEFINE @"<<vm->fieldName(i)<<".sumsquare "<<tostring(sc[i].sumsquare())<<endl;
00913             sfile<<"DEFINE @"<<vm->fieldName(i)<<".variance "<<tostring(sc[i].variance())<<endl;
00914         }
00915         for(int i=0;i<bins.size();i++)
00916         {
00917             int b=bins[i];
00918             PPath f_name = vm->getMetaDataDir() / "bins"+tostring(b)+".def";
00919             PStream bfile = openFile(f_name, PStream::raw_ascii, "w");
00920             RealMapping rm;
00921             for(int j=0;j<sc.size();j++)
00922             {
00923                 rm = sc[j].getBinMapping(int(vm.length()/real(b)),int(vm.length()/real(b)));
00924                 bfile<<"DEFINE @"<<vm->fieldName(j)<<".ranges"+tostring(b)+" "<<rm<<endl;
00925                 bfile<<"DEFINE @"<<vm->fieldName(j)<<".ranges"+tostring(b)+".nbins "<<rm.size()<<endl;
00926                 bfile<<"DEFINE @"<<vm->fieldName(j)<<".ranges"+tostring(b)+".nbins_m1 "<<rm.size()-1<<endl;
00927                 bfile<<"DEFINE @"<<vm->fieldName(j)<<".binned"+tostring(b)+" @"<<vm->fieldName(j)<<" @"
00928                      <<vm->fieldName(j)<<".ranges"+tostring(b)<<endl;
00929                 bfile<<"DEFINE @"<<vm->fieldName(j)<<".onehot"+tostring(b)+" @"<<vm->fieldName(j)<<".binned"
00930                     +tostring(b)+" @"<<vm->fieldName(j)<<".ranges"+tostring(b)+".nbins onehot"<<endl;
00931 
00932             }
00933         }
00934     }
00935     else if(command=="genkfold")
00936     {
00937         if(argc<5)
00938         {
00939             cerr<<"usage vmat genkfold <source_dataset> <fileprefix> <kvalue>\n";
00940             exit(1);
00941         }
00942         string dbname = argv[2];
00943         string prefix = argv[3];
00944         int kval=toint(argv[4]);
00945         VMat vm = getVMat(dbname, indexf);
00946         for(int i=0;i<kval;i++)
00947         {
00948             ofstream out((prefix+"_train_"+tostring(i+1)+".vmat").c_str());
00949             out<<"<SOURCES>"<<endl;
00950             int ntest = vm.length()/kval;
00951             int ntrain_before_test = i*ntest;
00952             int ntrain_after_test = vm.length()-(i+1)*ntest;
00953             if(ntrain_before_test>0)
00954                 out<<dbname<<":0:"<<ntrain_before_test<<endl;
00955             if(ntrain_after_test>0)
00956                 out<<dbname<<":"<<ntest+ntrain_before_test<<":"<<ntrain_after_test<<endl;
00957             out<<"</SOURCES>"<<endl;
00958             ofstream out2((prefix+"_test_"+tostring(i+1)+".vmat").c_str());
00959             out2<<"<SOURCES>"<<endl;
00960             out2<<dbname<<":"<<ntrain_before_test<<":"<<ntest<<endl;
00961             out2<<"</SOURCES>"<<endl;
00962         }
00963     }
00964     else if(command=="genvmat")
00965     {
00966         if(argc<5)
00967         {
00968             cerr<<"usage vmat genvmat <source_dataset> <dest_vmat> (binned{num} | onehot{num} | normalized)\n";
00969             exit(1);
00970         }
00971         string dbname = argv[2];
00972         string destvmat = argv[3];
00973         string type=argv[4];
00974         int typen= 0;
00975         int bins= 0;
00976         if(type.find("binned")!=string::npos)
00977         {
00978             typen=0;
00979             bins=toint(type.substr(6));
00980         }
00981         else if(type.find("onehot")!=string::npos)
00982         {
00983             typen=1;
00984             bins=toint(type.substr(6));
00985         }
00986         else if(type.find("normalized")!=string::npos)
00987             typen=2;
00988         else PLERROR("Unknown operation: %s",type.c_str());
00989 
00990         VMat vm = getVMat(dbname, indexf);
00991         ofstream out(destvmat.c_str());
00992       
00993         out<<"<SOURCES>"<<endl;
00994         out<<dbname<<endl;
00995         out<<"</SOURCES>"<<endl;
00996         out<<"<PROCESSING>"<<endl;
00997         out<<"INCLUDE "<<dbname+".metadata/stats.def"<<endl;
00998         if(typen!=2)
00999             out<<"INCLUDE "<<dbname+".metadata/bins"<<bins<<".def"<<endl;
01000 
01001         for(int i=0;i<vm.width();i++)
01002         {
01003             switch(typen)
01004             {
01005             case 0:
01006                 out<<"@"<<vm->fieldName(i)<<".binned"<<bins<<endl;
01007                 out<<":"<<vm->fieldName(i)<<endl;
01008                 break;
01009             case 1:
01010                 out<<"@"<<vm->fieldName(i)<<".onehot"<<bins<<endl;
01011                 out<<":"<<vm->fieldName(i)<<".:0:@"<<vm->fieldName(i)<<".ranges"<<bins<<".nbins_m1"<<endl;
01012                 break;
01013             case 2:
01014                 out<<"@"<<vm->fieldName(i)<<".normalized"<<endl;
01015                 out<<":"<<vm->fieldName(i)<<endl;
01016                 break;
01017             }
01018 
01019         }
01020         out<<"</PROCESSING>"<<endl;
01021         out.close();
01022     }
01023     else if(command=="diststat")
01024     {
01025         VMat vm = getVMat(argv[2], indexf);
01026         int inputsize = atoi(argv[3]);
01027         printDistanceStatistics(vm, inputsize);      
01028     }
01029     else if(command=="diff")
01030     {
01031         if(argc < 4)
01032             PLERROR("'vmat diff' must be used that way : vmat diff <dataset1> <dataset2> [<tolerance> [<verbose>]]");
01033 
01034         VMat vm1 = getVMat(argv[2], indexf);
01035         VMat vm2 = getVMat(argv[3], indexf);
01036         double tol = 1e-6;
01037         int verb = 1;
01038         if(argc >= 5)
01039             tol = atof(argv[4]);
01040         if (argc >= 6)
01041             verb = atoi(argv[5]);
01042         print_diff(cout, vm1, vm2, tol, verb);      
01043     }
01044     else if(command=="cat")
01045     {
01046         if(argc < 3)
01047             PLERROR("'vmat cat' must be used that way : vmat cat FILE... [--precision=N] [vplFilteringCode]");
01048         string code;
01049         int nb_file=argc-2;
01050         int precision = -1;
01051         for (int i=argc-1 ; i >=3 && argv[i] ; i--) {
01052             string curopt = removeblanks(argv[i]);
01053             if(curopt.substr(0,12) == "--precision="){
01054                 precision = toint(curopt.substr(12));
01055                 nb_file--;
01056             }else if(!isfile(argv[argc-1])){
01057                 code=argv[argc-1];
01058                 nb_file--;
01059             }
01060         }
01061         if(precision>0){
01062             char tmpbuf[100];
01063             snprintf(tmpbuf,100,"%%#.%df",precision);
01064             pout.setDoubleFormat(tmpbuf);
01065             pout.setFloatFormat(tmpbuf);
01066         }
01067         for(int file=0;file<nb_file;file++)
01068         {
01069             string dbname=argv[file+2];
01070             if(nb_file>1)
01071                 pout<<dbname<<endl;
01072             VMat vm = getVMat(dbname, indexf);
01073             Vec tmp(vm.width());
01074             if(code.length()>0){
01075                 VMatLanguage vpl(vm);
01076                 vector<string> fn; 
01077                 for(int i=0;i<vm->width();i++)
01078                     fn.push_back(vm->fieldName(i));
01079                 vpl.compileString(code,fn);
01080                 Vec answer(1);
01081                 for(int i=0;i<vm.length();i++)
01082                 {
01083                     vpl.run(i,answer);
01084                     if(!fast_exact_is_equal(answer[0], 0)) {
01085                         vm->getRow(i, tmp);
01086                         pout<<tmp<<endl;
01087                     }
01088                 }
01089 
01090             }
01091             else
01092                 for(int i=0;i<vm.length();i++)
01093                 {
01094                     vm->getRow(i,tmp);      
01095                     pout<<tmp<<endl;
01096                 }
01097         }
01098     }
01099     else if(command=="catstr")
01100     {
01101         if(argc!=3 && argc != 4)
01102             PLERROR("'vmat catstr' must be used that way : vmat cat FILE [separator]");
01103         string dbname = argv[2];
01104         string sep = "\t";
01105         if(argc==4)
01106             sep = argv[3];
01107         VMat vm = getVMat(dbname, indexf);
01108         Vec tmp(vm.width());
01109         string out = "";
01110         for(int i=0;i<vm.length();i++)
01111         {
01112             vm->getRow(i,tmp);
01113             for(int j=0; j<vm.width(); j++)
01114             {
01115                 out = vm->getValString(j,tmp[j]);
01116                 if(out == "") out = tostring(tmp[j]);
01117                 pout << out << sep;
01118             }
01119             pout << endl;
01120         }
01121     }
01122     else if(command=="sascat")
01123     {
01124         if(argc!=4)
01125             PLERROR("'vmat sascat' must be used that way : vmat sascat <in-dataset> <out-filename.txt>");
01126         string dbname = argv[2];
01127         string outname = argv[3];
01128         string code;
01129         VMat vm = getVMat(dbname, indexf);
01130         ofstream out(outname.c_str());
01131         for (int i=0;i<vm.width();i++)
01132             out << vm->fieldName(i) << "\t";
01133         out << endl;
01134         for(int i=0;i<vm.length();i++)
01135         {
01136             for (int j=0;j<vm.width();j++)
01137                 out << vm->getString(i,j) << "\t";
01138             out<<endl;
01139         }
01140     }
01141     else if(command=="plot")
01142     {
01143         if(0 != argc%2)
01144             PLERROR("Bad number of arguments. Syntax for option plot:\n"
01145                     "%s plot <dbname0> <col0>[:<row0>:<nrows0>] {<dbnameN> <colN>[:<rowN>:<nrowsN>]}", argv[0]);
01146         plotVMats(argv+2, argc-2);
01147     }
01148     else if(command=="help")
01149     {
01150         pout << getDataSetHelp() << endl;
01151     }
01152     else if(command=="compare_stats")
01153     {
01154         if(!(argc==4||argc==5||argc==6))
01155             PLERROR("vmat compare_stats must be used that way: vmat compare_stats <dataset1> <dataset2> [[stderror threshold [missing threshold]]");
01156 
01157         VMat m1 = getVMat(argv[2], indexf);
01158         VMat m2 = getVMat(argv[3], indexf);
01159 
01160         m1->compatibleSizeError(m2);
01161 
01162         real stderror_threshold = 1;
01163         real missing_threshold = 10;
01164         if(argc>4)
01165             stderror_threshold=toreal(argv[4]);
01166         if(argc>5)
01167             missing_threshold=toreal(argv[5]);
01168         Vec missing(m1->width());
01169         Vec stderror(m1->width());
01170 
01171         pout << "Test of difference that suppose gaussiane variable"<<endl;
01172         m1->compareStats(m2, stderror_threshold, missing_threshold,
01173                          stderror, missing);
01174 
01175         Mat score(m1->width(),3);
01176 
01177         for(int col = 0;col<m1->width();col++)
01178         {
01179             score(col,0)=col;
01180             score(col,1)=stderror[col];
01181             score(col,2)=missing[col];
01182         }
01183         
01184         int nbdiff = 0;
01185 
01186         pout<<"Print the field that do not pass the threshold sorted by the stderror"<<endl;
01187         sortRows(score,1,false);
01188         for(int i=0;i<score.length();i++)
01189         {
01190             if(score(i,1)>stderror_threshold)
01191             {
01192                 const StatsCollector tstats = m1->getStats(i);
01193                 const StatsCollector lstats = m2->getStats(i);
01194                 real tmean = tstats.mean();
01195                 real lmean = lstats.mean();
01196                 real tstderror = sqrt(pow(tstats.stderror(), 2) + 
01197                                       pow(lstats.stderror(), 2));
01198 
01199                 pout<<i<<"("<<m1->fieldName(int(round(score(i,0))))<<")"
01200                     <<" differ by "<<score(i,1)<<" stderror."
01201                     <<" The mean is "<<lmean<<" while the target mean is "<<tmean
01202                     <<" and the used stderror is "<<tstderror<<endl;
01203                 nbdiff++;
01204             }
01205         }
01206 
01207         cout<<"Print the field that do not pass the threshold sorted by the missing error"<<endl;
01208         sortRows(score,2,false);
01209         for(int i=0;i<score.length();i++)
01210         {
01211             if(score(i,2)>missing_threshold)
01212             {
01213                 const StatsCollector tstats = m1->getStats(i);
01214                 const StatsCollector lstats = m2->getStats(i);
01215                 real tmissing = tstats.nmissing()/tstats.n();
01216                 real lmissing = lstats.nmissing()/lstats.n();
01217                 pout<<i<<"("<<m1->fieldName(int(round(score(i,0))))<<")"
01218                     <<" The missing stats difference is "<< score(i,2)
01219                     <<". There are "<<lmissing<<" missing while target has "
01220                     <<tmissing<<" missing."<<endl;
01221                 nbdiff++;
01222             }
01223         }
01224 
01225         pout<<"There are "<<nbdiff<<"/"<<m1.width()
01226             <<" fields that have different stats"<<endl;
01227 
01228     }
01229     else if(command=="compare_stats_ks")
01230     {
01231         bool err = false;
01232         real threshold = REAL_MAX;
01233         bool mat_to_mem = false;
01234         if(argc<4||argc>6)
01235             err = true;
01236         if(argc==5)
01237         {
01238             if(argv[4]==string("--mat_to_mem"))
01239                 mat_to_mem=true;
01240             else if(!pl_isnumber(string(argv[4]),&threshold))
01241                 err = true;
01242         }
01243         else if(argc==6)
01244         {
01245              if(argv[5]!=string("--mat_to_mem"))
01246                  err = true;
01247              else if(!pl_isnumber(string(argv[4]),&threshold))
01248                  err = true;
01249         }
01250         if(err)
01251             PLERROR("vmat compare_stats_ks must be used that way:"
01252                     " vmat compare_stats_ks <dataset1> <dataset2> [threshold]"
01253                     " [--mat_to_mem]");
01254 
01255         VMat m1 = getVMat(argv[2], indexf);
01256         VMat m2 = getVMat(argv[3], indexf);
01257         if(mat_to_mem)
01258         {
01259             m1.precompute();
01260             m2.precompute();
01261         }
01262 
01263         m1->compatibleSizeError(m2);
01264         int pc_value_99=0;
01265         int pc_value_95=0;
01266         int pc_value_90=0;
01267         int pc_value_0=0;
01268 
01269         uint size_fieldnames=m1->maxFieldNamesSize();
01270 
01271         Vec Ds(m1->width());
01272         Vec p_values(m1->width());
01273         KS_test(m1,m2,10,Ds,p_values,true);
01274         Mat score(m1->width(),3);
01275             
01276         for(int col = 0;col<m1->width();col++)
01277         {
01278             score(col,0)=col;
01279             score(col,1)=Ds[col];
01280             real p_value = p_values[col];
01281             score(col,2)=p_value;
01282             if(p_value>0.99)
01283                 pc_value_99++;
01284             if(p_value>0.95)
01285                 pc_value_95++;
01286             if(p_value>0.90)
01287                 pc_value_90++;
01288             else
01289                 pc_value_0++;
01290         }
01291 
01292         sortRows(score,2,false);
01293         pout <<"Kolmogorov Smirnov two sample test"<<endl<<endl;
01294         if(threshold<REAL_MAX)
01295             pout<<"Variables that are under the threshold"<<endl;
01296         pout<<"Sorted by p_value"<<endl;
01297         cout << std::left << setw(8) << "# "
01298              << setw(size_fieldnames) << " fieldname " << std::right
01299              << setw(15) << " D"
01300              << setw(15) << " p_value"
01301              <<endl;
01302         int threshold_fail=0;
01303         for(int col=0;col<score.length();col++)
01304         {
01305             if(threshold>=score(col,2))
01306             {
01307                 cout << std::left << setw(8) << tostring(col)+"/"+tostring(score(col,0))
01308                      << setw(size_fieldnames) << m1->fieldName(int(round(score(col,0))))
01309                      << std::right
01310                      << setw(15) << score(col,1)
01311                      << setw(15) << score(col,2)
01312                      <<endl;
01313                 threshold_fail++;
01314             }
01315         }
01316         if(threshold<REAL_MAX)
01317             pout << "There are "<<threshold_fail<<" variables that are under the threshold"<<endl;
01318         if(threshold==REAL_MAX)
01319         {
01320             pout << "99% cutoff: "<<pc_value_99<<endl;
01321             pout << "95% cutoff: "<<pc_value_95<<endl;
01322             pout << "90% cutoff: "<<pc_value_90<<endl;
01323             pout << "0-90% cutoff: "<<pc_value_0<<endl;
01324         }
01325         pout <<"Kolmogorov Smirnov two sample test end"<<endl<<endl;
01326     }
01327     else if(command=="compare_stats_desjardins")
01328     {      
01329         bool err=false;
01330         bool mat_to_mem=false;
01331         if(!(argc==8||argc==9))
01332             err=true;
01333         if(argc==9)
01334         {
01335             if(argv[8]!=string("--mat_to_mem"))
01336                  err = true;
01337             else
01338                 mat_to_mem=true;
01339         }
01340         if(err)
01341             PLERROR("vmat compare_stats_desjardins must be used that way:"
01342                     " vmat compare_stats_desjardins <orig dataset1> <orig dataset2> <new dataset3> <ks_threshold> <stderror_threshold> <missing_threshold> [--mat_to_mem]");
01343 
01344         VMat m1 = getVMat(argv[2], indexf);
01345         VMat m2 = getVMat(argv[3], indexf);
01346         VMat m3 = getVMat(argv[4], indexf);
01347         real ks_threshold = toreal(argv[5]);
01348 
01349         m3->compatibleSizeError(m1);
01350         m3->compatibleSizeError(m2);
01351 
01352         Vec Ds(m1->width());
01353         Vec p_values(m1->width());
01354         Mat score(m1->width(),3);
01355         uint size_fieldnames=m1->maxFieldNamesSize();
01356         if(mat_to_mem==true)
01357         {
01358             m1.precompute();
01359             m2.precompute();
01360             m3.precompute();
01361         }
01362 
01363         KS_test(m1,m3,10,Ds,p_values,true);
01364         for(int col = 0;col<m1->width();col++)
01365         {
01366             score(col,0)=col;
01367             score(col,1)=Ds[col];
01368             real p_value = p_values[col];
01369             score(col,2)=p_value;
01370         }
01371 
01372         KS_test(m2,m3,10,Ds,p_values,true);
01373         for(int col = 0;col<m1->width();col++)
01374         {
01375             if(p_values[col]>score(col,2))
01376             {
01377                 score(col,1)=Ds[col];
01378                 score(col,2)=p_values[col];
01379             }
01380         }
01381 
01382         sortRows(score,2,false);
01383         pout <<"Kolmogorov Smirnov two sample test"<<endl<<endl;
01384         pout<<"Variables that are under the ks_threshold"<<endl;
01385         pout<<"Sorted by p_value"<<endl;
01386         cout << std::left << setw(8) << "# "
01387              << setw(size_fieldnames) << " fieldname " << std::right
01388              << setw(15) << " D"
01389              << setw(15) << " p_value"
01390              <<endl;
01391         int threshold_fail = 0;
01392         for(int col=0;col<score.length();col++)
01393         {
01394             if(ks_threshold>=score(col,2))
01395             {
01396                 cout << std::left << setw(8) << tostring(col)+"/"+tostring(score(col,0))
01397                      << setw(size_fieldnames) << m1->fieldName(int(round(score(col,0))))
01398                      << std::right
01399                      << setw(15) << score(col,1)
01400                      << setw(15) << score(col,2)
01401                      <<endl;
01402                 threshold_fail++;
01403             }
01404         }
01405         pout << "There are "<<threshold_fail<<"/"<<m1->width()<<
01406             " variables that are under the threshold"<<endl<<
01407             " Kolmogorov Smirnov two sample test end"<<endl<<endl;
01408 
01409 
01410 //         real stderror_threshold = 1;
01411 //         real missing_threshold = 10;
01412 //         stderror_threshold=toreal(argv[6]);
01413 //         missing_threshold=toreal(argv[7]);
01414 
01415 //         pout << "Test of difference that suppose gaussiane variable"<<endl;
01416 //         pout << "Comparing with dataset1"<<endl;
01417 
01418 //         m3->compareStats(m1, stderror_threshold, missing_threshold,
01419 //                          stderr, missing);
01420 //         pout << "Comparing with dataset2"<<endl;
01421 //         m3->compareStats(m2, stderror_threshold, missing_threshold,
01422 //                          stderr, missing);
01423 //         pout<<"There are "<<diff<<"/"<<m1.width()
01424 //             <<" fields that have different stats"<<endl;
01425     }
01426     else if(command=="characterize")
01427     {
01428         if(argc!=3)
01429             PLERROR("The command 'vmat characterize' must be used that way: vmat caracterize <dataset1>");
01430         VMat m1 = getVMat(argv[2], indexf);
01431         TVec<StatsCollector> stats = 
01432             m1->getStats();//"stats_all.psave",-1,true);
01433         TVec<string> caracs;
01434         uint size_fieldnames=m1->maxFieldNamesSize();
01435 
01436         for(int i=0;i<stats.size();i++)
01437         {
01438             StatsCollector& stat=stats[i];
01439             string carac = tostring(i)+"\t"+left(m1->fieldName(i),size_fieldnames);
01440             if(stat.isbinary())
01441                 carac += "\tBinary";
01442             else if(stat.isinteger())
01443                 carac+="\tInteger";
01444             else
01445                 carac+="\tReal";
01446 
01447             //find is normal or not
01448             int m=min(100,int(round(sqrt(stat.nnonmissing()))));
01449             int nelem=int(stat.nnonmissing()/m);
01450             int row=0;
01451             real gsum=0;
01452             for(int bloc=0;bloc<m;bloc++)
01453             {
01454                 real bloc_sum=0;
01455                 for(int bloc_elem=0;bloc_elem<nelem;)
01456                 {
01457                     real v = m1->get(row,i);
01458                     if(is_missing(v))
01459                         continue;
01460                     else
01461                     {
01462                         bloc_elem++;
01463                         row++;
01464                         bloc_sum+=v;
01465                     }
01466                     
01467                 }
01468                 gsum+=pow((bloc_sum/nelem)-stat.mean(),2);
01469             }
01470             real s2=(stat.variance()/nelem);
01471             real mu_square = stat.mean()*stat.mean();
01472             real th = ((1./(m-1))*gsum)/s2;
01473             real th2= mu_square+s2-12*mu_square*s2+mu_square*mu_square
01474                 * s2*s2;
01475             bool b = th>(s2+2*th2)/s2;
01476                 carac+="\tnormal test value: "+tostring(th)+" "+tostring(th2)
01477                 + " "+tostring(b);
01478             caracs.append(carac);
01479             pout<<carac<<endl;
01480         }
01481 
01482     }
01483     else if(command=="mtime")
01484     {    
01485         if(argc!=3)
01486             PLERROR("The command 'vmat mtime' must be used that way: vmat mtime <dataset>");
01487         VMat m1 = getVMat(argv[2], indexf);
01488         pout<<m1->getMtime()<<endl;
01489     }
01490     else
01491         PLERROR("Unknown command : %s",command.c_str());
01492     return 0;
01493 }
01494 
01495 } // end of namespace PLearn
01496 
01497 
01498 /*
01499   Local Variables:
01500   mode:c++
01501   c-basic-offset:4
01502   c-file-style:"stroustrup"
01503   c-file-offsets:((innamespace . 0)(inline-open . 0))
01504   indent-tabs-mode:nil
01505   fill-column:79
01506   End:
01507 */
01508 // 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