enumeffknocks.cpp

Go to the documentation of this file.

#include <iostream>
#include <fstream>   // file I/O
#include <sstream>
#include <stdlib.h>
#include <iomanip>
#include <algorithm>
#include <deque>
#include <list>
#include <queue>
#include <set>
#include <time.h>
#include <math.h>
#include <assert.h>
#include "enumeffknocks.h"
#include "ufes.h"
#include "rk.h"
#include "util.h"
#include "options.h"
#include "general_hash_functions.h"
using namespace std;



EnumEffKnocks::EnumEffKnocks() {}

EnumEffKnocks::EnumEffKnocks(MetabolicNetwork *metab_network) {

        this->metab_network=metab_network;
        this->erk_enum_time=0;
        this->ufes_enum_time=0;
        this->max_rk_length=this->metab_network->get_options().get_max_reac_subset_length();
        
}

void  EnumEffKnocks::enum_ufes() {

        vector<long>& eff_efm_indices=this->metab_network->get_eff_efm_indices();
        long first,last;


        int nump=1,proc_id=0;
        util::get_proc_id_nump(proc_id,nump);
        first=(proc_id)*(eff_efm_indices.size()/nump);
        if (proc_id<nump-1)
                last= (proc_id+1)*(eff_efm_indices.size()/nump);
        else
                last= eff_efm_indices.size();

        if (proc_id==0) {
           cout << endl << "\t"<<util::get_time_str()<<": # of modes explored: " << this->metab_network->get_num_modes();
           cout << endl << "\t"<<util::get_time_str() <<": enumerating UFES....";
        }

        for (long i=first;i<last;i++) {
                if (i%50==0 & proc_id==0)
                   cout << endl <<"\t"<<util::get_time_str() <<": crunching efm index "<<i<<" ....";
                enum_ufes_from_single(i);
        }
        
        if (proc_id==0) {
                cout << endl << "\t"<<util::get_time_str()  << ": # of UFES: "<<this->ufes_map.size();
                //cout << endl << "\t"<<util::get_time_str()  << ": # of UFES (hash): "<<this->ufes_hash.size();
        }

}
void EnumEffKnocks::enum_ufes_from_single(long eff_efm_index) {
        //enumerate ufes and add to this->ufes_set

        double total_time,start,end;
        this-> portion_time=0;
        queue<UFES> q;
        vector<long>& eff_efm_indices=this->metab_network->get_eff_efm_indices();
        
        util::current_time(&start);
        total_time=start;
        vector<unsigned char> bit_sets_table256;

        bit_sets_table256.resize(256,0);

        for (int i=0;i<=255;i++) 
           for (int j=0,mask=1;j<8;j++) {
              bit_sets_table256[i]+=((i & mask)>0);
              mask<<=1;            
           }

        //create UFES from single eff_efm_indices[ind] mode
        vector<long> sv;
        sv.push_back(eff_efm_indices[eff_efm_index]);
        UFES s(sv);
        //is s FES? if it is not return, else continue; but no need to check this. every single EFM is a FES

        //add UFES to q
        
        q.push(s);

        //warning! noo point in placing subsets that are not feasible to queue. make sure 

        while (!q.empty()) {
                //fetch UFES from queue
                UFES v=q.front();
                q.pop();

                if (v.get_efm_indices().size()==this->max_ufes_length)  {
                        //append v to ufes_set
                        this->ufes_map[v.get_key()]=v;
                        //this->ufes_hash[v.get_key().c_str()]=v;
                        continue;
                }

                //fetch index of last element in v
                vector<long> curr_efm_indices=v.get_efm_indices();

                vector<long>::reverse_iterator  ri=curr_efm_indices.rbegin();
                long last_efm_index=*ri;//curr_efm_indices[curr_efm_indices.size()-1];

                //vector<long>::iterator qi=find(eff_efm_indices.begin(),eff_efm_indices.end(), last_efm_index);
                vector<long>::iterator qi=lower_bound(eff_efm_indices.begin(),eff_efm_indices.end(), last_efm_index);
                long range_start=qi-(eff_efm_indices.begin());
                range_start++;
                long range_end=eff_efm_indices.end()-eff_efm_indices.begin();
        
                //rank modes in array this->eff_efm_indices[last_efm_index+1 ... this->eff_efm_indices->size()-1]

                curr_efm_indices.push_back(0); //add an extra element

                multimap<int,long>  ranked_efms;  //ranked modes, where key is equal to rank, and value to index of the mode

                compute_ranks_vector_bitwise(v, range_start, range_end, ranked_efms,bit_sets_table256);
                //util::current_time(&end);
                //portion_time+=end;
                int k=0;
                int not_expanding=1;
                
                int curr_max_efm_cands=this->metab_network->get_options().get_max_efm_cands();
                
                int curr_ufes_length=curr_efm_indices.size();
                if ((curr_max_efm_cands>>curr_ufes_length)>=1)
                        curr_max_efm_cands>>=curr_ufes_length;
                else
                        curr_max_efm_cands=1;
                

                for (multimap<int,long>::iterator mi=ranked_efms.begin();mi!=ranked_efms.end() && k<curr_max_efm_cands; mi++,k++)  {
                        curr_efm_indices[curr_efm_indices.size()-1]=mi->second;
                        if (!is_fes(curr_efm_indices))
                                continue;
                        not_expanding=0;
                        UFES v_new(curr_efm_indices);
                        q.push(v_new);                  
                }
                if (not_expanding) {
                        this->ufes_map[v.get_key()]=v;
                        //this->ufes_hash[v.get_key().c_str()]=v;
                }

        }

        util::current_time(&end);
        total_time=end-total_time;
        //cout << endl <<"\t total_time, portion_time: "<<total_time<<","<<this->portion_time;
}

