Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F85148231
ProbToModelApplication.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Fri, Sep 27, 02:48
Size
7 KB
Mime Type
text/x-c
Expires
Sun, Sep 29, 02:48 (2 d)
Engine
blob
Format
Raw Data
Handle
21135000
Attached To
R8820 scATAC-seq
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
Log In to Comment