mussa_nway_other.cc

   1 //  This file is part of the Mussa source distribution.
   2 //  http://mussa.caltech.edu/
   3 //  Contact author: Tristan  De Buysscher, tristan@caltech.edu
   4
   5 // This program and all associated source code files are Copyright (C) 2005
   6 // the California Institute of Technology, Pasadena, CA, 91125 USA.  It is
   7 // under the GNU Public License; please see the included LICENSE.txt
   8 // file for more information, or contact Tristan directly.
   9
  10
  11 #include "mussa_nway.hh"
  12 #include <iostream>
  13
  14 using namespace std;
  15
  16 void
  17 Nway_Paths::radiate_path_search(vector<vector<FLPs> > all_comparisons)
  18 {
  19   vector<int> path;
  20   int win_i, sp_i, sp_depth, window_num;
  21   bool some_matches, still_paths, not_advanced;
  22   list<int> new_nodes;
  23   vector<list<int> > all_matches;
  24   vector<list<int>::iterator> sp_nodes, sp_nodes_end;
  25   list<int>::iterator debug_i;
  26
  27   pathz.clear();
  28   window_num = all_comparisons[0][1].win_num();
  29   cout << window_num << endl;    // just a sanity check
  30
  31   sp_nodes.reserve(species_num);
  32   sp_nodes_end.reserve(species_num);
  33
  34   // loop thru all windows in first species
  35   for (win_i = 0; win_i < window_num; win_i++)
  36   {
  37     // first we see if the first seq has matches to all other seqs at this win
  38     some_matches = true;
  39     sp_i = 1;
  40     all_matches.clear();
  41     while ( (sp_i < species_num) && (some_matches) )
  42     {
  43       new_nodes.clear();
  44       new_nodes = all_comparisons[0][sp_i].matches(win_i);
  45       if (new_nodes.empty())
  46         some_matches = false;
  47
  48       all_matches.push_back(new_nodes);
  49       sp_i++;
  50     }
  51     //cout << "fee\n";
  52     /*
  53     if (some_matches)
  54     {
  55       cout << win_i << " plrf" << endl;
  56       for (sp_i = 0; sp_i < species_num-1; sp_i++)
  57       {
  58         debug_i = all_matches[sp_i].begin();
  59         while (debug_i != all_matches[sp_i].end())
  60         {
  61           cout << *debug_i << " ";
  62           debug_i++;
  63         }
  64         cout << "plrfff"<< endl;
  65       }
  66     }
  67     */
  68     // if 1st seq does match to all others, make all possible paths
  69     // out of all these matches
  70     if (some_matches)
  71     {
  72       sp_nodes.clear();
  73       sp_nodes_end.clear();
  74       // set each species list of matches to beginning
  75       for (sp_i = 0; sp_i < species_num-1; sp_i++)
  76       {
  77         sp_nodes.push_back(all_matches[sp_i].begin());
  78         sp_nodes_end.push_back(all_matches[sp_i].end());
  79       }
  80
  81       still_paths = true;
  82       sp_depth = species_num - 2;       // this roves up & down species list
  83       //cout << "fie\n";
  84       while (still_paths)
  85       {
  86         // add path that each species iterator is pointing to
  87         path.clear();
  88         path.push_back(win_i);
  89
  90         //cout << win_i;
  91         for (sp_i = 0; sp_i < species_num-1; sp_i++)
  92         {
  93           //cout  << ", " << *(sp_nodes[sp_i]);
  94           path.push_back(*(sp_nodes[sp_i]));
  95         }
  96         //cout << endl;
  97
  98         pathz.push_back(path);
  99
 100         // now advance the right iterator
 101         not_advanced = true;
 102         while ( (not_advanced) && (sp_depth != -1) )
 103         {
 104           //cout << "foe\n";
 105           //cout << sp_depth << ".." << *(sp_nodes[sp_depth]) << endl;
 106           (sp_nodes[sp_depth])++;   // advance iter at current depth
 107           //cout << sp_depth << ".." << *(sp_nodes[sp_depth]) << endl;
 108
 109           // if we reached the end of the list, reset iter & go up one depth
 110           if (sp_nodes[sp_depth] == sp_nodes_end[sp_depth])
 111           {
 112             //cout << "fum\n";
 113             sp_nodes[sp_depth] = all_matches[sp_depth].begin();
 114             sp_depth--;
 115             //cout << "depth = " << sp_depth << endl;
 116           }
 117           else
 118             not_advanced = false;
 119         }
 120
 121         if (sp_depth == -1)   // jumped up to first species, all paths searched
 122           still_paths = false;
 123         else                 // otherwise just reset to max depth and continue
 124           sp_depth = species_num - 2;
 125       }
 126     }
 127   }
 128 }
 129
 130
 131
 132
 133 void
 134 Nway_Paths::trans_path_search(vector<vector<FLPs> > all_comparisons)
 135 {
 136   vector<int> path;
 137   int win_i, sp_i, sp_depth, window_num, i, cur_node;
 138   bool some_matches, still_paths, not_advanced, trans_good;
 139   list<int> new_nodes, trans_check_nodes;
 140   vector<list<int> > all_matches;
 141   vector<list<int>::iterator> sp_nodes, sp_nodes_end;
 142   list<int>::iterator debug_i;
 143
 144   cout << "trans: softhres = " << soft_thres;
 145   cout << ", window = " << win_size << ", ";
 146
 147   pathz.clear();
 148   window_num = all_comparisons[0][1].win_num();
 149   cout << "window number = " << window_num << endl;   // just a sanity check
 150   cout << "trans: species_num = " << species_num << endl;
 151   sp_nodes.reserve(species_num);
 152   cout << "trans: fie\n";
 153   sp_nodes_end.reserve(species_num);
 154   cout << "trans: foe\n";
 155   // loop thru all windows in first species
 156   for (win_i = 0; win_i < window_num; win_i++)
 157   {
 158     // first we see if the first seq has matches to all other seqs at this win
 159     some_matches = true;
 160     sp_i = 1;
 161     all_matches.clear();
 162     while ( (sp_i < species_num) && (some_matches) )
 163     {
 164       new_nodes.clear();
 165       // --thres
 166       //new_nodes = all_comparisons[0][sp_i].matches(win_i);
 167       new_nodes = all_comparisons[0][sp_i].thres_matches(win_i, soft_thres);
 168       if (new_nodes.empty())
 169         some_matches = false;
 170
 171       all_matches.push_back(new_nodes);
 172       sp_i++;
 173     }
 174     //cout << "fee\n";
 175     /*
 176     if (some_matches)
 177     {
 178       cout << win_i << " plrf" << endl;
 179       for (sp_i = 0; sp_i < species_num-1; sp_i++)
 180       {
 181         debug_i = all_matches[sp_i].begin();
 182         while (debug_i != all_matches[sp_i].end())
 183         {
 184           cout << *debug_i << " ";
 185           debug_i++;
 186         }
 187         cout << "plrfff"<< endl;
 188       }
 189     }
 190     */
 191     // if 1st seq does match to all others, make all possible paths
 192     // out of all these matches
 193     if (some_matches)
 194     {
 195       sp_nodes.clear();
 196       sp_nodes_end.clear();
 197       // set each species list of matches to beginning
 198       for (sp_i = 0; sp_i < species_num-1; sp_i++)
 199       {
 200         sp_nodes.push_back(all_matches[sp_i].begin());
 201         sp_nodes_end.push_back(all_matches[sp_i].end());
 202       }
 203
 204       still_paths = true;
 205       //cout << "fie\n";
 206       while (still_paths)
 207       {
 208         // add path that each species iterator is pointing to
 209         path.clear();
 210         path.push_back(win_i);
 211
 212         //cout << win_i;
 213         for (sp_i = 0; sp_i < species_num-1; sp_i++)
 214         {
 215           //cout  << ", " << *(sp_nodes[sp_i]);
 216           path.push_back(*(sp_nodes[sp_i]));
 217         }
 218         //cout << endl;
 219
 220         // check transitivity
 221         sp_depth = 1;
 222         trans_good = true;
 223         while ( (sp_depth != species_num-1) && (trans_good) )
 224         {
 225           cur_node = path[sp_depth];
 226           for(i = sp_depth+1; i < species_num; i++)
 227           {
 228             // --thres
 229             //trans_check_nodes = all_comparisons[sp_depth][i].matches(cur_node);
 230             trans_check_nodes =
 231               all_comparisons[sp_depth][i].thres_matches(cur_node, soft_thres);
 232
 233             if ( (trans_check_nodes.end() == find(trans_check_nodes.begin(),
 234                                                   trans_check_nodes.end(),
 235                                                   path[i]) ) &&
 236                  (trans_check_nodes.end() == find(trans_check_nodes.begin(),
 237                                                   trans_check_nodes.end(),
 238                                                   path[i] * -1) ) )
 239               trans_good = false;
 240           }
 241
 242           sp_depth++;
 243         }
 244
 245         if (trans_good)
 246           pathz.push_back(path);
 247
 248         // now advance the right iterator
 249         not_advanced = true;
 250         sp_depth = species_num - 2;       // this roves up & down species list
 251         while ( (not_advanced) && (sp_depth != -1) )
 252         {
 253           //cout << "foe\n";
 254           //cout << sp_depth << ".." << *(sp_nodes[sp_depth]) << endl;
 255           (sp_nodes[sp_depth])++;   // advance iter at current depth
 256           //cout << sp_depth << ".." << *(sp_nodes[sp_depth]) << endl;
 257
 258           // if we reached the end of the list, reset iter & go up one depth
 259           if (sp_nodes[sp_depth] == sp_nodes_end[sp_depth])
 260           {
 261             //cout << "fum\n";
 262             sp_nodes[sp_depth] = all_matches[sp_depth].begin();
 263             sp_depth--;
 264             //cout << "depth = " << sp_depth << endl;
 265           }
 266           else
 267             not_advanced = false;
 268         }
 269
 270         if (sp_depth == -1)   // jumped up to first species, all paths searched
 271           still_paths = false;
 272         else                 // otherwise just reset to max depth and continue
 273           sp_depth = species_num - 2;
 274       }
 275     }
 276   }
 277 }
 278
 279
 280
 281
 282 /*
 283   // loop thru all windows in first species
 284   for (win_i = 0; win_i < window_num; win_i++)
 285   {
 286     path.clear();
 287     path.push_back(win_i);
 288     cur_nodes = all_comparisons[0][1].matches(path[0]);
 289     cur_nodes_i = cur_nodes.begin();
 290     cur_nodes_end = cur_nodes.end();
 291     path.push_back(0);
 292     for (sp_i = 2; sp_i < species_num; sp_i++)
 293     {
 294       new_nodes = all_comparisons[0][sp_i].matches(path[0]);
 295       new_nodes_i = new_nodes.begin();
 296       new_nodes_end = new_nodes.end();
 297
 298       while(cur_nodes_i != cur_nodes_end)
 299       {
 300         intersect_bad = false;
 301         if ( (new_nodes.end() == find(new_nodes.begin(), new_nodes.end(),
 302                                       *cur_nodes_i) ) &&
 303              (new_nodes.end() == find(new_nodes.begin(), new_nodes.end(),
 304                                       *cur_nodes_i * -1) ) )
 305           intersect_bad = true;
 306
 307         // if no matches to this sequence, we don't need to check this node
 308         // against any of the remaining ones
 309         if (intersect_bad)
 310           cur_nodes_i = cur_nodes.erase(cur_nodes_i); // this advances the iter
 311         else
 312           ++cur_nodes_i;  // just advance normally
 313
 314       }
 315     }
 316
 317   }
 318 */
 319
 320
 321 /*
 322     while(new_nodes_i != new_nodes_end)
 323
 324
 325       path[1] = *new_nodes_i;
 326       path_search(all_comparisons, path, 2);
 327
 328 */
 329     //if (new_nodes_i != new_nodes_end)
 330     //cout << "* species 0 node: " << win_i << endl;
 331     //  cout << "foookin hell\n";
 332
 333      //cout << "  * species 1 node: " << *new_nodes_i << endl;
 334
 335     //cout << "    * species " << depth << " node: " << *new_nodes_i << endl;
 336     // check transitivity with previous nodes in path
 337
 338 /*
 339 void
 340 Nway_Paths::path_search(vector<vector<FLPs> > all_comparisons, vector<int> path, int depth)
 341 {
 342   list<int> new_nodes, trans_check_nodes;
 343   list<int>::iterator new_nodes_i, new_nodes_end;
 344   int i;
 345
 346   new_nodes = all_comparisons[depth - 1][depth].matches(path[depth-1]);
 347   new_nodes_i = new_nodes.begin();
 348   new_nodes_end = new_nodes.end();
 349
 350     if (trans_check_good)
 351     {
 352       // this makes sure path nodes are recorded with RC status relative to
 353       // the base species
 354       if ( path[depth-1] >= 0)
 355         path.push_back(*new_nodes_i);
 356       else
 357         path.push_back(*new_nodes_i * -1);
 358
 359       if (depth < species_num - 1)
 360         path_search(all_comparisons, path, depth + 1);
 361       else
 362         pathz.push_back(path);
 363       path.pop_back();
 364     }
 365     ++new_nodes_i;
 366   }
 367 }
 368 */
 369 //          cout << "    ****I have the power...\n";
 370
 371 /* use this if I ever get the friggin seqcomp match lists to sort...
 372       if (binary_search(trans_check_nodes.begin(), trans_check_nodes.end(),
 373                         *new_nodes_i))
 374 */