int EnumEffKnocks::communicate_ufes_map() {

   //serialize this->ufes_map into vector<long>. keep two datastructures. ufes_map_srl, ufes_map_srl_lengths
   

   vector<unsigned long> ufes_cnts;
   unsigned long local_ufes_cnt=this->ufes_map.size();
   //unsigned long local_ufes_cnt=this->ufes_hash.size();

   int nump=1,proc_id=0;
   #ifdef PARALLEL
   util::get_proc_id_nump(proc_id,nump);

   //communicate size of ufes_map to every compute node
   ufes_cnts.resize(nump);
   
   MPI_Allgather(&local_ufes_cnt, 1, MPI_UNSIGNED_LONG,  &ufes_cnts[0],  1,  MPI_UNSIGNED_LONG,  MPI_COMM_WORLD );
   
   vector<int> rcv_count, rcv_displs;

   rcv_count.resize(nump);
   rcv_displs.resize(nump);

   unsigned long total_ufes_cnt=0;
   for (int i=0;i<nump;i++) {

        rcv_displs[i]=total_ufes_cnt;
        rcv_count[i]=ufes_cnts[i];
        total_ufes_cnt+=ufes_cnts[i];

   }
        
   //allocate array to store the sizes of ufes_map at each compute node
   vector<unsigned long> all_ufes_lengths, local_ufes_lengths;

   all_ufes_lengths.resize(total_ufes_cnt);


   unsigned long total_local_ufes_lengths=0;

   for (map<string,UFES>::iterator it=this->ufes_map.begin();it!=this->ufes_map.end();it++) {
        UFES u=it->second;
        local_ufes_lengths.push_back(u.get_efm_indices().size());
        total_local_ufes_lengths+=u.get_efm_indices().size();
   }
                
   MPI_Allgatherv(&local_ufes_lengths[0],local_ufes_lengths.size(),MPI_UNSIGNED_LONG, &all_ufes_lengths[0], &rcv_count[0], &rcv_displs[0],MPI_UNSIGNED_LONG,MPI_COMM_WORLD);

   unsigned long all_ufes_length=0;
   for (unsigned long i=0;i<total_ufes_cnt;i++)
        all_ufes_length+=all_ufes_lengths[i];

   vector<unsigned long> all_ufes,local_ufes;
   all_ufes.resize(all_ufes_length);
   local_ufes.resize(total_local_ufes_lengths);

   unsigned long k=0;
   for (map<string,UFES>::iterator it=this->ufes_map.begin();it!=this->ufes_map.end();it++) {
        UFES u=it->second;
        vector<long>& efm_indices=u.get_efm_indices();
        for (int j=0;j<efm_indices.size();j++)
                local_ufes[k++]=efm_indices[j]; 
   }
   rcv_count.clear();
   rcv_displs.clear();
   rcv_count.resize(nump,0);
   rcv_displs.resize(nump,0);

   unsigned long curr_index=0, total_length=0;
   for (int i=0;i<nump;i++) {

        rcv_displs[i]=total_length;
        for (unsigned long j=0;j<ufes_cnts[i];j++) {
                rcv_count[i]+=all_ufes_lengths[curr_index];
                curr_index++;
        }
        total_length+=rcv_count[i];

   }

   //cout << endl << "\t"<<util::get_time_str()  << ": total length: "<<total_length;
   //cout << endl << "\t"<<util::get_time_str()  << ": total local ufes length: "<<total_local_ufes_lengths;

   MPI_Allgatherv(&local_ufes[0],local_ufes.size(),MPI_UNSIGNED_LONG, &all_ufes[0], &rcv_count[0], &rcv_displs[0],MPI_UNSIGNED_LONG,MPI_COMM_WORLD);
  
   //communicate arrays of lengths of ufes_map to each communicate node
   curr_index=0;
   unsigned long ind_all_ufes=0;
   for (int i=0;i<nump;i++) {
        for (unsigned long j=0;j<ufes_cnts[i];j++) {
                vector<long> ufes_indices;
                for (unsigned long k=0;k<all_ufes_lengths[curr_index];k++) 
                        ufes_indices.push_back(all_ufes[ind_all_ufes++]);

                UFES q(ufes_indices);
                this->ufes_map[q.get_key()]=q;
                curr_index++;
        }
   }
  
   //communicate serialized ufes_map to each compute node
   #endif
   return 0;
}

int EnumEffKnocks::communicate_erk_map() {

   //serialize this->ufes_map into vector<long>. keep two datastructures. ufes_map_srl, ufes_map_srl_lengths
   
   vector<string> reacs_names=this->metab_network->get_reacs_names();
   vector<unsigned long> erk_cnts;
   unsigned long local_erk_cnt=this->erk_map.size();

   int nump=1,proc_id=0;
   #ifdef PARALLEL
   util::get_proc_id_nump(proc_id,nump);
 
   //communicate size of ufes_map to every compute node
   erk_cnts.resize(nump);
   MPI_Allgather(&local_erk_cnt, 1, MPI_UNSIGNED_LONG,  &erk_cnts[0],  1,  MPI_UNSIGNED_LONG,  MPI_COMM_WORLD );
   
   vector<int> rcv_count, rcv_displs;

   rcv_count.resize(nump);
   rcv_displs.resize(nump);

   unsigned long total_erk_cnt=0;
   for (int i=0;i<nump;i++) {

        rcv_displs[i]=total_erk_cnt;
        rcv_count[i]=erk_cnts[i];
        total_erk_cnt+=erk_cnts[i];

   }


   //allocate array to store the sizes of ufes_map at each compute node
   vector<unsigned long> all_erk_lengths, local_erk_lengths;

   all_erk_lengths.resize(total_erk_cnt);


   unsigned long total_local_erk_lengths=0;

   for (map<string,RK>::iterator it=this->erk_map.begin();it!=this->erk_map.end();it++) {
        RK u=it->second;
        local_erk_lengths.push_back(u.get_knock_reacs().size());
        total_local_erk_lengths+=u.get_knock_reacs().size();
   }     
      
   MPI_Allgatherv(&local_erk_lengths[0],local_erk_lengths.size(),MPI_UNSIGNED_LONG, &all_erk_lengths[0], &rcv_count[0], &rcv_displs[0],MPI_UNSIGNED_LONG,MPI_COMM_WORLD);

   unsigned long all_erk_length=0;
   for (unsigned long i=0;i<total_erk_cnt;i++)
        all_erk_length+=all_erk_lengths[i];

   vector<int> all_erk,local_erk;
   all_erk.resize(all_erk_length);
   local_erk.resize(total_local_erk_lengths);

   unsigned long k=0;
   for (map<string,RK>::iterator it=this->erk_map.begin();it!=this->erk_map.end();it++) {
        RK u=it->second;
        map<int,string>& knock_reacs=u.get_knock_reacs();
        for (map<int,string>::iterator mit=knock_reacs.begin();mit!=knock_reacs.end();mit++)
                local_erk[k++]=mit->first;      
   }
   
   rcv_count.clear();
   rcv_displs.clear();
   rcv_count.resize(nump,0);
   rcv_displs.resize(nump,0);

   unsigned long curr_index=0, total_length=0;
   for (int i=0;i<nump;i++) {

        rcv_displs[i]=total_length;
        for (unsigned long j=0;j<erk_cnts[i];j++) {
                rcv_count[i]+=all_erk_lengths[curr_index];
                curr_index++;
        }
        total_length+=rcv_count[i];

   }

   
   MPI_Allgatherv(&local_erk[0],local_erk.size(),MPI_INT, &all_erk[0], &rcv_count[0],&rcv_displs[0],MPI_INT,MPI_COMM_WORLD);
  
   //communicate arrays of lengths of erk_map to each communicate node
   curr_index=0;

   unsigned long ind_all_erk=0;
   for (int i=0;i<nump;i++) {
        for (unsigned long j=0;j<erk_cnts[i];j++) {

                map<int,string> knock_reacs;
                for (unsigned long k=0;k<all_erk_lengths[curr_index];k++)  {
                        knock_reacs[all_erk[ind_all_erk]]=reacs_names[all_erk[ind_all_erk]];
                        ind_all_erk++;
                }
                RK rk(knock_reacs,this->metab_network);

                this->erk_map[rk.get_key()]=rk;
                curr_index++;
        }
   }

   //communicate serialized ufes_map to each compute node
   #endif
   return 0;
}


