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

Work in progress. Going nicely!

parent b8642c21
......@@ -87,7 +87,10 @@ struct ostream_manager {
class HeaderInserter : public std::streambuf {
std::streambuf* dest;
bool start_of_line;
int indent;
/*int indent;*/
std::vector<std::string> prefix;
bool prefix_insert_mode;
std::stringstream prefix_ss;
protected:
int overflow(int ch) override;
......@@ -95,7 +98,10 @@ struct ostream_manager {
HeaderInserter(std::streambuf* dest)
: dest(dest)
, start_of_line(true)
, indent(0)
/*, indent(0)*/
, prefix()
, prefix_insert_mode(false)
, prefix_ss()
{}
};
......@@ -276,14 +282,34 @@ int ostream_manager::HeaderInserter::overflow(int ch)
{
int retval = 0;
if (ch == 1) {
indent += 3;
/*indent += 3;*/
prefix_insert_mode = !prefix_insert_mode;
if (prefix_insert_mode) {
prefix_ss.clear();
prefix_ss.str(std::string());
} else {
prefix.push_back(prefix_ss.str());
}
} else if (ch == 2) {
indent -= 3 * (indent > 0);
if (prefix_insert_mode) {
overflow('\x1');
}
if (prefix.size()) {
prefix.pop_back();
}
/*indent -= 3 * (indent > 0);*/
} else if (prefix_insert_mode) {
prefix_ss << (char) ch;
} else if (ch != traits_type::eof()) {
if (start_of_line) {
for (int i = 0; i < indent; ++i) {
dest->sputc(' ');
for (const std::string& s: prefix) {
for (const char c: s) {
dest->sputc(c);
}
}
/*for (int i = 0; i < indent; ++i) {*/
/*dest->sputc(' ');*/
/*}*/
}
retval = dest->sputc( ch );
start_of_line = ch == '\n';
......@@ -389,10 +415,17 @@ void message_queue::run()
} while(0)
#endif
#define MSG_DEBUG_INDENT CREATE_MESSAGE(msg_channel::Log, "\x1")
#define MSG_DEBUG_INDENT CREATE_MESSAGE(msg_channel::Log, "\x1 \x1")
#define MSG_DEBUG_INDENT_EXPR(_str_) CREATE_MESSAGE(msg_channel::Log, MESSAGE('\x1' << _str_ << '\x1'))
#define MSG_DEBUG_DEDENT CREATE_MESSAGE(msg_channel::Log, "\x2")
struct scoped_indent { scoped_indent() { MSG_DEBUG_INDENT; } ~scoped_indent() { MSG_DEBUG_DEDENT; } };
struct scoped_indent {
scoped_indent() { MSG_DEBUG_INDENT; }
scoped_indent(const std::string& str) { MSG_DEBUG_INDENT_EXPR(str); }
~scoped_indent() { MSG_DEBUG_DEDENT; }
};
inline void msg_handler_t::state_t::check(bool fatal)
{
......
......@@ -18,6 +18,7 @@
/*typedef std::set<int> subset;*/
typedef std::vector<int> subset;
#if 1
#define _EPSILON (1.e-10)
namespace std {
......@@ -29,6 +30,7 @@ namespace std {
}
};
}
#endif
#if 0
template <typename K, typename V, typename H=std::hash<K>, typename A=std::allocator<std::pair<const K, V>>>
......@@ -199,7 +201,7 @@ struct lumper {
stack.push_back(&x);
}
while (stack.size()) {
const subset& C = *stack.back();
subset C = *stack.back();
/*const subset& C = *stack.front();*/
/*std::cout << "** " << C << std::endl;*/
/*std::cout << '.' << std::flush;*/
......@@ -215,6 +217,7 @@ struct lumper {
/*}*/
key_map subsets = compute_keys(C, B);
/*MSG_DEBUG("KEYS: " << subsets);*/
/*MSG_DEBUG("KEYS " << labels[subsets.begin()->second.front()] << ": " << subsets);*/
/*size_t size = 0;*/
/*while (size != subsets.size()) {*/
/*size = subsets.size();*/
......
This diff is collapsed.
This diff is collapsed.
......@@ -25,7 +25,7 @@ COMPUTATIONS_SRC=basic_data.cc probabilities.cc model.cc frontends.cc
SRC=$(addprefix input/,$(INPUT_SRC)) $(addprefix computations/,$(COMPUTATIONS_SRC)) $(MAIN_SRC)
OBJ=$(subst .cc,.o,$(SRC))
DEP=$(subst .cc,.d,$(SRC)) test_script.d matrixexp.d matrixexp3.d
DEP=$(subst .cc,.d,$(SRC)) test_script.d matrixexp4.d
COV_OBJ=$(subst .cc,.cov.o,$(SRC))
ALL_BUT_MAIN_OBJ=$(subst main.o ,,$(OBJ))
......@@ -164,5 +164,5 @@ test_bayes_dai: ../sandbox/bayes.cc bayes.d $(ALL_BUT_MAIN_OBJ)
test_sib: pedigree_analysis.cc static_data.o
$C $< static_data.o -o $@
matrixexp: matrixexp3.cc ../include/lumping2.h
matrixexp: matrixexp4.cc ../include/lumping2.h
$C -DORDER=$(ORDER) $< -o $@
for (i in (2:30)*2) {
x <- read.table(paste('min_angles_', i, '.txt', sep=''), header=T)
cat("Best gen for", i, "states:", x$gen[which.max(x$max - x$min)], "\n")
}
#values-AU.F2_inherit.txt
#values-AU.F2_over.txt
#values-AU.F2_normal.txt
#values-AU.F3_inherit.txt
#values-AU.F3_over.txt
#values-AU.F3_normal.txt
#values-AU.F4_inherit.txt
#values-AU.F4_over.txt
#values-AU.F5_normal.txt
#values-AU.F4_normal.txt
#values-AU.F5_inherit.txt
#values-AU.F5_over.txt
#values-AU.F6_normal.txt
#values-AU.F6_inherit.txt
#values-AU.F6_over.txt
#values-AU.F7_normal.txt
#values-AU.F7_inherit.txt
#values-AU.F7_over.txt
#
A.col <- "red"
B.col <- "blue"
H.col <- "green"
normal.type <- "p"
over.type <- "l"
inherit.type <- "l"
normal.lty <- "blank"
over.lty <- "solid"
inherit.lty <- "dashed"
all_cols <- function(df, prefix) colnames(df)[which(substr(colnames(df), 1, nchar(prefix)) == prefix)]
sum_cols <- function(df, prefix) {
ac <- unlist(lapply(c(prefix), function(p) all_cols(df, p)))
if (length(ac) > 1) rowSums(df[, ac]) else df[, ac]
}
f <- function(gen, ...) {
g.normal <- read.table(paste("values-AU.F", gen, "_normal.txt", sep=""), header=T)
g.over <- read.table(paste("values-AU.F", gen, "_over.txt", sep=""), header=T)
g.inherit <- read.table(paste("values-AU.F", gen, "_inherit.txt", sep=""), header=T)
plot(sum_cols(g.normal, 'bb'), ..., type=normal.type, lty=normal.lty, col=B.col)
lines(sum_cols(g.over, 'bb'), type=over.type, lty=over.lty, col=B.col)
lines(sum_cols(g.inherit, 'bb'), type=inherit.type, lty=inherit.lty, col=B.col)
lines(sum_cols(g.normal, 'aa'), type=normal.type, lty=normal.lty, col=A.col)
lines(sum_cols(g.over, 'aa'), type=over.type, lty=over.lty, col=A.col)
lines(sum_cols(g.inherit, 'aa'), type=inherit.type, lty=inherit.lty, col=A.col)
lines(sum_cols(g.normal, c('ab', 'ba')), type=normal.type, lty=normal.lty, col=H.col)
lines(sum_cols(g.over, c('ab', 'ba')), type=over.type, lty=over.lty, col=H.col)
lines(sum_cols(g.inherit, c('ab', 'ba')), type=inherit.type, lty=inherit.lty, col=H.col)
}
This diff is collapsed.
source('make_plots_RIL.R')
make.plot(5, 5, 'AAA')
make.plot(5, 1, 'AAB')
dim(all.data)
all.data[,, 1,1,1, 6, 6)
all.data[,, 1,1,1, 6, 6]
all.data[1:26,, 1,1,1, 6, 6]
26:1
all.data[1:26,, 1,1,1, 6, 6] - all.data[26:1,, 1,1,1, 6,6]
plot(all.data[1:26,, 1,1,1, 6, 6] - all.data[26:1,1, 1,1,1, 6,6])
plot(all.data[1:26,1, 1,1,1, 6, 6] - all.data[26:1,1, 1,1,1, 6,6])
all.data[1:26,, 1,3,1, 6, 6]
all.data[1:26,, 1,3,1, 6, 1]
all.data[1:26,, 1,3,1, 6, 2]
ABA.62.A <- all.data[1:26,1, 1,3,1, 6, 2]
ABA.62.B <- all.data[1:26,3, 1,3,1, 6, 2]
rev(1:10)
ABA.62.A - rev(ABA.62.B)
plot(ABA.62.A - rev(ABA.62.B))
AAA.62.A <- all.data[1:26,1, 1,1,1, 6, 2]
AAB.62.A <- all.data[1:26,1, 1,1,3, 6, 2]
AAA.62.A - AAB.62.A
plot(AAA.62.A - AAB.62.A)
source('make_plots_RIL.R')
aba.62 <- get.data('ABA', 6, 2)
abb.62 <- get.data('ABB', 6, 2)
aba.62$A
aba.62[,1]
aba.62 - abb.62
plot((aba.62 - abb.62)[,1])
aba.11.2 <- get.data('ABA', 11, 2)
aba.10.1 <- get.data('ABA', 10, 1)
abb.10.1 <- get.data('ABB', 10, 1)
plot((aba.10.1 - abb.10.1)[,1])
plot(aba.10.1[,1])
aaa.10.1 <- get.data('AAA', 10, 1)
aab.10.1 <- get.data('AAB', 10, 1)
plot(aaa.10.1[,1])
plot((aaa.10.1 - aab.10.1)[,1])
make.plot(10, 1, 'AAA')
make.plot(10, 1, 'AAB')
make.plot(10, 1, 'AAA')
make.plot(10, 1, 'AAB')
make.plot
dev.off()
source('make_plots_RIL.R')
make.plot(10, 1, 'AAB')
make.plot(10, 1, 'AAA')
make.plot(10, 1, 'AAB')
make.plot(10, 1, 'AAA')
make.plot(10, 1, 'AAB')
make.plot(10, 1, 'AAA')
make.plot(10, 1, 'AAB')
make.plot(10, 1, 'AAA')
make.plot(10, 1, 'AAB')
make.plot(10, 1, 'AAA')
make.plot(10, 1, 'AAB')
make.plot(10, 1, 'AAA')
make.plot(10, 1, 'AAB')
make.plot(10, 1, 'ABB')
make.plot(10, 1, 'ABA')
make.plot(10, 1, 'ABB')
make.plot(10, 1, 'ABA')
make.plot(10, 1, 'ABB')
make.plot(10, 1, 'ABA')
make.plot(3, 1, 'ABA')
make.plot(3, 1, 'ABB')
make.plot(3, 1, 'ABA')
make.plot(3, 1, 'ABB')
make.plot(3, 1, 'ABA')
make.plot(3, 1, 'ABB')
make.plot(3, 1, 'ABA')
make.plot(3, 1, 'ABB')
make.plot(3, 1, 'ABA')
make.plot(3, 1, 'ABB')
make.plot(3, 1, 'ABA')
make.plot(3, 1, 'ABB')
make.plot(1, 1, 'ABB')
make.plot(1, 1, 'ABA')
make.plot(1, 1, 'ABB')
make.plot(1, 1, 'ABA')
make.plot(1, 1, 'ABB')
make.plot(1, 1, 'ABA')
make.plot(1, 1, 'ABB')
sapply(0:10, function(r1) sapply(0:10, function(r2) max(abs(get.data('AAA', r1, r2) - get.data('AAB', r1, r2)))
sapply(0:10, function(r1) sapply(0:10, function(r2) max(abs(get.data('AAA', r1, r2) - get.data('AAB', r1, r2))))
sapply(0:10, function(r1) sapply(0:10, function(r2) max(abs(get.data('AAA', r1, r2) - get.data('AAB', r1, r2)))))
which.max(sapply(0:10, function(r1) sapply(0:10, function(r2) max(abs(get.data('AAA', r1, r2) - get.data('AAB', r1, r2))))), arr.ind=T)
which.max(sapply(0:10, function(r1) sapply(0:10, function(r2) max(abs(get.data('AAA', r1, r2) - get.data('AAB', r1, r2))))))
aa <- sapply(0:10, function(r1) sapply(0:10, function(r2) max(abs(get.data('AAA', r1, r2) - get.data('AAB', r1, r2)))))
ab <- sapply(0:10, function(r1) sapply(0:10, function(r2) max(abs(get.data('ABA', r1, r2) - get.data('ABB', r1, r2)))))
which(aa == max(aa), arr.ind=T)
which(ab == max(ab), arr.ind=T)
ab
aa <- sapply(0:10, function(r1) sapply(0:10, function(r2) max(abs(get.data('AAA', r1, r2) - get.data('AAB', r1, r2)))))[-1,-1]
ab <- sapply(0:10, function(r1) sapply(0:10, function(r2) max(abs(get.data('ABA', r1, r2) - get.data('ABB', r1, r2)))))[-1,-1]
which(aa == max(aa), arr.ind=T)
which(ab == max(ab), arr.ind=T)
make.plot(10, 1, 'AAB')
make.plot(10, 1, 'AAA')
make.plot(10, 1, 'AAA')
make.plot(10, 1, 'AAA')
make.plot(10, 1, 'AAA')
make.plot(10, 1, 'AAB')
make.plot(10, 1, 'AAA')
make.plot(10, 1, 'AAB')
make.plot(10, 1, 'AAA')
make.plot(10, 1, 'AAB')
make.plot(10, 1, 'AAA')
make.plot(10, 1, 'AAB')
make.plot(10, 1, 'AAA')
sleep(.1)
sleep
for (i in 1:1000000);
for (i in 1:1000000) 0
for (i in 1:1000000000) 0
for (i in 1:100000000) 0
for (i in 1:10000000) 0
for (i in 1:50000000) 0
for (i in 1:20000000) 0
for (k in 1:100) { for (i in 1:20000000) 0; make.plot(10, 1, 'AAB'); for (i in 1:20000000) 0; make.plot(10, 1, 'AAA'); }
for (k in 1:100) { for (i in 1:20000000) 0; make.plot(10, 1, 'AAB'); for (i in 1:20000000) 0; make.plot(10, 1, 'AAA'); }
for (k in 1:100) { for (i in 1:20000000) 0; make.plot(10, 1, 'AAB'); for (i in 1:20000000) 0; make.plot(10, 1, 'AAA'); }
for (k in 1:100) { for (i in 1:20000000) 0; make.plot(10, 1, 'AAB') }
for (k in 1:100) { make.plot(10, 1, 'AAB') }
traceback()
source('make_plots_RIL.R')
for (k in 1:100) { for (i in 1:20000000) 0; make.plot(10, 1, 'AAB'); for (i in 1:20000000) 0; make.plot(10, 1, 'AAA'); }
for (k in 1:10000) { for (i in 1:20000000) 0; make.plot(10, 1, 'AAB'); for (i in 1:20000000) 0; make.plot(10, 1, 'AAA'); }
dev.off()
aas <- sapply(0:10, function(r1) sapply(0:10, function(r2) sum(abs(get.data('AAA', r1, r2) - get.data('AAB', r1, r2)))))[-1,-1]
aas <- sapply(1:10, function(r1) sapply(1:10, function(r2) sum(abs(get.data('AAA', r1, r2) - get.data('AAB', r1, r2)))))
aa.s <- sapply(0:10, function(r1) sapply(0:10, function(r2) sum(abs(get.data('AAA', r1, r2) - get.data('AAB', r1, r2)))))[-1,-1]
ab.s <- sapply(1:10, function(r1) sapply(1:10, function(r2) sum(abs(get.data('ABA', r1, r2) - get.data('ABB', r1, r2)))))
which(ab.s == max(ab.s), arr.ind=T)
which(aa.s == max(aa.s), arr.ind=T)
for (k in 1:10000) { for (i in 1:20000000) 0; make.plot(8, 1, 'AAB'); for (i in 1:20000000) 0; make.plot(8, 1, 'AAA'); }
for (k in 1:10000) { for (i in 1:20000000) 0; make.plot(8, 1, 'ABB'); for (i in 1:20000000) 0; make.plot(8, 1, 'ABA'); }
for (k in 1:10000) { for (i in 1:20000000) 0; make.plot(8, 1, 'ABB'); for (i in 1:20000000) 0; make.plot(8, 1, 'ABA'); }
source('make_plots_RIL.R')
make.plot(5, 5, 'AAA')
make.plot(5, 1, 'AAB')
dim(all.data)
all.data[,, 1,1,1, 6, 6)
all.data[,, 1,1,1, 6, 6]
all.data[1:26,, 1,1,1, 6, 6]
26:1
all.data[1:26,, 1,1,1, 6, 6] - all.data[26:1,, 1,1,1, 6,6]
plot(all.data[1:26,, 1,1,1, 6, 6] - all.data[26:1,1, 1,1,1, 6,6])
plot(all.data[1:26,1, 1,1,1, 6, 6] - all.data[26:1,1, 1,1,1, 6,6])
all.data[1:26,, 1,3,1, 6, 6]
all.data[1:26,, 1,3,1, 6, 1]
all.data[1:26,, 1,3,1, 6, 2]
ABA.62.A <- all.data[1:26,1, 1,3,1, 6, 2]
ABA.62.B <- all.data[1:26,3, 1,3,1, 6, 2]
rev(1:10)
ABA.62.A - rev(ABA.62.B)
plot(ABA.62.A - rev(ABA.62.B))
AAA.62.A <- all.data[1:26,1, 1,1,1, 6, 2]
AAB.62.A <- all.data[1:26,1, 1,1,3, 6, 2]
AAA.62.A - AAB.62.A
plot(AAA.62.A - AAB.62.A)
source('make_plots_RIL.R')
aba.62 <- get.data('ABA', 6, 2)
abb.62 <- get.data('ABB', 6, 2)
aba.62$A
aba.62[,1]
aba.62 - abb.62
plot((aba.62 - abb.62)[,1])
aba.11.2 <- get.data('ABA', 11, 2)
aba.10.1 <- get.data('ABA', 10, 1)
abb.10.1 <- get.data('ABB', 10, 1)
plot((aba.10.1 - abb.10.1)[,1])
plot(aba.10.1[,1])
aaa.10.1 <- get.data('AAA', 10, 1)
aab.10.1 <- get.data('AAB', 10, 1)
plot(aaa.10.1[,1])
plot((aaa.10.1 - aab.10.1)[,1])
make.plot(10, 1, 'AAA')
make.plot(10, 1, 'AAB')
make.plot(10, 1, 'AAA')
make.plot(10, 1, 'AAB')
make.plot
dev.off()
source('make_plots_RIL.R')
make.plot(10, 1, 'AAB')
make.plot(10, 1, 'AAA')
make.plot(10, 1, 'AAB')
make.plot(10, 1, 'AAA')
make.plot(10, 1, 'AAB')
make.plot(10, 1, 'AAA')
make.plot(10, 1, 'AAB')
make.plot(10, 1, 'AAA')
make.plot(10, 1, 'AAB')
make.plot(10, 1, 'AAA')
make.plot(10, 1, 'AAB')
make.plot(10, 1, 'AAA')
make.plot(10, 1, 'AAB')
make.plot(10, 1, 'ABB')
make.plot(10, 1, 'ABA')
make.plot(10, 1, 'ABB')
make.plot(10, 1, 'ABA')
make.plot(10, 1, 'ABB')
make.plot(10, 1, 'ABA')
make.plot(3, 1, 'ABA')
make.plot(3, 1, 'ABB')
make.plot(3, 1, 'ABA')
make.plot(3, 1, 'ABB')
make.plot(3, 1, 'ABA')
make.plot(3, 1, 'ABB')
make.plot(3, 1, 'ABA')
make.plot(3, 1, 'ABB')
make.plot(3, 1, 'ABA')
make.plot(3, 1, 'ABB')
make.plot(3, 1, 'ABA')
make.plot(3, 1, 'ABB')
make.plot(1, 1, 'ABB')
make.plot(1, 1, 'ABA')
make.plot(1, 1, 'ABB')
make.plot(1, 1, 'ABA')
make.plot(1, 1, 'ABB')
make.plot(1, 1, 'ABA')
make.plot(1, 1, 'ABB')
sapply(0:10, function(r1) sapply(0:10, function(r2) max(abs(get.data('AAA', r1, r2) - get.data('AAB', r1, r2)))
sapply(0:10, function(r1) sapply(0:10, function(r2) max(abs(get.data('AAA', r1, r2) - get.data('AAB', r1, r2))))
sapply(0:10, function(r1) sapply(0:10, function(r2) max(abs(get.data('AAA', r1, r2) - get.data('AAB', r1, r2)))))
which.max(sapply(0:10, function(r1) sapply(0:10, function(r2) max(abs(get.data('AAA', r1, r2) - get.data('AAB', r1, r2))))), arr.ind=T)
which.max(sapply(0:10, function(r1) sapply(0:10, function(r2) max(abs(get.data('AAA', r1, r2) - get.data('AAB', r1, r2))))))
aa <- sapply(0:10, function(r1) sapply(0:10, function(r2) max(abs(get.data('AAA', r1, r2) - get.data('AAB', r1, r2)))))
ab <- sapply(0:10, function(r1) sapply(0:10, function(r2) max(abs(get.data('ABA', r1, r2) - get.data('ABB', r1, r2)))))
which(aa == max(aa), arr.ind=T)
which(ab == max(ab), arr.ind=T)
ab
aa <- sapply(0:10, function(r1) sapply(0:10, function(r2) max(abs(get.data('AAA', r1, r2) - get.data('AAB', r1, r2)))))[-1,-1]
ab <- sapply(0:10, function(r1) sapply(0:10, function(r2) max(abs(get.data('ABA', r1, r2) - get.data('ABB', r1, r2)))))[-1,-1]
which(aa == max(aa), arr.ind=T)
which(ab == max(ab), arr.ind=T)
make.plot(10, 1, 'AAB')
make.plot(10, 1, 'AAA')
make.plot(10, 1, 'AAA')
make.plot(10, 1, 'AAA')
make.plot(10, 1, 'AAA')
make.plot(10, 1, 'AAB')
make.plot(10, 1, 'AAA')
make.plot(10, 1, 'AAB')
make.plot(10, 1, 'AAA')
make.plot(10, 1, 'AAB')
make.plot(10, 1, 'AAA')
make.plot(10, 1, 'AAB')
make.plot(10, 1, 'AAA')
sleep(.1)
sleep
for (i in 1:1000000);
for (i in 1:1000000) 0
for (i in 1:1000000000) 0
for (i in 1:100000000) 0
for (i in 1:10000000) 0
for (i in 1:50000000) 0
for (i in 1:20000000) 0
for (k in 1:100) { for (i in 1:20000000) 0; make.plot(10, 1, 'AAB'); for (i in 1:20000000) 0; make.plot(10, 1, 'AAA'); }
for (k in 1:100) { for (i in 1:20000000) 0; make.plot(10, 1, 'AAB'); for (i in 1:20000000) 0; make.plot(10, 1, 'AAA'); }
for (k in 1:100) { for (i in 1:20000000) 0; make.plot(10, 1, 'AAB'); for (i in 1:20000000) 0; make.plot(10, 1, 'AAA'); }
for (k in 1:100) { for (i in 1:20000000) 0; make.plot(10, 1, 'AAB') }
for (k in 1:100) { make.plot(10, 1, 'AAB') }
traceback()
source('make_plots_RIL.R')
for (k in 1:100) { for (i in 1:20000000) 0; make.plot(10, 1, 'AAB'); for (i in 1:20000000) 0; make.plot(10, 1, 'AAA'); }
for (k in 1:10000) { for (i in 1:20000000) 0; make.plot(10, 1, 'AAB'); for (i in 1:20000000) 0; make.plot(10, 1, 'AAA'); }
dev.off()
aas <- sapply(0:10, function(r1) sapply(0:10, function(r2) sum(abs(get.data('AAA', r1, r2) - get.data('AAB', r1, r2)))))[-1,-1]
aas <- sapply(1:10, function(r1) sapply(1:10, function(r2) sum(abs(get.data('AAA', r1, r2) - get.data('AAB', r1, r2)))))
aa.s <- sapply(0:10, function(r1) sapply(0:10, function(r2) sum(abs(get.data('AAA', r1, r2) - get.data('AAB', r1, r2)))))[-1,-1]
ab.s <- sapply(1:10, function(r1) sapply(1:10, function(r2) sum(abs(get.data('ABA', r1, r2) - get.data('ABB', r1, r2)))))
which(ab.s == max(ab.s), arr.ind=T)
which(aa.s == max(aa.s), arr.ind=T)
for (k in 1:10000) { for (i in 1:20000000) 0; make.plot(8, 1, 'AAB'); for (i in 1:20000000) 0; make.plot(8, 1, 'AAA'); }
for (k in 1:10000) { for (i in 1:20000000) 0; make.plot(8, 1, 'ABB'); for (i in 1:20000000) 0; make.plot(8, 1, 'ABA'); }
for (k in 1:10000) { for (i in 1:20000000) 0; make.plot(8, 1, 'ABB'); for (i in 1:20000000) 0; make.plot(8, 1, 'ABA'); }
This diff is collapsed.
This diff is collapsed.
all.data <- array(dim=c(27, 3, 3, 3, 3, 11, 11, 19))
all.counts <- array(dim=c(27, 3, 3, 3, 3, 11, 11, 19))
L <- c('A', 'H', 'B')
for (g in 2:20)
for (d1 in 0:10)
for (d2 in 0:10) {
cat("reading gen_", g, " ", d1, "/", d2, " \r", sep="")
for (m1 in 1:3)
for (m2 in 1:3)
for (m3 in 1:3) {
path <- paste('gen_', g, '/counts_10/', d1, '/', d2, '/', L[m1], L[m2], L[m3], '.txt', sep='')
t <- read.table(path, sep="\t", header=T)
if (nrow(t) > 0) {
m <- as.matrix(t[, -1])
all.counts[,, m1,m2,m3, d1 + 1, d2 + 1, g - 1] <- m
m <- m / rowSums(m)
all.data[,, m1,m2,m3, d1 + 1, d2 + 1, g - 1] <- m
}
}
}
cat(" \n")
all.data[is.nan(all.data)] <- 0
all.data[is.na(all.data)] <- 0
all.counts[is.nan(all.counts)] <- 0
all.counts[is.na(all.counts)] <- 0
inv.h <- function(r) if (r < .5) round(50 * log(1 / (1 - 2 * r)), 3) else "+∞"
make.plot <- function(g, d1, d2, m123) {
m1 <- which(L == substr(m123, 1, 1))
m2 <- which(L == substr(m123, 2, 2))
m3 <- which(L == substr(m123, 3, 3))
tab <- all.data[1:26,, m1, m2, m3, 1 + d1, 1 + d2, g]
if (all(is.nan(tab))) {
cat("No data.\n")
return(invisible())
}
titolo <- paste("Genotype frequencies between M1 and M2 in generation ", g, "\nM1[", L[m1], "]--", inv.h(d1 / 20), "cM--M2[", L[m2], "]--", inv.h(d2 / 20), "cM--M3[", L[m3], "]", sep="")
plot(tab[, 1], ylab="Frequency", ylim=c(0, 1), col="blue", type="p", main=titolo)
lines(tab[, 2], col="green", type="p"); lines(tab[, 3], col="red", type="p")
}
all.data <- array(dim=c(28, 2, 2, 2, 2, 11, 11))
all.counts <- array(dim=c(28, 2, 2, 2, 2, 11, 11))
L <- c('A', 'H')
for (d1 in 0:10)
for (d2 in 0:10) {
cat("reading gen_BC2 ", d1, "/", d2, " \r", sep="")
for (m1 in 1:2)
for (m2 in 1:2)
for (m3 in 1:2) {
path <- paste('gen_BC2/counts_10/', d1, '/', d2, '/', L[m1], L[m2], L[m3], '.txt', sep='')
t <- read.table(path, sep="\t", header=T)
if (nrow(t) > 0) {
m <- as.matrix(t[, c(-1, -4)])
all.counts[,, m1,m2,m3, d1 + 1, d2 + 1] <- m
m <- m / rowSums(m)
all.data[,, m1,m2,m3, d1 + 1, d2 + 1] <- m
}
}
}
cat(" \n")
all.data[is.nan(all.data)] <- 0
all.data[is.na(all.data)] <- 0
all.counts[is.nan(all.counts)] <- 0
all.counts[is.na(all.counts)] <- 0
inv.h <- function(r) if (r < .5) round(50 * log(1 / (1 - 2 * r)), 3) else "+∞"
.d1 <- 0
make.plot <- function(d1, d2, m123) {
m1 <- which(L == substr(m123, 1, 1))
m2 <- which(L == substr(m123, 2, 2))
m3 <- which(L == substr(m123, 3, 3))
tab <- all.data[1:27,, m1, m2, m3, 1 + d1, 1 + d2]
if (all(is.nan(tab))) {
cat("No data.\n")
return(invisible())
}
.d1 <<- inv.h(d1 / 20)
titolo <- paste("Genotype frequencies between M1 and M2 in generation BC2\nM1[", L[m1], "]--", .d1, "cM--M2[", L[m2], "]--", inv.h(d2 / 20), "cM--M3[", L[m3], "]", sep="")
x <- seq(0, .d1, length.out=27)
plot(x, tab[, 1], ylab="Frequency", ylim=c(0, 1), col="blue", type="p", main=titolo)
lines(x, tab[, 2], col="green", type="p")