denovodp.cpp

来自「MS-Clustering is designed to rapidly clu」· C++ 代码 · 共 1,254 行 · 第 1/3 页

		if (edge_idxs.size()==0)
		{
			cout << "Error: path with 0 edges!" << endl;
			exit(0);
		}

		if (debug_mode)
		{
			cout << "MP " << i << "\t";
			multi_paths[i].print(config);
		}
	}


	if (debug_mode)
	{

		exit(0);
	}
}



/*****************************************************************************
Same as the above function, but only returns paths that have good start/end idxs
******************************************************************************/
void  DeNovoDp::get_top_scoring_antisymetric_paths_with_specified_start_end_idxs(
		const vector<bool>& ind_allowed_start_end_idxs,
		vector<MultiPath>& multi_paths, 
		int required_num_paths, 
		int min_length, 
		int max_length, 
		score_t sym_penalty,
		score_t min_score_needed,
		bool try_complete_sequences,
		bool only_complete_sequences,
		double half_life_time) const
{
	const vector<Node>& nodes = prm->nodes;
	const vector<MultiEdge>& multi_edges = prm->multi_edges;
	const int num_nodes = nodes.size();
	
	int last_heap_pos = required_num_paths - 1;
	int last_alt_heap_pos = required_num_paths - 1;

	const score_t non_complete_penalty = 0; // this is a penalty added for each end of the peptide
											  // that does not reach the terminal. Used to bias 										     // against sequences that do not read ends.
							
	vector<edge_idx_set> heap,		// heap is used to store the paths. If running time takes too long,
									// we start reducing its size so the  prunning becomes more efficient
									// in such a case, we start putting removed paths into the alt_heap
									// so we don't lose good paths
						 alt_heap;
	
	vector< vector<score_t> > max_gains_for_length; // the maximal score attainable from each node
										            // using a given number of amino acids.
												    // length, node_idx

	vector< int > node_ordering; // holds the optimal order of nodes for the BB search

	vector<score_t> added_scores;  // holds for each depth in the tree, the score that was
								   // added by using the edges in the current path          

	vector<int> out_idx_counters; // for each node, what branch are we going down

	vector<bool> used_nodes;  // indicators for each node if it was used in the current path


	int i;
	if (max_length>31)
		max_length=31;
	
	added_scores.resize(num_nodes,NEG_INF);
	out_idx_counters.resize(num_nodes,0);
	used_nodes.resize(num_nodes,false);
	heap.clear();
	alt_heap.clear();
	heap.resize(required_num_paths);

	find_max_gains_per_length(max_length, max_gains_for_length);
	prm->sort_outgoing_edges_according_to_max_gains(max_gains_for_length);
	
	clock_t start_t,end_t;
	start_t = clock();

	if (only_complete_sequences)
	{
		node_ordering.clear();
		node_ordering.push_back(0);
	}
	else
		prm->get_node_ordering_according_to_max_gains(max_gains_for_length, node_ordering);

	int ns;
	for (ns=0; ns<node_ordering.size(); ns++)
	{
		const int start_idx = node_ordering[ns];

		if (! ind_allowed_start_end_idxs[start_idx])
			continue;

		if (max_gains_for_length[start_idx][max_length]<heap[0].score)
			break;

		const int num_first_out_edges = nodes[start_idx].out_edge_idxs.size();
		edge_idx_set current_path;
		int          current_node_idx;
	
		current_node_idx=start_idx;
		current_path.score = nodes[start_idx].score;
		used_nodes[start_idx] = true;

		if (try_complete_sequences && start_idx>0)
			current_path.score += non_complete_penalty;

		while (1)
		{
			
			// check if the search is running too long, if so decrease the heap size
			// and send the excess paths to the alternate heap
			end_t = clock();
			const double iteration_time = (end_t - start_t)/(double)CLOCKS_PER_SEC;
			if (0 && iteration_time>half_life_time)
			{
				if (alt_heap.size()==0)
					alt_heap.resize(heap.size());

				// reduce heap only if large enough
				if (heap.size()>7)
				{
					int half_size = heap.size()/2;
					while (heap.size()>half_size)
					{
						pop_heap(heap.begin(),heap.end());
						const edge_idx_set& removed_path = heap[heap.size()-1];
						if (removed_path.score>min_score_needed && removed_path.score>alt_heap[0].score)
						{
							pop_heap(alt_heap.begin(),alt_heap.end());
							alt_heap[last_alt_heap_pos] = removed_path;
							push_heap(alt_heap.begin(),alt_heap.end());
						}
						heap.pop_back();
					}
					last_heap_pos = heap.size()-1;
				//	cout << "Reduced heap to : " << heap.size() << endl;
				}
				start_t = end_t;
			}

			if (out_idx_counters[current_node_idx] >= nodes[current_node_idx].out_edge_idxs.size())
			{
				if (current_node_idx == start_idx)
					break; // we've returned all the way back, exhausted this tree
				// store path if necessary
				if (ind_allowed_start_end_idxs[current_node_idx] &&
					nodes[current_node_idx].out_edge_idxs.size() == 0 &&
					current_path.num_aa >= min_length &&
					current_path.num_aa <= max_length &&
					current_path.score > heap[0].score)
				{
					if (only_complete_sequences)
					{
						if (nodes[current_node_idx].type == NODE_C_TERM)
						{
							pop_heap(heap.begin(),heap.end());

							const edge_idx_set& removed_path = heap[last_heap_pos];
							if (alt_heap.size()>0 && 
								removed_path.score>min_score_needed &&
								removed_path.score>alt_heap[0].score)
							{
								pop_heap(alt_heap.begin(),alt_heap.end());
								alt_heap[last_alt_heap_pos] = removed_path;
								push_heap(alt_heap.begin(),alt_heap.end());
							}

							heap[last_heap_pos] = current_path;
							push_heap(heap.begin(),heap.end());
						}
					}
					else if (! try_complete_sequences || nodes[current_node_idx].type == NODE_C_TERM)
					{
						pop_heap(heap.begin(),heap.end());
						
						const edge_idx_set& removed_path = heap[last_heap_pos];
						if (removed_path.score>min_score_needed &&
							alt_heap.size()>0 &&
							removed_path.score>alt_heap[0].score)
						{
							pop_heap(alt_heap.begin(),alt_heap.end());
							alt_heap[last_alt_heap_pos] = removed_path;
							push_heap(alt_heap.begin(),alt_heap.end());
						}

						heap[last_heap_pos] = current_path;
						push_heap(heap.begin(),heap.end());
					}
					else
					{
						const score_t score_with_penalty = current_path.score + non_complete_penalty;
						if (score_with_penalty > min_score_needed &&
							score_with_penalty > heap[0].score)
						{
							current_path.score += non_complete_penalty;
							pop_heap(heap.begin(),heap.end());
							const edge_idx_set& removed_path = heap[last_heap_pos];
							if (alt_heap.size()>0 && removed_path.score>alt_heap[0].score)
							{
								pop_heap(alt_heap.begin(),alt_heap.end());
								alt_heap[last_alt_heap_pos] = removed_path;
								push_heap(alt_heap.begin(),alt_heap.end());
							}

							heap[last_heap_pos] = current_path;
							push_heap(heap.begin(),heap.end());
							current_path.score -= non_complete_penalty;
						}
					}
				}

				// backtrack
				out_idx_counters[current_node_idx] =0;
				used_nodes[current_node_idx]=false;
				current_path.length--;

				const int& path_length     = current_path.length;
				const MultiEdge& back_edge = multi_edges[current_path.edge_idxs[path_length]];
				current_path.num_aa		 -= back_edge.num_aa;
				current_node_idx          = back_edge.n_idx;
				current_path.score       -= added_scores[path_length];
				continue;
			}


			// discard this path if we are using too many edges or 
			// the score will not be able to improve enough
			// since the edges are sorted according to the gain they can bring,
			// none of the rest can help so we skip the rest

			const int remaining_aas = max_length - current_path.num_aa;
			const score_t threshold_score = (min_score_needed > heap[0].score ?  min_score_needed : heap[0].score);
			const score_t maximal_achievable_score = (remaining_aas>=0 ? current_path.score + max_gains_for_length[current_node_idx][remaining_aas] : NEG_INF);
			if (current_path.num_aa > max_length || maximal_achievable_score<threshold_score)
			{
				out_idx_counters[current_node_idx] = nodes[current_node_idx].out_edge_idxs.size();
				continue; 
			}
		
			// advance on the edge
			const int edge_idx = nodes[current_node_idx].out_edge_idxs[out_idx_counters[current_node_idx]];
			const MultiEdge& e = multi_edges[edge_idx];

			out_idx_counters[current_node_idx]++;
			current_node_idx = e.c_idx;
			out_idx_counters[current_node_idx] =0;
			used_nodes[current_node_idx]=true;
			
			added_scores[current_path.length] = e.max_variant_score + nodes[e.c_idx].score;

			// check if forbidden pair is used..
			if (forbidden_idxs[current_node_idx]>=0 && used_nodes[forbidden_idxs[current_node_idx]])
				added_scores[current_path.length] -= sym_penalty; 

			current_path.edge_idxs[current_path.length] = edge_idx;
			current_path.num_aa += multi_edges[edge_idx].num_aa;
			current_path.score += added_scores[current_path.length];
			current_path.length++;

	
			// 
			score_t heap_score = heap[0].score;
			score_t added_score = added_scores[current_path.length-1];

			// check if the path should be stored at this stage, and if we can mark
			if (ind_allowed_start_end_idxs[current_node_idx] &&
				added_scores[current_path.length-1]>-5 && 
				nodes[current_node_idx].out_edge_idxs.size()>0 &&
				current_path.num_aa <= max_length &&
				current_path.num_aa >= min_length &&
				current_path.score> heap[0].score )
			{
				if (only_complete_sequences)
				{
					if (nodes[current_node_idx].type == NODE_C_TERM)
					{
						pop_heap(heap.begin(),heap.end());
						const edge_idx_set& removed_path = heap[last_heap_pos];
						if (alt_heap.size()>0 &&
							removed_path.score > min_score_needed &&
							removed_path.score > alt_heap[0].score)
						{
							pop_heap(alt_heap.begin(),alt_heap.end());
							alt_heap[last_alt_heap_pos] = removed_path;
							push_heap(alt_heap.begin(),alt_heap.end());
						}

						heap[last_heap_pos] = current_path;
						push_heap(heap.begin(),heap.end());
					}
				}
				else if (! try_complete_sequences || nodes[current_node_idx].type == NODE_C_TERM)
				{
					pop_heap(heap.begin(),heap.end());
					const edge_idx_set& removed_path = heap[last_heap_pos];
					if (alt_heap.size()>0 && 
						removed_path.score > min_score_needed &&
						removed_path.score>alt_heap[0].score)
					{
						pop_heap(alt_heap.begin(),alt_heap.end());
						alt_heap[last_alt_heap_pos] = removed_path;
						push_heap(alt_heap.begin(),alt_heap.end());
					}

					heap[last_heap_pos] = current_path;
				}
				else
				{

					if (current_path.score + non_complete_penalty > heap[0].score)
					{
						current_path.score += non_complete_penalty;
						pop_heap(heap.begin(),heap.end());
						const edge_idx_set& removed_path = heap[last_heap_pos];
						if (alt_heap.size()>0 &&
							removed_path.score> min_score_needed &&
							removed_path.score>alt_heap[0].score)
						{
							pop_heap(alt_heap.begin(),alt_heap.end());
							alt_heap[last_alt_heap_pos] = removed_path;
							push_heap(alt_heap.begin(),alt_heap.end());
						}

						heap[last_heap_pos] = current_path;
						push_heap(heap.begin(),heap.end());
						current_path.score -= non_complete_penalty;
					}
				}
			}
		}

		used_nodes[start_idx] = false;
	}


	if (heap.size()==0)
		return;

	if (alt_heap.size()>0) 	// transfer all paths that are in current heap
	{
		int i;
		for (i=0; i<heap.size(); i++)
			if (heap[i].score>alt_heap[0].score)
			{
				pop_heap(alt_heap.begin(),alt_heap.end());
				alt_heap[last_alt_heap_pos]=heap[i];
				push_heap(alt_heap.begin(),alt_heap.end());
			}
	}

	
	// work on the alt_heap if necessary
	vector<edge_idx_set>& final_heap = (alt_heap.size()>0 ? alt_heap : heap);
	
	sort(final_heap.begin(),final_heap.end());

	while (final_heap.size()>0 && final_heap[final_heap.size()-1].score<-40)
		final_heap.pop_back();

	int actual_num_paths = required_num_paths;
	if (final_heap.size()<actual_num_paths)
		actual_num_paths = final_heap.size();

	multi_paths.resize(actual_num_paths);

	for (i=0; i<actual_num_paths; i++)
	{
		int j;
		vector<int> edge_idxs;
		edge_idxs.resize(final_heap[i].length);
		for (j=0; j<final_heap[i].length; j++)
			edge_idxs[j]=final_heap[i].edge_idxs[j];

		prm->create_path_from_edges(edge_idxs, multi_paths[i]);

		multi_paths[i].path_score = final_heap[i].score;
		multi_paths[i].original_rank = i;

		const vector<int>& node_idxs = multi_paths[i].node_idxs;
		const int max_idx = node_idxs.size()-1;
		int num_forbidden =0;
		for (j=0; j<max_idx; j++)
		{
			const int n_idx = node_idxs[j];
			if (forbidden_idxs[n_idx]<0)
				continue;
			int k;
			for (k=j+1; k<node_idxs.size(); k++)
				if (node_idxs[k]==forbidden_idxs[n_idx])
					break;
			if (k<node_idxs.size())
				num_forbidden++;
		}
		multi_paths[i].num_forbidden_nodes = num_forbidden;
	}
}