void EnumEffKnocks::enum_erk() {
        map<string,UFES>::iterator pi;
        map<unsigned int,UFES>::iterator qi;
        //hash_map<const char*, UFES, hash<const char*>, eqstr>::iterator pi;
        int nump=1,proc_id=0;
        long first, last;
        util::get_proc_id_nump(proc_id,nump);

        map<unsigned int,UFES>  unique_ufes_map;

        for (pi=this->ufes_map.begin();pi!=this->ufes_map.end();pi++) {
                //get available reactions from pi->second
                vector<int> avail_reacs;
                stringstream reac_str;
                
                //cout << endl << "\t"<<util::get_time_str()<<": reacs_for_knock(begin) " << erk_map.size()<<" ERK sets";
                reacs_for_knock(pi->second, avail_reacs);
                for (int i=0;i<avail_reacs.size();i++)
                        reac_str<<avail_reacs[i];
                unique_ufes_map[JSHash(reac_str.str())]=pi->second;
        }
        cout << endl << "\t"<<util::get_time_str()<<": size of ufes_map: "<<this->ufes_map.size();
        cout << endl << "\t"<<util::get_time_str()<<": size of unique_ufes_map: "<<unique_ufes_map.size();
        long ufes_length=unique_ufes_map.size();
        //long ufes_length=ufes_hash.size();
        //proc_id=0;nump=1;
        first=proc_id*(ufes_length/nump);
        if (proc_id==nump-1)
                last=ufes_length;
        else
                last=(proc_id+1)*(ufes_length/nump);    

        if (proc_id==0)
                cout << endl << "\t"<<util::get_time_str()<<": (pid="<<proc_id<<") Exploring "<< (last-first) << " UFES ..."<<endl;
        long j=0;


        for (qi=unique_ufes_map.begin();qi!=unique_ufes_map.end();qi++,j++)
                if (j>=first && j<last) {
                        enum_erk_from_single_ufes(qi->second);          
                        if (j%20==0)
                                   cout << endl <<"\t"<<util::get_time_str() <<": crunching ufes no. "<<j<<" ....";
                }

        if (proc_id==0)
                cout << endl << "\t"<<util::get_time_str()<<": (pid="<<proc_id<<")  Enumerated " << erk_map.size()<<" ERK sets";

}




