bitefm.cpp

Go to the documentation of this file.

#include <iostream>
#include <fstream>   // file I/O
#include <sstream>
#include <stdlib.h>
#include <vector>
#include <algorithm>
#include <math.h>
#include <iomanip>
#include "options.h"
#include "bitefm.h"
#include "util.h"
#include "mpi_utils.h"


BitEfm::BitEfm() { 

}
int BitEfm::load_efm_bit(string efm_bit_fname, t_UnsignedInt  num_ems, vector< vector<double> > &stoich) {

        ifstream fp_in;  // file for reading
        int nump=1,proc_id=0;
        util::get_proc_id_nump(proc_id,nump);

        this->efm_bit_fname=util::trim(efm_bit_fname);
        this->num_ems=num_ems;
        this->stoich=stoich;
        this->num_reacs=stoich[0].size();
        this->stoich.resize(stoich.size(),vector<double>(stoich[0].size(),0));
        for (int i=0;i<stoich.size();i++)
                for (int j=0;j<stoich[0].size();j++)
                        this->stoich[i][j]=stoich[i][j];


        //initialization finished

        //load bit efm matrix into buffer
        if (efm_bit_fname=="") {

                return 0;
        }
        
        long file_length;
        if (proc_id==0) {
                fp_in.open(efm_bit_fname.c_str(), ios::in|ios::binary); 

                if (fp_in.fail()) {
                      cerr << endl << "  unable to open file "<<efm_bit_fname<<"for reading" << endl;
                      this->num_ems=0;
                      return 0;
                }

                // get length of file:
                fp_in.seekg (0, ios::end);
                this->file_length = fp_in.tellg();
                file_length=this->file_length;
                fp_in.seekg (0, ios::beg);
                //this->bytes_per_efm=this->file_length/this->num_ems;  

                //load entire bit compressed matrix 
                // allocate memory:
                this->buffer = new  char [this->file_length];

                // read data as a block:
                fp_in.read(this->buffer,this->file_length);
                fp_in.close();
        }
        mpi_bcast_scalar(file_length);

        this->bytes_per_efm=file_length/this->num_ems;  
        /*
        #ifdef PARALLEL
                MPI_Bcast(&file_length,1,MPI_LONG,0,MPI_COMM_WORLD);
        #endif

        this->bytes_per_efm=file_length/this->num_ems;  
        */
        if (proc_id!=0) 
                this->buffer=new char[file_length];

        #ifdef PARALLEL
        MPI_Bcast(this->buffer,file_length,MPI_CHAR,0,MPI_COMM_WORLD);
        #endif
        
        if (proc_id==0) {
                cout << endl << "\t"<<util::get_time_str()<<": loaded bit efm matrix."<<endl;
                cout << endl << "\t"<<util::get_time_str()<<": num modes."<<this->num_ems;
        }
        return 0;
}
int BitEfm::fast_is_fes(vector<long> &ineff_efm_indices, set<int> &reacs) {
        vector<unsigned char> reac_map;
        reac_map.resize((this->num_reacs/8)+1,0);
        for (set<int>::iterator sit=reacs.begin();sit!=reacs.end();sit++) {     
                reac_map[(*sit)/8]|= 1<<((*sit)%8);
        }
        for (vector<long>::iterator it=ineff_efm_indices.begin();it!=ineff_efm_indices.end();it++) {
                int contained=1;
                for (int i=0;i<this->bytes_per_efm;i++) {
                //is this->buffer[(*it)*this->bytes_per_efm+i] equal to this->buffer[index*this->bytes_per_efm+i] || reac_map[i]
                        unsigned char w=this->buffer[(*it)*this->bytes_per_efm+i];
                        if (reac_map[i] != (w | reac_map[i])) {
                                contained=0;
                                break;
                        }
                }
                //
                if (contained)
                        return 0;

        }
        return 1;
        
}

int BitEfm::is_efm_elem_zero(unsigned long efm_index, int reac, int &is_zero) {

        
        unsigned char w=this->buffer[efm_index*this->bytes_per_efm+reac/8];
        unsigned char mask=1<<(reac%8);
        if (mask & w) 
                is_zero=0;
        else            
                is_zero=1;
        
        return 0;

}

int BitEfm::read_single_efm(unsigned long efm_index, vector<double> &efm) {

        vector<int> indices;
        indices.resize(this->stoich[0].size(),0);

        int num_nnz=0;
        for(int j=0;j<stoich[0].size()/8+1;j++){
            //indices((j-1)*8+1:j*8)=bitget(results.bit_ems(em_no(i)*bytes_per_em+j),1:8);
            //explore byte buffer[index*this->bytes_per_em+j]   
            unsigned char w=buffer[efm_index*this->bytes_per_efm+j];
            unsigned char mask=1;
            for (int i=0;i<8 && j*8+i<stoich[0].size();i++) {
                        if (mask & w) {
                                 indices[j*8+i]=1;
                                 num_nnz++;
                        }
                        mask<<=1;
            }
        }
        //find kernel of the matrix 
        vector< vector<double> > nnz_cols;
        nnz_cols.resize(stoich.size(),vector<double>(num_nnz));
        for (int i=0;i<this->stoich.size();i++) {
            int k=0;
            for (int j=0;j<this->stoich[0].size();j++)
                if (indices[j]) {
                        nnz_cols[i][k++]=this->stoich[i][j];
                }
        }


        vector< vector<double> > nsp;
        util::right_nullspace_nonint(nnz_cols, nsp);

        efm.resize(this->num_reacs,0);
        for (int i=0;i<this->num_reacs;i++)
                efm[i]=0;       
        if (nsp.size()==1) {


        for (int i=0,k=0;i<efm.size();i++)
                if (indices[i])
                        efm[i]=nsp[0][k++];
        } else {
                cerr << endl <<"error: nullity equal to "<<nsp.size();
                //exit(1);
        }       
        return 0;
        
        
}
int BitEfm::print_efm(unsigned long index, vector<string> &reacs_names) {
        vector<double> efm;
        read_single_efm(index, efm);
        //vector<string> reacs_names=this->metab_network->get_reacs_names();
        for (int i=0;i<efm.size();i++) 
                cout << "   "<<setprecision(2)<<efm[i]<<" "<<reacs_names[i];    
}

