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
198b00c0
Commit
198b00c0
authored
Dec 23, 2017
by
damien
Browse files
First implementation of spell-map.
parent
27ec4c15
Changes
8
Hide whitespace changes
Inline
Side-by-side
CMakeLists.txt
View file @
198b00c0
...
...
@@ -100,6 +100,12 @@ set(SPELL_PEDIGREE_SRC
src/pedigree/cli.cc
)
set
(
SPELL_MAP_SRC
src/static_data.cc
src/map-likelihood/main.cc
src/map-likelihood/cli.cc
)
set
(
SPELL_MARKER_SRC
src/input/read_mark.cc src/input/read_map.cc
# src/input/design.cc src/input/xml/xml_design.cc src/input/xml/xml_format.cc
...
...
@@ -166,6 +172,8 @@ else()
add_executable
(
spell-pedigree
${
SPELL_PEDIGREE_SRC
}
)
add_executable
(
spell-marker
${
SPELL_MARKER_SRC
}
)
add_executable
(
spell-qtl
${
SPELL_QTL_SRC
}
)
# TODO add spell-map to BUILD_FOR_DEPLOYMENT
add_executable
(
spell-map
${
SPELL_MAP_SRC
}
)
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
}
"
)
...
...
include/bayes/generalized_product.h
View file @
198b00c0
...
...
@@ -526,13 +526,13 @@ struct joint_variable_product_type {
bool
move_cursor_if_state_is_valid
(
table_descr
&
t
,
state_index_type
cursor
)
{
MSG_DEBUG
(
"table @"
<<
t
.
data
<<
" cursor="
<<
dump
(
cursor
));
//
MSG_DEBUG("table @" << t.data << " cursor=" << dump(cursor));
auto
it
=
t
.
domain_to_offset
.
find
(
t
.
local
(
cursor
));
if
(
it
==
t
.
domain_to_offset
.
end
())
{
MSG_DEBUG
(
"not found"
);
//
MSG_DEBUG("not found");
return
false
;
}
MSG_DEBUG
(
"found"
);
//
MSG_DEBUG("found");
t
.
offset
=
it
->
second
;
return
true
;
}
...
...
@@ -584,7 +584,7 @@ struct joint_variable_product_type {
d
.
bit_width
=
1
;
while
((
1ULL
<<
d
.
bit_width
)
<
d
.
size
)
{
++
d
.
bit_width
;
}
total_bits
+=
d
.
bit_width
;
MSG_DEBUG
(
"Have domain for variable #"
<<
d
.
variable_name
<<
" size "
<<
d
.
size
<<
" bit_width "
<<
d
.
bit_width
);
//
MSG_DEBUG("Have domain for variable #" << d.variable_name << " size " << d.size << " bit_width " << d.bit_width);
/*MSG_QUEUE_FLUSH();*/
d
.
index_to_label
.
reserve
(
domain
.
size
());
for
(
size_t
i
=
0
;
i
<
domain
.
size
();
++
i
)
{
...
...
@@ -603,9 +603,9 @@ struct joint_variable_product_type {
shift
+=
d
.
bit_width
;
d
.
assign
(
max_cursor
,
d
.
max
);
/*MSG_QUEUE_FLUSH();*/
MSG_DEBUG
(
"Have domain for variable #"
<<
d
.
variable_name
<<
" size "
<<
d
.
size
<<
" bit_width "
<<
d
.
bit_width
<<
" mask "
<<
boost
::
dynamic_bitset
<
uint64_t
>
(
total_bits
,
d
.
mask
)
<<
" max "
<<
boost
::
dynamic_bitset
<
uint64_t
>
(
total_bits
,
d
.
max
));
//
MSG_DEBUG("Have domain for variable #" << d.variable_name << " size " << d.size << " bit_width " << d.bit_width << " mask " << boost::dynamic_bitset<uint64_t>(total_bits, d.mask) << " max " << boost::dynamic_bitset<uint64_t>(total_bits, d.max));
/*MSG_QUEUE_FLUSH();*/
if
(
1
)
{
if
(
0
)
{
std
::
stringstream
ss
;
for
(
const
auto
&
kv
:
d
.
label_to_index
)
{
ss
<<
' '
<<
kv
.
first
<<
':'
<<
kv
.
second
;
...
...
@@ -638,9 +638,9 @@ struct joint_variable_product_type {
tab
.
variables_left
[
v
]
=
true
;
}
/*MSG_QUEUE_FLUSH();*/
MSG_DEBUG
(
"Table "
<<
*
tab
.
data
);
//
MSG_DEBUG("Table " << *tab.data);
/*MSG_QUEUE_FLUSH();*/
MSG_DEBUG
(
"Table variables "
<<
tab
.
variable_names
<<
" mask "
<<
tab
.
mask
);
//
MSG_DEBUG("Table variables " << tab.variable_names << " mask " << tab.mask);
/*MSG_QUEUE_FLUSH();*/
tab
.
offset
=
0
;
coordinates
.
resize
(
tab
.
variable_names
.
size
());
...
...
@@ -654,7 +654,7 @@ struct joint_variable_product_type {
tab
.
domain_to_offset
[
cursor
]
=
i
;
tab
.
unif
=
tab
.
unif
&&
std
::
abs
(
coef
-
tab
.
data
->
m_combination
[
i
].
coef
)
>
_EPSILON
;
}
if
(
1
)
{
if
(
0
)
{
std
::
stringstream
ss
;
ss
<<
"Table"
;
for
(
const
auto
&
x
:
tab
.
domain_to_offset
)
{
...
...
@@ -675,26 +675,26 @@ struct joint_variable_product_type {
compile_output
()
{
scoped_indent
_
(
"[output] "
);
MSG_DEBUG
(
"all output variables "
<<
output_variable_names
);
//
MSG_DEBUG("all output variables " << output_variable_names);
output_variable_names
=
output_variable_names
%
all_variable_names
;
MSG_DEBUG
(
"all_variable_names "
<<
all_variable_names
<<
" output_variable_names "
<<
output_variable_names
);
//
MSG_DEBUG("all_variable_names " << all_variable_names << " output_variable_names " << output_variable_names);
output_variables_in_use
.
resize
(
all_variable_names
.
size
(),
false
);
for
(
int
v
:
output_variable_names
)
{
size_t
index
=
std
::
find
(
all_variable_names
.
begin
(),
all_variable_names
.
end
(),
v
)
-
all_variable_names
.
begin
();
output_variables_in_use
[
index
]
=
true
;
MSG_DEBUG
(
"variable "
<<
v
<<
" has index "
<<
index
);
//
MSG_DEBUG("variable " << v << " has index " << index);
for
(
size_t
t
=
0
;
t
<
tables
.
size
();
++
t
)
{
auto
iter
=
std
::
find
(
tables
[
t
].
variable_names
.
begin
(),
tables
[
t
].
variable_names
.
end
(),
v
);
if
(
iter
!=
tables
[
t
].
variable_names
.
end
())
{
MSG_DEBUG
(
"variable "
<<
v
<<
" found in table "
<<
t
<<
" at index "
<<
(
iter
-
tables
[
t
].
variable_names
.
begin
()));
//
MSG_DEBUG("variable " << v << " found in table " << t << " at index " << (iter - tables[t].variable_names.begin()));
key_extractors
.
emplace_back
(
t
,
iter
-
tables
[
t
].
variable_names
.
begin
());
break
;
}
else
{
MSG_DEBUG
(
"variable "
<<
v
<<
" not found in table "
<<
t
);
//
} else {
//
MSG_DEBUG("variable " << v << " not found in table " << t);
}
}
}
if
(
1
)
{
if
(
0
)
{
std
::
stringstream
ss
;
ss
<<
"key extractors"
;
for
(
const
auto
&
t_ki
:
key_extractors
)
{
...
...
@@ -843,11 +843,11 @@ struct joint_variable_product_type {
operator
()
(
const
joint_variable_product_type
&
jvp
,
genotype_comb_type
::
key_list
&
key
)
const
{
auto
ki
=
key
.
begin
();
MSG_DEBUG
(
"extractors "
<<
jvp
.
key_extractors
);
//
MSG_DEBUG("extractors " << jvp.key_extractors);
for
(
const
auto
&
t_ki
:
jvp
.
key_extractors
)
{
*
ki
++
=
jvp
.
tables
[
t_ki
.
first
].
current_element
().
keys
.
keys
[
t_ki
.
second
];
}
MSG_DEBUG
(
"key: "
<<
key
);
//
MSG_DEBUG("key: " << key);
}
};
...
...
@@ -860,18 +860,18 @@ struct joint_variable_product_type {
auto
ki
=
key
.
begin
();
auto
ti
=
jvp
.
tables
.
begin
();
bn_label_type
G
=
ti
->
current_element
().
keys
.
keys
.
front
().
state
;
MSG_DEBUG
(
"Have gamete label "
<<
G
);
MSG_QUEUE_FLUSH
();
MSG_DEBUG
(
"Have output variable names "
<<
jvp
.
output_variable_names
);
//
MSG_DEBUG("Have gamete label " << G); MSG_QUEUE_FLUSH();
//
MSG_DEBUG("Have output variable names " << jvp.output_variable_names);
for
(
const
auto
&
t_ki
:
jvp
.
key_extractors
)
{
*
ki
++
=
jvp
.
tables
[
t_ki
.
first
].
current_element
().
keys
.
keys
[
t_ki
.
second
];
/**ki++ = ti->current_element().keys.keys.front();*/
}
MSG_DEBUG
(
"Have temp key "
<<
key
);
MSG_QUEUE_FLUSH
();
//
MSG_DEBUG("Have temp key " << key); MSG_QUEUE_FLUSH();
const
auto
&
p1
=
key
.
keys
[
G
.
first_allele
];
const
auto
&
p2
=
key
.
keys
[
G
.
second_allele
];
bool
f1
=
G
.
first
==
GAMETE_L
;
bool
f2
=
G
.
second
==
GAMETE_L
;
MSG_DEBUG
(
"Have p1 "
<<
p1
<<
" p2 "
<<
p2
<<
" f1 "
<<
f1
<<
" f2 "
<<
f2
);
MSG_QUEUE_FLUSH
();
//
MSG_DEBUG("Have p1 " << p1 << " p2 " << p2 << " f1 " << f1 << " f2 " << f2); MSG_QUEUE_FLUSH();
--
ki
;
*
ki
=
{
jvp
.
output_variable_names
.
back
(),
{
f1
?
p1
.
state
.
first
:
p1
.
state
.
second
,
...
...
@@ -879,7 +879,7 @@ struct joint_variable_product_type {
f1
?
p1
.
state
.
first_allele
:
p1
.
state
.
second_allele
,
f2
?
p2
.
state
.
first_allele
:
p2
.
state
.
second_allele
}
};
MSG_DEBUG
(
"Have complete key "
<<
key
);
MSG_QUEUE_FLUSH
();
//
MSG_DEBUG("Have complete key " << key); MSG_QUEUE_FLUSH();
}
};
...
...
@@ -892,25 +892,25 @@ struct joint_variable_product_type {
auto
ki
=
key
.
begin
();
auto
ti
=
jvp
.
tables
.
begin
();
bn_label_type
G
=
ti
->
current_element
().
keys
.
keys
.
front
().
state
;
MSG_DEBUG
(
"Have gamete label "
<<
G
);
MSG_QUEUE_FLUSH
();
MSG_DEBUG
(
"Have output variable names "
<<
jvp
.
output_variable_names
);
//
MSG_DEBUG("Have gamete label " << G); MSG_QUEUE_FLUSH();
//
MSG_DEBUG("Have output variable names " << jvp.output_variable_names);
for
(
const
auto
&
t_ki
:
jvp
.
key_extractors
)
{
*
ki
++
=
jvp
.
tables
[
t_ki
.
first
].
current_element
().
keys
.
keys
[
t_ki
.
second
];
/**ki++ = ti->current_element().keys.keys.front();*/
}
MSG_DEBUG
(
"Have temp key "
<<
key
);
MSG_QUEUE_FLUSH
();
//
MSG_DEBUG("Have temp key " << key); MSG_QUEUE_FLUSH();
const
auto
&
p1
=
key
.
keys
[
1
+
G
.
first_allele
];
const
auto
&
p2
=
key
.
keys
[
1
+
G
.
second_allele
];
bool
f1
=
G
.
first
==
GAMETE_L
;
bool
f2
=
G
.
second
==
GAMETE_L
;
MSG_DEBUG
(
"Have p1 "
<<
p1
<<
" p2 "
<<
p2
<<
" f1 "
<<
f1
<<
" f2 "
<<
f2
);
MSG_QUEUE_FLUSH
();
//
MSG_DEBUG("Have p1 " << p1 << " p2 " << p2 << " f1 " << f1 << " f2 " << f2); MSG_QUEUE_FLUSH();
*
ki
=
{
-
jvp
.
output_variable_names
.
front
(),
{
f1
?
p1
.
state
.
first
:
p1
.
state
.
second
,
f2
?
p2
.
state
.
first
:
p2
.
state
.
second
,
f1
?
p1
.
state
.
first_allele
:
p1
.
state
.
second_allele
,
f2
?
p2
.
state
.
first_allele
:
p2
.
state
.
second_allele
}
};
MSG_DEBUG
(
"Have complete key "
<<
key
);
MSG_QUEUE_FLUSH
();
//
MSG_DEBUG("Have complete key " << key); MSG_QUEUE_FLUSH();
}
};
...
...
@@ -1076,7 +1076,7 @@ struct joint_variable_product_type {
/*MSG_DEBUG(key.size()); MSG_QUEUE_FLUSH();*/
state_index_type
z
;
z
.
reset
();
MSG_DEBUG
(
"Initializing cursor"
);
//
MSG_DEBUG("Initializing cursor");
state_index_type
cursor
=
next_state
(
z
),
last
=
cursor
;
...
...
include/bayes/graphnode2.h
View file @
198b00c0
...
...
@@ -1070,17 +1070,31 @@ struct factor_graph_type : public recursive_graph_type<factor_graph_type> {
MSG_DEBUG
(
"keep_gametes "
<<
std
::
boolalpha
<<
keep_gametes
);
if
(
keep_gametes
)
{
gametes
.
emplace_back
(
-
spawnling
);
io
[
-
spawnling
]
=
Output
;
add_variables_to_node
(
n
,
{
-
spawnling
});
}
/* Add domain for spawnling */
MSG_DEBUG
(
"node_domains[n] "
<<
node_domains
[
n
]);
MSG_DEBUG
(
"spawnling vec "
<<
var_vec
{
spawnling
});
auto
msg
=
(
node_domains
[
n
]
%
var_vec
{
spawnling
});
auto
dom
=
msg
.
front
();
for
(
auto
&
e
:
dom
)
{
e
.
coef
=
1
;
{
auto
msg
=
(
node_domains
[
n
]
%
var_vec
{
spawnling
});
auto
dom
=
msg
.
front
();
for
(
auto
&
e
:
dom
)
{
e
.
coef
=
1
;
}
propagate_spawnling_domain
({
spawnling
},
dom
);
}
if
(
keep_gametes
)
{
auto
msg
=
(
node_domains
[
n
]
%
var_vec
{
-
spawnling
});
auto
dom
=
msg
.
front
();
for
(
auto
&
e
:
dom
)
{
e
.
coef
=
1
;
}
propagate_spawnling_domain
({
-
spawnling
},
dom
);
}
propagate_spawnling_
domain
(
{
spawnling
}
,
dom
);
/*
MSG_DEBUG("
Domain for spawnling #" << spawnling
<< "
:
" <<
domains[{spawnling}]);*/
MSG_DEBUG
(
"Domain for spawnling #"
<<
spawnling
<<
": "
<<
domain
s
[
{
spawnling
}
]
);
MSG_DEBUG
(
"
Variables for corresponding node #"
<<
n
<<
"
are
"
<<
variables_of
(
n
));
}
/*MSG_DEBUG("COMPUTED DOMAIN " << node_domains[n]);*/
dump_node
(
n
);
...
...
@@ -1138,6 +1152,8 @@ struct factor_graph_type : public recursive_graph_type<factor_graph_type> {
compute_node_domain
(
n
);
}
}
std
::
sort
(
gametes
.
begin
(),
gametes
.
end
());
/*MSG_DEBUG("domains computes " << node_domain_computed);*/
...
...
@@ -1622,24 +1638,47 @@ struct instance_type {
compile_incoming_message_domains
(
output_operation
);
}
else
{
for
(
node_index_type
output
:
g
->
nodes
())
{
/*MSG_DEBUG("active node " << output << " has io " << g->node_io(output));*/
if
(
g
->
node_is_interface
(
output
)
&&
g
->
node_io
(
output
)
&
Output
)
{
V
.
output_operations
.
emplace_back
();
auto
&
output_operation
=
V
.
output_operations
.
back
();
/*output_operation.output = g->variables_of(output);*/
for
(
variable_index_type
v
:
g
->
variables_of
(
output
))
{
if
(
g
->
var_is_output
(
v
))
{
output_operation
.
output
.
push_back
(
v
);
MSG_DEBUG
(
"active node "
<<
output
<<
" has io "
<<
g
->
node_io
(
output
)
<<
" and variables "
<<
g
->
variables_of
(
output
));
if
(
g
->
node_is_interface
(
output
))
{
if
(
g
->
node_io
(
output
)
&
Output
)
{
V
.
output_operations
.
emplace_back
();
auto
&
output_operation
=
V
.
output_operations
.
back
();
/*output_operation.output = g->variables_of(output);*/
for
(
variable_index_type
v
:
g
->
variables_of
(
output
))
{
MSG_DEBUG
(
"Output variable "
<<
v
<<
"?"
);
if
(
g
->
var_is_output
(
v
))
{
output_operation
.
output
.
push_back
(
v
);
}
}
output_operation
.
receiver
=
0
;
output_operation
.
emitter
=
output
;
for
(
node_index_type
i
:
g
->
all_nei
(
output
))
{
size_t
m
=
get_message_index
(
i
,
output
);
get_message_order_rec
(
m
,
message_visited
,
V
.
operations
);
output_operation
.
incoming_messages
.
push_back
(
m
);
}
compile_incoming_message_domains
(
output_operation
);
}
output_operation
.
receiver
=
0
;
output_operation
.
emitter
=
output
;
for
(
node_index_type
i
:
g
->
all_nei
(
output
))
{
size_t
m
=
get_message_index
(
i
,
output
);
get_message_order_rec
(
m
,
message_visited
,
V
.
operations
);
output_operation
.
incoming_messages
.
push_back
(
m
);
}
else
{
MSG_DEBUG
(
"node variables "
<<
g
->
variables_of
(
output
));
MSG_DEBUG
(
"gametes "
<<
g
->
gametes
);
auto
gam
=
g
->
variables_of
(
output
)
%
g
->
gametes
;
if
(
gam
.
size
())
{
MSG_DEBUG
(
"node has gametes "
<<
gam
);
for
(
variable_index_type
v
:
gam
)
{
V
.
output_operations
.
emplace_back
();
auto
&
output_operation
=
V
.
output_operations
.
back
();
output_operation
.
output
.
push_back
(
v
);
output_operation
.
receiver
=
0
;
output_operation
.
emitter
=
output
;
for
(
node_index_type
i
:
g
->
all_nei
(
output
))
{
size_t
m
=
get_message_index
(
i
,
output
);
get_message_order_rec
(
m
,
message_visited
,
V
.
operations
);
output_operation
.
incoming_messages
.
push_back
(
m
);
}
compile_incoming_message_domains
(
output_operation
);
}
}
compile_incoming_message_domains
(
output_operation
);
}
}
}
...
...
include/bayes/output.h
View file @
198b00c0
...
...
@@ -714,7 +714,7 @@ struct rw_base : public rw_any<rw_base> {
struct
gamete_LV_database
{
/* IND.ped / MARK => 1 or 2 gametes (.second will be empty when cross type is DH) */
std
::
map
<
int
,
std
::
map
<
std
::
string
,
std
::
pair
<
Eigen
::
Vector2d
,
Eigen
::
Vector2d
>
>>
data
;
std
::
map
<
int
,
std
::
map
<
std
::
string
,
Eigen
::
Vector2d
>>
data
;
std
::
map
<
double
,
Eigen
::
Matrix2d
>
cache
;
gamete_LV_database
()
:
data
()
{}
...
...
@@ -723,7 +723,7 @@ struct gamete_LV_database {
get_TR
(
double
dist
)
{
auto
it
=
cache
.
find
(
dist
);
if
(
it
!
=
cache
.
end
())
{
if
(
it
=
=
cache
.
end
())
{
double
s
=
.5
+
exp
(
-
dist
*
.02
)
*
.5
;
double
r
=
1.
-
s
;
Eigen
::
Matrix2d
ret
;
...
...
@@ -737,17 +737,19 @@ struct gamete_LV_database {
void
add_gametes
(
const
std
::
string
&
mark
,
int
ind
,
const
VectorXd
&
lv
,
bool
dh
)
{
std
::
pair
<
Eigen
::
Vector2d
,
Eigen
::
Vector2d
>
g
;
Eigen
::
Vector2d
g1
,
g2
;
MSG_DEBUG
(
"add_gametes mark="
<<
mark
<<
" ind="
<<
ind
<<
" lv="
<<
lv
.
transpose
()
<<
" dh="
<<
std
::
boolalpha
<<
dh
);
if
(
!
dh
)
{
g
.
first
(
0
)
=
lv
(
0
)
+
lv
(
1
);
g
.
first
(
1
)
=
lv
(
2
)
+
lv
(
3
);
g
.
second
(
0
)
=
lv
(
0
)
+
lv
(
2
);
g
.
second
(
1
)
=
lv
(
1
)
+
lv
(
3
);
g
1
(
0
)
=
lv
(
0
)
+
lv
(
1
);
g
1
(
1
)
=
lv
(
2
)
+
lv
(
3
);
g
2
(
0
)
=
lv
(
0
)
+
lv
(
2
);
g
2
(
1
)
=
lv
(
1
)
+
lv
(
3
);
}
else
{
g
.
first
=
lv
;
g
.
second
(
0
)
=
g
.
second
(
1
)
=
0
;
g
1
=
lv
;
g
2
(
0
)
=
g2
(
1
)
=
0
;
}
data
[
ind
][
mark
]
=
g
;
data
[
ind
][
mark
]
=
g1
;
data
[
-
ind
][
mark
]
=
g2
;
}
double
...
...
@@ -766,20 +768,25 @@ struct gamete_LV_database {
TR
.
push_back
(
get_TR
(
d
));
}
double
ret
=
0
;
MSG_DEBUG
(
"DATA"
<<
std
::
endl
<<
data
);
for
(
auto
&
kv
:
data
)
{
int
ind
=
kv
.
first
;
// int ind = kv.first;
MSG_DEBUG
(
"ON GAMETES #"
<<
kv
.
first
);
auto
&
LV
=
kv
.
second
;
auto
name_i
=
marker_names
.
begin
();
auto
&
gam0
=
LV
[
*
name_i
++
];
Eigen
::
Vector2d
accum
;
accum
<<
1
,
1
;
Eigen
::
Vector2d
accum
=
gam0
;
//
accum << 1, 1;
std
::
vector
<
double
>
norms
;
double
norm_accum
=
0
;
norms
.
reserve
(
distances
.
size
());
norms
.
push_back
(
.
5
);
norms
.
push_back
(
1
.
);
norm_accum
-=
log
(
norms
.
back
());
for
(
size_t
i
=
0
;
i
<
TR
.
size
();
++
i
,
++
name_i
)
{
accum
=
accum
.
transpose
()
*
TR
[
i
]
*
LV
[
*
name_i
].
first
.
asDiagonal
();
MSG_DEBUG
(
"accum"
<<
std
::
endl
<<
accum
);
MSG_DEBUG
(
"TR["
<<
i
<<
']'
<<
std
::
endl
<<
TR
[
i
]);
MSG_DEBUG
(
"LV["
<<
(
*
name_i
)
<<
"] = "
<<
LV
[
*
name_i
].
transpose
());
accum
=
accum
.
transpose
()
*
TR
[
i
]
*
LV
[
*
name_i
].
asDiagonal
();
norms
.
push_back
(
accum
.
lpNorm
<
2
>
());
if
(
norms
.
back
()
==
0
)
{
MSG_ERROR
(
"Have a zero vector, aborting."
,
""
);
...
...
include/map-likelihood/cli.h
0 → 100644
View file @
198b00c0
/* Spell-QTL Software suite for the QTL analysis of modern datasets.
* Copyright (C) 2016,2017 Damien Leroux <damien.leroux@inra.fr>, Sylvain Jasson <sylvain.jasson@inra.fr>
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef _SPELL_MAPQA_CLI_H_
#define _SPELL_MAPQA_CLI_H_
#include
"error.h"
// #include <map>
// #include <vector>
/*#include "generation_rs.h"*/
#include
"commandline.h"
/*#include "bn.h"*/
/*#include "factor_var4.h"*/
// #include "file.h"
#include
"data/chromosome.h"
struct
mapqa_settings_t
{
/* Configuration */
std
::
string
name
;
std
::
string
work_directory
;
chromosome
group
;
mapqa_settings_t
()
:
name
()
,
work_directory
(
"."
)
{}
static
mapqa_settings_t
*
from_args
(
int
argc
,
const
char
**
argv
);
};
void
print_usage
();
#endif
src/bayes/jobs.cc
View file @
198b00c0
...
...
@@ -385,10 +385,21 @@ job_registry = {
gen
->
name
)
!=
settings
->
output_generations
.
end
())
{
output_prob
[
node
]
=
state_prob
[
node
];
}
if
(
settings
->
output_mode
&
bn_settings_t
::
OutputGametes
)
{
if
(
(
settings
->
output_mode
&
bn_settings_t
::
OutputGametes
)
&&
!
settings
->
pedigree
.
items
[
ind
-
1
].
is_ancestor
())
{
auto
&
mat
=
marginals
[
-
variable
];
MSG_DEBUG
(
"marginals for gamete "
<<
mat
);
//gamete_prob[ind] = mat.col(0);
std
::
map
<
label_type
,
double
>
prob
;
for
(
const
auto
&
kv
:
mat
)
{
prob
[{
kv
.
first
.
first
(),
kv
.
first
.
second
()}]
=
kv
.
second
;
}
if
(
settings
->
pedigree
.
items
[
ind
-
1
].
is_dh
())
{
gamete_prob
[
ind
].
resize
(
2
);
gamete_prob
[
ind
]
<<
prob
[{
GAMETE_L
,
0
}],
prob
[{
GAMETE_R
,
0
}];
}
else
{
gamete_prob
[
ind
].
resize
(
4
);
gamete_prob
[
ind
]
<<
prob
[{
GAMETE_L
,
GAMETE_L
}],
prob
[{
GAMETE_L
,
GAMETE_R
}],
prob
[{
GAMETE_R
,
GAMETE_L
}],
prob
[{
GAMETE_R
,
GAMETE_R
}];
}
MSG_DEBUG
(
"GAMETE #"
<<
ind
<<
" prob "
<<
gamete_prob
[
ind
-
1
].
transpose
());
}
}
...
...
@@ -402,176 +413,13 @@ job_registry = {
rw_base
()(
ofs
,
output_prob
);
if
(
settings
->
output_mode
&
bn_settings_t
::
OutputGametes
)
{
ofile
og
(
MESSAGE
(
"gametes."
<<
filename
));
ofile
og
(
settings
->
job_filename
(
"gametes-LV"
,
mark
));
rw_base
()(
og
,
gamete_prob
);
}
}
return
true
;
/*for (const auto& pop_obs: settings->observed_mark) {*/
/*MSG_DEBUG("Obs in generation " << pop_obs.first);*/
/*size_t n = 0;*/
/*for (size_t i: settings->pedigree.individuals_by_generation_name.find(pop_obs.first)->second) {*/
/*size_t node = settings->pedigree.ind(i);*/
/*size_t igen = settings->pedigree.node_generations[node];*/
/*const auto& gen = settings->pedigree.generations[igen];*/
/*gen->labels;*/
/*}*/
/*}*/
/**/
/*std::map<size_t, std::map<label_type, double>> prob_by_ind;*/
/*for (const auto& kv: marginals) {*/
/*genotype_comb_type::key_type key = kv.first.keys[0];*/
/*prob_by_ind[key.parent][{key.state.first, key.state.second}] += kv.second;*/
/*}*/
/*std::map<size_t, std::vector<double>> state_prob;*/
#if 0
MSG_DEBUG("MARGINALS:");
for (const auto& t: output) {
var_vec varz = get_parents(t);
for (variable_index_type v: varz) {
marginals[v] = project(t, {v}, {});
}
}
#if 0
std::vector<std::vector<double>> gen_dispatch;
for (const auto& gen: settings->pedigree.generations) {
gen_dispatch.emplace_back();
if (!gen) { continue; }
std::map<label_type, double> label_weights;
for (size_t i = 0; i < gen->labels.size(); ++i) {
label_weights[gen->labels[i]] += gen->stat_dist[i];
}
gen_dispatch.back().resize(gen->labels.size());
for (size_t i = 0; i < gen->labels.size(); ++i) {
gen_dispatch.back()[i] = gen->stat_dist[i] / label_weights[gen->labels[i]];
/*gen_dispatch.back()[i] = 1 / label_weights[gen->labels[i]];*/
}
}
std::map<size_t, std::map<label_type, double>> prob_by_ind;
for (const auto& kv: marginals) {
genotype_comb_type::key_type key = kv.first.keys[0];
prob_by_ind[key.parent][{key.state.first, key.state.second}] += kv.second;
}
std::map<size_t, std::vector<double>> state_prob;
auto
compute_geno_prob = [&] (const std::vector<gencomb_type>& lc, const std::vector<double>& weights, std::function<double(size_t, size_t)> get_state_prob)
{
std::vector<double> ret(lc.size());
double norm = 0;
for (size_t i = 0; i < ret.size(); ++i) {
ret[i] = 0;
for (const auto& e: lc[i]) {
double coef = e.coef;
for (const auto& k: e.keys) {
/*coef *= state_prob[k.parent][k.state];*/
coef *= get_state_prob(k.parent, k.state);
}
ret[i] += coef;
}
MSG_DEBUG("unweighted_state[" << i << "] = " << ret[i]);
ret[i] *= weights[i];
norm += ret[i];
}
if (norm > 0) {
norm = 1. / norm;
} else {
CREATE_MESSAGE(msg_channel::Out, MESSAGE(std::endl));
MSG_ERROR("The observations for marker " << marker_name << " are inconsistent. Aborting computation.", "Check your genotype data files.");
ret.clear();
}
for (double& d: ret) { d *= norm; }
return ret;
};
#endif
std::map<size_t, std::map<label_type, double>> prob_by_ind;
for (const auto& v_prob: marginals) {
int v = v_prob.first;
auto& pbi = prob_by_ind[v];
/*MSG_DEBUG("### v=" << v << " pbi " << pbi);*/
/*MSG_QUEUE_FLUSH();*/
for (const auto& e: v_prob.second) {
/*MSG_DEBUG("[element] " << e);*/
/*MSG_QUEUE_FLUSH();*/
const genotype_comb_type::key_type& key = e.keys.keys[0];
/*MSG_DEBUG("[geno prob for #" << v << "] " << key << " += " << e.coef);*/
/*MSG_QUEUE_FLUSH();*/
pbi[{key.state.first, key.state.second}] += e.coef;
/*MSG_DEBUG("--");*/
/*MSG_QUEUE_FLUSH();*/
}
}
#if 0
for (size_t i = 0; i < state.size(); ++i) {
variable_index_type v = marginals[i].back().variables.front();
/*auto dom = fg.get_domain(v);*/
for (const auto& e: state[i]->front()) {
const genotype_comb_type::key_type& key = e.keys.keys[0];
MSG_DEBUG("[geno prob for #" << v << "] " << key << " += " << e.coef);
prob_by_ind[v][{key.state.first, key.state.second}] += e.coef;
}
}
#endif
/*for (const auto& kv: state) {*/
/*genotype_comb_type::key_type key = kv.first.keys[0];*/
/*prob_by_ind[key.parent][{key.state.first, key.state.second}] += kv.second;*/
/*}*/
/*for (const auto& kv: prob_by_ind) {*/
/*std::stringstream msg;*/
/*msg << "IND PROB #" << kv.first;*/
/*for (const auto& lp: kv.second) {*/
/*msg << ' ' << lp.first << '=' << lp.second;*/
/*}*/
/*MSG_DEBUG(msg.str());*/
/*}*/
std::map<size_t, std::vector<double>> state_prob;
/*for (size_t ind: settings->pedigree.tree.m_ind_number_to_node_number) {*/
const auto& in3 = settings->pedigree.tree.m_ind_number_to_node_number;
for (const auto& ind_probs: prob_by_ind) {
size_t ind = ind_probs.first;
size_t I = in3[settings->pedigree.m_id[ind]];