void EnumEffKnocks::enum_erk_direct() {

        queue<RK> q;
        //all reactions in the network
        vector<string> reacs_names=this->metab_network->get_reacs_names();
        //vector<long>& ems_to_collapse=this->metab_network->get_ineff_efm_indices();

        //vector<long>& eff_ems_to_collapse=this->metab_network->get_eff_efm_indices();
        //vector<long>& ineff_ems_to_collapse=this->metab_network->get_ineff_efm_indices();

        //reactions availabel for knockout for the given ufes
        vector<int> avail_reacs;
        //cout << endl << "\t"<<util::get_time_str()<<": reacs_for_knock(begin) " << erk_map.size()<<" ERK sets";

        UFES empty_ufes;
        reacs_for_knock(empty_ufes, avail_reacs);
        //cout << endl << "\t"<<util::get_time_str()<<": reacs_for_knock(end) " << erk_map.size()<<" ERK sets";

        //initialize emtpy RK object
        map<int,string> knock_reacs;
        //initialize empty RK object;
        RK rk(knock_reacs,this->metab_network);

        q.push(rk);

        while (!q.empty()) {
                
                //fetch UFES from queue
                RK curr_rk=q.front();
                q.pop();

                //cout << endl <<"\t"<<util::get_time_str() <<": crunching..";

                //if (curr_rk.get_knock_reacs().size()<this->metab_network->get_options().get_l1_height()) 
                //      this->max_rk_cands=this->metab_network->get_options().get_max_rk_cands_l1();
                //else
                //      this->max_rk_cands=this->metab_network->get_options().get_max_rk_cands_l2();

                int curr_max_rk_cands=this->metab_network->get_options().get_max_rk_cands();
                if ((curr_max_rk_cands>>curr_rk.get_knock_reacs().size())>=1)
                        curr_max_rk_cands>>=curr_rk.get_knock_reacs().size();
                else
                        curr_max_rk_cands=1;

                //if core EFMs are all collapsed, then continue loop
                int cnt_not_collapsed;
                curr_rk.is_core_collapsed(cnt_not_collapsed);

                if (cnt_not_collapsed==0) {
                        continue;
                }

                //find how many non-core and inefficient EFMs are collapsed.
                //vector<long> uncollapsed_ems;
                vector<long> &uncollapsed_eff_ems=curr_rk.get_uncollapsed_eff_efm_indices();
                vector<long> &uncollapsed_ineff_ems=curr_rk.get_uncollapsed_ineff_efm_indices();
                //uncollapsed_ems.resize(0);
                //collapse_ems(curr_rk, ems_to_collapse, uncollapsed_ems);
                //collapse_ems(curr_rk, eff_ems_to_collapse, uncollapsed_eff_ems);
                //collapse_ems(curr_rk, ineff_ems_to_collapse, uncollapsed_ineff_ems);



                //if curr_rk collapses all modes in ems_to_collapse then add it to table of ERKs


                //have all the inefficient modes been collapsed. if so, stop the iteration and add the reaction subset to collection
                if (uncollapsed_ineff_ems.empty()) {
                        if  (!uncollapsed_eff_ems.empty())
                                erk_map[curr_rk.get_key()]=curr_rk;                                             
                        continue;
                }
                //if curr_rk reached maximal allowed length then continue loop
                if (curr_rk.get_knock_reacs().size()==this->max_rk_length) {
                        continue;
                }

                //find index in reacs_names of v[v.size()-1]
                map<int,string>& curr_knock_reacs=curr_rk.get_knock_reacs();

                //explore reactions in curr_avail_reacs=(avail_reacs-curr_knock_reacs)
                //find  
                vector<int> curr_avail_reacs;
                for (vector<int>::iterator ai=avail_reacs.begin();ai!=avail_reacs.end();ai++) {
                        // is avail_reacs[i] in curr_knock_reacs
                        map<int,string>::iterator mi=curr_knock_reacs.find(*ai);
                        if (mi==curr_knock_reacs.end()) //not found
                                curr_avail_reacs.push_back(*ai);
                }               
                
                //multimap<long,int>  ranked_reacs;  //ranked modes, where key is equal to rank, and value to index of the mode         
                //reac_buckets stores reactions in buckets with equal rank value
                map<double, vector<int> > reac_buckets;
                for (vector<int>::iterator vi=curr_avail_reacs.begin();vi!=curr_avail_reacs.end();vi++) {

                        //for the given subset of efficient modes, determine the rank, and if it is FES 

                        //compute the effect of reaction 'i' when added to poolset v.get_knock_reacs() on ufes.
                        //rank is found as number of inefficient modes that are collapsed
                        
                        //long rank=compute_reac_rank(curr_knock_reacs,*vi,uncollapsed_ems);    
                        //long rank=rank_reacs_for_direct_enumeration(curr_knock_reacs,*vi,uncollapsed_ems);    
                        double rank=rank_reac_for_direct_enumeration(curr_knock_reacs,*vi, uncollapsed_ineff_ems,uncollapsed_eff_ems);
                        //if (rank>=0)
                        //cout << endl <<"\t"<<util::get_time_str() <<": rank "<<rank;
                        //long rank=0;
                        //don't consider reactions with rank=0 for knockout. no purpose.
                        if (rank==0)
                                continue;
                        //ranked_reacs.insert(pair<long,int>(rank,*vi));                        
                        reac_buckets[rank].push_back(*vi);
                }
                /*
                //define iterator to traverse map of ranked reacs (i.e. pairs (rank,reaction_cand_name))
                multimap<long,int>::reverse_iterator pi;
                //explore only the top this->max_rk_cands reaction candidates
                int k=0;                
                
                //pi->first: rank value; pi->second: reaction #
                for (pi=ranked_reacs.rbegin();pi!=ranked_reacs.rend() && k<this->max_rk_cands;pi++,k++) {
                        RK new_rk(curr_rk);
                        long rank;
                        new_rk.add_reac(pi->second,rank);
                        q.push(new_rk);                         
                }
                */
                map<double, vector<int> >::reverse_iterator rmi;
                int kk=0;
                //for (rmi=reac_buckets.rbegin();rmi!=reac_buckets.rend() && kk<this->max_rk_cands;rmi++) {
                for (rmi=reac_buckets.rbegin();rmi!=reac_buckets.rend() && kk<curr_max_rk_cands;rmi++) {
                        vector<int> vv=rmi->second;
                        sort(vv.begin(),vv.end());
                        //for (vector<int>::iterator vvi=vv.begin();vvi!=vv.end() && kk<this->max_rk_cands;vvi++,kk++) {
                        for (vector<int>::iterator vvi=vv.begin();vvi!=vv.end() && kk<curr_max_rk_cands;vvi++,kk++) {
                                RK new_rk(curr_rk);
                                double rank;
                                new_rk.add_reac_direct(*vvi,rank);
                                q.push(new_rk);                         
                        } 
                }

        }


}

