Page MenuHomec4science
Files F121968904

ProbToModelApplication.cpp
No OneTemporary
Actions

File Metadata

Created
Tue, Jul 15, 03:16

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 <matrices.hpp>
namespace po = boost::program_options ;
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_i data = read_matrix2d_i(file_data) ;
matrix4d_d prob = read_matrix4d_d(this->file_prob) ;
if(data.size() != prob.size())
{ char msg[4096] ;
sprintf(msg, "Error! data and prob matrices have unequal "
"row numbers (%zu / %zu)!",
data.size(), prob.size()) ;
throw std::runtime_error(msg) ;
}
else if(data[0].size() < prob[0][0].size())
{ char msg[4096] ;
sprintf(msg, "Error! too many shift states for the data!"
"%zu shift states and %zu columns in data)!",
prob[0][0].size(), data[0].size()) ;
throw std::runtime_error(msg) ;
}
// get the data model
ModelComputer* ptr = nullptr ;
if(read_data)
{ ptr = new ReadModelComputer(data, prob) ; }
else if(seq_data)
{ ptr = new SequenceModelComputer(data, prob) ; }
matrix2d_d model = ptr->get_model() ;
delete ptr ;
ptr = nullptr ;
// compute the class prob
size_t n_row = prob.size() ;
size_t n_class = prob[0].size() ;
size_t n_shift = prob[0][0].size() ;
size_t n_flip = prob[0][0][0].size() ;
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_d model_final(model.size(),
vector_d(model[0].size() + 1)) ;
// 1st column contain the class prob
if(read_data)
{ for(size_t i=0; i<model_final.size(); 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.size(); i++)
{ if(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.size(); i++)
{ for(size_t j=0; j<model[0].size(); 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_parallel_msg = "The number of threads dedicated to the \n"
"computations, by default 1." ;
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())
("parallel", po::value<std::size_t>(&(this->n_threads)), opt_parallel_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) ;
}
// threads
if(this->n_threads == 0)
{ this->n_threads = 1 ; }
// 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