kruskalwalliscommand.cpp

   1 /*
   2  * File:   kruskalwalliscommand.cpp
   3  * Author: kiverson
   4  *
   5  * Created on June 26, 2012, 11:06 AM
   6  */
   7 #include "kruskalwalliscommand.h"
   8
   9 //**********************************************************************************************************************
  10 vector<string> KruskalWallisCommand::setParameters(){
  11         try {
  12                 CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir);
  13                 CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir);
  14
  15                 vector<string> myArray;
  16                 for (int i = 0; i < parameters.size(); i++) {   myArray.push_back(parameters[i].name);          }
  17                 return myArray;
  18         }
  19         catch(exception& e) {
  20                 m->errorOut(e, "KruskalWallisCommand", "setParameters");
  21                 exit(1);
  22         }
  23 }
  24 //**********************************************************************************************************************
  25 string KruskalWallisCommand::getHelpString(){
  26         try {
  27                 string helpString = "";
  28                 helpString += "The kruskalwallis command parameter options are \n";
  29         helpString += "Kruskal–Wallis one-way analysis of variance is a non-parametric method for testing whether samples originate from the same distribution.";
  30                 return helpString;
  31         }
  32         catch(exception& e) {
  33                 m->errorOut(e, "KruskalWallisCommand", "getHelpString");
  34                 exit(1);
  35         }
  36 }
  37 //**********************************************************************************************************************
  38 string KruskalWallisCommand::getOutputFileNameTag(string type, string inputName=""){
  39         try {
  40         string outputFileName = "";
  41                 map<string, vector<string> >::iterator it;
  42
  43         //is this a type this command creates
  44         it = outputTypes.find(type);
  45         if (it == outputTypes.end()) {  m->mothurOut("[ERROR]: this command doesn't create a " + type + " output file.\n"); }
  46         else {
  47             if (type == "summary") {  outputFileName =  "cooccurence.summary"; }
  48             else { m->mothurOut("[ERROR]: No definition for type " + type + " output file tag.\n"); m->control_pressed = true;  }
  49         }
  50         return outputFileName;
  51         }
  52         catch(exception& e) {
  53                 m->errorOut(e, "KruskalWallisCommand", "getOutputFileNameTag");
  54                 exit(1);
  55         }
  56 }
  57 //**********************************************************************************************************************
  58 KruskalWallisCommand::KruskalWallisCommand(){
  59         try {
  60                 abort = true; calledHelp = true;
  61                 setParameters();
  62         vector<string> tempOutNames;
  63                 outputTypes["summary"] = tempOutNames;
  64
  65         }
  66         catch(exception& e) {
  67                 m->errorOut(e, "KruskalWallisCommand", "KruskalWallisCommand");
  68                 exit(1);
  69         }
  70 }
  71 //**********************************************************************************************************************
  72 KruskalWallisCommand::KruskalWallisCommand(string option) {
  73         try {
  74                 abort = false; calledHelp = false;
  75
  76                 //allow user to run help
  77                 if(option == "help") { help(); abort = true; calledHelp = true; }
  78                 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
  79
  80                 else {
  81                         vector<string> myArray = setParameters();
  82
  83                         OptionParser parser(option);
  84                         map<string,string> parameters = parser.getParameters();
  85                         map<string,string>::iterator it;
  86
  87                         ValidParameters validParameter;
  88
  89                         //check to make sure all parameters are valid for command
  90                         for (it = parameters.begin(); it != parameters.end(); it++) {
  91                                 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
  92                         }
  93
  94
  95                         //if the user changes the input directory command factory will send this info to us in the output parameter
  96                         string inputDir = validParameter.validFile(parameters, "inputdir", false);
  97                         if (inputDir == "not found"){   inputDir = "";          }
  98                         else {
  99                                 string path;
 100                                 it = parameters.find("shared");
 101                                 //user has given a template file
 102                                 if(it != parameters.end()){
 103                                         path = m->hasPath(it->second);
 104                                         //if the user has not given a path then, add inputdir. else leave path alone.
 105                                         if (path == "") {       parameters["shared"] = inputDir + it->second;           }
 106                                 }
 107                         }
 108
 109             vector<string> tempOutNames;
 110             outputTypes["summary"] = tempOutNames;
 111
 112
 113                 }
 114
 115         }
 116         catch(exception& e) {
 117                 m->errorOut(e, "KruskalWallisCommand", "KruskalWallisCommand");
 118                 exit(1);
 119         }
 120 }
 121 //**********************************************************************************************************************
 122 int KruskalWallisCommand::execute(){
 123         try {
 124                 if (abort == true) { if (calledHelp) { return 0; }  return 2;   }
 125
 126         //math goes here
 127
 128         int N; //= thisLookUp.size();
 129         double H;
 130         double tmp = 0.0;
 131         vector<groupRank> vec;
 132         string group;
 133         int count;
 134         double sum;
 135
 136         //merge all groups into a vector
 137         //rank function here
 138         assignRank(vec);
 139
 140         //populate counts and ranSums vectors
 141         for (int i=0;i<N;i++) {
 142             count = 0;
 143             sum = 0;
 144             //group = next group
 145             for(int j;j<vec.size();j++) {
 146                 if (vec[j].group == group) {
 147                     count++;
 148                     sum = sum + vec[j].rank;
 149                 }
 150             }
 151             counts[i] = count;
 152             rankSums[i] = sum;
 153         }
 154
 155         //test statistic
 156         for (int i=0;i<N;i++) { tmp = tmp + (pow(rankSums[i],2) / counts[i]); }
 157
 158         H = (12 / (N*(N+1))) * tmp - (3*(N+1));
 159
 160         //ss = tmp - pow(accumulate(rankSums.begin(), rankSums.end(), 0), 2);
 161
 162         //H = ss / ( (N * (N + 1))/12 );
 163
 164         //correction for ties?
 165
 166         //p-value calculation
 167
 168                 return 0;
 169         }
 170         catch(exception& e) {
 171                 m->errorOut(e, "KruskalWallisCommand", "execute");
 172                 exit(1);
 173         }
 174 }
 175 //**********************************************************************************************************************
 176 void KruskalWallisCommand::assignRank(vector<groupRank> &vec) {
 177     try {
 178         double rank = 1;
 179         double numRanks, avgRank, j;
 180         vector<groupRank>::iterator it, oldit;
 181
 182         sort (vec.begin(), vec.end(), comparevalue);
 183
 184         it = vec.begin();
 185
 186         while ( it != vec.end() ) {
 187             j = rank;
 188             oldit = it;
 189             if (!equalvalue(*it, *it+1)) { *it->rank = rank; rank=rank+1; it++; }
 190             else {
 191                 while(equalrank(*it, *it+1)) {
 192                     j = j + (j+1.0);
 193                     rank++;
 194                     it++;
 195                 }
 196                 numRanks = double (distance(oldit, it));
 197                 avgRank = j / numRanks;
 198                 while(oldit != it) {
 199                     *oldit->rank = avgRank;
 200                     oldit++;
 201                 }
 202             }
 203
 204         }
 205
 206     }
 207     catch(exception& e) {
 208                 m->errorOut(e, "KruskalWallisCommand", "getRank");
 209                 exit(1);
 210         }
 211
 212 }
 213 //**********************************************************************************************************************
 214
 215 //**********************************************************************************************************************
 216 //**********************************************************************************************************************
 217 //**********************************************************************************************************************