void EnumEffKnocks::enum_erk_from_single_ufes(UFES ufes) {

        queue<RK> q;
        //all reactions in the network
        vector<string> reacs_names=this->metab_network->get_reacs_names();
        //vector<long>& ems_to_collapse=this->metab_network->get_ineff_efm_indices();
        //reactions availabel for knockout for the given ufes
        vector<int> avail_reacs;
        //cout << endl << "\t"<<util::get_time_str()<<": reacs_for_knock(begin) " << erk_map.size()<<" ERK sets";
        reacs_for_knock(ufes, avail_reacs);
        //cout << endl << "\t"<<util::get_time_str()<<": reacs_for_knock(end) " << erk_map.size()<<" ERK sets";

        //initialize emtpy RK object
        map<int,string> knock_reacs;
        //initialize empty RK object;
        RK rk(knock_reacs,this->metab_network);
        assert(knock_reacs.size()==0);

        double start, end, portion_time=0, portion_start, portion_end;
        util::current_time(&start);
        q.push(rk);

        while (!q.empty()) {
                //fetch UFES from queue
                RK curr_rk=q.front();
                q.pop();

                //if (curr_rk.get_knock_reacs().size()<this->metab_network->get_options().get_l1_height()) 
                //      this->max_rk_cands=this->metab_network->get_options().get_max_rk_cands_l1();
                //else
                //      this->max_rk_cands=this->metab_network->get_options().get_max_rk_cands_l2();

                int curr_max_rk_cands=this->metab_network->get_options().get_max_rk_cands();
                if ((curr_max_rk_cands>>curr_rk.get_knock_reacs().size())>=1)
                        curr_max_rk_cands>>=curr_rk.get_knock_reacs().size();
                else
                        curr_max_rk_cands=1;

                //if core EFMs are all collapsed, then continue loop
                int cnt_not_collapsed;
                curr_rk.is_core_collapsed(cnt_not_collapsed);

                if (cnt_not_collapsed==0) {
                        continue;
                }

                //find how many non-core and inefficient EFMs are collapsed.
                vector<long> &uncollapsed_ems=curr_rk.get_uncollapsed_ineff_efm_indices(); // zzzzP

                //collapse_ems(curr_rk, ems_to_collapse, uncollapsed_ems);
                assert(uncollapsed_ems.size()==curr_rk.get_uncollapsed_ineff_efm_indices().size());
                //cout << endl << uncollapsed_ems.size()<<","<<curr_rk.get_uncollapsed_ineff_efm_indices().size();

                //if curr_rk collapses all modes in ems_to_collapse then add it to table of ERKs
                //if (uncollapsed_ems.size()==0) {
                if (uncollapsed_ems.empty()) {
                        erk_map[curr_rk.get_key()]=curr_rk;                                             
                        continue;
                }
                //if curr_rk reached maximal allowed length then continue loop
                if (curr_rk.get_knock_reacs().size()==this->max_rk_length) {
                        continue;
                }

                //find index in reacs_names of v[v.size()-1]
                map<int,string>& curr_knock_reacs=curr_rk.get_knock_reacs();

                //explore reactions in curr_avail_reacs=(avail_reacs-curr_knock_reacs)
                //find  
                vector<int> curr_avail_reacs;
                for (vector<int>::iterator ai=avail_reacs.begin();ai!=avail_reacs.end();ai++) {
                        // is avail_reacs[i] in curr_knock_reacs
                        map<int,string>::iterator mi=curr_knock_reacs.find(*ai);
                        if (mi==curr_knock_reacs.end()) //not found
                                curr_avail_reacs.push_back(*ai);
                }               
                
                //multimap<long,int>  ranked_reacs;  //ranked modes, where key is equal to rank, and value to index of the mode         
                //reac_buckets stores reactions in buckets with equal rank value
                map<long, vector<int> > reac_buckets;
                for (vector<int>::iterator vi=curr_avail_reacs.begin();vi!=curr_avail_reacs.end();vi++) {

                        //for the given subset of efficient modes, determine the rank, and if it is FES 

                        //compute the effect of reaction 'i' when added to poolset v.get_knock_reacs() on ufes.
                        //rank is found as number of inefficient modes that are collapsed
                        long rank=compute_reac_rank(curr_knock_reacs,*vi,uncollapsed_ems);      
                        //don't consider reactions with rank=0 for knockout. no purpose.
                        if (rank==0)
                                continue;
                        //ranked_reacs.insert(pair<long,int>(rank,*vi));                        
                        reac_buckets[rank].push_back(*vi);
                }
                /*
                //define iterator to traverse map of ranked reacs (i.e. pairs (rank,reaction_cand_name))
                multimap<long,int>::reverse_iterator pi;
                //explore only the top this->max_rk_cands reaction candidates
                int k=0;                
                
                //pi->first: rank value; pi->second: reaction #
                for (pi=ranked_reacs.rbegin();pi!=ranked_reacs.rend() && k<this->max_rk_cands;pi++,k++) {
                        RK new_rk(curr_rk);
                        long rank;
                        new_rk.add_reac(pi->second,rank);
                        q.push(new_rk);                         
                }
                */
        
                map<long, vector<int> >::reverse_iterator rmi;
                int kk=0;
                //for (rmi=reac_buckets.rbegin();rmi!=reac_buckets.rend() && kk<this->max_rk_cands;rmi++) {
                for (rmi=reac_buckets.rbegin();rmi!=reac_buckets.rend() && kk<curr_max_rk_cands;rmi++) {
                        vector<int> vv=rmi->second;
                        sort(vv.begin(),vv.end());
                        //for (vector<int>::iterator vvi=vv.begin();vvi!=vv.end() && kk<this->max_rk_cands;vvi++,kk++) {
                        util::current_time(&portion_start);
                        for (vector<int>::iterator vvi=vv.begin();vvi!=vv.end() && kk<curr_max_rk_cands;vvi++,kk++) {
                
                                RK new_rk(curr_rk);
                                long rank;
                                new_rk.add_reac(*vvi,rank);
                
                                q.push(new_rk);                         
                        } 
                        util::current_time(&portion_end);
                }

                portion_time+=portion_end-portion_start;

        }
        
        util::current_time(&end);
        double total_time=end-start;
        //cout << endl <<"\t total_time, portion_time: "<<total_time<<","<<portion_time;

}
bool EnumEffKnocks::can_collapse_ineff_modes(vector<int> reac_knockout) {

        return true;
}
/*
   Description: Find out the rank of the knockout if mode with index cand_efm_index is appended to efm_indices

        return: if it is not FES, return -1, else return some value [0,1]
*/
int EnumEffKnocks::compute_rank(vector<long> efm_indices,long cand_efm_index) {
        //which reactions cannot be collapsed  to preserve efms[efm_indices]
        //see how many reactions are added to that pool with addition to efms[cand_efm_index]

        vector<string>& reacs_names=this->metab_network->get_reacs_names();
        //vector< vector<double> >& efms=this->metab_network->get_efms();

        set<int> cand_reacs;//, used_reacs;

        for (int i=0;i<reacs_names.size();i++) {
                int reac_used=0;
                for (long j=0;j<efm_indices.size();j++)  {
                        int is_zero;
                        this->metab_network->is_efm_elem_zero(efm_indices[j],i, is_zero); 
                        if (!is_zero)  {
                                reac_used=1;
                                break;
                        }               
                }
                if (reac_used) 
                        ;//used_reacs.insert(i);
                else
                        cand_reacs.insert(i);
        }       

        set<int> cand_efm_nnz(cand_reacs);
        for (int i=0;i<reacs_names.size();i++) {
                int is_zero;
                this->metab_network->is_efm_elem_zero(cand_efm_index,i, is_zero); 
                if (!is_zero)  {
                        cand_efm_nnz.insert(i);
                }
        }
        int rank=cand_efm_nnz.size()-cand_reacs.size();

        return rank;
}

