options.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 <sys/types.h>
#include <iomanip>
#include "options.h"
#include "util.h"
#include "mpi_utils.h"


string trimlr(const string &t)
{
    string str = t;
    size_t found;
    found = str.find_last_not_of(" \t");
    if (found != string::npos)
        str.erase(found+1);
    else
        str.clear();            // str is all whitespace

    return str;
}

Options::Options(string config_fname) {

        this->config_fname=config_fname;

}
void Options::load(string config_fname) {



        ifstream fp_in;  // file for reading
        int num_lines=0;
        string line;
        int nump=1, proc_id=0;
        string ko_fname;


        util::get_proc_id_nump(proc_id,nump);

        this->config_fname=config_fname;

        //parse the config_fname file
        //open config_fname

        //parse each line in file to reaad metatool_fname, targets, substrates, yields, core_yields, init knocks
                //load efm from  matrix (initialize num_mods, num_reacs, efms

        if (proc_id==0) {
                fp_in.open(this->config_fname.c_str(), ios::in); 
                if (fp_in.fail()) {
                      cerr << endl << "  unable to open file "<<this->config_fname<<"for reading" << endl;
                      exit(1);
                }

                while(!getline(fp_in, line, '\n').eof()) {
                
                        //parse line to read left and right side from : character
                        vector<string> result;
                        //skip comments
                        line=util::trim(line);
                        if (line.at(0)=='#')
                                continue;

                        parse_string(line, ':', result);
                        if (result.size()!=2)
                                continue;
                        string lhs,rhs;
                        lhs=result[0]; rhs=result[1];

                        std::stringstream trimmer;
                        trimmer << lhs;
                        lhs.clear();
                        trimmer >> lhs;

                        rhs=util::trim(rhs);

                        if (!lhs.compare("algorithm")) {
                                if (rhs.compare("direct")!=0 && rhs.compare("indirect")!=0) {
                                        cout << endl << "\t"<<util::get_time_str()<<": algorithm (direct, indirect) is incorrectly specified in config file "<< config_fname<<endl<<endl;
                                        exit(1);
                                }
                                this->algorithm=rhs;
                        }

                        if (!lhs.compare("network_capacity")) {
                                if (rhs.compare("low")!=0 && rhs.compare("high")!=0 ) {
                                        cout << endl << "\t"<<util::get_time_str()<<": network_capacity (low, high) is incorrectly specified in config file "<< config_fname<<endl<<endl;
                                        exit(1);
                                }
                                this->network_capacity=rhs;
                        }
                        
                        if (!lhs.compare("efm_compressed_file")) {

                                this->efm_bit_fname=rhs;
                        }
                        if (!lhs.compare("efm_meta_file")) {

                                this->efm_meta_fname=rhs;
                                bool ret_val=util::file_exists(this->efm_meta_fname);
                                if (ret_val==false) {
                                        cout << endl << "\t"<<util::get_time_str()<<": meta file '"<< this->efm_meta_fname<<"' does not exist"<<endl<<endl;
                                        exit(1);
                                }
                                         
                        }
                        if (!lhs.compare("use_xchg_reacs")) {
                                if (util::strcmp_ignore_case(rhs,"yes")==false && util::strcmp_ignore_case(rhs,"no")==false) {
                                        cout << endl << "\t"<<util::get_time_str()<<": 'use of exchange reactions' field value is incorrect. yes/no required "<< config_fname<<endl<<endl;
                                        exit(1);
                                }
                                this->use_xchg_reacs=rhs;
                        }

                        if (!lhs.compare("output_file")) {

                                this->output_fname=rhs;
                        }

                        if (!lhs.compare("max_efm_cands")) {

                                int max_efm_cands;
                                stringstream  ss;
                                ss << rhs;
                                ss >> max_efm_cands; 
                                this->max_efm_cands=max_efm_cands;
                        }

                        if (!lhs.compare("max_rk_cands")) {
                                int max_rk_cands;
                                stringstream  ss;
                                ss << rhs;
                                ss >> max_rk_cands; 
                                this->max_rk_cands=max_rk_cands;
                        }
                        if (!lhs.compare("max_reaction_knockout_subset_length")) {
                                int max_reac_subset_length;
                                stringstream  ss;
                                ss << rhs;
                                ss >> max_reac_subset_length; 
                                this->max_reac_subset_length=max_reac_subset_length;
                        }
        
                        if (!lhs.compare("init_knocks"))  {
                                //parse rhs
                                parse_string(rhs, ' ', this->init_knock_reacs);                 
                        }

                        if (!lhs.compare("substrates"))  {
                                //parse rhs
                                parse_string(rhs, ' ', this->substrate_reacs);                  
                        }

                        if (!lhs.compare("targets"))  {
                                //parse rhs
                                parse_string(rhs, ' ', this->target_reacs);                     
                        }

                        if (!lhs.compare("growth_eff_min_yield"))  {
                                //parse rhs
                                vector<string> yields_str;
                                parse_string(rhs, ' ', yields_str);
                                                
                                for (vector<string>::iterator it=yields_str.begin();it!=yields_str.end();it++) {
                                        string s_tmp=trimlr(*it);
                                        if (s_tmp.size()==0)
                                                continue;

                                        double dd;
                                        stringstream  ss;
                                        ss << *it;
                                        ss >> dd; 
                                        this->min_yields.push_back(dd);
                        
                                }
                        }

                        if (!lhs.compare("transport_reacs")) {
                                parse_string(rhs, ' ', this->transport_reacs);
                                                                        
                        }
                        if (!lhs.compare("production_eff_min_yield"))  {
                                //parse rhs
                                vector<string> yields_str;
                                parse_string(rhs, ' ', yields_str);
                                                
                                for (vector<string>::iterator it=yields_str.begin();it!=yields_str.end();it++) {
                                        string s_tmp=trimlr(*it);
                                        if (s_tmp.size()==0)
                                                continue;

                                        double dd;
                                        stringstream  ss;
                                        ss << *it;
                                        ss >> dd; 
                                        this->min_core_yields.push_back(dd);
                        
                                }
                        }
                
                        if (!lhs.compare("yields_file")) {

                                this->yields_fname=rhs;
                        }

                }

                fp_in.close();

        }

        mpi_bcast_string(this->algorithm);
        mpi_bcast_string(this->network_capacity);

        mpi_bcast_string(this->efm_bit_fname);
        mpi_bcast_string(this->efm_meta_fname);

        mpi_bcast_string(this->use_xchg_reacs);
        mpi_bcast_string(this->output_fname);
        mpi_bcast_scalar(this->max_efm_cands);

        mpi_bcast_scalar(this->max_rk_cands);

        mpi_bcast_vector_of_strings(this->init_knock_reacs) ;
        mpi_bcast_vector_of_strings(this->substrate_reacs);
        mpi_bcast_vector_of_strings(this->target_reacs); 
        mpi_bcast_vector(this->min_yields);
        mpi_bcast_vector_of_strings(this->transport_reacs);
        mpi_bcast_string(this->yields_fname);
        mpi_bcast_vector(this->min_core_yields);

        //broadcast
        //close  config_fname
        //check consistency 
        //test 1: length of arrays min_core_yields, min_yields, targets has to be same. value also has to be >0 and <num_reacs

        //test 2: min_core_yield and min_yields has to be (0,1]

        //test 3: use_xchg_reacs has to be yes/no and case insensitive

        //test 4: meta file has to be given and it has to exist (be able to open it)

        //test 5: at least one of two files has to be given with valid path: efm_full_fname, efm_bit_fname
        
}


