PLearn 0.1
|
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 :