Page MenuHomec4science
Files F121099533

ProbToModelApplication.cpp
No OneTemporary
Actions

File Metadata

Created
Tue, Jul 8, 16:25

ProbToModelApplication.cpp
View Options

#include <ProbToModelApplication.hpp>
#include <boost/program_options.hpp>
#include <string>
#include <iostream>
#include <stdexcept> // std::invalid_argument, std::runtime_error
#include <ModelComputer.hpp>
#include <ReadModelComputer.hpp>
#include <SequenceModelComputer.hpp>
#include <Matrix2D.hpp>
#include <Matrix4D.hpp>
namespace po = boost::program_options ;
typedef std::vector<double> vector_d ;
ProbToModelApplication::ProbToModelApplication(int argn, char** argv)
: file_read(""), file_seq(""), file_prob(""),
n_threads(0), runnable(false)
{ this->parseOptions(argn, argv) ; }
ProbToModelApplication::~ProbToModelApplication()
{}
int ProbToModelApplication::run()
{ if(this->runnable)
{
// load data
std::string file_data ;
bool read_data = false ;
bool seq_data = false ;
if(this->file_read != "")
{ file_data = this->file_read ;
read_data = true ; seq_data = false ;
}
else if(this->file_seq != "")
{ file_data = this->file_seq ;
read_data = false ; seq_data = true ;
}
else
{ std::string msg("Error! Could not determine the type of the data!") ;
throw std::runtime_error(msg) ;
}
Matrix2D<int> data(file_data) ;
Matrix4D<double> prob(this->file_prob) ;
if(data.get_nrow() != prob.get_dim()[0])
{ char msg[4096] ;
sprintf(msg, "Error! data and prob matrices have unequal "
"row numbers (%zu / %zu)!",
data.get_nrow(), prob.get_dim()[0]) ;
throw std::runtime_error(msg) ;
}
else if(data.get_ncol() < prob.get_dim()[2])
{ char msg[4096] ;
sprintf(msg, "Error! too many shift states for the data!"
"%zu shift states and %zu columns in data)!",
prob.get_dim()[2], data.get_ncol()) ;
throw std::runtime_error(msg) ;
}
// get the data model
ModelComputer* ptr = nullptr ;
if(read_data)
{ ptr = new ReadModelComputer(data, prob, this->n_threads) ; }
else if(seq_data)
{ ptr = new SequenceModelComputer(data, prob, this->n_threads) ; }
Matrix2D<double> model = ptr->get_model() ;
delete ptr ;
ptr = nullptr ;
// compute the class prob
size_t n_row = prob.get_dim()[0] ;
size_t n_class = prob.get_dim()[1] ;
size_t n_shift = prob.get_dim()[2] ;
size_t n_flip = prob.get_dim()[3] ;
vector_d class_prob(n_class, 0.) ;
double p_tot = 0. ;
for(size_t i=0; i<n_row; i++)
{ for(size_t j=0; j<n_class; j++)
{ for(size_t k=0; k<n_shift; k++)
{ for(size_t l=0; l<n_flip; l++)
{ class_prob[j] += prob(i,j,k,l) ;
p_tot += prob(i,j,k,l) ;
}
}
}
}
for(auto& prob : class_prob)
{ prob /= p_tot ; }
// create a matrix containing the class prob in the 1st
// column and the model in the remaining columns
Matrix2D<double> model_final(model.get_nrow(),
model.get_ncol() + 1) ;
// 1st column contain the class prob
if(read_data)
{ for(size_t i=0; i<model_final.get_nrow(); i++)
{ model_final(i,0) = class_prob[i] ; }
}
else if(seq_data)
{ size_t i_class = 0 ;
for(size_t i=0; i<model_final.get_nrow(); i++)
{ if((i != 0) and (i % 4 == 0))
{ i_class++ ; }
model_final(i,0) = class_prob[i_class] ;
}
}
// fill the remaining with the model parameters
for(size_t i=0; i<model.get_nrow(); i++)
{ for(size_t j=0; j<model.get_ncol(); j++)
{ model_final(i,j+1) = model(i,j) ; }
}
std::cout << model_final << std::endl ;
return 0 ;
}
else
{ return 1 ; }
}
void ProbToModelApplication::parseOptions(int argn, char** argv)
{ // no option to parse
if(argv == nullptr)
{ std::string message = "no options to parse!" ;
throw std::invalid_argument(message) ;
}
// help messages
std::string desc_msg = "\n"
"ProbToRef reconstructs the class models \n"
"from the classification results (posterior \n"
"probabilities computed by ChIPPartitioning \n"
"and the data).\n\n" ;
std::string opt_help_msg = "Produces this help message." ;
std::string opt_thread_msg = "The number of threads dedicated to parallelize the computations,\n "
"by default 0 (no parallelization)." ;
std::string opt_read_msg = "The path to the file containing the data, if the data are read counts.";
std::string opt_seq_msg = "The path to the file containing the data, if the data are sequences.";
std::string opt_prob_msg = "The path to the file containing the posterior \n"
"probabilities." ;
// option parser
boost::program_options::variables_map vm ;
boost::program_options::options_description desc(desc_msg) ;
std::string seeding_tmp ;
desc.add_options()
("help,h", opt_help_msg.c_str())
("read,", po::value<std::string>(&(this->file_read)), opt_read_msg.c_str())
("seq,", po::value<std::string>(&(this->file_seq)), opt_seq_msg.c_str())
("prob,", po::value<std::string>(&(this->file_prob)), opt_prob_msg.c_str())
("thread", po::value<std::size_t>(&(this->n_threads)), opt_thread_msg.c_str()) ;
// parse
try
{ po::store(po::parse_command_line(argn, argv, desc), vm) ;
po::notify(vm) ;
}
catch(std::invalid_argument& e)
{ std::string msg = std::string("Error! Invalid option given!\n") + std::string(e.what()) ;
throw std::invalid_argument(msg) ;
}
catch(...)
{ throw std::invalid_argument("An unknown error occured while parsing the options") ; }
bool help = vm.count("help") ;
// checks unproper option settings
if((this->file_read == "") and
(this->file_seq == "") and
(not help))
{ std::string msg("Error! No data file was given (--read or --seq)!") ;
throw std::invalid_argument(msg) ;
}
else if((this->file_read != "") and
(this->file_seq != "") and
(not help))
{ std::string msg("Error! --read and --seq are mutually exclusive!") ;
throw std::invalid_argument(msg) ;
}
else if(this->file_prob == "" and (not help))
{ std::string msg("Error! No posterior probabily file was given (--prob)!") ;
throw std::invalid_argument(msg) ;
}
// help invoked, run() cannot be invoked
if(help)
{ std::cout << desc << std::endl ;
this->runnable = false ;
return ;
}
// everything fine, run() can be called
else
{ this->runnable = true ;
return ;
}
}
int main(int argn, char** argv)
{ ProbToModelApplication app(argn, argv) ;
return app.run() ;
}

Event Timeline

c4science · Help