void EnumEffKnocks::print_ufes(string ufes_fname) {


        ofstream fp_out;


        fp_out.open(ufes_fname.c_str(), ios::out); 
        if (fp_out.fail()) {
              cerr << endl << "\tunable to open file "<<ufes_fname<<" for reading" << endl;
              exit(1);
        }
        vector< vector<double> >& efm_yields=this->metab_network->get_efm_yields();
        vector<string> reacs_names=this->metab_network->get_reacs_names();


        for (map<string,UFES>::iterator it=ufes_map.begin();it!=ufes_map.end();it++) {
        //for (hash_map<const char*, UFES, hash<const char*>, eqstr>::iterator it=ufes_hash.begin();it!=ufes_hash.end();it++) {
                UFES ufes=it->second;
                //map<int,string>& knock_reacs rk.get_knock_reacs();
                vector<long>& efm_indices=ufes.get_efm_indices();
                fp_out<< efm_indices.size()<<" | ";
                for (int i=0;i<efm_indices.size();i++) {
                        
                        fp_out << efm_indices[i]<<"(";
                        for (int j=0;j<efm_yields.size();j++) {
                                fp_out<< efm_yields[j][efm_indices[i]];
                                if (j<efm_yields.size()-1)
                                        fp_out<<",";
                        }
                        fp_out<<") ";
                }
                fp_out << endl;
        }       
        fp_out.close();

}

//find reactions availabel for knockout (not used for ufes, not substrate, not targets, not init knocks
int EnumEffKnocks::reacs_for_knock(UFES ufes, vector<int> &avail_reacs) {
        //find which reactions are non-zero in subset ufes

        //vector< vector<double> >& efms=this->metab_network->get_efms();
        vector<long> efm_indices=ufes.get_efm_indices();
        vector<string>& reacs_names=this->metab_network->get_reacs_names();
        vector<long>& knock_cand_reacs=this->metab_network->get_knock_cand_reacs();
        for (vector<long>::iterator it=knock_cand_reacs.begin(); it!=knock_cand_reacs.end();it++) {             
                int not_used=1;
                for (vector<long>::iterator vi=efm_indices.begin();vi!=efm_indices.end();vi++) {
                        int is_zero;
                        //this->metab_network->is_efm_elem_zero(*vi,i, is_zero); 
                        this->metab_network->is_efm_elem_zero(*vi,*it, is_zero);
                        if (!is_zero)  {
                                not_used=0;
                                break;
                        }
                }               
                if (not_used)
                        avail_reacs.push_back(*it);
        }
        return 0;
}

double EnumEffKnocks::rank_reac_for_direct_enumeration(map<int,string> reac_knocks, int reac, vector<long>& remain_ineff_efm_indices,vector<long>& remain_eff_efm_indices) {

        char *buffer=this->metab_network->get_bit_efm().get_buffer();
        int bytes_per_efm=this->metab_network->get_bit_efm().get_bytes_per_efm(), word_index=reac/8;
        char bit_mask=1<<(reac%8);

        long num_ineff_collapsed=0;
        for (vector<long>::iterator vi=remain_ineff_efm_indices.begin();vi!=remain_ineff_efm_indices.end();vi++){
                
                if (buffer[(*vi)*bytes_per_efm+word_index] & bit_mask)
                        num_ineff_collapsed++;
        }


        long num_eff_collapsed=0;
        for (vector<long>::iterator vi=remain_eff_efm_indices.begin();vi!=remain_eff_efm_indices.end();vi++){
                
                if (buffer[(*vi)*bytes_per_efm+word_index] & bit_mask)
                        num_eff_collapsed++;
        }
        
        //if (num_eff_collapsed==0)
        //      return -1.0;
        

        map<int,string> lc_reac_knocks;
        lc_reac_knocks=reac_knocks;
        lc_reac_knocks[reac]="";
        int cnt_not_collapsed;
        int cores_collapsed=this->metab_network->is_core_efms_collapsed(lc_reac_knocks,cnt_not_collapsed);
        double rank;
        if (this->metab_network->get_options().get_network_capacity().compare("high")==0) {
                // compute score for reaction 
                double frac_ineff_collapsed=num_ineff_collapsed*1.0/remain_ineff_efm_indices.size();
                double frac_eff_collapsed=num_eff_collapsed*1.0/remain_eff_efm_indices.size();
                rank=frac_ineff_collapsed*(1-frac_eff_collapsed)*(!cores_collapsed);  //or frac_ineff_collapsed*frac_eff_collapsed

        } else if (this->metab_network->get_options().get_network_capacity().compare("low")==0) {

                int not_all_eff_collapsed=(remain_eff_efm_indices.size()!=num_eff_collapsed);
                rank=((num_ineff_collapsed+num_eff_collapsed)*1.0/(remain_eff_efm_indices.size()+remain_ineff_efm_indices.size()))*not_all_eff_collapsed*(!cores_collapsed);  //or frac_ineff_collapsed*frac_eff_collapsed

        }
        return rank;
}

