Skip to content

Commit

Permalink
Merge pull request #130 from stevenhwu/feature_update_peeling
Browse files Browse the repository at this point in the history
Update peeling algorithm and test cases
reedacartwright committed Mar 26, 2016

Verified

This commit was created on GitHub.com and signed with GitHub’s verified signature.
2 parents 00209fb + deb35fe commit b4256ad
Showing 3 changed files with 292 additions and 32 deletions.
20 changes: 10 additions & 10 deletions src/lib/peeling.cc
Original file line number Diff line number Diff line change
@@ -70,7 +70,7 @@ void dng::peel::to_father(workspace_t &work, const family_members_t &family,
}
// Include Mom
work.paired_buffer.resize(10, 10);
work.lower[dad] *= (work.paired_buffer.matrix() * (work.upper[mom] *
work.lower[dad] *= (work.paired_buffer.matrix().transpose() * (work.upper[mom] *
work.lower[mom]).matrix()).array();
work.paired_buffer.resize(100, 1);
}
@@ -88,7 +88,7 @@ void dng::peel::to_father_fast(workspace_t &work,
}
// Include Mom
work.paired_buffer.resize(10, 10);
work.lower[dad] = (work.paired_buffer.matrix() * (work.upper[mom] *
work.lower[dad] = (work.paired_buffer.matrix().transpose() * (work.upper[mom] *
work.lower[mom]).matrix()).array();
work.paired_buffer.resize(100, 1);
}
@@ -106,7 +106,7 @@ void dng::peel::to_mother(workspace_t &work, const family_members_t &family,
}
// Include Dad
work.paired_buffer.resize(10, 10);
work.lower[mom] *= (work.paired_buffer.matrix().transpose() * (work.upper[dad] *
work.lower[mom] *= (work.paired_buffer.matrix() * (work.upper[dad] *
work.lower[dad]).matrix()).array();
work.paired_buffer.resize(100, 1);
}
@@ -124,7 +124,7 @@ void dng::peel::to_mother_fast(workspace_t &work,
}
// Include Dad
work.paired_buffer.resize(10, 10);
work.lower[mom] = (work.paired_buffer.matrix().transpose() * (work.upper[dad] *
work.lower[mom] = (work.paired_buffer.matrix() * (work.upper[dad] *
work.lower[dad]).matrix()).array();
work.paired_buffer.resize(100, 1);
}
@@ -204,14 +204,14 @@ void dng::peel::to_father_reverse(workspace_t &work,
IndividualVector::value_type mom_v = work.upper[mom] * work.lower[mom];

// Calculate P(dad-only data & dad = g)
work.paired_buffer.resize(10, 10);
work.paired_buffer.resize(10, 10); //TODO: FIXME: add transpose()
IndividualVector::value_type dad_v = work.upper[dad] * (work.lower[dad] /
((work.paired_buffer.matrix() * mom_v.matrix()).array() +
DNG_INDIVIDUAL_BUFFER_MIN));

// Calculate P(dependent data | mom = g)
work.lower[mom] *= (work.paired_buffer.matrix().transpose() *
dad_v.matrix()).array();
dad_v.matrix()).array(); //TODO: FIXME: remove transpose()
work.paired_buffer.resize(100, 1);

// Calculate P(data & dad & mom)
@@ -238,11 +238,11 @@ void dng::peel::to_mother_reverse(workspace_t &work,
}
IndividualVector::value_type dad_v = work.upper[dad] * work.lower[dad];

work.paired_buffer.resize(10, 10);
work.paired_buffer.resize(10, 10); //TODO: FIXME: remove transpose()
IndividualVector::value_type mom_v = work.upper[mom] * (work.lower[mom] /
((work.paired_buffer.matrix().transpose() * dad_v.matrix()).array() +
DNG_INDIVIDUAL_BUFFER_MIN));
work.lower[dad] *= (work.paired_buffer.matrix() * mom_v.matrix()).array();
work.lower[dad] *= (work.paired_buffer.matrix() * mom_v.matrix()).array(); //TODO: FIXME: add transpose()
work.paired_buffer.resize(100, 1);

work.paired_buffer *= kroneckerProduct(dad_v.matrix(), mom_v.matrix()).array();
@@ -271,9 +271,9 @@ void dng::peel::to_child_reverse(workspace_t &work,
IndividualVector::value_type dad_v = work.upper[dad] * work.lower[dad];
IndividualVector::value_type mom_v = work.upper[mom] * work.lower[mom];
work.paired_buffer.resize(10, 10);
work.lower[dad] *= (work.paired_buffer.matrix() * mom_v.matrix()).array();
work.lower[dad] *= (work.paired_buffer.matrix() * mom_v.matrix()).array(); //TODO: FIXME: add transpose()
work.lower[mom] *= (work.paired_buffer.matrix().transpose() *
dad_v.matrix()).array();
dad_v.matrix()).array(); //TODO: FIXME: remove transpose()
work.paired_buffer.resize(100, 1);

// Update Siblings
81 changes: 81 additions & 0 deletions test/script/R/generate_family_data.R
Original file line number Diff line number Diff line change
@@ -0,0 +1,81 @@
options(digits=20)

set.seed(0)
dip_count<- 10

# family order: dad, mom, child1, child2
lower<- sapply(1:4, function(x){runif(dip_count)}, simplify=F)
upper<- sapply(1:4, function(x){runif(dip_count)}, simplify=F)

transition<- sapply(1:2, function(x){ matrix(runif(dip_count^2), nrow=dip_count, byrow=T) }, simplify=F)
transition[[3]]<- matrix(runif(dip_count^3), nrow=dip_count*dip_count, byrow=T)
transition[[4]]<- matrix(runif(dip_count^3), nrow=dip_count*dip_count, byrow=T)

## up
up_fast <- transition[[2]] %*% lower[[2]]
up<- lower[[1]] * up_fast

## to father
buffer <- transition[[3]] %*% lower[[3]]
buffer_matrix <- matrix(buffer,nrow=dip_count, byrow=T)
mom_lower_upper <- lower[[2]] * upper[[2]]
father_fast <- buffer_matrix %*% mom_lower_upper
father <- lower[[1]] * father_fast

## to mother
buffer <- transition[[3]] %*% lower[[3]]
buffer_matrix <- matrix(buffer,nrow=dip_count, byrow=F)##by column here
dad_lower_upper <- lower[[1]] * upper[[1]]
mother_fast <- buffer_matrix %*% dad_lower_upper
mother <- lower[[2]] * mother_fast

## to child
buffer<- (lower[[1]]*upper[[1]]) %x% ((lower[[2]]*upper[[2]]))
other_child <- transition[[4]] %*% lower[[4]]
updated_buffer <- buffer * as.vector(other_child)

child_fast<- t(transition[[3]]) %*% buffer
child<- t(transition[[3]]) %*% updated_buffer

## output to c++
formatList<-function(name,v){
paste(name, "[", 0:(length(v)-1), "] << ",
sapply(v,function(x){
paste(x,collapse=", ")
}),
";", sep="")
}

formatVector<-function(name,v){
paste(name, " << ",
paste(v,collapse=", "),
";", sep="")
}

outputToCPP<- function(filename=""){
cat("", fill=T, file=filename)
cat(formatList("lower_array", lower), fill=T, file=filename, append=T)
cat(formatList("upper_array", upper), fill=T, file=filename, append=T)

cat(paste("trans_matrix[",0:(length(transition)-1),"] << ",
sapply(transition,function(x){
paste(t(x),collapse=", ") ## eigen << init take it by row
}),
";", sep=""),
fill=T, file=filename, append=T)

cat(formatVector("expected_up_fast", up_fast), fill=T, file=filename, append=T)
cat(formatVector("expected_up", up), fill=T, file=filename, append=T)

cat(formatVector("expected_father_fast", father_fast), fill=T, file=filename, append=T)
cat(formatVector("expected_father", father), fill=T, file=filename, append=T)

cat(formatVector("expected_mother_fast", mother_fast), fill=T, file=filename, append=T)
cat(formatVector("expected_mother", mother), fill=T, file=filename, append=T)

cat(formatVector("expected_child_fast", child_fast), fill=T, file=filename, append=T)
cat(formatVector("expected_child", child), fill=T, file=filename, append=T)
}

outputToCPP("generateFamilyData")

223 changes: 201 additions & 22 deletions test/src/lib/test2_peeling.cpp

Large diffs are not rendered by default.

0 comments on commit b4256ad

Please sign in to comment.