int BitEfm::get_num_reacs() {
        return this->num_reacs;
}
t_UnsignedInt BitEfm::get_num_modes() {
        return this->num_ems;
}


int BitEfm::retain_efm(vector<long> efm_indices2retain, vector< vector<double> > &yields) {

        char* new_buffer;

        //this->bytes_per_efm=this->file_length/this->num_ems;  
        new_buffer=new char[(efm_indices2retain.size())*this->bytes_per_efm];

        vector< vector<double> > new_yields;
        new_yields.resize(yields.size(), vector<double>((efm_indices2retain.size()),0));


        //scan efms to find 
        unsigned long k=0;
        //for (long i=0;i<this->num_ems;i++) {
        for (vector<long>::iterator it=efm_indices2retain.begin();it!=efm_indices2retain.end();it++) {

                //if i is not found put in efm_indices
                for (int j=0;j<this->bytes_per_efm;j++)
                                new_buffer[k*this->bytes_per_efm+j]=this->buffer[(*it)*this->bytes_per_efm+j];
                for (int j=0;j<yields.size();j++)
                                new_yields[j][k]=yields[j][(*it)];
                k++;
        }
        delete [] this->buffer;
        this->buffer=new_buffer;
        for (int j=0;j<yields.size();j++)
                yields[j].assign(new_yields[j].begin(),new_yields[j].end());

        this->num_ems=efm_indices2retain.size();

        return 0;       
}

int BitEfm::remove_efm(vector<long> efm_indices2remove, vector< vector<double> > &yields) {

        char* new_buffer;

        //this->bytes_per_efm=this->file_length/this->num_ems;  
        new_buffer=new char[(this->num_ems-efm_indices2remove.size())*this->bytes_per_efm];
        vector< vector<double> > new_yields;
        new_yields.resize(yields.size(), vector<double>((this->num_ems-efm_indices2remove.size()),0));

        //scan efms to find 
        unsigned long k=0;
        for (long i=0;i<this->num_ems;i++) {
                //is i in efm_indices2remove
                vector<long>::iterator p=find(efm_indices2remove.begin(),efm_indices2remove.end(),i);
                //if i is not found put in efm_indices
                if (p==efm_indices2remove.end()) {      
                        for (int j=0;j<this->bytes_per_efm;j++)
                                new_buffer[k*this->bytes_per_efm+j]=this->buffer[i*this->bytes_per_efm+j];

                        for (int j=0;j<yields.size();j++)
                                new_yields[j][k]=yields[j][i];
                        k++;
        
                }

        }
        for (int j=0;j<yields.size();j++)
                yields[j].assign(new_yields[j].begin(),new_yields[j].end());

        delete [] this->buffer;
        this->buffer=new_buffer;
        this->num_ems=this->num_ems-efm_indices2remove.size();

}

int BitEfm::knock_efms(map<string,int> knock_reacs, vector< vector<double> > &yields) {

        //vector< vector<double> > &efms=efms();

        //for reaction p->second remove elementary modes which have non-zero value in respective row
        vector<long> indices2remove;
        //vector<double> efm;
        //efm.resize(this->num_reacs);
        /*      
        for (long i=0;i<this->num_ems;i++)  {
                //read_single_efm(i,efm);
                if (i%1000==0)
                        cout << endl<<"\t examining index i :"<<i;                      
                map<string,int>::iterator p;
                for (p=knock_reacs.begin();p!=knock_reacs.end();p++)  {

                        int is_zero;
                        is_efm_elem_zero(i,p->second, is_zero); 
                        if (!is_zero) {
                        //if (efm[p->second]!=0)  {
                                indices2remove.push_back(i);
                                break;
                        }
                }
        }


        remove_efm(indices2remove);
        */

        vector<long> indices2retain;
        
        for (long i=0;i<this->num_ems;i++)  {
                map<string,int>::iterator p;
                int all_zero=1;
                for (p=knock_reacs.begin();p!=knock_reacs.end();p++)  {
                        int is_zero;
                        is_efm_elem_zero(i,p->second, is_zero); 
                        if (!is_zero) {
                                all_zero=0;
                                break;
                        }
                }
                if (all_zero)
                        indices2retain.push_back(i);
        }

        retain_efm(indices2retain, yields);

        return 0; 
}