Skip to content
GitLab
Projects
Groups
Snippets
/
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Menu
Maintenance - Mise à jour mensuelle Lundi 6 Février entre 7h00 et 9h00
Open sidebar
QTL
spell-qtl
Commits
fcb85d47
Commit
fcb85d47
authored
Nov 21, 2017
by
Damien Leroux
Browse files
various small fixes
parent
2d5d90de
Changes
17
Expand all
Hide whitespace changes
Inline
Side-by-side
CMakeLists.txt
View file @
fcb85d47
...
...
@@ -23,15 +23,15 @@ project(spell_qtl)
find_package
(
Boost 1.56.0 REQUIRED
)
#set(CMAKE_CONFIGURATION_TYPES Debug Release CACHE INTERNAL FORCE)
FIND_PACKAGE
(
PythonInterp 3
)
FIND_PACKAGE
(
PythonLibs 3
)
if
(
${
PYTHON_VERSION_MAJOR
}
EQUAL 3
)
FIND_PACKAGE
(
Boost COMPONENTS python3
)
SET
(
Boost_PYTHON_LIBRARY boost_python-py35
)
else
()
FIND_PACKAGE
(
Boost COMPONENTS python
)
endif
()
include_directories
(
AFTER
${
PYTHON_INCLUDE_DIR
}
)
#
FIND_PACKAGE(PythonInterp 3)
#
FIND_PACKAGE(PythonLibs 3)
#
if (${PYTHON_VERSION_MAJOR} EQUAL 3)
#
FIND_PACKAGE(Boost COMPONENTS python3)
#
SET(Boost_PYTHON_LIBRARY boost_python-py35)
#
else()
#
FIND_PACKAGE(Boost COMPONENTS python)
#
endif()
#
include_directories(AFTER ${PYTHON_INCLUDE_DIR})
set
(
CMAKE_VERBOSE_MAKEFILE ON
)
...
...
@@ -68,7 +68,7 @@ set(CMAKE_CXX_STANDARD 14)
SET
(
CMAKE_CXX_FLAGS
"
${
CMAKE_CXX_FLAGS
}
-Wextra -Wall -Wno-unused-parameter -pthread -fPIC"
)
SET
(
CMAKE_CXX_FLAGS_RELEASE
"
${
CMAKE_CXX_FLAGS
}
-O2 -DNDEBUG -DEIGEN_NO_DEBUG"
)
SET
(
CMAKE_CXX_FLAGS_DEBUG
"
${
CMAKE_CXX_FLAGS
}
-O
g
-ggdb"
)
SET
(
CMAKE_CXX_FLAGS_DEBUG
"
${
CMAKE_CXX_FLAGS
}
-O
0
-ggdb"
)
add_definitions
(
-DEIGEN_NO_DEPRECATED_WARNING -DVERSION_MAJOR=\"
${
VERSION_MAJOR
}
\" -DVERSION_MINOR=\"
${
VERSION_MINOR
}
\" -DVERSION_PATCH=\"
${
VERSION_PATCH
}
\"
)
...
...
@@ -86,11 +86,11 @@ include_directories(AFTER ${CMAKE_SOURCE_DIR}/include/ ${CMAKE_SOURCE_DIR}/inclu
include_directories
(
SYSTEM /usr/include
${
Boost_INCLUDE_DIRS
}
${
EXPAT_INCLUDE_DIR
}
)
include_directories
(
AFTER /usr/include
${
Boost_INCLUDE_DIRS
}
${
EXPAT_INCLUDE_DIR
}
)
if
(
${
BUILD_FOR_DEPLOYMENT
}
)
set
(
libstdcpp
${
CMAKE_BINARY_DIR
}
/libstdc++.a
)
else
(
${
BUILD_FOR_DEPLOYMENT
}
)
#
if(${BUILD_FOR_DEPLOYMENT})
#
set(libstdcpp ${CMAKE_BINARY_DIR}/libstdc++.a)
#
else(${BUILD_FOR_DEPLOYMENT})
set
(
libstdcpp
""
)
endif
(
${
BUILD_FOR_DEPLOYMENT
}
)
#
endif(${BUILD_FOR_DEPLOYMENT})
MESSAGE
(
STATUS
"libstdcpp =
${
libstdcpp
}
"
)
...
...
@@ -157,9 +157,7 @@ if (${BUILD_FOR_DEPLOYMENT})
add_dependencies
(
spell-marker glibc.h pouet
)
add_dependencies
(
spell-qtl glibc.h
)
MESSAGE
(
STATUS
" PROUT
${
libstdcpp
}
"
)
set
(
CMAKE_EXE_LINKER_FLAGS
"-include
${
CMAKE_BINARY_DIR
}
/glibc.h -rdynamic -static-libgcc"
)
set
(
CMAKE_EXE_LINKER_FLAGS
"-include
${
CMAKE_BINARY_DIR
}
/glibc.h -rdynamic -static-libgcc -static-libstdc++"
)
target_link_libraries
(
spell-marker
${
libstdcpp
}
${
CMAKE_BINARY_DIR
}
/libgmp.a
)
target_link_libraries
(
spell-qtl
${
libstdcpp
}
)
...
...
@@ -171,12 +169,14 @@ else()
SET_SOURCE_FILES_PROPERTIES
(
${
SPELL_PEDIGREE_SRC
}
${
SPELL_MARKER_SRC
}
${
SPELL_QTL_SRC
}
PROPERTIES
COMPILE_FLAGS
"-Wno-int-in-bool-context -U_FORTIFY_SOURCE -D_FORTIFY_SOURCE=0
${
SANITIZER_OPT
}
"
)
target_compile_definitions
(
spell-qtl PUBLIC SPELL_USE_XLNT
)
set
(
CMAKE_EXE_LINKER_FLAGS
"-rdynamic
${
SANITIZER_OPT
}
"
)
endif
()
target_compile_definitions
(
spell-qtl PUBLIC SPELL_USE_XLNT
)
target_link_libraries
(
spell-marker dl gmp
)
target_link_libraries
(
spell-qtl
${
Boost_PYTHON_LIBRARY
}
${
PYTHON_LIBRARY
}
xlnt dl rt
)
# target_link_libraries(spell-qtl ${Boost_PYTHON_LIBRARY} ${PYTHON_LIBRARY} xlnt dl rt)
target_link_libraries
(
spell-qtl xlnt dl rt
)
SET
(
EXECUTABLE_OUTPUT_PATH
${
PROJECT_BINARY_DIR
}
/bin
)
...
...
include/bayes/cli.h
View file @
fcb85d47
...
...
@@ -157,6 +157,13 @@ struct bn_settings_t {
*/
/*const_cast<pedigree_type&>(pedigree).load(pedigree_filename, false);*/
pedigree
.
load
(
pedigree_filename
,
false
);
for
(
const
auto
&
pg
:
pedigree
.
generations
)
{
if
(
pg
)
{
MSG_DEBUG
(
"Generation "
<<
pg
->
name
<<
" has "
<<
pg
->
labels
);
}
else
{
MSG_DEBUG
(
"NULL pg"
);
}
}
}
};
...
...
include/bayes/output.h
View file @
fcb85d47
...
...
@@ -378,7 +378,7 @@ void write_vector(ofile& ofs, const std::vector<V>& vec, WRITE_ELEM_FUNC write_e
inline
label_type
read_label
(
ifile
&
ifs
)
{
char
f
,
s
;
char
f
=
0
,
s
=
0
;
ifs
>>
f
>>
s
;
return
{
f
,
s
};
}
...
...
@@ -427,6 +427,7 @@ void read_geno_matrix(ifile& ifs, geno_matrix& mat)
read_matrix
(
ifs
,
mat
.
collect
);
read_matrix
(
ifs
,
mat
.
dispatch
);
/* also skip symmetries for now. */
MSG_DEBUG
(
"read geno_matrix "
<<
MATRIX_SIZE
(
mat
.
inf_mat
)
<<
" labels "
<<
mat
.
labels
);
}
...
...
include/cache2.h
View file @
fcb85d47
...
...
@@ -120,8 +120,7 @@ template <typename Ret, typename... Args>
*/
std
::
string
func_name
=
get_func_name
(
func
);
chrono_trace
_
(
func_name
);
/*
if (0) {
if
(
1
)
{
// msg_handler_t::cout() << "INVOKING " << func_name << "(…)" << std::endl;
std
::
stringstream
ss
;
ss
<<
"INVOKING "
<<
func_name
<<
"("
<<
std
::
endl
;
...
...
@@ -130,13 +129,14 @@ template <typename Ret, typename... Args>
for
(
size_t
i
=
0
;
i
<
avec
.
size
()
-
1
;
++
i
)
{
ss
<<
avec
[
i
]
<<
", "
;
}
msg_handler_t::cout() << ss.str() << ')';
ss
<<
')'
<<
std
::
endl
;
msg_handler_t
::
cout
()
<<
ss
.
str
();
}
// std::thread::id this_id = std::this_thread::get_id();
// active_settings->thread_stacks[this_id].push_back(func_name);
// MSG_DEBUG("[" << this_id << ',' << (this_id == active_settings->main_thread) << "] ENTER " << func_name);
msg_handler_t::run_hooks();
*/
//
msg_handler_t::run_hooks();
auto
ret
=
func
(
*
args
...);
// active_settings->thread_stacks[this_id].pop_back();
// msg_handler_t::run_hooks();
...
...
@@ -196,24 +196,23 @@ template <typename Ret, typename... Args>
*/
std
::
string
func_name
=
get_func_name
(
func
);
chrono_trace
_
(
func_name
);
/*
if (0) {
// msg_handler_t::cout() << "INVOKING " << func_name << "(" << std::endl;
if
(
1
)
{
// msg_handler_t::cout() << "INVOKING " << func_name << "(…)" << std::endl;
std
::
stringstream
ss
;
ss << "INVOKING " << func_name << "(
…)
" << std::endl;
ss
<<
"INVOKING "
<<
func_name
<<
"("
<<
std
::
endl
;
std
::
vector
<
std
::
string
>
avec
;
do_with_arg_pack
(
avec
.
push_back
(
MESSAGE
(
args
)));
for
(
size_t
i
=
0
;
i
<
avec
.
size
()
-
1
;
++
i
)
{
ss
<<
avec
[
i
]
<<
", "
;
}
msg_handler_t::cout() << ss.str() << ')';
ss
<<
')'
<<
std
::
endl
;
msg_handler_t
::
cout
()
<<
ss
.
str
();
}
// std::thread::id this_id = std::this_thread::get_id();
// active_settings->thread_stacks[this_id].push_back(func_name);
// MSG_DEBUG("[" << this_id << ',' << (this_id == active_settings->main_thread) << "] p ENTER " << func_name);
// msg_handler_t::run_hooks();
// m_promise.set_value(proxy(*args...));
*/
// mutex.lock();
auto
ret
=
proxy
(
*
args
...);
// mutex.unlock();
...
...
include/computations/report_output.h
View file @
fcb85d47
...
...
@@ -359,6 +359,9 @@ struct report_stream_text : public report_stream_filesystem {
void
matrix_impl
(
std
::
string
title
,
const
mws_result
&
mat
)
override
{
header
(
title
);
for
(
const
auto
&
sec
:
mat
.
m_column_sections
)
{
MSG_DEBUG
(
"[matrix_impl keys] "
<<
sec
.
field_labels
());
}
m_file
<<
mat
<<
std
::
endl
;
}
void
matrix_impl
(
std
::
string
title
,
const
mws_result_by_block
&
mat
)
override
...
...
@@ -480,7 +483,7 @@ make_full_map_item_list(const point_of_interest_type& poi, const region_of_inter
auto
ri
=
R
.
begin
(),
rj
=
R
.
end
();
auto
pi
=
P
.
begin
(),
pj
=
P
.
end
();
auto
li
=
chrom
.
raw
.
marker_locus
.
begin
(),
lj
=
chrom
.
raw
.
marker_locus
.
end
();
auto
ni
=
chrom
.
raw
.
marker_name
.
begin
(),
nj
=
chrom
.
raw
.
marker_name
.
end
();
auto
ni
=
chrom
.
raw
.
marker_name
.
begin
()
/*
, nj = chrom.raw.marker_name.end()
*/
;
for
(;
li
!=
lj
;
++
li
,
++
ni
)
{
double
l
=
*
li
;
while
(
pi
!=
pj
&&
pi
->
first
<
l
)
{
...
...
include/data/genoprob_computer2.h
View file @
fcb85d47
...
...
@@ -100,7 +100,22 @@ struct geno_prob_computer {
{
auto
it
=
tr_cache
.
find
(
dist
);
if
(
it
==
tr_cache
.
end
())
{
return
tr_cache
[
dist
]
=
this_gen
->
exp
(
dist
);
auto
&
TR
=
tr_cache
[
dist
]
=
this_gen
->
exp
(
dist
);
// if ((TR.array() > 1.).any() || (TR.array() < -1.e-13).any()) {
// double max = TR.array().abs().maxCoeff();
// double min = TR.minCoeff();
// MSG_DEBUG("BAD TR @" << dist << "cM min=" << min << " max=" << max << std::endl << TR);
// abort();
// }
double
min
=
TR
.
minCoeff
();
if
(
min
<
0
)
{
TR
.
array
()
+=
min
;
// double max = TR.array().abs().maxCoeff();
// MSG_DEBUG("BAD TR @" << dist << "cM min=" << min << " max=" << max << std::endl << TR);
// abort();
}
normalize_by_col
(
TR
);
return
tr_cache
[
dist
];
}
return
it
->
second
;
}
...
...
@@ -193,11 +208,19 @@ struct geno_prob_computer {
/*MSG_DEBUG("0 " << lstring[0]);*/
/*MSG_DEBUG(left[0]);*/
for
(
++
l
;
l
<=
r
;
++
l
)
{
MatrixXd
old_kernel
=
kernel
;
kernel
=
LV
->
col
(
l
).
asDiagonal
()
*
TR
[
l
-
1
]
*
kernel
;
double
s
=
kernel
.
sum
();
if
(
s
==
0
)
{
return
false
;
}
if
(
s
!=
s
)
{
MSG_ERROR
(
"#"
<<
l
<<
" LEFT KERNEL IS NaN"
,
""
);
MSG_DEBUG
(
"LV->col("
<<
l
<<
") ="
<<
std
::
endl
<<
LV
->
col
(
l
).
transpose
());
MSG_DEBUG
(
"TR["
<<
(
l
-
1
)
<<
"] ="
<<
std
::
endl
<<
TR
[
l
-
1
]);
MSG_DEBUG
(
"old kernel ="
<<
std
::
endl
<<
old_kernel
);
abort
();
}
kernel
/=
s
;
left
[
l
]
=
kernel
;
std
::
stringstream
ss
;
...
...
@@ -213,11 +236,19 @@ struct geno_prob_computer {
right
[
r
-
1
]
=
right
[
r
]
=
kernel
=
(
this_gen
->
stat_dist
.
array
()
*
LV
->
col
(
r
).
array
()).
matrix
();
//LV->col(r) / LV->col(r).sum();
{
std
::
stringstream
s
;
s
<<
".LV"
<<
r
;
rstring
[
r
-
1
]
=
rstring
[
r
]
=
s
.
str
();
}
for
(
--
r
;
r
>
0
;
--
r
)
{
MatrixXd
old_kernel
=
kernel
;
kernel
=
LV
->
col
(
r
).
asDiagonal
()
*
TR
[
r
]
*
kernel
;
double
s
=
kernel
.
sum
();
if
(
s
==
0
)
{
return
false
;
}
if
(
s
!=
s
)
{
MSG_ERROR
(
"#"
<<
r
<<
" RIGHT KERNEL IS NaN"
,
""
);
MSG_DEBUG
(
"LV->col("
<<
l
<<
") ="
<<
std
::
endl
<<
LV
->
col
(
l
).
transpose
());
MSG_DEBUG
(
"TR["
<<
(
l
-
1
)
<<
"] ="
<<
std
::
endl
<<
TR
[
l
-
1
]);
MSG_DEBUG
(
"old kernel ="
<<
std
::
endl
<<
old_kernel
);
abort
();
}
kernel
/=
s
;
right
[
r
-
1
]
=
kernel
;
std
::
stringstream
ss
;
...
...
@@ -272,11 +303,16 @@ struct geno_prob_computer {
MSG_DEBUG(L);
MSG_DEBUG(MATRIX_SIZE(R));
MSG_DEBUG(R);
MSG_DEBUG(MATRIX_SIZE(
this_gen->get_
redux
()
));
MSG_DEBUG(
this_gen->get_
redux
()
);
MSG_DEBUG(MATRIX_SIZE(redux));
MSG_DEBUG(redux);
#endif
VectorXd
vv
=
redux
*
(
L
*
R
*
inv_stat_dist
.
array
()).
matrix
();
double
s
=
vv
.
sum
();
if
(
s
!=
s
)
{
MSG_ERROR
(
"geno_prob_computer spotted a nan!"
,
""
);
MSG_DEBUG
(
vv
);
abort
();
}
if
(
s
==
0
)
{
return
vv
;
}
...
...
@@ -287,7 +323,7 @@ struct geno_prob_computer {
fast_compute_over_segment_raw
(
const
std
::
vector
<
double
>&
steps
)
{
MatrixXd
ret
(
unique_labels
.
size
(),
steps
.
size
());
/*
MSG_DEBUG(MATRIX_SIZE(ret));
*/
MSG_DEBUG
(
MATRIX_SIZE
(
ret
));
size_t
l
=
0
,
last
=
steps
.
size
()
-
1
,
r
=
!!
(
LV
->
outerSize
()
-
1
);
for
(
size_t
i
=
0
;
i
<
steps
.
size
();
++
i
)
{
/* FIXME: optimize the CORRECT initialization of l & r */
...
...
@@ -300,7 +336,7 @@ struct geno_prob_computer {
}
if
((
ret
.
col
(
i
).
array
()
<
0
).
any
())
{
MSG_DEBUG
(
"Bad value at column #"
<<
i
<<
std
::
endl
<<
ret
.
col
(
i
).
transpose
());
MSG_DEBUG
((
*
this_gen
));
//
MSG_DEBUG((*this_gen));
#if 0
MSG_DEBUG("steps: " << steps);
MSG_DEBUG("LV:" << std::endl << (*LV));
...
...
@@ -329,7 +365,7 @@ struct geno_prob_computer {
/*r += i < last;*/
}
}
/*
MSG_DEBUG("STATE PROBABILITIES" << std::endl << ret);
*/
//
MSG_DEBUG("STATE PROBABILITIES" << std::endl << ret);
return
ret
;
}
...
...
include/lumping2.h
View file @
fcb85d47
...
...
@@ -222,7 +222,7 @@ struct lumper {
/*const subset& C = *stack.front();*/
/*std::cout << "** " << C << std::endl;*/
/*std::cout << '.' << std::flush;*/
/*
MSG_DEBUG("TACKLING " << C);
*/
MSG_DEBUG
(
"TACKLING "
<<
C
);
stack
.
pop_back
();
/*stack.pop_front();*/
if
(
C
.
size
()
<=
1
)
{
...
...
@@ -233,20 +233,19 @@ struct lumper {
/*continue;*/
/*}*/
key_map
subsets
=
compute_keys
(
C
,
B
);
/*MSG_DEBUG("KEYS: " << subsets);*/
/*MSG_DEBUG("KEYS " << labels[subsets.begin()->second.front()] << ": " << subsets);*/
// MSG_DEBUG("KEYS: " << subsets);
/*size_t size = 0;*/
/*while (size != subsets.size()) {*/
/*size = subsets.size();*/
/*}*/
/* MARK */
if
(
subsets
.
size
()
>
1
)
{
/*
MSG_DEBUG("
KEYS " << labels[subsets.begin()->second
.front()] << ": " << subsets);
*/
MSG_DEBUG
(
"
WE HAVE A SPLIT "
<<
labels
[
C
.
front
()]
<<
'/'
<<
labels
[
B
.
front
()]
<<
": "
<<
subsets
);
P
.
erase
(
C
);
for
(
auto
&
s
:
subsets
)
{
auto
io
=
P
.
insert
(
s
.
second
);
if
(
!
io
.
second
)
{
/*
MSG_DEBUG("PROUT SUBSET ALREADY IN P " << (*io.first));
*/
MSG_DEBUG
(
"PROUT SUBSET ALREADY IN P "
<<
(
*
io
.
first
));
}
if
(
io
.
first
->
size
()
>
1
)
{
stack
.
push_back
(
&*
io
.
first
);
...
...
include/model/model.h
View file @
fcb85d47
...
...
@@ -916,7 +916,9 @@ struct model {
:
m_Y
(),
m_blocks
(),
m_X
(),
m_rank
(),
m_rss
(),
m_coefficients
(),
m_solver_type
(),
m_computed
(
false
)
,
m_with_constraints
(
false
),
m_ghost_constraint
(),
m_all_pops
(),
m_ancestor_names
(),
m_max_order
(
1
)
,
m_threshold
(
0
),
m_rank_fixer
()
{}
{
MSG_DEBUG
(
"[NEW MODEL] "
<<
keys
());
}
model
(
const
value
<
MatrixXd
>&
y
,
double
threshold
,
const
collection
<
const
qtl_pop_type
*>&
pops
,
const
std
::
map
<
char
,
std
::
string
>&
anam
,
size_t
mo
=
1
,
SolverType
st
=
SolverType
::
QR
)
:
m_Y
(
y
)
...
...
@@ -932,7 +934,9 @@ struct model {
,
m_threshold
(
threshold
)
,
m_rank_fixer
()
/*{ MSG_DEBUG("new model " << __LINE__ << " with Y(" << y.innerSize() << ',' << y.outerSize() << ')'); }*/
{}
{
MSG_DEBUG
(
"[NEW MODEL] "
<<
keys
());
}
model
(
const
value
<
MatrixXd
>&
y
,
double
threshold
,
const
collection
<
const
qtl_pop_type
*>&
pops
,
size_t
mo
=
1
,
SolverType
st
=
SolverType
::
QR
)
:
m_Y
(
y
)
...
...
@@ -948,7 +952,9 @@ struct model {
,
m_threshold
(
threshold
)
,
m_rank_fixer
()
/*{ MSG_DEBUG("new model " << __LINE__ << " with Y(" << y.innerSize() << ',' << y.outerSize() << ')'); }*/
{}
{
MSG_DEBUG
(
"[NEW MODEL] "
<<
keys
());
}
model
(
const
model
&
mo
)
:
m_Y
(
mo
.
m_Y
)
...
...
@@ -970,6 +976,7 @@ struct model {
model_block_key
mbk
=
std
::
make_shared
<
model_block_key_struc
>
(
*
kv
.
first
);
m_blocks
.
emplace
(
mbk
,
kv
.
second
);
}
MSG_DEBUG
(
"[NEW MODEL] "
<<
keys
());
}
model
(
const
value
<
MatrixXd
>&
y
,
const
model
&
mo
)
...
...
@@ -992,6 +999,7 @@ struct model {
model_block_key
mbk
=
std
::
make_shared
<
model_block_key_struc
>
(
*
kv
.
first
);
m_blocks
.
emplace
(
mbk
,
kv
.
second
);
}
MSG_DEBUG
(
"[NEW MODEL] "
<<
keys
());
}
model
&
...
...
@@ -1011,6 +1019,7 @@ struct model {
m_ancestor_names
=
mo
.
m_ancestor_names
;
m_max_order
=
mo
.
m_max_order
;
m_threshold
=
mo
.
m_threshold
;
MSG_DEBUG
(
"[ASSIGN MODEL] "
<<
keys
());
return
*
this
;
}
...
...
@@ -1031,6 +1040,7 @@ struct model {
m_ancestor_names
.
swap
(
mo
.
m_ancestor_names
);
m_max_order
=
mo
.
m_max_order
;
m_threshold
=
mo
.
m_threshold
;
MSG_DEBUG
(
"[ASSIGN MODEL] "
<<
keys
());
return
*
this
;
}
...
...
@@ -2022,11 +2032,13 @@ public:
{
for
(
const
auto
&
kv
:
m_blocks
)
{
if
(
kv
.
first
->
type
==
mbk_CI
)
{
std
::
vector
<
std
::
vector
<
char
>>
labels
;
for
(
const
auto
&
kp
:
m_all_pops
)
{
const
auto
&
qgn
=
(
*
kp
)
->
qtl_generation_name
;
labels
.
emplace_back
(
qgn
.
begin
(),
qgn
.
end
());
}
// std::vector<std::vector<char>> labels;
// for (const auto& kp: m_all_pops) {
// const auto& qgn = (*kp)->qtl_generation_name;
// labels.emplace_back(qgn.begin(), qgn.end());
// }
const
std
::
vector
<
std
::
vector
<
char
>>&
labels
=
kv
.
second
->
column_labels
;
MSG_DEBUG
(
"mbk_CI to print… Labels are "
<<
labels
);
mws
.
add_column_section
(
MESSAGE
(
""
<<
kv
.
first
),
labels
);
}
else
{
mws
.
add_column_section
(
MESSAGE
(
""
<<
kv
.
first
),
labels_to_ancestor_names
(
*
kv
.
second
));
...
...
@@ -2039,11 +2051,12 @@ public:
{
for
(
const
auto
&
kv
:
m_blocks
)
{
if
(
kv
.
first
->
type
==
mbk_CI
)
{
std
::
vector
<
std
::
vector
<
char
>>
labels
;
for
(
const
auto
&
kp
:
m_all_pops
)
{
const
auto
&
qgn
=
(
*
kp
)
->
qtl_generation_name
;
labels
.
emplace_back
(
qgn
.
begin
(),
qgn
.
end
());
}
const
std
::
vector
<
std
::
vector
<
char
>>&
labels
=
kv
.
second
->
column_labels
;
MSG_DEBUG
(
"mbk_CI to print… Labels are "
<<
labels
);
// for (const auto& kp: m_all_pops) {
// const auto& qgn = (*kp)->qtl_generation_name;
// labels.emplace_back(qgn.begin(), qgn.end());
// }
mws
.
add_column_section
(
kv
.
first
,
labels
);
mws
.
add_row_section
(
kv
.
first
,
labels
);
}
else
{
...
...
@@ -2123,14 +2136,30 @@ public:
std
::
ostream
&
output_keys
(
std
::
ostream
&
os
)
const
{
auto
output_key
=
[
&
]
(
model_block_key
k
,
value
<
model_block_type
>
b
)
{
os
<<
k
;
os
<<
'{'
;
bool
first
=
true
;
for
(
const
auto
&
label
:
b
->
column_labels
)
{
if
(
first
)
{
first
=
false
;
}
else
{
os
<<
','
<<
' '
;
}
os
<<
std
::
string
{
label
.
begin
(),
label
.
end
()};
}
os
<<
"} ("
<<
b
->
data
.
rows
()
<<
','
<<
b
->
data
.
cols
()
<<
')'
;
};
auto
i
=
m_blocks
.
begin
();
auto
j
=
m_blocks
.
end
();
if
(
i
==
j
)
{
return
os
<<
"[empty]"
;
}
o
s
<<
i
->
first
;
o
utput_key
(
i
->
first
,
i
->
second
)
;
for
(
++
i
;
i
!=
j
;
++
i
)
{
os
<<
", "
<<
i
->
first
;
os
<<
", "
;
output_key
(
i
->
first
,
i
->
second
);
}
return
os
;
}
...
...
include/python.h
0 → 100644
View file @
fcb85d47
This diff is collapsed.
Click to expand it.
include/settings.h
View file @
fcb85d47
...
...
@@ -608,6 +608,13 @@ settings_t::finalize()
MSG_DEBUG
(
"COVAR struct "
<<
covar_per_pop
.
back
());
pop_data
->
extract_subpops
(
populations
,
gt
.
second
,
observed_traits
,
linked_pops
,
covar_per_pop
.
back
());
}
for
(
const
auto
&
n_p
:
linked_pops
)
{
for
(
const
auto
&
vp
:
n_p
.
second
)
{
for
(
const
auto
&
p
:
vp
)
{
MSG_DEBUG
(
"trait "
<<
n_p
.
first
<<
" pop "
<<
p
->
name
<<
' '
<<
MATRIX_SIZE
(
p
->
gen
->
inf_mat
)
<<
" labels "
<<
p
->
gen
->
labels
);
}
}
}
for
(
const
auto
&
pleio
:
pleiotropic_traits_descr
)
{
auto
&
pops
=
linked_pops
[
pleio
.
name
];
MatrixXd
xtx
=
MatrixXd
::
Zero
(
pleio
.
traits
.
size
(),
pleio
.
traits
.
size
());
...
...
src/bayes/jobs.cc
View file @
fcb85d47
...
...
@@ -135,6 +135,12 @@ job_registry = {
[](
bn_settings_t
*
settings
)
{
return
settings
->
count_markers
();
},
[](
bn_settings_t
*
settings
,
size_t
mark
)
{
std
::
string
filename
=
settings
->
job_filename
(
"compute-alleles-per-marker"
,
mark
);
if
(
check_file
(
filename
,
false
,
true
,
false
))
{
return
true
;
}
/*const std::string& marker_name = settings->marker_names[mark];*/
auto
get_obs
=
[
&
]
(
const
std
::
string
&
gen
,
size_t
i
)
...
...
@@ -148,7 +154,7 @@ job_registry = {
return
ret
;
};
ofile
dump
(
settings
->
job_filename
(
"compute-alleles-per-marker"
,
mark
)
);
ofile
dump
(
filename
);
std
::
set
<
char
>
used_alleles
;
for
(
const
auto
&
kv
:
settings
->
observed_mark
)
{
...
...
@@ -181,6 +187,10 @@ job_registry = {
ifile
unafs
(
settings
->
job_filename
(
"unique_n_alleles"
,
0
));
rw_base
()(
unafs
,
unique_n_alleles
);
}
std
::
string
filename
=
settings
->
job_filename
(
"factor-graph"
,
unique_n_alleles
[
n
]);
if
(
check_file
(
filename
,
false
,
true
,
false
))
{
return
true
;
}
/*MSG_DEBUG("Pedigree" << std::endl << settings->pedigree);*/
/*MSG_DEBUG("COMPUTING FACTOR GRAPH FOR " << unique_n_alleles[n] << " ALLELES.");*/
/*factor_graph fg(settings->pedigree, unique_n_alleles[n], settings->noise);*/
...
...
@@ -191,7 +201,7 @@ job_registry = {
/*fg->dump();*/
//fg->to_image(MESSAGE("factor-graph-" << unique_n_alleles[n]), "png");
auto
I
=
instance
(
fg
);
ofile
of
(
settings
->
job_filename
(
"factor-graph"
,
unique_n_alleles
[
n
])
);
ofile
of
(
filename
);
I
->
file_io
(
of
);
rw_comb
<
int
,
bn_label_type
>
()(
of
,
fg
->
domains
);
return
true
;
...
...
@@ -201,6 +211,11 @@ job_registry = {
[](
bn_settings_t
*
settings
)
{
return
settings
->
count_markers
();
},
[](
bn_settings_t
*
settings
,
size_t
mark
)
{
std
::
string
filename
=
settings
->
job_filename
(
"compute-LV"
,
mark
);
if
(
check_file
(
filename
,
false
,
true
,
false
))
{
return
true
;
}
const
std
::
string
&
marker_name
=
settings
->
marker_names
[
mark
];
std
::
set
<
char
>
used_alleles
;
std
::
map
<
char
,
char
>
allele_obs_to_idx
;
...
...
@@ -377,7 +392,7 @@ job_registry = {
MSG_DEBUG
(
" * "
<<
kv
.
first
<<
' '
<<
kv
.
second
);
}
ofile
ofs
(
settings
->
job_filename
(
"compute-LV"
,
mark
)
);
ofile
ofs
(
filename
);
rw_base
()(
ofs
,
output_prob
);
}
return
true
;
...
...
@@ -581,7 +596,7 @@ job_registry = {
pop_data
.
marker_observation_filenames
[
om
.
first
]
=
om
.
second
.
filename
;
}
settings
->
load_full_pedigree_data
();
//
settings->load_full_pedigree_data();
/*pop_data.genetic_map_filename = settings->map_filename;*/
pop_data
.
pedigree_filename
=
settings
->
pedigree_filename
;
...
...
src/computations/frontends.cc
View file @
fcb85d47
...
...
@@ -312,7 +312,7 @@ model_manager&
qtl_detect_iqtlm
(
model_manager
&
mm
)
{
report_algo_phase
(
"QTL detection: iQTLm"
);
MSG_INFO
(
"iQTLm starts with: "
<<
mm
.
get_selection
());
//
MSG_INFO("iQTLm starts with: " << mm.get_selection());
std
::
set
<
std
::
map
<
chromosome_value
,
locus_key
>>
history
;
std
::
map
<
chromosome_value
,
locus_key
>
last_sel
,
new_sel
;
new_sel
=
mm
.
get_selection
();
...
...
src/computations/model.cc
View file @
fcb85d47
...
...
@@ -190,10 +190,13 @@ model_block_type
cross_indicator
(
const
std
::
string
&
trait
)
{
const
auto
&
pops
=
active_settings
->
linked_pops
[
trait
];
// MSG_DEBUG("[Building Cross Indicator] Have " << pops.size() << " distinct pops");
int
n_cols
=
pops
.
size
();
int
n_rows
=
0
;
for
(
const
auto
&
pop_vec
:
pops
)
{
// MSG_DEBUG("[Building Cross Indicator] Pop has " << pop_vec.size() << " subpops");
for
(
const
auto
&
pop
:
pop_vec
)
{
// MSG_DEBUG("[Building Cross Indicator] Subpop has " << pop->size() << " individuals");
n_rows
+=
pop
->
size
();
}
}
...
...
@@ -217,9 +220,9 @@ cross_indicator(const std::string& trait)
val
.
column_labels
.
clear
();
val
.
column_labels
.
reserve
(
val
.
data
.
outerSize
());
for
(
const
auto
&
pop_vec
:
pops
)
{
for
(
const
auto
&
pop
:
pop_vec
)
{
val
.
column_labels
.
emplace_back
(
pop
->
qtl_generation_name
.
begin
(),
pop
->
qtl_generation_name
.
end
());
}
//
for (const auto& pop: pop_vec) {
val
.
column_labels
.
emplace_back
(
pop
_vec
.
front
()
->
qtl_generation_name
.
begin
(),
pop
_vec
.
front
()
->
qtl_generation_name
.
end
());
//
}
}
/*for (const auto& pop: pops) {*/
...
...
@@ -230,7 +233,7 @@ cross_indicator(const std::string& trait)
/*n_rows += sz;*/
/*++p;*/
/*}*/
/*
MSG_DEBUG("cross indicator " << MATRIX_SIZE(val) << std::endl << val);
*/
//
MSG_DEBUG("
[Building Cross Indicator]
cross indicator " << MATRIX_SIZE(val) << std::endl << val);
/*MSG_QUEUE_FLUSH();*/
return
val
;
}
...
...
@@ -492,7 +495,9 @@ init_model(const value<MatrixXd>& Y, const value<model_block_type>& indicator, c
// const char* decomp_method = getenv("DECOMP_METHOD"); /* FIXME get rid of this sh\bt. */
// std::string s(decomp_method ? decomp_method : "");
M
.
m_max_order
=
1
;
MSG_DEBUG
(
"[NEW MODEL] indicator labels: "
<<
indicator
->
column_labels
);
M
.
add_block
(
model_block_key_struc
::
cross_indicator
(),
indicator
);
MSG_DEBUG
(
"[NEW MODEL] init_model with CI "
<<
M
.
keys
());
// MSG_DEBUG("covariables " << active_settings->with_covariables << " covars.cols=" << covars->data.cols());
if
(
covars
&&
covars
->
data
.
cols
()
>
0
)
{
M
.
add_block
(
model_block_key_struc
::
covariables
(),
covars
);
...
...
src/computations/report.cc
View file @
fcb85d47
...
...
@@ -311,9 +311,12 @@ analysis_report::report_algo_selection(model& Mbase, std::string chrom, double s
{
init_selection_report
();
Mbase
.
compute
();
mm
.
vMcurrent
->
compute
();
MSG_DEBUG
(
"[report_algo_selection] keys: "
<<
mm
.
vMcurrent
->
keys
());
++
m_selection_index
;