void Options::print() {

        cout << "\t\t options" << endl;

        cout << "\talgorithm: "<<this->algorithm<<endl;

        cout << "\tnetwork_capacity: "<<this->network_capacity<<endl;


        cout << "\tEFM compressed file name:"<<this->efm_bit_fname<<endl;
        cout << "\tEFM meta file name:"<<this->efm_meta_fname<<endl;
        cout << "\tyields file: "<<this->yields_fname<< endl;
        cout << "\tinitial reaction knockouts: ";
        for (vector<string>::iterator it=this->init_knock_reacs.begin(); it!=this->init_knock_reacs.end();it++) 
                cout << *it <<" ";
        cout <<endl;
        cout << "\tsubstrate reactions: ";
        for (vector<string>::iterator it=this->substrate_reacs.begin(); it!=this->substrate_reacs.end();it++) 
                cout << *it <<" ";
        cout <<endl;
        cout << "\ttarget reactions: ";
        for (vector<string>::iterator it=this->target_reacs.begin(); it!=this->target_reacs.end();it++) 
                cout << *it <<" ";
        cout <<endl;
        cout << "\ttransport reactions: ";
        for (vector<string>::iterator it=this->transport_reacs.begin(); it!=this->transport_reacs.end();it++) 
                cout << *it <<" ";
        cout <<endl;

        cout << endl<<"\tmax efm candidates:"<<this->max_efm_cands;
        cout << endl<<"\tmax erk candidates:"<<this->max_rk_cands;
        cout << endl<<"\tmax reaction subset length:"<<this->max_reac_subset_length;

        cout << endl<< "\tminimum growth-efficient mode yields [0-1]: ";
        for (vector<double>::iterator it=this->min_yields.begin(); it!=this->min_yields.end();it++) 
                cout << *it <<" ";
        cout <<endl;
        cout << endl<< "\tminimum production-efficient mode yields [0-1]: ";
        for (vector<double>::iterator it=this->min_core_yields.begin(); it!=this->min_core_yields.end();it++) 
                cout << *it <<" ";
        cout <<endl;
        cout << "\tuse exchange reactions in knockout subset enumeration: "<<this->use_xchg_reacs <<endl;
        cout << "\toutput file: "<<this->output_fname<< endl;
        //cout << "\tbit-compressed EFM files: "<<this->ems_bit_fname<< "/"<<this->ems_bit_meta_fname<<endl;

}
void parse_string( const string& str, const char separator, vector<string>& result ) {

    if( str.size() ) {
        string::const_iterator last = str.begin();
        for( string::const_iterator i=str.begin(); i!=str.end(); ++i ) {
            if( *i == separator ) {
                const string str_tmp(last,i);
                //string token = str_tmp;
                result.push_back(str_tmp);
                last = i;
                ++ last;
            }
        }
        if( last != str.end() ) 
                result.push_back(&*last);
    }
}

