Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F121099533
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
Tue, Jul 8, 16:25
Size
7 KB
Mime Type
text/x-c
Expires
Thu, Jul 10, 16:25 (2 d)
Engine
blob
Format
Raw Data
Handle
27252228
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