Permalink
Browse files

Seems to be working in it's entirety.. Need to check that probabiliti…

…es are being correctly calculated.
  • Loading branch information...
1 parent 3944ab0 commit 5f930de9e6f70fe904633e1eefd3b485abf64a70 @cstrope committed Oct 5, 2012
Showing with 50 additions and 73 deletions.
  1. BIN data/Markov/isg-test-one
  2. +9 −9 src/dependency.cpp
  3. +12 −28 src/motif.cpp
  4. +0 −4 src/propose_path.cpp
  5. +0 −4 src/rate_type.cpp
  6. +29 −28 src/treefile.cpp
View
Binary file not shown.
View
@@ -930,21 +930,21 @@ contextDependence::setOffset()
}
}
- cerr << "Offset setup: " << endl;
- int p = 0;
- for (at = index_offset.begin(); at != index_offset.end(); ++at, ++p) {
- cerr << "Level " << p << ": " << endl;
- for (vector<int>::iterator bt = (*at).begin(); bt != (*at).end(); ++bt) {
- cerr << (*bt) << " ";
- }
- cerr << endl;
+ //cerr << "Offset setup: " << endl;
+ //int p = 0;
+ //for (at = index_offset.begin(); at != index_offset.end(); ++at, ++p) {
+ // cerr << "Level " << p << ": " << endl;
+ // for (vector<int>::iterator bt = (*at).begin(); bt != (*at).end(); ++bt) {
+ // cerr << (*bt) << " ";
+ // }
+ // cerr << endl;
// for (int j = 0; j < numStates; ++j) { // current == sequence i
// for (int k = 0; k < numStates; ++k) { // proposed == sequence j
// cerr << (*at).at(j*numStates+k) << " ";
// }
// cerr << endl;
// }
- }
+ //}
//exit(0);
}
View
@@ -145,26 +145,10 @@ Sequence::zeroStopCodons (
// If we end up zeroing some possible transition, this needs to be reflected in the overall sequence probability.
// For this, we return the amount that the sequence needs to be adjusted by.
double probability_adjustment = 0;
-/* // First, set the event site to the nearest codon boundary.
- cerr << "???event_site(" << event_site << ") -= 3 + " << (event_site%3) << endl;
- event_site -= 3 + (event_site%3);
- if (event_site < 0) event_site = 0;
-
- // and the end site similarly.
- cerr << "???end_site(" << end_site << ") += 3 + " << (3-(end_site%3)) << " ----> ";
- end_site += 3 + (3-(end_site%3)); // ? End site at the end of a codon should be a 2... but if less than end site?
-
-// cerr << "???now event_site=" << event_site << " and end_site=" << end_site << endl;
-*/
-
vector<Site>::iterator end;
- if (end_site > evolutionaryAttributes.size()) {
- end = evolutionaryAttributes.end();
- cerr << "end" << endl;
- } else {
- end = evolutionaryAttributes.begin()+end_site;
- cerr << end_site << endl;
- }
+
+ if (end_site > evolutionaryAttributes.size()) end = evolutionaryAttributes.end();
+ else end = evolutionaryAttributes.begin()+end_site;
event_site = 0;
end = evolutionaryAttributes.end();
@@ -175,7 +159,7 @@ Sequence::zeroStopCodons (
codon[1] = (*(it+1)).returnState();
codon[2] = (*(it+2)).returnState();
- cerr << "site " << i << ": " << stateCharacters.at(codon[0]) << stateCharacters.at(codon[1]) << stateCharacters.at(codon[2]) << " ";
+ //cerr << "site " << i << ": " << stateCharacters.at(codon[0]) << stateCharacters.at(codon[1]) << stateCharacters.at(codon[2]) << " ";
// All stop codons start with T, this checks for that.
if (codon[0] == 3) { // T
@@ -185,39 +169,39 @@ Sequence::zeroStopCodons (
if (codon[1] == 0) { // TAX
probability_adjustment += (*(it+2)).site_rate_away.at(0);
(*(it+2)).site_rate_away.at(0) = 0;
- cerr << "G2->0 ";
+ //cerr << "G2->0 ";
} else if (codon[1] == 2) { // TGX
probability_adjustment += (*(it+2)).site_rate_away.at(2);
(*(it+2)).site_rate_away.at(2) = 0;
- cerr << "A2->0 ";
+ //cerr << "A2->0 ";
} else if (codon[2] == 0) { // TXA -> can change to TAA or TGA
probability_adjustment
+= (*(it+1)).site_rate_away.at(0)
+ (*(it+1)).site_rate_away.at(2);
(*(it+1)).site_rate_away.at(0) = (*(it+1)).site_rate_away.at(2) = 0;
- cerr << "{A1,G1}->0 ";
+ //cerr << "{A1,G1}->0 ";
} else if (codon[2] == 2) { // TXG
probability_adjustment += (*(it+1)).site_rate_away.at(2);
(*(it+1)).site_rate_away.at(2) = 0;
- cerr << "A1->0 ";
+ //cerr << "A1->0 ";
} // else this codon is fine.
// Now check to see if second and third codon positions are stoppers.
} else if (codon[1] == 0) {
if (codon[2] == 2 || codon[2] == 0) {
probability_adjustment += (*it).site_rate_away.at(3);
(*it).site_rate_away.at(3) = 0;
- cerr << "T0->0 ";
+ //cerr << "T0->0 ";
}
} else if (codon[1] == 2) {
if (codon[2] == 0) {
probability_adjustment += (*it).site_rate_away.at(3);
(*it).site_rate_away.at(3) = 0;
- cerr << "T0->0 ";
+ //cerr << "T0->0 ";
}
}
- cerr << endl;
+ //cerr << endl;
}
- cerr << "EXITING ZEROSTOPCODONS" << endl;
+ //cerr << "EXITING ZEROSTOPCODONS" << endl;
return probability_adjustment;
}
View
@@ -508,8 +508,6 @@ PathProposal::EvolveStep(
else i_z->evolvingSequence->Qidot = k_0->seq_evo.size();
num_diff = k_0->evolvingSequence->compare_sequence(i_z->evolvingSequence);
- cerr << "Number of differenced between sequences: " << num_diff << endl;
-
// For independent sites, branch needs to be scaled from substitutions per site to an approximation
// of the dependent number of substitutions per site. S is the value that we use.
if (!Pc && !nij) i_z->branch->S = i_z_rate_away / i_z->seq_evo.size();
@@ -529,8 +527,6 @@ PathProposal::EvolveStep(
}
}
- cerr << "Point-> PathProposal::EvolveStep() preprocessing done." << endl;
-
lambda_T = i_z->calculateEndpointRateAwayFromSequence(tree, k_0, T, t_0, -1);
begin_event = branch_terminal_event(k_0->mytipNo, BRANCH_BEGIN, k_0->DistanceFromRoot, k_0->bipartition, t_0, T);
begin_event->assign_Q(i_z, i_z->seq_evo.at(0), SUBSTITUTION, num_diff);
View
@@ -135,11 +135,7 @@ iQdPc(
/// affected by a change for subsequent calculation of the site-specific transition matrices.
//////////
i_z->set_site_window(tree->dep.front()->context.return_order(), &event_site, &end_site);
- cerr << "Point-> iQdPc() event_site " << event_site << " end_site " << end_site << endl;
i_z->site_specific_Qmat(tree, event_site, end_site);
- cerr << "Point-> iQdPc() returning rate matrix." << endl;
-
- exit(0);
return initial_rates;
}
View
@@ -726,7 +726,7 @@ TNode::set_site_window(
}
}
- cerr << "Point-> TNode::set_site_window(order, *start, *end) EXIT" << endl;
+ //cerr << "Point-> TNode::set_site_window(order, *start, *end) EXIT" << endl;
}
double
@@ -775,18 +775,18 @@ TNode::calculateForwardRateAwayFromSequence__order3Markov(
}
}
- if (Human_Data_simulation) {
- evolvingSequence->printSequenceRateAway();
- cerr << "Rate away from sequence: " << R_D << endl;
+ //if (Human_Data_simulation) {
+ //evolvingSequence->printSequenceRateAway();
+ //cerr << "Rate away from sequence: " << R_D << endl;
- int i = 0;
- for (vector<Site>::iterator seqi = seq_evo.begin(); seqi != seq_evo.end(); ++seqi, i++) {
- cerr << i << " " << stateCharacters.at((*seqi).returnState()) << " "
- << (*seqi).return_lookup_table_environment_index() << " "
- << (*seqi).return_lookup_table_sequence_index() << endl;
- }
+ //int i = 0;
+ //for (vector<Site>::iterator seqi = seq_evo.begin(); seqi != seq_evo.end(); ++seqi, i++) {
+ // cerr << i << " " << stateCharacters.at((*seqi).returnState()) << " "
+ // << (*seqi).return_lookup_table_environment_index() << " "
+ // << (*seqi).return_lookup_table_sequence_index() << endl;
+ //}
//exit(0);
- }
+ //}
//cerr << "Point-> TNode::calculateForwardRateAwayFromSequence__order3Markov(tree, event_site) RATE AWAY = " << R_D << endl;
return R_D;
@@ -1165,20 +1165,22 @@ TNode::calculateEndpointRateAwayFromSequence(
unsigned int end_site;
vector<Site>::iterator target_it = k_0->seq_evo.begin();
- //RateMatrix temp_rates = initialize_rate_matrices(tree, k_0, T, at_dt, event_site);
-//cerr << "Temp_rates being made: " << endl;
+ //cerr << "Temp_rates being made: " << endl;
+ // This makes temporary rates (set to global model of rates), but more importantly, allocates
+ // e_QijDt (from setRatesAway) and sets each position in the sequence to the appropriate, context
+ // dependent, e_QijDt matrix (ONLY 4x4, even when using codons; why? only one site changes at a time,
+ // and each nucleotide in the codon stores the proper rate away according to the entirety of the codon).
RateMatrix temp_rates = (*ptr2init)(tree, this, k_0, T, at_dt, event_site);
-//cerr << "Temp rates done." << endl;
-// cerr << "XXXXXXX" << temp_rates.Pij.at(0).at(0) << " " << temp_rates.Pij.at(0).at(1) << endl;
+ //cerr << "Temp rates done." << endl;
+ // cerr << "XXXXXXX" << temp_rates.Pij.at(0).at(0) << " " << temp_rates.Pij.at(0).at(1) << endl;
// Checking the cumulative sum with Pij versus Nij.
int seq_pos = 0;
for (vector<Site>::iterator it = seq_evo.begin(); it != seq_evo.end(); ++it, ++target_it, ++seq_pos) {
- //update_rate_matrices(&temp_rates, it, T, at_dt);
- (*ptr2update)(&temp_rates, this, it, T, at_dt); // Points to rate matrices in rate_type.cpp. //
-
-// cerr << "pos " << seq_pos << ": " << temp_rates.Pij.at(0).at(0) << " " << temp_rates.Pij.at(0).at(1) << endl;
-
+ // Sets rate matrix; exponentiates e_QijDt for the site to make the transition probabilities for the site.
+ // ptr2update points to rate matrix updates (uXxXx) in rate_type.cpp
+ (*ptr2update)(&temp_rates, this, it, T, at_dt);
+ //cerr << "pos " << seq_pos << ": " << temp_rates.Pij.at(0).at(0) << " " << temp_rates.Pij.at(0).at(1) << endl;
calculateEndpointRateAwayFromSite(it, target_it, &cumulative_site_sum_away, &forward_site_sum_away, &temp_rates);
}
@@ -1312,7 +1314,6 @@ TNode::site_specific_Qmat(
unsigned int end // End position
)
{
- cerr << "Point-> TNode::site_specific_Qmat(tree, start, end) ENTERED" << endl;
int seq_pos = 0;
int i,j;
vector<unsigned int> power_1 (numStates, 0);
@@ -1391,13 +1392,13 @@ TNode::site_specific_Qmat(
// ATGCCCAAGGGGGAT (0.268073 -0.832671 0.344024 0.220574 --> 0.268072779911734 0.268072779911734 0.612096796818472 0.832671112727411 )
// ATGGCCAAGGGGGAT (0.189035 0.174901 -0.516156 0.15222 --> 0.189035482893738 0.363936077856475 0.363936077856475 0.516156453888057 )
// ATGTCCAAGGGGGAT (0.301444 0.281933 0.382703 -0.96608 --> 0.301444108560112 0.583377210066327 0.966079793639185 0.966079793639185)
- cerr << "Position " << seq_pos << " Qmat: " << endl;
- int endline = 0;
- for (jt = (*it).e_QijDt.begin(); jt != (*it).e_QijDt.end(); ++jt, ++endline) {
- if (endline == 4) { cerr << endl; endline = 0; }
- cerr << (*jt) << " ";
- }
- cerr << endl;
+ //cerr << "Position " << seq_pos << " Qmat: " << endl;
+ //int endline = 0;
+ //for (jt = (*it).e_QijDt.begin(); jt != (*it).e_QijDt.end(); ++jt, ++endline) {
+ // if (endline == 4) { cerr << endl; endline = 0; }
+ // cerr << (*jt) << " ";
+ //}
+ //cerr << endl;
}
}

0 comments on commit 5f930de

Please sign in to comment.