Mercurial > repos > youngkim > ezbamqc
view ezBAMQC/src/ezBAMQC/IntervalTree.cpp @ 19:caed0168f704
wrapper
author | youngkim |
---|---|
date | Wed, 30 Mar 2016 13:46:41 -0400 |
parents | dfa3745e5fd8 |
children |
line wrap: on
line source
// // IntervalTree.cpp // BAMQC-0.5 // // Created by Ying Jin on 9/15/15. // Copyright (c) 2015 Ying Jin. All rights reserved. // #include "IntervalTree.h" #include <iostream> Interval::Interval(int g, int e, int st, int end, int t) { gene = g; start = st; stop = end; type = t; exon = e; } Interval::Interval(){} Interval::~Interval() { } IntervalTree::IntervalTree(){} IntervalTree::IntervalTree(std::vector<Interval> intervals, int depth, int minbucket, int extent_st , int extent_end, int maxbucket){ std::vector<Interval> lefts ; std::vector<Interval> rights ; int left_ext = 0; int right_ext = 0; int max_stop = 0; float center_pos = 0.0 ; Interval itv ; size_t i ; left = nullptr; right = nullptr; depth -= 1; if ((depth == 0 || intervals.size() < (size_t)minbucket) && intervals.size() < (size_t)maxbucket){ itvlist = intervals; left = right = nullptr; center = -1; return ; } //if (extent_st == -1){ // sorting the first time through allows it to get // better performance in searching later. // quick_sort(intervals,0,intervals.size()-1); //} for(i=0;i<intervals.size(); i++){ itv = intervals[i]; if (itv.stop > max_stop) { max_stop = itv.stop; } } if(extent_st != -1){ left_ext = extent_st; right_ext = extent_end; } else { left_ext = intervals[0].start; right_ext = max_stop; } //left, right = _extent or (intervals[0].start, max_stop) //center = intervals[len(intervals)/ 2].stop center_pos = (left_ext + right_ext) / 2.0; //std::cout << "left_ext " << left_ext << std::endl; //std::cout << "right_ext " << right_ext << std::endl; //std::cout << "center " << center_pos << std::endl; for(i=0;i<intervals.size();i++){ itv = intervals[i]; if (itv.stop < center_pos) { lefts.push_back(itv); } else { if (itv.start > center_pos) { rights.push_back(itv);} else { itvlist.push_back(itv) ;} } } if (lefts.empty()) { left = nullptr; } else{ left = new IntervalTree(lefts,depth,minbucket,intervals[0].start,center_pos,maxbucket); } if (rights.empty()) { right = nullptr; } else { right = new IntervalTree(rights, depth, minbucket, center_pos,right_ext,maxbucket); } center = center_pos; } void IntervalTree::clear(IntervalTree * node){ // std::cout << "clear node" << std::endl; //if (node != nullptr) { // std::cout << node->center << std::endl; //} if (node != nullptr) { clear(node->left); clear(node->right); // delete node; } } IntervalTree::~IntervalTree() { clear(left); clear(right); /* if (left != nullptr){ if( left->left == nullptr && left->right == nullptr) { delete left; } else { //if (left->left != nullptr) { (*left).~IntervalTree(); // } //if (left->right != nullptr) { // (*left->right).~IntervalTree(); //} } } if (right != nullptr ) { if(right->left == nullptr && right->right == nullptr) { delete right; } else{ // if (right->left != nullptr) { // (*right->left).~IntervalTree(); //} //if (right->right != nullptr) { // (*right->right).~IntervalTree(); //} (*right).~IntervalTree(); } } delete this;*/ } std::vector<Interval> IntervalTree::find(int start, int stop){ //"""find all elements between (or overlapping) start and stop""" std::vector<Interval> overlapping ; std::vector<Interval> temp ; size_t i; if (!itvlist.empty() && stop >= itvlist[0].start) { for (i=0;i< itvlist.size(); i++) { if (itvlist[i].stop >= start && itvlist[i].start <= stop) { overlapping.push_back(itvlist[i]); } } } if( left != nullptr && start <= center) { temp = (*left).find(start, stop); overlapping.insert(overlapping.end(), temp.begin(),temp.end()); } if(right != nullptr && stop >= center){ temp = (*right).find(start, stop); overlapping.insert(overlapping.end(),temp.begin(),temp.end()); } return overlapping ; } std::vector<int> IntervalTree::find_gene(int start, int stop){ std::vector<int> overlapping ; std::vector<int> temp ; size_t i; if (!itvlist.empty() && stop >= itvlist[0].start) { for (i=0;i< itvlist.size(); i++) { if (itvlist[i].stop >= start && itvlist[i].start <= stop) { overlapping.push_back(itvlist[i].gene); } } } if( left != nullptr && start <= center) { temp = (*left).find_gene(start, stop); overlapping.insert(overlapping.end(), temp.begin(),temp.end()); } if(right != nullptr && stop >= center){ temp = (*right).find_gene(start, stop); overlapping.insert(overlapping.end(),temp.begin(),temp.end()); } return overlapping ; } /* int main(){ std::vector<Interval> ll; IntervalTree *Idx; std::vector<Interval> res; std::vector<std::string> res_gene; int i; ll.push_back(Interval("gid1",100,300,"cds")); ll.push_back(Interval("gid2",10,100,"cds")); ll.push_back(Interval("gid3",1000,3000,"cds")); ll.push_back(Interval("gid4",500,800,"cds")); ll.push_back(Interval("gid5",200,700,"cds")); std::cout << "before create idx" << std::endl; Idx = new IntervalTree(ll,16, 2, -1, -1, 3); std::cout << "after create idx" << std::endl; res = (*Idx).find(20,600); res_gene = (*Idx).find_gene(20,600); for (i = 0; i<res.size(); i++) { std::cout << res[i].gene << std::endl; } for (i = 0; i<res_gene.size(); i++) { std::cout << res[i].gene << std::endl; } } */