Commit c58471bf authored by Damien Leroux's avatar Damien Leroux
Browse files

New tool to prepare for spell-marker and spell-qtl.

Name is spell-pedigree.
parent 19f11ef9
......@@ -335,6 +335,7 @@ void write_vector(std::ofstream& ofs, const std::vector<V>& vec, WRITE_ELEM_FUNC
/** LABEL_TYPE **/
inline
label_type read_label(std::ifstream& ifs)
{
char f, s;
......@@ -342,6 +343,7 @@ label_type read_label(std::ifstream& ifs)
return {f, s};
}
inline
void write_label(std::ofstream& ofs, const label_type& l)
{
ofs << l.first() << l.second();
......@@ -350,7 +352,8 @@ void write_label(std::ofstream& ofs, const label_type& l)
/** GENO_MATRIX **/
inline void write_geno_matrix(std::ofstream& ofs, const geno_matrix& mat)
inline
void write_geno_matrix(std::ofstream& ofs, const geno_matrix& mat)
{
write_fourcc(ofs, "SGEM");
write_str(ofs, mat.name);
......@@ -387,11 +390,13 @@ void read_geno_matrix(std::ifstream& ifs, geno_matrix& mat)
}
inline void write_geno_matrix(std::ofstream& ofs, const geno_matrix* ptr)
inline
void write_geno_matrix(std::ofstream& ofs, const geno_matrix* ptr)
{
write_geno_matrix(ofs, *ptr);
}
inline
geno_matrix* read_geno_matrix(std::ifstream& ifs)
{
geno_matrix* ret = new geno_matrix();
......
......@@ -234,6 +234,7 @@ struct msg_handler_t {
std::set<std::string> workarounds;
int count;
int debug_indent;
bool debug_enabled;
std::vector<std::function<void()>> hooks;
message_queue queue;
......@@ -243,7 +244,7 @@ struct msg_handler_t {
const char* normal() { return color ? _NORMAL : ""; }
state_t()
: color(!!isatty(fileno(stdout))), workarounds(), count(0), debug_indent(0), hooks()
: color(!!isatty(fileno(stdout))), workarounds(), count(0), debug_indent(0), debug_enabled(true), hooks()
, queue()
{
/*std::cout << "Message handler instance created." << std::endl;*/
......@@ -285,6 +286,8 @@ struct msg_handler_t {
return w.ws_col;
}
static bool& debug_enabled() { return instance().debug_enabled; }
static lock_type mutex;
};
......@@ -386,8 +389,10 @@ void message_queue::run()
#define MSG_DEBUG(_msg_expr_) \
do {\
CREATE_MESSAGE(msg_channel::Log, MESSAGE(_msg_expr_ << std::endl));\
} while (0)
if (msg_handler_t::debug_enabled()) {\
CREATE_MESSAGE(msg_channel::Log, MESSAGE(_msg_expr_ << std::endl));\
}\
} while (0);
#define MSG_QUEUE_FLUSH() do { \
CREATE_MESSAGE(msg_channel::Log, _DBG_SYNC_MARK_S); \
......@@ -431,6 +436,9 @@ void message_queue::run()
#define MSG_DEBUG_DEDENT CREATE_MESSAGE(msg_channel::Log, _DBG_DEDENT_MARK_S)
/*#ifndef _SPELL_ERROR_H_MSGQ_CATCH_SIGSEG_*/
/*#define _SPELL_ERROR_H_MSGQ_CATCH_SIGSEG_*/
inline
void message_queue::catch_SIGSEG(int)
{
MSG_DEBUG("**SEGFAULT DETECTED**");
......@@ -438,7 +446,7 @@ void message_queue::catch_SIGSEG(int)
signal(SIGSEGV, old_sighandler());
raise(SIGSEGV);
}
/*#endif*/
struct scoped_indent {
scoped_indent() { MSG_DEBUG_INDENT; }
......
......@@ -1759,6 +1759,7 @@ struct experimental_lumper {
}
inline
std::map<label_type, VectorXd>
geno_matrix::LV_map() const
{
......
......@@ -51,11 +51,13 @@ bool check_file(const std::string& path, bool req_directory, bool req_writable,
return false;
}
if (req_directory && !fs.is_dir) {
if (display) {
MSG_ERROR(path << " is not a directory", "Specify the path to directory, not a file [" << path << ']');
if (req_directory) {
if (!fs.is_dir) {
if (display) {
MSG_ERROR(path << " is not a directory", "Specify the path to directory, not a file [" << path << ']');
}
return false;
}
return false;
} else if (!fs.is_file) {
if (display) {
MSG_ERROR(path << " is not a regular file or a symbolic link to a regular file", "Specify the path to a file, not a directory [" << path << ']');
......
......@@ -133,6 +133,7 @@ read_csv(const std::string& pedigree_file, char field_sep=';');
typedef std::map<size_t, size_t> ancestor_list_type;
inline
ancestor_list_type reentrants(const ancestor_list_type& a)
{
ancestor_list_type ret;
......@@ -145,6 +146,7 @@ ancestor_list_type reentrants(const ancestor_list_type& a)
}
inline
ancestor_list_type operator + (const ancestor_list_type& a1, const ancestor_list_type& a2)
{
ancestor_list_type ret(a1);
......@@ -155,6 +157,7 @@ ancestor_list_type operator + (const ancestor_list_type& a1, const ancestor_list
}
inline
ancestor_list_type operator / (const ancestor_list_type& a, const ancestor_list_type& restr)
{
ancestor_list_type ret;
......@@ -168,6 +171,7 @@ ancestor_list_type operator / (const ancestor_list_type& a, const ancestor_list_
}
inline
ancestor_list_type operator % (const ancestor_list_type& a, const ancestor_list_type& restr)
{
ancestor_list_type ret;
......@@ -180,6 +184,7 @@ ancestor_list_type operator % (const ancestor_list_type& a, const ancestor_list_
}
inline
ancestor_list_type operator - (const ancestor_list_type& a, const ancestor_list_type& restr)
{
ancestor_list_type ret;
......@@ -195,6 +200,7 @@ ancestor_list_type operator - (const ancestor_list_type& a, const ancestor_list_
}
inline
ancestor_list_type operator * (const ancestor_list_type& a, size_t weight)
{
ancestor_list_type ret;
......@@ -205,6 +211,7 @@ ancestor_list_type operator * (const ancestor_list_type& a, size_t weight)
}
inline
std::ostream& operator << (std::ostream& os, const ancestor_list_type& a)
{
auto i = a.begin();
......@@ -242,6 +249,7 @@ label_type operator * (label_type a, label_type b)
#define SELECT_A(__p, __b) ((__b) == GAMETE_R ? (__p).second_allele : (__p).first_allele)
inline
bn_label_type operator * (bn_label_type a, bn_label_type b)
{
bn_label_type ret;
......@@ -650,7 +658,7 @@ struct pedigree_type {
int n = tree.add_node();
MSG_DEBUG("node=" << n << " ind=" << ind);
compute_generation(generation_name, n);
/*compute_LC(n);*/
compute_LC(n);
MSG_DEBUG_DEDENT;
}
break;
......@@ -662,7 +670,8 @@ struct pedigree_type {
int g1 = spawn_gamete("M", tree.ind2node(p1));
int n = tree.add_node(g1, g1);
MSG_DEBUG("node=" << n << " ind=" << ind);
/*compute_generation(generation_name, n);*/
compute_generation(generation_name, n);
compute_LC(n);
/*compute_data_for_bn(n);*/
MSG_DEBUG_DEDENT;
}
......@@ -681,7 +690,8 @@ struct pedigree_type {
int g2 = spawn_gamete("F", n2);
int n = tree.add_node(g1, g2);
MSG_DEBUG("node=" << n << " ind=" << ind);
/*compute_generation(generation_name, n);*/
compute_generation(generation_name, n);
compute_LC(n);
/*compute_data_for_bn(n);*/
MSG_DEBUG_DEDENT;
}
......@@ -1449,32 +1459,37 @@ struct pedigree_type {
template <typename STREAM_TYPE>
void load_save_impl(STREAM_TYPE& fs)
void load_save_impl(STREAM_TYPE& fs, bool skip_gen)
{
rw_comb<size_t, size_t> comb_rw;
rw_comb<size_t, bn_label_type> bn_comb_rw;
rw_tree tree_rw;
rw_base rw;
if (rw.fourcc(fs, "PDAT")) {
return;
}
tree_rw(fs, tree);
rw(fs, node_generations);
rw(fs, generations);
rw(fs, ancestor_letters);
rw(fs, generation_names);
rw(fs, geno_matrix_by_generation_name);
rw(fs, individuals_by_generation_name);
rw(fs, max_states);
rw(fs, n_alleles);
rw(fs, filename);
comb_rw(fs, LC);
bn_comb_rw(fs, factor_messages);
if (!skip_gen) {
if (rw.fourcc(fs, "PGEN")) {
return;
}
rw(fs, node_generations);
rw(fs, generations);
}
}
void load(const std::string& filename)
void load(const std::string& filename, bool skip_gen)
{
std::ifstream ifs(filename);
load_save_impl(ifs);
load_save_impl(ifs, skip_gen);
generation_name_by_individual.clear();
for (int i = 1; i < tree.next_ind_idx(); ++i) {
generation_name_by_individual[i] = &get_gen(i)->name;
......@@ -1484,19 +1499,22 @@ struct pedigree_type {
void save(const std::string& filename)
{
std::ofstream ofs(filename);
load_save_impl(ofs);
load_save_impl(ofs, false);
}
};
inline
pedigree_type
read_pedigree(const std::string& filename, char sep=';')
read_pedigree(const std::string& filename, char sep=';', size_t max_states=std::numeric_limits<size_t>::max())
{
auto items = read_csv(filename, sep);
pedigree_type ret;
ret.max_states = max_states;
std::map<int, int> id;
size_t cur = 0;
for (const auto& i: items) {
CREATE_MESSAGE(msg_channel::Out, MESSAGE("\rProcessing entry " << (++cur) << "/" << items.size() << "..."));
if (i.is_ancestor()) {
id[i.id] = ret.ancestor(i.gen_name);
} else if (i.is_dh()) {
......@@ -1510,6 +1528,7 @@ read_pedigree(const std::string& filename, char sep=';')
return {};
}
}
CREATE_MESSAGE(msg_channel::Out, MESSAGE(std::endl));
return ret;
}
......
#ifndef _SPELL_PEDIGREE_SETTINGS_H_
#define _SPELL_PEDIGREE_SETTINGS_H_
#include "error.h"
#include "commandline.h"
#include "pedigree.h"
struct ped_settings_t {
std::string prg_name;
std::vector<std::string> command_line;
std::string pedigree_filename;
char csv_sep;
std::string work_directory;
size_t max_states;
ped_settings_t()
: prg_name()
, command_line()
, pedigree_filename()
, csv_sep(';')
, work_directory()
, max_states(std::numeric_limits<decltype(max_states)>::max())
{}
static ped_settings_t* from_args(int argc, const char** argv);
};
#endif
......@@ -10,6 +10,7 @@
typedef std::map<int, int> ancestor_node_list_type;
inline
ancestor_node_list_type reentrants(const ancestor_node_list_type& a)
{
ancestor_node_list_type ret;
......@@ -22,6 +23,7 @@ ancestor_node_list_type reentrants(const ancestor_node_list_type& a)
}
inline
ancestor_node_list_type operator + (const ancestor_node_list_type& a1, const ancestor_node_list_type& a2)
{
ancestor_node_list_type ret(a1);
......@@ -32,6 +34,7 @@ ancestor_node_list_type operator + (const ancestor_node_list_type& a1, const anc
}
inline
ancestor_node_list_type operator / (const ancestor_node_list_type& a, const ancestor_node_list_type& restr)
{
ancestor_node_list_type ret;
......@@ -45,6 +48,7 @@ ancestor_node_list_type operator / (const ancestor_node_list_type& a, const ance
}
inline
ancestor_node_list_type operator % (const ancestor_node_list_type& a, const ancestor_node_list_type& restr)
{
ancestor_node_list_type ret;
......@@ -57,6 +61,7 @@ ancestor_node_list_type operator % (const ancestor_node_list_type& a, const ance
}
inline
ancestor_node_list_type operator - (const ancestor_node_list_type& a, const ancestor_node_list_type& restr)
{
ancestor_node_list_type ret;
......@@ -72,6 +77,7 @@ ancestor_node_list_type operator - (const ancestor_node_list_type& a, const ance
}
inline
ancestor_node_list_type operator * (const ancestor_node_list_type& a, int weight)
{
ancestor_node_list_type ret;
......@@ -82,6 +88,7 @@ ancestor_node_list_type operator * (const ancestor_node_list_type& a, int weight
}
inline
std::ostream& operator << (std::ostream& os, const ancestor_node_list_type& a)
{
auto i = a.begin();
......
......@@ -333,16 +333,19 @@ struct transposed_permutation_type : public permutation_base<transposed_permutat
const std::vector<int>& transposed_map() const { return ptr->m_map; }
};
inline
permutation_type& permutation_type::operator = (const transposed_permutation_type& p) { m_map = p.map(); m_transposed = p.transposed_map(); return *this; }
/*static*/
inline
typename transposition_result<permutation_type>::return_type
transposition_result<permutation_type>::
result(const permutation_type* p) { return {p}; }
/*static*/
inline
typename transposition_result<transposed_permutation_type>::return_type
transposition_result<transposed_permutation_type>::
result(const transposed_permutation_type* p) { return *p->ptr; }
......
......@@ -996,8 +996,15 @@ int main(int argc, char** argv)
MSG_DEBUG("G");
size_t G = ped.crossing("G", F, E);
MSG_DEBUG("prout");
MSG_DEBUG(ped.get_gen(G)->lim_inf());
/*MSG_DEBUG("prout");*/
/*MSG_DEBUG(ped.get_gen(G)->lim_inf());*/
factor_graph fg(ped, 1);
MSG_DEBUG("" << fg);
bn_message_type marginals;
double delta = fg.run(marginals); (void) delta;
MSG_DEBUG("Marginal probabilities:");
MSG_DEBUG("" << marginals);
}
/* 16 */
......
......@@ -27,6 +27,11 @@ current_all_files = None
all_include_files = set()
try:
os.mkdir('include-graphs')
except OSError, ose:
pass
for f in os.popen(all_includes_cmd).xreadlines():
all_include_files.add(os.path.relpath(os.path.realpath(f.strip())))
......@@ -111,7 +116,7 @@ def print_cluster(out, name, contents):
for curfile in all_files:
safe = curfile.replace('.', '').replace('/', '-')
with open('includes-' + safe + '.dot', 'w') as out:
with open('include-graphs/includes-' + safe + '.dot', 'w') as out:
print >> out, 'digraph "includes-' + safe + '" {'
pathitems = [clean_path(f) for f in all_files[curfile]]
......@@ -141,7 +146,8 @@ for curfile in all_files:
print >> out, "}"
# now invoke dot
os.system('for i in includes-*dot; do dot -Tpng $i -o ${i%%dot}png; done')
os.system('for i in include-graphs/includes-*dot; do '
+ 'dot -Tpng $i -o ${i%%dot}png; done')
used_inc = set()
......@@ -151,6 +157,7 @@ for k, v in all_files.iteritems():
used_inc.add(f)
#print used_inc
print "UNUSED INCLUDES"
for unused in (all_include_files - set(used_inc)):
print unused
with open('include-graphs/unused.txt', 'w') as out:
print >> out, "UNUSED INCLUDES"
for unused in (all_include_files - set(used_inc)):
print >> out, unused
include Makefile.conf
#TARGET=spell-qtl.experimental
TARGET=spell-pedigree
#LD=g++-4.9
#CXX=clang++
COV=gcov$(GCC_VERSION)
LD=$(CXX)
LIBS=-rdynamic -lexpat -ldl
VERSION_PATCH=$(shell git rev-list --count `git rev-list --tags --max-count=1`..HEAD 2>/dev/null || echo 0)
VERSION_DEFINES=$(shell git tag | tail -1 | sed 's/\([0-9]*\)[.]\([0-9]*\)\([.][0-9]*\)\?/-DVERSION_MAJOR=\\\"\1\\\" -DVERSION_MINOR=\\\"\2\\\" -DVERSION_PATCH=\\\"$(VERSION_PATCH)\\\"/')
CXXARGS=-Wextra -Wall -pthread -fPIC -DEIGEN_NO_DEPRECATED_WARNING $(VERSION_DEFINES)
#CXXARGS=-Wextra -Wall -pthread -fPIC -DEIGEN_NO_DEPRECATED_WARNING $(VERSION_DEFINES)
INPUT_SRC=
MAIN_SRC= #static_data.cc
PEDIGREE_SRC=main.cc cli.cc
#SRC=malloc.cc $(addprefix input/,$(INPUT_SRC)) $(addprefix computations/,$(COMPUTATIONS_SRC)) $(MAIN_SRC)
SRC=$(addprefix ../input/,$(INPUT_SRC)) $(addprefix ../,$(MAIN_SRC)) $(PEDIGREE_SRC)
OBJ=$(subst .cc,.o,$(SRC))
DEP=$(subst .cc,.d,$(SRC))
COV_OBJ=$(subst .cc,.cov.o,$(SRC))
#DEBUG_OPTS=-ggdb
#DEBUG_OPTS=-ggdb -DNDEBUG
OPT_OPTS=-O3 -DEIGEN_NO_DEBUG -DNDEBUG
#PROF_OPTS=-O2 -DEIGEN_NO_DEBUG -DNDEBUG -g -pg
#DEBUG_OPTS=-ggdb
#OPT_OPTS=-O3 -DNDEBUG
#OPT_OPTS=-O -DEIGEN_NO_DEBUG -DNDEBUG
#OPT_OPTS=-DEIGEN_NO_DEBUG -DNDEBUG
#OPT_OPTS=-O -DNDEBUG
#OPT_OPTS=-DNDEBUG
#OPT_OPTS=-O3 -DNDEBUG
#OPT_OPTS=-O3 -DNDEBUG
#DEBUG_OPTS=-Winvalid-pch -ggdb #-H
#DEBUG_OPTS=-ggdb -DNDEBUG
#EXTRA_OPTS=-g
#EXTRA_OPTS=-fomit-frame-pointer -g
#OPT_OPTS=-O0 -ggdb -DEIGEN_NO_DEBUG -DNDEBUG
#OPT_OPTS=-O2 -DEIGEN_NO_DEBUG -DNDEBUG $(EXTRA_OPTS)
#OPT_OPTS=-O2 -DEIGEN_NO_DEBUG -DNDEBUG -g -fPIE -pie -fsanitize=thread
#OPT_OPTS=-O0 -ggdb -DEIGEN_NO_DEBUG -DNDEBUG -g -fPIE -pie -fsanitize=thread
#COV_OPTS=$(DEBUG_OPTS) --coverage
.header: .FORCE
@echo "*********************************************************"
@echo "** MAKE ***********************************************"
@echo "*********************************************************"
-include .header
INFO_FILE=test-coverage.info
C=$(CXX) $(CXXARGS) $(INC) $(OPT_OPTS) $(DEBUG_OPTS) $(PROF_OPTS) $(EXTRA_OPTS)
L=$(LD) $(CXXARGS) $(INC) $(OPT_OPTS) $(DEBUG_OPTS) $(PROF_OPTS) $(EXTRA_OPTS)
all: $(TARGET)
include .depend
test: $(OBJ)
$C $(OBJ) -o $@ $(LIBS)
test-coverage: test_fast_polynom
lcov --zerocounters --directory .
./$< | tail -3
lcov --gcov-tool $(COV) --no-checksum --directory . --capture --output-file $(INFO_FILE)
lcov --remove $(INFO_FILE) "/usr*" -o $(INFO_FILE)
lcov --remove $(INFO_FILE) "catch.hpp" -o $(INFO_FILE)
genhtml --highlight --legend --output-directory TestCodeCoverage $(INFO_FILE)
#.depend: $(SRC) Makefile
# $C --depend $(SRC) > $@
.depend: $(DEP)
cat $(DEP) > $@
$(DEP):%.d: %.cc
$C -MM $< -MT $@\ $(subst .d,.o,$@) -o $@
clean:
rm -f $(DEP) .depend
rm -f $(OBJ)
rm -f $(TARGET)
.FORCE:
$(TARGET): $(OBJ)
$L $(OBJ) $(LIBS) -o $@
$(OBJ):%.o: %.cc
$C -c $< -o $@
GCC_VERSION=-4.9
CXX=g++$(GCC_VERSION) --std=gnu++0x -Wno-unused-local-typedefs -fPIC
#CXX=clang++ -std=c++11 -stdlib=libc++
INC=-I../../include -I../../include/input -I/home/daleroux/include/eigen3.2 -I/home/daleroux/include -I../../include/bayes
COV=gcov #llvm-cov-3.4 #gcov #$(GCC_VERSION)
LD=$(CXX)
#include "pedigree_settings.h"
#define PREDICATE CONSTRAINT_PREDICATE(ped_settings_t)
#define CALLBACK_ARGS ARGUMENT_CALLBACK_ARGS(ped_settings_t)
ped_settings_t* ensure(ped_settings_t*& t)
{
if (t == NULL) {
t = new ped_settings_t();
}
return t;
}
static
argument_section_list_t<ped_settings_t>
arguments = {
{"Advanced (NOT INTENDED FOR USERS!)", "", true, {
{{"--full-help", "--advanced-help"},
{},
"Display usage.",
false,
{false},
[](CALLBACK_ARGS)
{
print_usage_impl(true, arguments);
exit(0);
SAFE_IGNORE_CALLBACK_ARGS;
}},
{{"-D", "--debug"},
{},
"Display debug messages",
false,
{false},
[](CALLBACK_ARGS)
{
msg_handler_t::debug_enabled() = true;
SAFE_IGNORE_CALLBACK_ARGS;
}},
}},
{"Miscellaneous", "", false, {
{{"-h", "--help"},
{},
"Display usage.",
false,
{false},
[](CALLBACK_ARGS)
{
print_usage_impl(false, arguments);
exit(0);
SAFE_IGNORE_CALLBACK_ARGS;
}},
{{"-wd", "--work-directory"},
{"path"},
"Path to directory for output files",
false,
{&ped_settings_t::work_directory},
[](CALLBACK_ARGS) {
ensure(target)->work_directory = *++ai;
check_file(ensure(target)->work_directory, true, true, true);
SAFE_IGNORE_CALLBACK_ARGS;