int EnumEffKnocks::compute_reac_rank(map<int,string> reac_knocks, int reac, vector<long>& remain_ineff_efm_indices) {

        //long rank=0;  
        //vector< vector<double> >& efms=this->metab_network->get_efms();

        //store in remain_ineff_efm_indices the indices of those modes which are collapsed when reactions from reac_knocks are deleted
        
        //now find what fraction of modes indexed by remain_ineff_efm_indices are collapsed if  'reac' is deleted
        long num_collapsed=0;
        char *buffer=this->metab_network->get_bit_efm().get_buffer();
        int bytes_per_efm=this->metab_network->get_bit_efm().get_bytes_per_efm(), word_index=reac/8;
        char bit_mask=1<<(reac%8);
        
        for (vector<long>::iterator vi=remain_ineff_efm_indices.begin();vi!=remain_ineff_efm_indices.end();vi++){
                /*
                int is_zero;
                this->metab_network->is_efm_elem_zero(*vi,reac, is_zero); 
                if (!is_zero)           
                        num_collapsed++;
                */
                if (buffer[(*vi)*bytes_per_efm+word_index] & bit_mask)
                        num_collapsed++;
                
        }

        
        //rank=num_collapsed; //*1.0/remain_ineff_efm_indices.size();
        return num_collapsed;
}
/*
        Description: Rank the reaction 'reac' by how many remaining inefficient modes it collapses if deleted.

        Return: value in the range (0,1) representing the ratio of collapsed ineff. modes vs. total


*/
void EnumEffKnocks::print_erk(string erk_fname) {


        this->erk_fname=erk_fname;

        map<string,int>& target_names=this->metab_network->get_target_names();
        map<string,double>& max_yields=this->metab_network->get_max_yields();

        ofstream fp_out;
        fp_out.open(this->erk_fname.c_str(), ios::out); 
        if (fp_out.fail()) {
              cerr << endl << "\tunable to open file "<<this->erk_fname<<" for reading" << endl;
              exit(1);
        }

        //fp_out << this->metab_network->get_options();
        this->metab_network->get_options().print_to_file(fp_out);


        fp_out << "\tenumerated "<<this->ufes_map.size()<<" MFES subsets"<<endl;
        fp_out << "\tenumerated "<<this->erk_map.size()<<" ERK subsets"<<endl;
        fp_out << "\t# of efficient modes: "<<this->metab_network->get_eff_efm_indices().size()<<endl;
        fp_out << "\t# of inefficient modes: "<<this->metab_network->get_ineff_efm_indices().size()<<endl;
        fp_out << "\t# of core modes: "<<this->metab_network->get_core_efm_indices().size()<<endl;
        fp_out << "\t# of all modes (w/o core modes): "<<this->metab_network-> get_num_modes()<<endl;




        fp_out << "LEGEND:"<<endl;
        fp_out << "\tKO_size - number of reactions in knockout subset"<<endl;
        fp_out << "\t# eff - number of remaining efficient modes (with positive target fluxes) after knockout"<<endl;
        fp_out << "\t# ineff - number of remaining inefficient modes after knockout"<<endl;
        fp_out << "\tmin path length - minimal number of active reactions in remaining modes"<<endl;
        fp_out << "\tmax path length - maximal number of active reactions in remaining modes"<<endl;
        fp_out << "\t# core EFMs - number of remaining core elementary modes after knockout"<<endl;
        fp_out << "\tIR - reactions that are implicitly deactivated after deletion of knockout subset"<<endl;

        fp_out << endl << endl; 

        map<int, vector<RK> > rk_by_length;

        for (map<string,RK>::iterator it=erk_map.begin();it!=erk_map.end();it++) {
                RK rk=it->second;

                rk_by_length[rk.get_knock_reacs().size()].push_back(rk);
        }
        for (map<int, vector<RK> >::iterator rbli=rk_by_length.begin();rbli!=rk_by_length.end();rbli++) {
                vector<RK> rks=rbli->second;
                fp_out << endl<<"\tsubset length:"<<rbli->first << "\t# of subsets:"<<rks.size();
        }

        fp_out << endl << endl; 
        fp_out << "\tenumerating MFES time: "<<this->ufes_enum_time;
        fp_out << "\tenumerating ERK time: "<<this->erk_enum_time;

        fp_out << endl << endl; 

        fp_out << "\t\t\t\t# eff.";
        fp_out << "\t# ineff.";
        fp_out << "\tmin path length";
        fp_out << "\tmax path length";
        fp_out << "\t# core EFMs";


        for (map<string,int>::iterator mpi=target_names.begin();mpi!=target_names.end();mpi++) 
                fp_out << "\tmin yield("<<mpi->first<<")";

        fp_out << "\tknockout subset(s)";
        fp_out << "\tinactive reactions"<<endl;


        for (map<int, vector<RK> >::iterator rbli=rk_by_length.begin();rbli!=rk_by_length.end();rbli++) {
                vector<RK> rks=rbli->second;
                fp_out << endl<<"\t\tsubset length:"<<rbli->first << "\t# of subsets:"<<rks.size()<<endl;
                for (vector<RK>::iterator rksi=rks.begin();rksi!=rks.end();rksi++) {
                        
                        vector<double> min_yields;
                        long num_eff_efm_rem, num_ineff_efm_rem;
                        rksi->test_knock_reacs2(rksi->get_knock_reacs(), min_yields, num_eff_efm_rem, num_ineff_efm_rem);       
                        fp_out<< "\t\t\t\t"<<num_eff_efm_rem<<"\t"<<num_ineff_efm_rem;

                        fp_out<< "\t"<<rksi->get_min_path_length()<<"\t"<<rksi->get_max_path_length();          
        
                        //how many core modes remain
                        int cnt_not_collapsed;
                        
                        rksi->is_core_collapsed(cnt_not_collapsed);
                        fp_out << "\t" << cnt_not_collapsed;
                        int i=0;
                        fp_out << "";
                        map<string,double>::iterator iit;
                        
                        for (iit=max_yields.begin();iit!=max_yields.end();iit++,i++) 
                                fp_out << "\t"<<setprecision(2)<<min_yields[i]*100/(iit->second) <<"%";
                        
                        fp_out << "\t"<<*rksi;
                        fp_out << endl;
                }
                //int num_shared_efms;
                //find_common_efm(rks, num_shared_efms);
                //fp_out << endl<<"   # shared efms "<<num_shared_efms <<endl;  

        }

        fp_out.close();
}
/*

*/
int EnumEffKnocks::is_fes(vector<long> efm_indices) {

        //is this->efm_indices a feasible set. UFES or FES
        //find reactions that participate in modes indexed by this->efm_indices
        vector<long>& ineff_efm_indices=this->metab_network->get_ineff_efm_indices();
        double start, end;
        
        set<int> reacs_used;  //reactions that are not used in efms[efm_indices] modes  
        util::current_time(&start);
        //knockoutReactionCandidates=find(~any(net.ems(:,newSubsetEms),2))';
        for (int i=0;i<this->metab_network->get_num_reacs();i++) {
                //if reaction i is not used in effms[efm_indices] add it to ko_cand_reacs
                int used=0;
                for (vector<long>::iterator it=efm_indices.begin();it!=efm_indices.end();it++) {

                        int is_zero;
                        this->metab_network->is_efm_elem_zero(*it,i,is_zero);
                        
                        if (!is_zero) {
                                used=1;
                                break;
                        }
                }
                if (used)
                        reacs_used.insert(i);
                
        }



        int ret_val=this->metab_network->get_bit_efm().fast_is_fes(ineff_efm_indices, reacs_used);
        return ret_val;

        return 1;
}
int EnumEffKnocks::is_all_ems_collapsed(RK rk, vector<long> ems_to_collapse) {


        map<int,string>& knock_reacs=rk.get_knock_reacs();
        long num_collapsed=0;

        //find those modes indexed by ems_to_collapse that are not deleted with knock_reacs
        for (vector<long>::iterator it=ems_to_collapse.begin();it!=ems_to_collapse.end();it++) {
                int not_collapsed=1; //indicator for current mode that it is not collapsed
                for (map<int,string>::iterator mi=knock_reacs.begin();mi!=knock_reacs.end();mi++) {

                        int is_zero;
                        this->metab_network->is_efm_elem_zero(*it,mi->first, is_zero); 

                        if (!is_zero)  {
                                not_collapsed=0;
                                break;
                        }
                }
                if (!not_collapsed)
                        num_collapsed++;
        }
        if (num_collapsed==ems_to_collapse.size())
                return 1;
        
        return 0;

}