void Options::load_efm_meta_file(vector< vector<double> > &stoich, vector<string> &reacs_names, t_UnsignedInt &num_ems) {


        ifstream fp_in;  // file for reading

        int num_lines=0, num_metabs,num_reacs;

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

        string line;

        num_ems=0;

        if (proc_id==0) {

                fp_in.open(this->efm_meta_fname.c_str(), ios::in); 

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

                while(!getline(fp_in, line, '\n').eof()) {

                        if (num_lines==0) {

                                vector<string> words;
                                parse_string(line, ' ', words);                 
                                num_metabs=atoi(words[0].c_str());
                                num_reacs=atoi(words[1].c_str());               
                                stoich.resize(atoi(words[0].c_str()),vector<double>(atoi(words[1].c_str())));                   
                                num_ems=atol(words[2].c_str()); //numems
                        } 

                
                        if (num_lines>0 && num_lines<=stoich.size()) {

                                //parse line
                                vector<string> words;

                                parse_string(line, ' ', words);                                         
                                for (int i=0;i<words.size();i++) {
                                        double val=atof(words[i].c_str());
                                        stoich[num_lines-1][i]=val;
                                }
                        }
        
                        if (num_lines>stoich.size()) {

                                vector<string> words;

                                parse_string(line, ' ', words);                                         
                                for (int i=0;i<words.size();i++) {
                                        reacs_names.push_back(words[i]);
                                }

                        }                       
                        num_lines++;    

                
                }               
                fp_in.close();
                cout << endl << "\t"<<util::get_time_str()<<": Loaded EFM meta file..."<<endl;
        }
        mpi_bcast_scalar(num_ems) ;
        mpi_bcast_scalar(num_metabs); 
        mpi_bcast_scalar(num_reacs);
        
        mpi_bcast_vector_of_strings(reacs_names);
        mpi_bcast_matrix(stoich);

}

