mussa_nway_other.cc

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