long EnumEffKnocks::collapse_ems(RK rk, vector<long> ems_to_collapse, vector<long>& uncollapsed_ems) {

        map<int,string>& knock_reacs=rk.get_knock_reacs();

        uncollapsed_ems.reserve(ems_to_collapse.size());
        //find those modes indexed by ems_to_collapse that are not deleted with knock_reacs
        for (vector<long>::iterator it=ems_to_collapse.begin();it!=ems_to_collapse.end();it++) {
                int not_collapsed=1;
                for (map<int,string>::iterator mi=knock_reacs.begin();mi!=knock_reacs.end();mi++) {
                        int is_zero;
                        this->metab_network->is_efm_elem_zero(*it,mi->first, is_zero); 

                        if (!is_zero)  {
                                not_collapsed=0;
                                break;
                        }
                }
                if (not_collapsed)
                        uncollapsed_ems.push_back(*it);
        }
        return 0;                       
}

int EnumEffKnocks::compute_ranks_vector(UFES v, long range_start, long range_end, multimap<int,long>&  ranked_efms) {

                vector<long>& eff_efm_indices=this->metab_network->get_eff_efm_indices();
                vector<string>& reacs_names=this->metab_network->get_reacs_names();
                vector<long> efm_indices=v.get_efm_indices();

                set<int> cand_reacs;

                for (int i=0;i<reacs_names.size();i++) {
                        int reac_used=0;
                        for (long j=0;j<efm_indices.size();j++)  {
                                int is_zero;
                                this->metab_network->is_efm_elem_zero(efm_indices[j],i, is_zero); 
                                if (!is_zero)  {                
                                        reac_used=1;
                                        break;
                                }               
                        }
                        if (!reac_used)                                         
                                cand_reacs.insert(i);
                }       
                for (long ii=range_start;ii<range_end;ii++) {
                                
                        long cand_efm_index=eff_efm_indices[ii];                

                        //which reactions cannot be collapsed  to preserve efms[efm_indices]
                        //see how many reactions are added to that pool with addition to efms[cand_efm_index]
                        set<int> cand_efm_nnz(cand_reacs);
                        for (int i=0;i<reacs_names.size();i++) {
                                int is_zero;
                                this->metab_network->is_efm_elem_zero(cand_efm_index,i, is_zero); 
                                if (!is_zero)  
                                        cand_efm_nnz.insert(i);                         
                        }
                        int rank=cand_efm_nnz.size()-cand_reacs.size();
                        ranked_efms.insert(pair<int,long>(rank,cand_efm_index));                        
                }
}

int EnumEffKnocks::compute_ranks_vector_bitwise(UFES v, long range_start, long range_end, multimap<int,long>&  ranked_efms,vector<unsigned char> &bit_sets_table256) {

                vector<long>& eff_efm_indices=this->metab_network->get_eff_efm_indices();
                int num_reacs=this->metab_network->get_num_reacs();
                vector<long> efm_indices=v.get_efm_indices();

                vector<unsigned char> used_reacs_bits;
                used_reacs_bits.resize(num_reacs/8+1,0);
                for (int i=0;i<num_reacs;i++) {
                        int reac_used=0;
                        for (long j=0;j<efm_indices.size();j++)  {
                                int is_zero;
                                this->metab_network->is_efm_elem_zero(efm_indices[j],i, is_zero); 
                                if (!is_zero)  {                
                                        reac_used=1;
                                        break;
                                }               
                        }
                        if (!reac_used)                                         
                                used_reacs_bits[i/8] |= 1<<(i%8);
                }       

                char *buffer=this->metab_network->get_bit_efm().get_buffer();
                int bytes_per_efm=this->metab_network->get_bit_efm().get_bytes_per_efm();
                int rank;
                for (long ii=range_start;ii<range_end;ii++) {
        
                        //int rank=compute_rank(v.get_efm_indices(),eff_efm_indices[i]);
                        long cand_efm_index=eff_efm_indices[ii];
                        rank=0;
                        for (int i=0;i<num_reacs/8+1;i++) {
                                unsigned char w=(used_reacs_bits[i]|buffer[cand_efm_index*bytes_per_efm+i])^(used_reacs_bits[i]);
                                rank+=bit_sets_table256[w];
                        }               
                        ranked_efms.insert(pair<int,long>(rank,cand_efm_index));                                                
                }
}