Skip to content
GitLab
Commit 63f9dadd authored by Damien Leroux's avatar Damien Leroux
Browse files

Generating reports in XLSX format. Much redesign to do.

parent 6304b230
Hide whitespace changes
Inline Side-by-side
CMakeLists.txt
View file @ 63f9dadd
......@@ -67,6 +67,8 @@ find_path(EXPAT_INCLUDE_DIR expat.h HINTS /usr/include /usr/local/include /usr/i
# find_path(X2C_INCLUDE_DIR x2c/x2c.h HINTS /usr/include/ /usr/local/include/ /home/daleroux/include/)
find_path(EIGEN_INCLUDE_DIR Eigen/Eigen HINTS /usr/include /usr/local/include /usr/include/eigen3/ /usr/local/include/eigen3/ /home/daleroux/include/eigen3/)
find_library(XLNT_LIBRARY_NAMES xlnt)
find_path(XLNT_INCLUDE_DIR xlnt/xlnt.hpp HINTS /usr/include /usr/local/include)
include_directories(AFTER 3rd-party/ThreadPool)
include_directories(AFTER ${CMAKE_SOURCE_DIR}/include/ ${CMAKE_SOURCE_DIR}/include/input/ ${CMAKE_SOURCE_DIR}/include/bayes/ ${EIGEN_INCLUDE_DIR})
include_directories(SYSTEM /usr/include ${Boost_INCLUDE_DIRS} ${EXPAT_INCLUDE_DIR})
......@@ -128,7 +130,7 @@ if(${BUILD_FOR_DEPLOYMENT})
SET_SOURCE_FILES_PROPERTIES(${SPELL_PEDIGREE_SRC} ${SPELL_MARKER_SRC} ${SPELL_QTL_SRC}
PROPERTIES
OBJECT_DEPENDS glibc.h
COMPILE_FLAGS "-include ${CMAKE_BINARY_DIR}/glibc.h -flto -U_FORTIFY_SOURCE -D_FORTIFY_SOURCE=0")
COMPILE_FLAGS "-include ${CMAKE_BINARY_DIR}/glibc.h -U_FORTIFY_SOURCE -D_FORTIFY_SOURCE=0")
#add_dependencies(stdc++ libstdc++.a)
#add_dependencies(gmp libgmp.a)
......@@ -138,7 +140,7 @@ if(${BUILD_FOR_DEPLOYMENT})
MESSAGE(STATUS " PROUT ${libstdcpp}")
set(CMAKE_EXE_LINKER_FLAGS "-include ${CMAKE_BINARY_DIR}/glibc.h -flto -rdynamic -static-libgcc")
set(CMAKE_EXE_LINKER_FLAGS "-include ${CMAKE_BINARY_DIR}/glibc.h -rdynamic -static-libgcc")
target_link_libraries(spell-marker ${libstdcpp} ${CMAKE_BINARY_DIR}/libgmp.a)
target_link_libraries(spell-qtl ${libstdcpp})
......@@ -149,12 +151,13 @@ else()
add_executable(spell-qtl ${SPELL_QTL_SRC})
SET_SOURCE_FILES_PROPERTIES(${SPELL_PEDIGREE_SRC} ${SPELL_MARKER_SRC} ${SPELL_QTL_SRC}
PROPERTIES
COMPILE_FLAGS "-flto -U_FORTIFY_SOURCE -D_FORTIFY_SOURCE=0")
set(CMAKE_EXE_LINKER_FLAGS "-flto -rdynamic")
COMPILE_FLAGS "-U_FORTIFY_SOURCE -D_FORTIFY_SOURCE=0")
target_compile_definitions(spell-qtl PUBLIC SPELL_USE_XLNT)
set(CMAKE_EXE_LINKER_FLAGS "-rdynamic")
endif()
target_link_libraries(spell-marker dl gmp)
target_link_libraries(spell-qtl dl rt)
target_link_libraries(spell-qtl xlnt dl rt)
SET(EXECUTABLE_OUTPUT_PATH ${PROJECT_BINARY_DIR}/bin)
......
include/computations/excel_report.h 0 → 100644
View file @ 63f9dadd
#ifndef _SPELL_EXCEL_REPORT_H_
#define _SPELL_EXCEL_REPORT_H_
#include "labelled_matrix.h"
#include "excel.h"
#include "model/print.h"
#define SGNF_0 "sgnf0"
#define SGNF_1 "sgnf1"
#define SGNF_2 "sgnf2"
#define SGNF_3 "sgnf3"
#define SGNF_4 "sgnf4"
namespace excel_report {
inline
void
create_signif_styles(excel::stream& s)
{
xlnt::font w, b;
w.color(xlnt::color::white());
b.color(xlnt::color::black());
xlnt::color c0(xlnt::indexed_color(11));
xlnt::color c1(xlnt::indexed_color(13));
xlnt::color c2(xlnt::indexed_color(55));
xlnt::color c3(xlnt::indexed_color(53));
xlnt::color c4(xlnt::indexed_color(10));
// Significance codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
s << excel::style(SGNF_0, [&](xlnt::style& sty) { sty.fill(xlnt::fill().solid(c0)); })
<< excel::style(SGNF_1, [&](xlnt::style& sty) { sty.fill(xlnt::fill().solid(c1)); })
<< excel::style(SGNF_2, [&](xlnt::style& sty) { sty.fill(xlnt::fill().solid(c2)); })
<< excel::style(SGNF_3, [&](xlnt::style& sty) { sty.fill(xlnt::fill().solid(c3)); })
<< excel::style(SGNF_4, [&](xlnt::style& sty) { sty.fill(xlnt::fill().solid(c4)); })
;
}
inline
std::string
get_signif_style(double v)
{
if (v < .001) { return SGNF_0; }
if (v < .01) { return SGNF_1; }
if (v < .05) { return SGNF_2; }
if (v < .1) { return SGNF_3; }
return SGNF_4;
}
}
inline excel::stream& operator << (excel::stream& s, const std::vector<char>& vc) { std::string str(vc.begin(), vc.end()); return s << str; }
inline excel::stream& operator << (excel::stream& s, const model_block_key& mbk) { return s << MESSAGE(mbk); }
namespace excel {
template <typename RS, typename RF, typename CS, typename CF, typename M>
manipulator
matrix_with_sections(std::string title, const model_print::matrix_with_sections<RS, RF, CS, CF, M>& mws, std::function<void(int, int, xlnt::cell)> formatter = [] (int, int, xlnt::cell) {})
{
static auto set_mat_title = excel::merge(2, 2) << excel::align([] (xlnt::alignment& al) {
al.horizontal(xlnt::horizontal_alignment::center).vertical(xlnt::vertical_alignment::center);
})
<< excel::border_right(xlnt::color::black(), xlnt::border_style::thin)
<< excel::border_bottom(xlnt::color::black(), xlnt::border_style::thin)
;
static excel::manipulator vertical_label = excel::align([](xlnt::alignment& al) {
al.rotation(90).horizontal(xlnt::horizontal_alignment::center).vertical(xlnt::vertical_alignment::center);
});
static excel::manipulator horizontal_label = excel::align([](xlnt::alignment& al) {
al.horizontal(xlnt::horizontal_alignment::center).vertical(xlnt::vertical_alignment::center);
});
xlnt::color border_color = xlnt::color::black();
xlnt::border_style border_style = xlnt::border_style::thin;
static excel::manipulator BL = excel::border_left(border_color, border_style);
static excel::manipulator BT = excel::border_top(border_color, border_style);
static excel::manipulator BR = excel::border_right(border_color, border_style);
static excel::manipulator BB = excel::border_bottom(border_color, border_style);
return [&] (stream& s) {
xlnt::cell_reference origin = s.origin;
excel::push(s);
#if 0
s << excel::move(20, 2) << BT << BR << BL << BB;
s << excel::move(23, 2) << BT << BR << BB << BL;
s << excel::move(26, 2) << BT << BB << BR << BL;
s << excel::move(29, 2) << BB << BT << BR << BL;
s << excel::move(32, 2) << BT << BL << BB << BR;
s << excel::move(20, 5) << BT << BR << BL << BB;
s << excel::move(22, 5) << BT << BR << BB << BL;
s << excel::move(24, 5) << BT << BB << BR << BL;
s << excel::move(26, 5) << BB << BT << BR << BL;
s << excel::move(28, 5) << BT << BL << BB << BR;
s << excel::move(28, 8) << BT << BL << BB << BR;
s << excel::move(26, 8) << BB << BT << BR << BL;
s << excel::move(24, 8) << BT << BB << BR << BL;
s << excel::move(22, 8) << BT << BR << BB << BL;
s << excel::move(20, 8) << BT << BR << BL << BB;
#endif
// s << "ORIGIN";
/* column section labels */
s << excel::move(origin.make_offset(2, 0));
for (const auto& sec: mws.m_column_sections) {
s << sec.label() << horizontal_label << excel::merge(sec.n_fields(), 1) << BL << BT << excel::next_col;
}
s << BL;
/* column field labels */
s << excel::move(origin.make_offset(2, 1));
for (const auto& sec: mws.m_column_sections) {
s << BL;
for (const auto& l: sec.field_labels()) {
s << l << horizontal_label << BB << excel::next_col;
}
}
s << excel::prev_col << BR;
/* row section labels */
s << excel::move(origin.make_offset(0, 1));
for (const auto& sec: mws.m_row_sections) {
s << excel::next_row << excel::merge(1, sec.n_fields()) << sec.label() << vertical_label << BL << BT;
}
s << BB;
/* row field labels */
s << excel::move(origin.make_offset(1, 2));
int rborder = 1;
for (const auto& sec: mws.m_row_sections) {
s << BT;
for (const auto& l: sec.field_labels()) {
s << l << vertical_label << BR << excel::next_row;
}
rborder += sec.n_fields();
s << excel::push << excel::move(origin.make_offset(1, rborder)) << BB << BR << excel::pop;
}
/* data */
auto& mat = mws.m_matrix;
int r0 = 0, r1 = 0;
for (const auto& rs: mws.m_row_sections) {
r1 += rs.n_fields();
int c0 = 0, c1 = 0;
for (const auto& cs: mws.m_column_sections) {
c1 += cs.n_fields();
s << excel::move(origin.make_offset(2 + c0, 2 + r0));
int r, c;
for (r = r0; r < r1 - 1; ++r) {
for (c = c0; c < c1 - 1; ++c) {
MSG_DEBUG('(' << r << ", " << c << ')');
s << mat(r, c);
formatter(r, c, s.sheet[s.cursor]);
s << excel::next_col;
}
s << mat(r, c);
formatter(r, c, s.sheet[s.cursor]);
s << BR << excel::next_row;
}
for (c = c0; c < c1 - 1; ++c) {
// MSG_DEBUG('(' << r << ", " << c << ')');
s << mat(r, c) << BB;
formatter(r, c, s.sheet[s.cursor]);
s << excel::next_col;
// s << mat(r, c) << BB << excel::next_col;
}
s << mat(r, c) << BB << BR;
formatter(r, c, s.sheet[s.cursor]);
s << excel::next_row;
// s << mat(r, c) << BB << BR << excel::next_row;
c0 += cs.n_fields();
}
s << excel::next_row;
r0 += rs.n_fields();
}
// for (int r = 0; r < mat.rows(); ++r) {
// for (int c = 0; c < mat.cols(); ++c) {
// s << mat(r, c) << excel::next_col;
// }
// s << excel::next_row;
// }
excel::pop(s);
s.cursor = s.origin;
s << title << set_mat_title;
s.cursor = origin.make_offset(1 + mat.cols(), 1 + mat.rows());
return s;
};
}
inline
stream&
autosize(stream& s) {
auto dimensions = s.sheet.calculate_dimension();
for (xlnt::column_t::index_t c = dimensions.top_left().column_index(); c <= dimensions.top_right().column_index(); ++c) {
MSG_DEBUG("Computed width " << s.sheet.column_width(c));
}
for (xlnt::row_t r = dimensions.top_left().row(); r <= dimensions.bottom_left().row(); ++r) {
MSG_DEBUG("Computed height " << s.sheet.row_height(r));
}
return s;
}
} /* namespace excel */
#endif
include/computations/frontends4.h
View file @ 63f9dadd
......@@ -21,12 +21,14 @@
#include "cache2.h"
#include "basic_data.h"
#include "model.h"
#include "../../attic/frontends2.h"
// #include "../../attic/frontends2.h"
/*#include "model/tests.h"*/
#include <regex>
#include <boost/math/distributions/normal.hpp> // for normal_distribution
#include "excel_report.h"
#define SIGNAL_DISPLAY_ONELINER
#define CONSTRAINT_RESIDUAL_EPSILON 1.e-5
......@@ -468,8 +470,10 @@ struct analysis_report {
}
void report_final_model(model_manager& mm);
void report_qtls(std::vector<QTL> &qtls, ComputationType lod_test_type);
void report_final_model_excel(model_manager& mm);
void report_qtls_excel(excel::stream& XLS, std::vector<QTL> &qtls, ComputationType lod_test_type);
};
......@@ -2076,7 +2080,7 @@ struct model_manager {
model_manager sub(*this);
sub.vMcurrent = model{*vMcurrent};
sub.deselect(chr, loc);
sub.custom_test(sub.lod_test_type(), ComputationResults::Test, sub.vMcurrent, {}, chr).find(chr)->second;
sub.custom_test(sub.lod_test_type(), ComputationResults::Test, sub.vMcurrent, {}, chr).find(chr);
MSG_DEBUG("sub.last_test_positions " << sub.last_test_positions);
ret.emplace_back(chr->name, loc, sub.last_test_positions[chr], *sub.last_computation[chr]);
}
......@@ -2377,6 +2381,12 @@ void analysis_report::report_final_model(model_manager& mm)
/*header(report_file, "RSS ") << mm.vMcurrent->rss() << std::endl;*/
//*/
MSG_DEBUG("pouet");
xlnt::workbook wb;
MSG_DEBUG("coin!");
excel::stream XLS(&wb);
excel_report::create_signif_styles(XLS);
/* TODO ajouter R2 global */
MSG_DEBUG(__FILE__ << ':' << __LINE__ << " keys " << mm.vMcurrent->keys());
......@@ -2471,6 +2481,26 @@ void analysis_report::report_final_model(model_manager& mm)
MatrixXd vcov = mm.vMcurrent->XtX_pseudo_inverse() * pseudo_variances(i); // mm.vMcurrent->vcov(i);
model::xtx_printable vc(*mm.vMcurrent, vcov);
header(report_file, "Covariance matrix") << vc << std::endl;
// XLS << excel::move(2, 2) << excel::matrix_with_sections(vc.mws,
// [&] (int r, int c, xlnt::cell cell) {
// static int i = -1;
// static std::vector<std::string> sty = {SGNF_0, SGNF_1, SGNF_2, SGNF_3, SGNF_4};
// cell.style(sty[(++i) % sty.size()]);
// });
// XLS << excel::move(XLS.origin.make_offset(0, 10));
//
// for (int r = 0; r < 16; ++r) {
// for (int c = 0; c < 16; ++c) {
// XLS << excel::move(31 + r, 1 + c);
// int i = (r << 4) + c;
// XLS << excel::style(MESSAGE(r << '_' << c), [&] (xlnt::style& s) { s.fill(xlnt::fill().solid(xlnt::indexed_color(i))); });
// XLS << excel::style(MESSAGE(r << '_' << c)) << i;
// }
// }
//
// wb.save("prout.xlsx");
auto conmat = X.bottomRows(X.rows() - mm.vMcurrent->n_obs());
for (int r = 0; r < conmat.rows(); ++r) {
......@@ -2570,6 +2600,226 @@ void analysis_report::report_final_model(model_manager& mm)
section_header(report_file, "XtX^-1") << std::endl << xtx << std::endl;
}
inline
void analysis_report::report_final_model_excel(model_manager& mm)
{
xlnt::workbook wb;
excel::stream XLS(&wb);
excel_report::create_signif_styles(XLS);
XLS << excel::rename_worksheet("Metadata");
static std::string emptystr;
MSG_DEBUG(__FILE__ << ':' << __LINE__ << " keys " << mm.vMcurrent->keys());
mm.vMcurrent->compute();
mm.vMcurrent->output_X_to_file(full_path);
mm.vMcurrent->output_XtX_inv_to_file(full_path);
model::xtx_printable xtx(*mm.vMcurrent, true);
const auto& TM = active_settings->trait_metadata[mm.trait_name];
/*const auto& XtX_inv = mm.vMcurrent->XtX_pseudo_inverse();*/
MatrixXd X = mm.vMcurrent->X();
/*const auto& XtX = mm.vMcurrent->XtX();*/
/*size_t dof = mm.vMcurrent->rank();*/
/*size_t n = mm.vMcurrent->Y().rows();*/
auto qtls = mm.QTLs();
XLS << "Spell-QTL" << excel::next_row << excel::next_row;
XLS << "Report for " << (TM.dim_names.size() > 1 ? "pleiotropic" : "single") << " trait " << mm.trait_name << excel::next_row;
report_qtls_excel(XLS, qtls, mm.lod_test_type());
//*
/*header(report_file, "RSS ") << mm.vMcurrent->rss() << std::endl;*/
//*/
/* TODO ajouter R2 global */
MSG_DEBUG(__FILE__ << ':' << __LINE__ << " keys " << mm.vMcurrent->keys());
mm.vMcurrent->compute();
/* Checks */
{
auto X = mm.vMcurrent->X();
if (mm.vMcurrent->rank() < X.cols()) {
MSG_ERROR("Insufficient rank for model matrix.", "Report to the developers");
}
int n_constraints = X.rows();
for (const auto& p: mm.all_pops) {
n_constraints -= (*p)->size();
}
auto res = mm.vMcurrent->residuals().bottomRows(n_constraints);
if (res.lpNorm<1>() > CONSTRAINT_RESIDUAL_EPSILON) {
MSG_DEBUG(mm.vMcurrent->residuals().transpose());
MSG_DEBUG(res.transpose());
MSG_ERROR("Inconsistent contraints for model matrix. " << res.lpNorm<1>() << " > " << CONSTRAINT_RESIDUAL_EPSILON, "Report to the developers");
}
}
MSG_DEBUG(__FILE__ << ':' << __LINE__ << " keys " << mm.vMcurrent->keys());
if (mm.vMcurrent->n_obs() <= mm.vMcurrent->dof() + 1) {
XLS << std::string("* There are not enough observations to compute a variance/covariance matrix. Results not displayed.") << excel::next_row;
} else {
MatrixXd Pt = TM.P.transpose();
MatrixXd all_coefficients = mm.vMcurrent->coefficients() * Pt;
/* FIXME use actual full Y instead of un-projecting the traits? */
MatrixXd all_residuals = mm.vMcurrent->residuals() * Pt;
VectorXd all_rss = all_residuals.array().square().colwise().sum();
MSG_QUEUE_FLUSH();
VectorXd pseudo_variances = (mm.vMcurrent->rss().transpose() * Pt.array().square().matrix()).transpose() / (mm.vMcurrent->n_obs() - mm.vMcurrent->dof() - 1);
for (int i = 0; i < (int) TM.dim_names.size(); ++i) {
std::string title;
if (TM.dim_names.size() > 1) {
// section_header(report_file, MESSAGE(mm.trait_name << " # " << TM.dim_names[i]));
title = MESSAGE(mm.trait_name << " # " << TM.dim_names[i]);
} else {
// section_header(report_file, mm.trait_name);
title = mm.trait_name;
}
XLS << excel::new_worksheet(title);
XLS << title << excel::next_row << excel::next_row;
{
MatrixXd R2(1, mm.vMcurrent->m_blocks.size());
auto bi = mm.vMcurrent->m_blocks.begin(), bj = mm.vMcurrent->m_blocks.end();
int c = 0;
for (; bi != bj; ++bi, ++c) {
auto M0 = *mm.vMcurrent;
M0.remove_block(bi->first);
M0.compute();
MatrixXd M0residuals = M0.residuals() * Pt;
VectorXd M0rss = M0residuals.array().square().colwise().sum();
double r1 = all_rss(i);
double r0 = M0rss(i);
R2(0, c) = fabs(r0) < 1.e-6 ? 0 : (r0 - r1) / r0; /* FIXME use some real epsilon */
}
model_print::matrix_with_sections<std::string, void, model_block_key, std::vector<char>> mws(R2);
static std::vector<std::vector<char>> empty_label = {{' '}};
bi = mm.vMcurrent->m_blocks.begin();
for (/*++bi*/; bi != bj; ++bi) { mws.add_column_section(bi->first, empty_label); }
mws.add_row_section(emptystr, 1);
// header(report_file, "R2") << mws << std::endl;
auto orig = XLS.cursor;
XLS << excel::set_origin << excel::matrix_with_sections("R2", mws) << excel::next_row << excel::next_row << excel::next_row;
}
MatrixXd coef = all_coefficients.col(i).transpose();
/*MSG_DEBUG(MATRIX_SIZE(coef));*/
{
model_print::matrix_with_sections<std::string, void, std::string, std::vector<char>> mws(coef);
mm.vMcurrent->set_column_sections(mws);
mws.add_row_section(emptystr, 1);
// header(report_file, "Coefficients") << mws << std::endl;
XLS << excel::set_origin << excel::matrix_with_sections("Coefficients", mws) << excel::next_row << excel::next_row << excel::next_row;
}
MatrixXd vcov = mm.vMcurrent->XtX_pseudo_inverse() * pseudo_variances(i); // mm.vMcurrent->vcov(i);
model::xtx_printable vc(*mm.vMcurrent, vcov);
// header(report_file, "Covariance matrix") << vc << std::endl;
XLS << excel::set_origin << excel::matrix_with_sections("Covariance matrix", vc.mws) << excel::next_row << excel::next_row << excel::next_row;
// XLS << excel::move(2, 2) << excel::matrix_with_sections(vc.mws,
// [&] (int r, int c, xlnt::cell cell) {
// static int i = -1;
// static std::vector<std::string> sty = {SGNF_0, SGNF_1, SGNF_2, SGNF_3, SGNF_4};
// cell.style(sty[(++i) % sty.size()]);
// });
auto conmat = X.bottomRows(X.rows() - mm.vMcurrent->n_obs());
bool legend_once = false;
excel::manipulator legend = [&] (excel::stream& s) {
if (legend_once) {
return;
}
legend_once = true;
auto R = excel::align([] (xlnt::alignment& al) { al.horizontal(xlnt::horizontal_alignment::right); });
s << excel::push << excel::set_origin;
auto tl = s.origin;
s << "Significance codes" << excel::merge(3, 1) << excel::align([] (xlnt::alignment& al) { al.horizontal(xlnt::horizontal_alignment::center); }) << excel::next_row
<< "0 <=" << R << excel::next_col << "***" << excel::style(SGNF_0) << excel::next_col << "< .001" << excel::next_row
<< ".001 <=" << R << excel::next_col << "**" << excel::style(SGNF_1) << excel::next_col << "< .01" << excel::next_row
<< ".01 <=" << R << excel::next_col << "*" << excel::style(SGNF_2) << excel::next_col << "< .05" << excel::next_row
<< ".05 <=" << R << excel::next_col << "." << excel::style(SGNF_3) << excel::next_col << "< .1" << excel::next_row
<< ".1 <=" << R << excel::next_col << excel::style(SGNF_4) << excel::next_col << "< 1.0";
auto br = s.cursor;
s << excel::pop << excel::border_box(tl, br, {{xlnt::color::black(), xlnt::border_style::thin}});
};
for (int r = 0; r < conmat.rows(); ++r) {
auto subset = filter_columns(conmat.row(r), mm.vMcurrent->m_blocks);
const auto& columns = subset.first;
const auto& blocks = subset.second;
/*Eigen::Matrix<std::string, Eigen::Dynamic, Eigen::Dynamic> pvmat = Eigen::Matrix<std::string, Eigen::Dynamic, Eigen::Dynamic>::Constant(columns.size(), columns.size(), empty_cell());*/
Eigen::Matrix<std::string, Eigen::Dynamic, Eigen::Dynamic> pvmat(columns.size(), columns.size());
Eigen::Matrix<std::string, Eigen::Dynamic, Eigen::Dynamic> ctrmat(columns.size(), columns.size());
MatrixXd pvalue_mat = MatrixXd::Zero(pvmat.rows(), pvmat.cols());
/*MatrixXd ctrmat = MatrixXd::Constant(columns.size(), columns.size(), std::numeric_limits<double>::quiet_NaN());*/
for (int i1 = 0; i1 < pvmat.cols(); ++i1) {
for (int i2 = i1 + 1; i2 < pvmat.cols(); ++i2) {
int a1 = columns[i1];
int a2 = columns[i2];
double sigma = vcov(a1, a1) + vcov(a2, a2) - 2 * vcov(a1, a2); /* TODO vérifier auprès de NV² pour rassurer Sylvain DONE en fait c'est bon. */
double delta = coef(a1) - coef(a2);
normal s(0, sigma > 0 ? sqrt(sigma) : 1.e-9);
double pvalue = 2 * cdf(complement(s, fabs(delta)));
pvalue_mat(i1, i2) = pvalue_mat(i2, i1) = pvalue;
ctrmat(i1, i2) = MESSAGE(std::setprecision(3) << -delta << ' ' << std::setw(3) << std::left << significance_string(pvalue));
ctrmat(i2, i1) = MESSAGE(std::setprecision(3) << delta << ' ' << std::setw(3) << std::left << significance_string(pvalue));
pvmat(i1, i2) = pvmat(i2, i1) = format_pvalue(pvalue);
}
}
model_print::matrix_with_sections<void, std::vector<char>, model_block_key, std::vector<char>, Eigen::Matrix<std::string, Eigen::Dynamic, Eigen::Dynamic>> ctrsig(pvmat);
model_print::matrix_with_sections<void, std::vector<char>, model_block_key, std::vector<char>, Eigen::Matrix<std::string, Eigen::Dynamic, Eigen::Dynamic>> ctr(ctrmat);
for (const auto& kl: blocks) {
std::vector<std::vector<char>> names;
if (kl.first->type == mbk_CI) {
for (const auto& kp: mm.vMcurrent->m_all_pops) {
const auto& qgn = (*kp)->qtl_generation_name;
names.emplace_back(qgn.begin(), qgn.end());
}
} else {
names = mm.vMcurrent->labels_to_ancestor_names(kl.second);
}
ctr.add_row_section(names);
ctrsig.add_row_section(names);
for (auto& n: names) { n.push_back(' '); n.push_back(' '); n.push_back(' '); n.push_back(' '); }
ctr.add_column_section(kl.first, names);
ctrsig.add_column_section(kl.first, names);
}
// header(report_file, "Contrasts") << ctr << std::endl << significance_legend() << std::endl;
// header(report_file, "Contrast significance") << ctrsig << std::endl << significance_legend() << std::endl;
auto cell_background = [&] (int r, int c, xlnt::cell cell) {
double v = pvalue_mat(r, c);
if (r == c) {
return;
}
if (v < .001) {
cell.style(SGNF_0);
} else if (v < .01) {
cell.style(SGNF_1);
} else if (v < .05) {
cell.style(SGNF_2);
} else if (v < .1) {
cell.style(SGNF_3);
} else {
cell.style(SGNF_4);
}
};
{
// int r = XLS.cursor.row();
XLS << excel::set_origin << excel::matrix_with_sections("Contrasts", ctr, cell_background);
XLS << excel::push << excel::moveby(-ctr.m_matrix.rows() - 1, 2) << legend << excel::next_row << excel::pop << excel::next_row << excel::next_row << excel::next_row;
XLS << excel::set_origin << excel::matrix_with_sections("Contrast significance", ctrsig, cell_background) << excel::next_row << excel::next_row << excel::next_row;
}
}
}
}
// section_header(report_file, "Final model") << std::endl << (*mm.vMcurrent) << std::endl;
// section_header(report_file, "XtX^-1") << std::endl << xtx << std::endl;
XLS << excel::new_worksheet("Final model") << excel::matrix_with_sections("Final model", mm.vMcurrent->printable_X());
XLS << excel::new_worksheet("XtX^-1") << excel::matrix_with_sections("XtX^-1", xtx.mws);
wb.save(MESSAGE(full_path << '/' << trait_name << "_report.xlsx"));
}
inline
void analysis_report::report_qtls(std::vector<QTL> &qtls, ComputationType lod_test_type)
{
......@@ -2586,6 +2836,27 @@ void analysis_report::report_qtls(std::vector<QTL> &qtls, ComputationType lod_te
}
}
inline
void analysis_report::report_qtls_excel(excel::stream& XLS, std::vector<QTL> &qtls, ComputationType lod_test_type)
{
auto M = excel::align([] (xlnt::alignment& al) { al.horizontal(xlnt::horizontal_alignment::center); });
XLS << "Detected QTLs" << excel::next_row << excel::next_row;
XLS << "Chromosome" << M << excel::next_col << "Position" << M << excel::next_col << "Confidence interval" << excel::merge(2, 1) << M << excel::next_row;
for (auto& qtl: qtls) {
const auto& ci = roi[qtl.chromosome][qtl.locus] = qtl.confidence_interval();
report_lod(qtl);
std::string name = MESSAGE(trait_name << " @ " << qtl.locus << " [" << ci.first << ':' << ci.second << ']');
if (poi[qtl.chromosome][qtl.locus].size()) {
poi[qtl.chromosome][qtl.locus] = MESSAGE(poi[qtl.chromosome][qtl.locus] << ',' << name);
} else {
poi[qtl.chromosome][qtl.locus] = name;
}
XLS << qtl.chromosome << excel::next_col << qtl.locus << excel::next_col
<< roi[qtl.chromosome][qtl.locus].first << excel::next_col
<< roi[qtl.chromosome][qtl.locus].second << excel::next_row;
}
}