/*
ostream& operator<<(ostream& out, const Options& option){

        Options opt=option;
        out << "\t\t options" << endl;
        out << "\tMETATOOL file name: "<<opt.metatool_fname<<endl;
        out << "\tEFM file name:"<<opt.efm_fname<<endl;
        out << "\tinit reaction knockouts: ";
        for (vector<string>::iterator it=opt.init_knock_reacs.begin(); it!=opt.init_knock_reacs.end();it++) 
                out << *it <<" ";
        out <<endl;
        out << "\tsubstrate reactions: ";
        for (vector<string>::iterator it=opt.substrate_reacs.begin(); it!=opt.substrate_reacs.end();it++) 
                out << *it <<" ";
        out <<endl;
        out << "\ttarget reactions: ";
        for (vector<string>::iterator it=opt.target_reacs.begin(); it!=opt.target_reacs.end();it++) 
                out << *it <<" ";
        out <<endl;
        out << "\ttransport reactions: ";
        for (vector<string>::iterator it=opt.transport_reacs.begin(); it!=opt.transport_reacs.end();it++) 
                out << *it <<" ";
        out <<endl;
        //max_efm_cands,max_rk_cands_l1,max_rk_cands_l2,l1_height;
        out << endl<<"\tmax efm candidates:"<<opt.max_efm_cands;
        out << endl<<"\tmax erk candidaates (level 1):"<<opt.max_rk_cands_l1;
        out << endl<<"\tmax erk candidaates (level 2):"<<opt.max_rk_cands_l2;
        out << endl<<"\theight (level 1):"<<opt.l1_height<<endl;

        out << "\tyields: ";
        for (vector<double>::iterator it=opt.min_yields.begin(); it!=opt.min_yields.end();it++) 
                out << *it <<" ";
        out <<endl;
        out << "\tcore yields: ";
        for (vector<double>::iterator it=opt.min_core_yields.begin(); it!=opt.min_core_yields.end();it++) 
                out << *it <<" ";
        out <<endl;
        out << "\tuse exchange reactionsin knockout: "<<opt.use_xchg_reacs <<endl;
        out << "\toutput file: "<<opt.output_fname<< endl;
        out << "\tbit-compressed EFM files: "<<opt.ems_bit_fname<< "/"<<opt.ems_bit_meta_fname<<endl;
        return out;

}
*/

void  Options::print_to_file(ostream &out) {
        out << "PROGRAM OPTIONS" << endl;

        out << "\talgorithm: "<<this->algorithm<<endl;
        if (this->algorithm.compare("direct")==0) {
                out << "\tnetwork_capacity: "<<this->network_capacity<<endl;
        }

        //out << "\tEFM expanded file name:"<<this->efm_full_fname<<endl;;
        out << "\tEFM compressed file name:"<<this->efm_bit_fname<<endl;
        out << "\tEFM meta file name:"<<this->efm_meta_fname<<endl;
        out << "\tyields file: "<<this->yields_fname<< endl;
        out << "\tinitial reaction knockouts: ";
        for (vector<string>::iterator it=this->init_knock_reacs.begin(); it!=this->init_knock_reacs.end();it++) 
                out << *it <<" ";
        out <<endl;
        out << "\tsubstrate reactions: ";
        for (vector<string>::iterator it=this->substrate_reacs.begin(); it!=this->substrate_reacs.end();it++) 
                out << *it <<" ";
        out <<endl;
        out << "\ttarget reactions: ";
        for (vector<string>::iterator it=this->target_reacs.begin(); it!=this->target_reacs.end();it++) 
                out << *it <<" ";
        out <<endl;
        out << "\ttransport reactions: ";
        for (vector<string>::iterator it=this->transport_reacs.begin(); it!=this->transport_reacs.end();it++) 
                out << *it <<" ";


        out << endl<<"\tmax efm candidates:"<<this->max_efm_cands;
        out << endl<<"\tmax erk candidates:"<<this->max_rk_cands;
        out << endl<<"\tmax reaction subset length:"<<this->max_reac_subset_length;

        out << endl<<endl<<"\tminimum growth-efficient mode yields [0-1]: ";
        for (vector<double>::iterator it=this->min_yields.begin(); it!=this->min_yields.end();it++) 
                out << *it <<" ";
        out <<endl;
        out << "\tminimum production-efficient mode yields [0-1]: ";
        for (vector<double>::iterator it=this->min_core_yields.begin(); it!=this->min_core_yields.end();it++) 
                out << *it <<" ";
        out <<endl;
        out << "\tuse exchange reactions in knockout subset enumeration: "<<this->use_xchg_reacs <<endl;
        out << "\toutput file: "<<this->output_fname<< endl;


        out << endl;

}