pcoacommand.cpp

   1
   2 /*
   3  *  pcacommand.cpp
   4  *  Mothur
   5  *
   6  *  Created by westcott on 1/4/10.
   7  *  Copyright 2010 Schloss Lab. All rights reserved.
   8  *
   9  */
  10
  11 #include "pcoacommand.h"
  12
  13 //**********************************************************************************************************************
  14 vector<string> PCOACommand::getValidParameters(){
  15         try {
  16                 string Array[] =  {"phylip", "metric","outputdir","inputdir"};
  17                 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
  18                 return myArray;
  19         }
  20         catch(exception& e) {
  21                 m->errorOut(e, "PCOACommand", "getValidParameters");
  22                 exit(1);
  23         }
  24 }
  25 //**********************************************************************************************************************
  26 PCOACommand::PCOACommand(){
  27         try {
  28                 abort = true;
  29                 //initialize outputTypes
  30                 vector<string> tempOutNames;
  31                 outputTypes["pcoa"] = tempOutNames;
  32                 outputTypes["loadings"] = tempOutNames;
  33                 outputTypes["corr"] = tempOutNames;
  34         }
  35         catch(exception& e) {
  36                 m->errorOut(e, "PCOACommand", "PCOACommand");
  37                 exit(1);
  38         }
  39 }
  40 //**********************************************************************************************************************
  41 vector<string> PCOACommand::getRequiredParameters(){
  42         try {
  43                 string Array[] =  {"phylip"};
  44                 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
  45                 return myArray;
  46         }
  47         catch(exception& e) {
  48                 m->errorOut(e, "PCOACommand", "getRequiredParameters");
  49                 exit(1);
  50         }
  51 }
  52 //**********************************************************************************************************************
  53 vector<string> PCOACommand::getRequiredFiles(){
  54         try {
  55                 vector<string> myArray;
  56                 return myArray;
  57         }
  58         catch(exception& e) {
  59                 m->errorOut(e, "PCOACommand", "getRequiredFiles");
  60                 exit(1);
  61         }
  62 }
  63 //**********************************************************************************************************************
  64
  65 PCOACommand::PCOACommand(string option)  {
  66         try {
  67                 abort = false;
  68
  69                 //allow user to run help
  70                 if(option == "help") { help(); abort = true; }
  71
  72                 else {
  73                         //valid paramters for this command
  74                         string Array[] =  {"phylip","metric","outputdir", "inputdir"};
  75                         vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
  76
  77                         OptionParser parser(option);
  78                         map<string, string> parameters = parser. getParameters();
  79
  80                         ValidParameters validParameter;
  81                         map<string, string>::iterator it;
  82
  83                         //check to make sure all parameters are valid for command
  84                         for (it = parameters.begin(); it != parameters.end(); it++) {
  85                                 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
  86                         }
  87                         //if the user changes the input directory command factory will send this info to us in the output parameter
  88                         string inputDir = validParameter.validFile(parameters, "inputdir", false);
  89                         if (inputDir == "not found"){   inputDir = "";          }
  90                         else {
  91                                 string path;
  92                                 it = parameters.find("phylip");
  93                                 //user has given a template file
  94                                 if(it != parameters.end()){
  95                                         path = m->hasPath(it->second);
  96                                         //if the user has not given a path then, add inputdir. else leave path alone.
  97                                         if (path == "") {       parameters["phylip"] = inputDir + it->second;           }
  98                                 }
  99                         }
 100
 101                         //initialize outputTypes
 102                         vector<string> tempOutNames;
 103                         outputTypes["pcoa"] = tempOutNames;
 104                         outputTypes["loadings"] = tempOutNames;
 105                         outputTypes["corr"] = tempOutNames;
 106
 107                         //required parameters
 108                         phylipfile = validParameter.validFile(parameters, "phylip", true);
 109                         if (phylipfile == "not open") { abort = true; }
 110                         else if (phylipfile == "not found") { phylipfile = ""; abort = true; }
 111                         else {  filename = phylipfile;  }
 112
 113                         //if the user changes the output directory command factory will send this info to us in the output parameter
 114                         outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){
 115                                 outputDir = "";
 116                                 outputDir += m->hasPath(phylipfile); //if user entered a file with a path then preserve it
 117                         }
 118
 119                         //error checking on files
 120                         if (phylipfile == "")   { m->mothurOut("You must provide a distance file before running the pcoa command."); m->mothurOutEndLine(); abort = true; }
 121
 122                         string temp = validParameter.validFile(parameters, "metric", false);    if (temp == "not found"){       temp = "T";                             }
 123                         metric = m->isTrue(temp);
 124                 }
 125
 126         }
 127         catch(exception& e) {
 128                 m->errorOut(e, "PCOACommand", "PCOACommand");
 129                 exit(1);
 130         }
 131 }
 132 //**********************************************************************************************************************
 133 void PCOACommand::help(){
 134         try {
 135
 136                 m->mothurOut("The pcoa command parameters are phylip and metric"); m->mothurOutEndLine();
 137                 m->mothurOut("The phylip parameter allows you to enter your distance file."); m->mothurOutEndLine();
 138                 m->mothurOut("The metric parameter allows indicate you if would like the pearson correlation coefficient calculated. Default=True"); m->mothurOutEndLine();
 139                 m->mothurOut("Example pcoa(phylip=yourDistanceFile).\n");
 140                 m->mothurOut("Note: No spaces between parameter labels (i.e. phylip), '=' and parameters (i.e.yourDistanceFile).\n\n");
 141         }
 142         catch(exception& e) {
 143                 m->errorOut(e, "PCOACommand", "help");
 144                 exit(1);
 145         }
 146 }
 147 //**********************************************************************************************************************
 148 PCOACommand::~PCOACommand(){}
 149 //**********************************************************************************************************************
 150 int PCOACommand::execute(){
 151         try {
 152
 153                 if (abort == true) { return 0; }
 154
 155                 cout.setf(ios::fixed, ios::floatfield);
 156                 cout.setf(ios::showpoint);
 157                 cerr.setf(ios::fixed, ios::floatfield);
 158                 cerr.setf(ios::showpoint);
 159
 160                 vector<string> names;
 161                 vector<vector<double> > D;
 162
 163                 fbase = outputDir + m->getRootName(m->getSimpleName(filename));
 164
 165                 read(filename, names, D);
 166
 167                 if (m->control_pressed) { return 0; }
 168
 169                 double offset = 0.0000;
 170                 vector<double> d;
 171                 vector<double> e;
 172                 vector<vector<double> > G = D;
 173                 vector<vector<double> > copy_G;
 174                 //int rank = D.size();
 175
 176                 m->mothurOut("\nProcessing...\n\n");
 177
 178                 for(int count=0;count<2;count++){
 179                         recenter(offset, D, G);         if (m->control_pressed) { return 0; }
 180                         tred2(G, d, e);                         if (m->control_pressed) { return 0; }
 181                         qtli(d, e, G);                          if (m->control_pressed) { return 0; }
 182                         offset = d[d.size()-1];
 183                         if(offset > 0.0) break;
 184                 }
 185
 186                 if (m->control_pressed) { return 0; }
 187
 188                 output(fbase, names, G, d);
 189
 190                 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) {        remove(outputNames[i].c_str());  } return 0; }
 191
 192                 if (metric) {
 193
 194                         for (int i = 1; i < 4; i++) {
 195
 196                                 vector< vector<double> > EuclidDists = calculateEuclidianDistance(G, i); //G is the pcoa file
 197
 198                                 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) {        remove(outputNames[i].c_str());  } return 0; }
 199
 200                                 double corr = calcPearson(EuclidDists, D); //G is the pcoa file, D is the users distance matrix
 201
 202                                 m->mothurOut("Pearson's coefficient using " + toString(i) + " axis: " + toString(corr)); m->mothurOutEndLine();
 203
 204                                 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) {        remove(outputNames[i].c_str());  } return 0; }
 205                         }
 206                 }
 207
 208                 m->mothurOutEndLine();
 209                 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
 210                 for (int i = 0; i < outputNames.size(); i++) {  m->mothurOut(outputNames[i]); m->mothurOutEndLine();    }
 211                 m->mothurOutEndLine();
 212
 213                 return 0;
 214         }
 215         catch(exception& e) {
 216                 m->errorOut(e, "PCOACommand", "execute");
 217                 exit(1);
 218         }
 219 }
 220 /*********************************************************************************************************************************/
 221 vector< vector<double> > PCOACommand::calculateEuclidianDistance(vector< vector<double> >& axes, int dimensions){
 222         try {
 223                 //make square matrix
 224                 vector< vector<double> > dists; dists.resize(axes.size());
 225                 for (int i = 0; i < dists.size(); i++) {  dists[i].resize(axes.size(), 0.0); }
 226
 227                 if (dimensions == 1) { //one dimension calc = abs(x-y)
 228
 229                         for (int i = 0; i < dists.size(); i++) {
 230
 231                                 if (m->control_pressed) { return dists; }
 232
 233                                 for (int j = 0; j < i; j++) {
 234                                         dists[i][j] = abs(axes[i][0] - axes[j][0]);
 235                                         dists[j][i] = dists[i][j];
 236                                 }
 237                         }
 238
 239                 }else if (dimensions == 2) { //two dimension calc = sqrt ((x1 - y1)^2 + (x2 - y2)^2)
 240
 241                         for (int i = 0; i < dists.size(); i++) {
 242
 243                                 if (m->control_pressed) { return dists; }
 244
 245                                 for (int j = 0; j < i; j++) {
 246                                         double firstDim = ((axes[i][0] - axes[j][0]) * (axes[i][0] - axes[j][0]));
 247                                         double secondDim = ((axes[i][1] - axes[j][1]) * (axes[i][1] - axes[j][1]));
 248
 249                                         dists[i][j] = sqrt((firstDim + secondDim));
 250                                         dists[j][i] = dists[i][j];
 251                                 }
 252                         }
 253
 254                 }else if (dimensions == 3) { //two dimension calc = sqrt ((x1 - y1)^2 + (x2 - y2)^2 + (x3 - y3)^2)
 255
 256                         for (int i = 0; i < dists.size(); i++) {
 257
 258                                 if (m->control_pressed) { return dists; }
 259
 260                                 for (int j = 0; j < i; j++) {
 261                                         double firstDim = ((axes[i][0] - axes[j][0]) * (axes[i][0] - axes[j][0]));
 262                                         double secondDim = ((axes[i][1] - axes[j][1]) * (axes[i][1] - axes[j][1]));
 263                                         double thirdDim = ((axes[i][2] - axes[j][2]) * (axes[i][2] - axes[j][2]));
 264
 265                                         dists[i][j] = sqrt((firstDim + secondDim + thirdDim));
 266                                         dists[j][i] = dists[i][j];
 267                                 }
 268                         }
 269
 270                 }else { m->mothurOut("[ERROR]: too many dimensions, aborting."); m->mothurOutEndLine(); m->control_pressed = true; }
 271
 272                 return dists;
 273         }
 274         catch(exception& e) {
 275                 m->errorOut(e, "PCOACommand", "calculateEuclidianDistance");
 276                 exit(1);
 277         }
 278 }
 279 /*********************************************************************************************************************************/
 280 double PCOACommand::calcPearson(vector< vector<double> >& euclidDists, vector< vector<double> >& userDists){
 281         try {
 282
 283                 //find average for - X
 284                 vector<float> averageEuclid; averageEuclid.resize(euclidDists.size(), 0.0);
 285                 for (int i = 0; i < euclidDists.size(); i++) {
 286                         for (int j = 0; j < euclidDists[i].size(); j++) {
 287                                 averageEuclid[i] += euclidDists[i][j];
 288                         }
 289                 }
 290                 for (int i = 0; i < averageEuclid.size(); i++) {  averageEuclid[i] = averageEuclid[i] / (float) euclidDists.size();   }
 291
 292                 //find average for - Y
 293                 vector<float> averageUser; averageUser.resize(userDists.size(), 0.0);
 294                 for (int i = 0; i < userDists.size(); i++) {
 295                         for (int j = 0; j < userDists[i].size(); j++) {
 296                                 averageUser[i] += userDists[i][j];
 297                         }
 298                 }
 299                 for (int i = 0; i < averageUser.size(); i++) {  averageUser[i] = averageUser[i] / (float) userDists.size();  }
 300
 301                 double numerator = 0.0;
 302                 double denomTerm1 = 0.0;
 303                 double denomTerm2 = 0.0;
 304
 305                 for (int i = 0; i < euclidDists.size(); i++) {
 306
 307                         for (int k = 0; k < i; k++) {
 308
 309                                 float Yi = userDists[i][k];
 310                                 float Xi = euclidDists[i][k];
 311
 312                                 numerator += ((Xi - averageEuclid[k]) * (Yi - averageUser[k]));
 313                                 denomTerm1 += ((Xi - averageEuclid[k]) * (Xi - averageEuclid[k]));
 314                                 denomTerm2 += ((Yi - averageUser[k]) * (Yi - averageUser[k]));
 315                         }
 316                 }
 317
 318                 double denom = (sqrt(denomTerm1) * sqrt(denomTerm2));
 319                 double r = numerator / denom;
 320
 321                 return r;
 322         }
 323         catch(exception& e) {
 324                 m->errorOut(e, "PCOACommand", "calculateEuclidianDistance");
 325                 exit(1);
 326         }
 327 }
 328 /*********************************************************************************************************************************/
 329
 330 inline double SIGN(const double a, const double b)
 331 {
 332     return b>=0 ? (a>=0 ? a:-a) : (a>=0 ? -a:a);
 333 }
 334
 335 /*********************************************************************************************************************************/
 336
 337 void PCOACommand::get_comment(istream& f, char begin, char end){
 338         try {
 339                 char d=f.get();
 340                 while(d != end){        d = f.get();    }
 341                 d = f.peek();
 342         }
 343         catch(exception& e) {
 344                 m->errorOut(e, "PCOACommand", "get_comment");
 345                 exit(1);
 346         }
 347 }
 348
 349 /*********************************************************************************************************************************/
 350
 351 int PCOACommand::read_phylip(istream& f, int square_m, vector<string>& name_list, vector<vector<double> >& d){
 352         try {
 353                 //     int count1=0;
 354                 //     int count2=0;
 355
 356                 int rank;
 357                 f >> rank;
 358
 359                 name_list.resize(rank);
 360                 d.resize(rank);
 361                 if(square_m == 1){
 362                         for(int i=0;i<rank;i++)
 363                                 d[i].resize(rank);
 364                         for(int i=0;i<rank;i++) {
 365                                 f >> name_list[i];
 366                                 //                      cout << i << "\t" << name_list[i] << endl;
 367                                 for(int j=0;j<rank;j++) {
 368                                         if (m->control_pressed) { return 0; }
 369
 370                                         f >> d[i][j];
 371                                         if (d[i][j] == -0.0000)
 372                                                 d[i][j] = 0.0000;
 373                                 }
 374                         }
 375                 }
 376                 else if(square_m == 2){
 377                         for(int i=0;i<rank;i++){
 378                                 d[i].resize(rank);
 379                         }
 380                         d[0][0] = 0.0000;
 381                         f >> name_list[0];
 382                         for(int i=1;i<rank;i++){
 383                                 f >> name_list[i];
 384                                 d[i][i]=0.0000;
 385                                 for(int j=0;j<i;j++){
 386                                         if (m->control_pressed) { return 0; }
 387                                         f >> d[i][j];
 388                                         if (d[i][j] == -0.0000)
 389                                                 d[i][j] = 0.0000;
 390                                         d[j][i]=d[i][j];
 391                                 }
 392                         }
 393                 }
 394
 395                 return 0;
 396         }
 397         catch(exception& e) {
 398                 m->errorOut(e, "PCOACommand", "read_phylip");
 399                 exit(1);
 400         }
 401
 402 }
 403
 404 /*********************************************************************************************************************************/
 405
 406 void PCOACommand::read(string fname, vector<string>& names, vector<vector<double> >& D){
 407         try {
 408                 ifstream f;
 409                 m->openInputFile(fname, f);
 410
 411                 //check whether matrix is square
 412                 char d;
 413                 int q = 1;
 414                 int numSeqs;
 415                 string name;
 416
 417                 f >> numSeqs >> name;
 418
 419                 while((d=f.get()) != EOF){
 420
 421                         //is d a number meaning its square
 422                         if(isalnum(d)){
 423                                 q = 1;
 424                                 break;
 425                         }
 426
 427                         //is d a line return meaning its lower triangle
 428                         if(d == '\n'){
 429                                 q = 2;
 430                                 break;
 431                         }
 432                 }
 433                 f.close();
 434
 435                 //reopen to get back to beginning
 436                 m->openInputFile(fname, f);
 437                 read_phylip(f, q, names, D);
 438         }
 439                 catch(exception& e) {
 440                 m->errorOut(e, "PCOACommand", "read");
 441                 exit(1);
 442         }
 443 }
 444
 445 /*********************************************************************************************************************************/
 446
 447 double PCOACommand::pythag(double a, double b)  {       return(pow(a*a+b*b,0.5));       }
 448
 449 /*********************************************************************************************************************************/
 450
 451 void PCOACommand::matrix_mult(vector<vector<double> > first, vector<vector<double> > second, vector<vector<double> >& product){
 452         try {
 453                 int first_rows = first.size();
 454                 int first_cols = first[0].size();
 455                 int second_cols = second[0].size();
 456
 457                 product.resize(first_rows);
 458                 for(int i=0;i<first_rows;i++){
 459                         product[i].resize(second_cols);
 460                 }
 461
 462                 for(int i=0;i<first_rows;i++){
 463                         for(int j=0;j<second_cols;j++){
 464                                 product[i][j] = 0.0;
 465                                 for(int k=0;k<first_cols;k++){
 466                                         product[i][j] += first[i][k] * second[k][j];
 467                                 }
 468                         }
 469                 }
 470         }
 471         catch(exception& e) {
 472                 m->errorOut(e, "PCOACommand", "matrix_mult");
 473                 exit(1);
 474         }
 475
 476 }
 477
 478 /*********************************************************************************************************************************/
 479
 480 void PCOACommand::recenter(double offset, vector<vector<double> > D, vector<vector<double> >& G){
 481         try {
 482                 int rank = D.size();
 483
 484                 vector<vector<double> > A(rank);
 485                 vector<vector<double> > C(rank);
 486                 for(int i=0;i<rank;i++){
 487                         A[i].resize(rank);
 488                         C[i].resize(rank);
 489                 }
 490
 491                 double scale = -1.0000 / (double) rank;
 492
 493                 for(int i=0;i<rank;i++){
 494                         A[i][i] = 0.0000;
 495                         C[i][i] = 1.0000 + scale;
 496                         for(int j=i+1;j<rank;j++){
 497                                 A[i][j] = A[j][i] = -0.5 * D[i][j] * D[i][j] + offset;
 498                                 C[i][j] = C[j][i] = scale;
 499                         }
 500                 }
 501
 502                 matrix_mult(C,A,A);
 503                 matrix_mult(A,C,G);
 504         }
 505         catch(exception& e) {
 506                 m->errorOut(e, "PCOACommand", "recenter");
 507                 exit(1);
 508         }
 509
 510 }
 511
 512 /*********************************************************************************************************************************/
 513
 514 //  This function is taken from Numerical Recipes in C++ by Press et al., 2nd edition, pg. 479
 515
 516 void PCOACommand::tred2(vector<vector<double> >& a, vector<double>& d, vector<double>& e){
 517         try {
 518                 double scale, hh, h, g, f;
 519
 520                 int n = a.size();
 521
 522                 d.resize(n);
 523                 e.resize(n);
 524
 525                 for(int i=n-1;i>0;i--){
 526                         int l=i-1;
 527                         h = scale = 0.0000;
 528                         if(l>0){
 529                                 for(int k=0;k<l+1;k++){
 530                                         scale += fabs(a[i][k]);
 531                                 }
 532                                 if(scale == 0.0){
 533                                         e[i] = a[i][l];
 534                                 }
 535                                 else{
 536                                         for(int k=0;k<l+1;k++){
 537                                                 a[i][k] /= scale;
 538                                                 h += a[i][k] * a[i][k];
 539                                         }
 540                                         f = a[i][l];
 541                                         g = (f >= 0.0 ? -sqrt(h) : sqrt(h));
 542                                         e[i] = scale * g;
 543                                         h -= f * g;
 544                                         a[i][l] = f - g;
 545                                         f = 0.0;
 546                                         for(int j=0;j<l+1;j++){
 547                                                 a[j][i] = a[i][j] / h;
 548                                                 g = 0.0;
 549                                                 for(int k=0;k<j+1;k++){
 550                                                         g += a[j][k] * a[i][k];
 551                                                 }
 552                                                 for(int k=j+1;k<l+1;k++){
 553                                                         g += a[k][j] * a[i][k];
 554                                                 }
 555                                                 e[j] = g / h;
 556                                                 f += e[j] * a[i][j];
 557                                         }
 558                                         hh = f / (h + h);
 559                                         for(int j=0;j<l+1;j++){
 560                                                 f = a[i][j];
 561                                                 e[j] = g = e[j] - hh * f;
 562                                                 for(int k=0;k<j+1;k++){
 563                                                         a[j][k] -= (f * e[k] + g * a[i][k]);
 564                                                 }
 565                                         }
 566                                 }
 567                         }
 568                         else{
 569                                 e[i] = a[i][l];
 570                         }
 571
 572                         d[i] = h;
 573                 }
 574
 575                 d[0] = 0.0000;
 576                 e[0] = 0.0000;
 577
 578                 for(int i=0;i<n;i++){
 579                         int l = i;
 580                         if(d[i] != 0.0){
 581                                 for(int j=0;j<l;j++){
 582                                         g = 0.0000;
 583                                         for(int k=0;k<l;k++){
 584                                                 g += a[i][k] * a[k][j];
 585                                         }
 586                                         for(int k=0;k<l;k++){
 587                                                 a[k][j] -= g * a[k][i];
 588                                         }
 589                                 }
 590                         }
 591                         d[i] = a[i][i];
 592                         a[i][i] = 1.0000;
 593                         for(int j=0;j<l;j++){
 594                                 a[j][i] = a[i][j] = 0.0;
 595                         }
 596                 }
 597         }
 598         catch(exception& e) {
 599                 m->errorOut(e, "PCOACommand", "tred2");
 600                 exit(1);
 601         }
 602
 603 }
 604
 605 /*********************************************************************************************************************************/
 606
 607 //  This function is taken from Numerical Recipes in C++ by Press et al., 2nd edition, pg. 479
 608
 609 void PCOACommand::qtli(vector<double>& d, vector<double>& e, vector<vector<double> >& z) {
 610         try {
 611                 int m, i, iter;
 612                 double s, r, p, g, f, dd, c, b;
 613
 614                 int n = d.size();
 615                 for(int i=1;i<=n;i++){
 616                         e[i-1] = e[i];
 617                 }
 618                 e[n-1] = 0.0000;
 619
 620                 for(int l=0;l<n;l++){
 621                         iter = 0;
 622                         do {
 623                                 for(m=l;m<n-1;m++){
 624                                         dd = fabs(d[m]) + fabs(d[m+1]);
 625                                         if(fabs(e[m])+dd == dd) break;
 626                                 }
 627                                 if(m != l){
 628                                         if(iter++ == 30) cerr << "Too many iterations in tqli\n";
 629                                         g = (d[l+1]-d[l]) / (2.0 * e[l]);
 630                                         r = pythag(g, 1.0);
 631                                         g = d[m] - d[l] + e[l] / (g + SIGN(r,g));
 632                                         s = c = 1.0;
 633                                         p = 0.0000;
 634                                         for(i=m-1;i>=l;i--){
 635                                                 f = s * e[i];
 636                                                 b = c * e[i];
 637                                                 e[i+1] = (r=pythag(f,g));
 638                                                 if(r==0.0){
 639                                                         d[i+1] -= p;
 640                                                         e[m] = 0.0000;
 641                                                         break;
 642                                                 }
 643                                                 s = f / r;
 644                                                 c = g / r;
 645                                                 g = d[i+1] - p;
 646                                                 r = (d[i] - g) * s + 2.0 * c * b;
 647                                                 d[i+1] = g + ( p = s * r);
 648                                                 g = c * r - b;
 649                                                 for(int k=0;k<n;k++){
 650                                                         f = z[k][i+1];
 651                                                         z[k][i+1] = s * z[k][i] + c * f;
 652                                                         z[k][i] = c * z[k][i] - s * f;
 653                                                 }
 654                                         }
 655                                         if(r == 0.00 && i >= l) continue;
 656                                         d[l] -= p;
 657                                         e[l] = g;
 658                                         e[m] = 0.0;
 659                                 }
 660                         } while (m != l);
 661                 }
 662
 663                 int k;
 664                 for(int i=0;i<n;i++){
 665                         p=d[k=i];
 666                         for(int j=i;j<n;j++){
 667                                 if(d[j] >= p){
 668                                         p=d[k=j];
 669                                 }
 670                         }
 671                         if(k!=i){
 672                                 d[k]=d[i];
 673                                 d[i]=p;
 674                                 for(int j=0;j<n;j++){
 675                                         p=z[j][i];
 676                                         z[j][i] = z[j][k];
 677                                         z[j][k] = p;
 678                                 }
 679                         }
 680                 }
 681         }
 682         catch(exception& e) {
 683                 m->errorOut(e, "PCOACommand", "qtli");
 684                 exit(1);
 685         }
 686 }
 687
 688 /*********************************************************************************************************************************/
 689
 690 void PCOACommand::output(string fnameRoot, vector<string> name_list, vector<vector<double> >& G, vector<double> d) {
 691         try {
 692                 int rank = name_list.size();
 693                 double dsum = 0.0000;
 694                 for(int i=0;i<rank;i++){
 695                         dsum += d[i];
 696                         for(int j=0;j<rank;j++){
 697                                 if(d[j] >= 0)   {       G[i][j] *= pow(d[j],0.5);       }
 698                                 else                    {       G[i][j] = 0.00000;                      }
 699                         }
 700                 }
 701
 702                 ofstream pcaData((fnameRoot+"pcoa").c_str(), ios::trunc);
 703                 pcaData.setf(ios::fixed, ios::floatfield);
 704                 pcaData.setf(ios::showpoint);
 705                 outputNames.push_back(fnameRoot+"pcoa");
 706                 outputTypes["pcoa"].push_back(fnameRoot+"pcoa");
 707
 708                 ofstream pcaLoadings((fnameRoot+"pcoa.loadings").c_str(), ios::trunc);
 709                 pcaLoadings.setf(ios::fixed, ios::floatfield);
 710                 pcaLoadings.setf(ios::showpoint);
 711                 outputNames.push_back(fnameRoot+"pcoa.loadings");
 712                 outputTypes["loadings"].push_back(fnameRoot+"pcoa.loadings");
 713
 714                 pcaLoadings << "axis\tloading\n";
 715                 for(int i=0;i<rank;i++){
 716                         pcaLoadings << i+1 << '\t' << d[i] * 100.0 / dsum << endl;
 717                 }
 718
 719                 pcaData << "group";
 720                 for(int i=0;i<rank;i++){
 721                         pcaData << '\t' << "axis" << i+1;
 722                 }
 723                 pcaData << endl;
 724
 725                 for(int i=0;i<rank;i++){
 726                         pcaData << name_list[i] << '\t';
 727                         for(int j=0;j<rank;j++){
 728                                 pcaData << G[i][j] << '\t';
 729                         }
 730                         pcaData << endl;
 731                 }
 732         }
 733         catch(exception& e) {
 734                 m->errorOut(e, "PCOACommand", "output");
 735                 exit(1);
 736         }
 737 }
 738
 739 /*********************************************************************************************************************************/
 740