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