Commit 062029db authored by Damien Leroux's avatar Damien Leroux
Browse files

Added option to manually select the thresholds and made chrono output to stderr.

parent c0d8ada8
......@@ -9,7 +9,7 @@ struct chrono {
~chrono_map()
{
for (auto& kv: *this) {
std::cout << kv.first << ": " << kv.second.accum << " seconds." << std::endl;
std::cerr << kv.first << ": " << kv.second.accum << " seconds." << std::endl;
}
}
};
......
......@@ -158,7 +158,7 @@ struct settings_t {
CIM
};
double qtl_threshold;
std::map<std::string, double> qtl_thresholds;
int n_permutations;
double qtl_threshold_quantile;
double cofactor_threshold;
......@@ -194,7 +194,7 @@ struct settings_t {
, cofactor_marker_selection_filename()
, cofactor_marker_selection_distance(10)
, epistasis_qtl_selection_filename()
, qtl_threshold(0)
, qtl_thresholds()
, n_permutations(10000)
, qtl_threshold_quantile(0.05)
, cofactor_threshold(0)
......@@ -295,6 +295,8 @@ struct settings_t {
ss >> parallel;
}
}
double get_threshold(const std::string& trait);
};
std::ostream& operator << (std::ostream&, const settings_t&);
......
TARGET=spell-qtl.experimental
#TARGET=spell-qtl
#TARGET=spell-qtl.experimental
TARGET=spell-qtl
GCC_VERSION=-4.9
CXX=g++$(GCC_VERSION) -fPIC
......
......@@ -120,6 +120,34 @@ read_locus_list(std::string& s, settings_t* target)
}
void
read_threshold_values(const std::string& s, settings_t* target)
{
std::map<std::string, double>& ql = target->qtl_thresholds;
std::string token;
auto i = s.cbegin();
auto j = s.cend();
size_t locus_end, colon;
while (i != j) {
locus_end = s.find_first_of(',');
if (locus_end == std::string::npos) {
locus_end = j - i;
}
colon = s.find_first_of('=');
if (colon == std::string::npos) {
colon = j - i;
}
if (locus_end <= colon) {
MSG_ERROR("Invalid trait threshold around \"" << std::string(i, i + locus_end) << "\": Trait threshold format is trait_name=value", "Specify trait thresholds in the correct format");
}
std::string name(i, i + colon);
ql[name] = to<double>(std::string(i + colon + 1, i + locus_end));
i += locus_end;
i += (i != j); /* skip comma if not at end */
}
}
std::vector<std::string>
cut(const std::string& s, int width)
{
......@@ -402,21 +430,30 @@ argument_section_list_t arguments = {
SAFE_IGNORE_CALLBACK_ARGS;
}},
{{"qtl-threshold-permutations"},
{"value"},
"Set the number of permutations to compute the QTL threshold value in automatic mode",
[](CALLBACK_ARGS)
{
ensure(target)->n_permutations = to<int>(*++ai);
SAFE_IGNORE_CALLBACK_ARGS;
}},
{{"qtl-threshold-quantile"},
{"value"},
"Set the quantile to select the QTL threshold value",
"Set the quantile value in range [0:1] to select the QTL threshold value in automatic mode",
[](CALLBACK_ARGS)
{
std::cout << "TODO" << std::endl;
ensure(target)->qtl_threshold_quantile = to<double>(*++ai);
SAFE_IGNORE_CALLBACK_ARGS;
}},
{{"qtl-threshold-value"},
{"value"},
"Set the QTL threshold value manually",
{"value_list"},
"Set the QTL threshold value manually for some traits",
[](CALLBACK_ARGS)
{
std::cout << "TODO" << std::endl;
read_threshold_values(*++ai, ensure(target));
SAFE_IGNORE_CALLBACK_ARGS;
}},
......
......@@ -60,6 +60,10 @@ DTD_START(analysis_dtd, analysis, settings_t)
ELEMENT(threshold, settings_t);
ELEMENT(permutations, int);
ELEMENT(quantile, double);
typedef std::map<std::string, double> threshold_map_type;
typedef std::pair<std::string, double> threshold_map_value_type;
ELEMENT(values, settings_t);
ELEMENT(trait, threshold_map_type);
ELEMENT_WITH_NAME(wd, "work-directory", std::string);
ELEMENT(model, settings_t);
ELEMENT(connected, bool);
......@@ -121,7 +125,15 @@ DTD_START(analysis_dtd, analysis, settings_t)
= A("value");
threshold
= (E(quantile, &settings_t::qtl_threshold_quantile),
E(permutations, &settings_t::n_permutations));
E(permutations, &settings_t::n_permutations)
/*,*/
/*E(values)*/
);
/*values*/
/*= E(trait, &settings_t::qtl_thresholds);*/
/*trait*/
/*= (A("name", &std::pair<std::string, double>::first),*/
/*A("value", &std::pair<std::string, double>::second));*/
model
= (E(connected, &settings_t::connected),
E(epistasis, &settings_t::epistasis),
......
/*#include "all.h"*/
#include "malloc.h"
#include "computations.h"
#include "chrono.h"
double settings_t::get_threshold(const std::string& trait)
{
double& thres = qtl_thresholds[trait];
if (thres == 0) {
chrono::start("thresholds");
std::stringstream title;
title << "Computing threshold for trait '" << trait << "' with " << n_permutations << " permutations and quantile " << qtl_threshold_quantile;
set_title(title.str());
thres = *make_value<Disk|Sync>(qtl_threshold_all_chromosomes,
value<std::string>(trait),
value<double>(qtl_threshold_quantile),
value<int>(n_permutations));
chrono::stop("thresholds");
}
return thres;
}
......@@ -41,6 +61,7 @@ int main(int argc, const char** argv)
exit(0);
}
chrono::start("initialization");
active_settings = settings_t::from_args(argc, argv);
if (!active_settings) {
......@@ -53,6 +74,8 @@ int main(int argc, const char** argv)
active_settings->finalize();
chrono::stop("initialization");
/* some screen management */
if (msg_handler_t::color() && active_settings->parallel > 1) {
MSG_DEBUG("\x1b[2J\x1b[0;0H" << std::endl << std::endl << std::endl << std::endl);
......@@ -90,7 +113,7 @@ int main(int argc, const char** argv)
/*computations::f_test_along_chromosome ftac(M0, M0, pop, chr, gen, computations::selected_qtls_on_chromosome(qtl_chr));*/
collection<population_value> all_pop = all_populations();
active_settings->set_title("Computing thresholds");
/*active_settings->set_title("Computing thresholds");*/
/*collection<double> thresholds*/
/*= make_collection<Disk|Sync>(qtl_threshold,*/
/*all_traits(), qtl_chr,*/
......@@ -98,23 +121,26 @@ int main(int argc, const char** argv)
/*value<int>(active_settings->n_permutations));*/
auto atr = all_traits();
*atr[0];
/*MSG_DEBUG("Traits: " << atr);*/
/*MSG_INFO("Traits: " << atr);*/
/**thresholds[0];*/
int i = 0;
collection<double> thresholds
= make_collection<Disk|Sync>(qtl_threshold_all_chromosomes,
all_traits(),
value<double>(active_settings->qtl_threshold_quantile),
value<int>(active_settings->n_permutations));
/*collection<double> thresholds*/
/*= make_collection<Disk|Sync>(qtl_threshold_all_chromosomes,*/
/*all_traits(),*/
/*value<double>(active_settings->qtl_threshold_quantile),*/
/*value<int>(active_settings->n_permutations));*/
if (*thresholds[0] == 0) {
double threshold = active_settings->get_threshold(*atr[0]);
if (threshold == 0) {
std::cout << "FOIRURE" << std::endl;
return -1;
}
active_settings->set_title("Species Perscrutandis Enixe Locis Locabuntur");
chrono::start("spell along chromosomes");
for (auto& qc: active_settings->working_set) {
qtl_chr = &qc;
chr = (*qtl_chr)->chr;
......@@ -145,7 +171,7 @@ int main(int argc, const char** argv)
/*} else {*/
/*std::cout << " NA NA";*/
/*}*/
std::cout << ' ' << (*thresholds[0]) << ' ' << ftac(max);
std::cout << ' ' << threshold << ' ' << ftac(max);
std::cout << ' ' << r2(*M0, *M1);
auto& loci = active_settings->estimation_loci[*chr];
......@@ -202,6 +228,7 @@ int main(int argc, const char** argv)
++i;
}
chrono::stop("spell along chromosomes");
#if 0
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment