C++ IntervalSet类代码示例

IntervalSet IntervalSet::Or(const std::vector<IntervalSet> &sets) {
  IntervalSet result;
  for (auto &s : sets) {
    result.addAll(s);
  }
  return result;
}

TEST(IntervalTest, IntervalIteration)
{
  IntervalSet<int> set;

  set += (Bound<int>::closed(0), Bound<int>::closed(1));
  set += (Bound<int>::open(2), Bound<int>::open(4));
  set += (Bound<int>::closed(5), Bound<int>::open(7));
  set += (Bound<int>::open(7), Bound<int>::closed(9));

  EXPECT_EQ(4u, set.intervalCount());

  int index = 0;

  foreach (const Interval<int>& interval, set) {
    if (index == 0) {
      EXPECT_EQ(0, interval.lower());
      EXPECT_EQ(2, interval.upper());
    } else if (index == 1) {
      EXPECT_EQ(3, interval.lower());
      EXPECT_EQ(4, interval.upper());
    } else if (index == 2) {
      EXPECT_EQ(5, interval.lower());
      EXPECT_EQ(7, interval.upper());
    } else if (index == 3) {
      EXPECT_EQ(8, interval.lower());
      EXPECT_EQ(10, interval.upper());
    }
    index++;
  }
}

    /**
     * @brief clip a ray by a Cone.
     * @param ray_ the ray to clip
     * @param bounds ray sorted bounds of the resulting clipping segment.
     * @return  true if ray was clipped. bounds parameter is filld by this method
     *
     * For simple convex objects, there is two values in  bounds that represent in and out events.
     * An in event is whe the ray enters the geometry, an out is when the ray leaves the geometry.
     *
     * Adapted from http://www.geometrictools.com/LibMathematics/Intersection/Wm5IntrLine3Cone3.cpp
     *
     */
    bool clip(Ray& ray_, IntervalSet &bounds) const {
        Diff_Geom ipoints[2];
        float t[2];

        int numintersection = get_clip_points(ray_, ipoints,  t);
        if (numintersection == 0)
            return false;
        else if (numintersection == 2) {
            Diff_Geom *in = new Diff_Geom(ipoints[0].pos(), ipoints[0].normal(), t[0]);
            Diff_Geom *out = new Diff_Geom(ipoints[1].pos(), ipoints[1].normal(), t[1]);
            bounds.add(in, out);
            return true;
        } else {
            Diff_Geom *in, *out;
            IntervalSet capset;
            // check with cap plane
            Plane cap(m_vertex + m_axis * m_height, m_axis);
            if (cap.clip(ray_, capset) ) {
				/* Beuark mais ca facilite la gestion mémoire ... */
                Diff_Geom *tmp =capset.bounds()[0].data;
				in = new Diff_Geom(*tmp);
                if (in->t()< t[0]) {
                    out = new Diff_Geom(ipoints[0].pos(), ipoints[0].normal(), t[0]);
                } else {
                    out = in;
                    in = new Diff_Geom(ipoints[0].pos(), ipoints[0].normal(), t[0]);
                }

                bounds.add(in, out);
                return true;
            } else {
                return false;
            }
        }
    }

XHMM::xhmmInputManager::IntervalSet* XHMM::xhmmInputManager::readIntervalsFromFile(string intervalsFile, const set<char>& excludeLinesStartingWith) {
	HMM_PP::istreamLineReader* intervalsStream = HMM_PP::utils::getIstreamLineReaderFromFile(intervalsFile);
	if (intervalsStream == NULL)
		throw new Exception("Unable to read table from file '" + intervalsFile + "'");

	IntervalSet* intervSet = new IntervalSet();
	while (!intervalsStream->eof()) {
		string* line = new string();
		*intervalsStream >> *line;
		if (line->empty()) {
			delete line;
			continue;
		}
		char firstChar = (*line)[0];
		if (firstChar != NO_EXCLUDE_CHAR && excludeLinesStartingWith.find(firstChar) != excludeLinesStartingWith.end()) {
			delete line;
			continue;
		}
		stringstream* lineStream = new stringstream(*line);
		delete line;

		intervSet->insert(Interval(*lineStream));
		delete lineStream;
	}
	delete intervalsStream;

	return intervSet;
}

// Initialize intervals, one for each unique attribute value
void initialize_intervals(int dimIndex)
{
	g_intervals.clear();

	IntervalSet interval;

	TupleVec::iterator it = g_data.begin(), end = g_data.end(), next;
	int index = 0;
	for ( ; it != end; ++it, ++index) {
		next = it+1;

		// Add element (index) to interval
		interval.insert(index);

		bool insertInterval = (next == end) ||
			( (next->first)[dimIndex] != (it->first)[dimIndex] );

		if (insertInterval) {
			// Insert interval into list if
			//  (a) end of sequence or
			//  (b) next element is not the same as current element
			g_intervals.push_back(interval);
			interval.clear();
		}
	}

	// Debug
	print_all_intervals();
}

void print_set(const IntervalSet &s)
{
    for(IntervalSet::const_iterator i=s.begin(); i!=s.end(); ++i)
    {
        printf("%i-%i ", (int)i->start, (int)i->end);
    }
    printf("\n");
}

typename IntervalSet<DomainT,Interval,Compare,Alloc>::interval_type 
enclosure(const IntervalSet<DomainT,Interval,Compare,Alloc>& object)
{
    typedef IntervalSet<DomainT,Interval,Compare,Alloc> IntervalSetT;
    typedef typename IntervalSetT::interval_type interval_type;
    return 
        object.empty() ? neutron<interval_type>::value()
        : (*object.begin()).span(*object.rbegin());
}

XHMM::xhmmInputManager::LoadedReadDepths
XHMM::xhmmInputManager::loadReadDepthsFromFile(string readDepthFile, IntervalSet* excludeTargets, StringSet* excludeTargetChromosomes, StringSet* excludeSamples,
		const ullint minTargetSize, const ullint maxTargetSize) {

	HMM_PP::DoubleMat* rdMat = MatrixReader<double>::readMatrixFromFile(readDepthFile);
	StringSet* excludedSamples = new StringSet();
	IntervalSet* excludedTargets = new IntervalSet();

	set<ullint>* excludeSampleIndices = new set<ullint>();
	if (excludeSamples != NULL) {
		for (ullint row = 0; row < rdMat->nrow(); ++row) {
			string samp = rdMat->rowName(row);

			if (excludeSamples->find(samp) != excludeSamples->end()) {
				cerr << "Excluded sample " << samp << endl;

				excludedSamples->insert(samp);
				excludeSampleIndices->insert(row);
			}
		}
	}

	set<ullint>* excludeTargetIndices = new set<ullint>();
	if (excludeTargets != NULL || excludeTargetChromosomes != NULL || minTargetSize > 0 || maxTargetSize < ULLINT_INFINITY) {
		for (ullint j = 0; j < rdMat->ncol(); ++j) {
			const Interval targJ(rdMat->colName(j));
			const ullint targLen = targJ.span();
			bool targLenFails = (targLen < minTargetSize || targLen > maxTargetSize);

			if ((excludeTargets != NULL && excludeTargets->find(targJ) != excludeTargets->end()) ||
					(excludeTargetChromosomes != NULL && excludeTargetChromosomes->find(targJ.getChr()) != excludeTargetChromosomes->end()) ||
					targLenFails) {
				cerr << "Excluded target " << targJ;
				if (targLenFails)
					cerr << " of length " << targLen;
				cerr << endl;

				excludeTargetIndices->insert(j);
				excludedTargets->insert(targJ);
			}
		}
	}

	if (!excludeSampleIndices->empty() || !excludeTargetIndices->empty()) {
		HMM_PP::DoubleMat* newRdMat = rdMat->deleteRowsAndColumns(excludeSampleIndices, excludeTargetIndices);
		delete rdMat;
		rdMat = newRdMat;
	}
	delete excludeSampleIndices;
	delete excludeTargetIndices;

	return LoadedReadDepths(rdMat, excludedTargets, excludedSamples);
}

pair<DataType, DataType> get_range(IntervalSet& interval, int dimIndex)
{
	assert(!interval.empty());

	int beginIndex = *interval.begin();
	IntervalSet::iterator end = interval.end();
	--end;
	int endIndex = *end;

	assert(beginIndex >= 0 && beginIndex < (int)g_data.size());
	assert(endIndex >= 0 && endIndex < (int)g_data.size());

	return make_pair(
		g_data[beginIndex].first[dimIndex], 
		g_data[endIndex].first[dimIndex]);
}

bool is_disjoint
(
          interval_base_set<SubType,DomainT,Interval,Compare,Alloc>& object,
    const IntervalSet              <DomainT,Interval,Compare,Alloc>& operand
)
{
    typedef interval_base_set<SubType,DomainT,Interval,Compare,Alloc> object_type;
    typedef IntervalSet              <DomainT,Interval,Compare,Alloc> operand_type;
    object_type intersection;

    if(operand.empty())
        return true;

    typename operand_type::const_iterator common_lwb;
    typename operand_type::const_iterator common_upb;

    if(!Set::common_range(common_lwb, common_upb, operand, object))
        return true;

    typename operand_type::const_iterator it = common_lwb;
    while(it != common_upb)
    {
        object.add_intersection(intersection, *it++);
        if(!intersection.empty())
            return false;
    }

    return true; 
}

IntervalSet IntervalSet::And(const IntervalSet &other) const {
  IntervalSet intersection;
  size_t i = 0;
  size_t j = 0;

  // iterate down both interval lists looking for nondisjoint intervals
  while (i < _intervals.size() && j < other._intervals.size()) {
    Interval mine = _intervals[i];
    Interval theirs = other._intervals[j];

    if (mine.startsBeforeDisjoint(theirs)) {
      // move this iterator looking for interval that might overlap
      i++;
    } else if (theirs.startsBeforeDisjoint(mine)) {
      // move other iterator looking for interval that might overlap
      j++;
    } else if (mine.properlyContains(theirs)) {
      // overlap, add intersection, get next theirs
      intersection.add(mine.intersection(theirs));
      j++;
    } else if (theirs.properlyContains(mine)) {
      // overlap, add intersection, get next mine
      intersection.add(mine.intersection(theirs));
      i++;
    } else if (!mine.disjoint(theirs)) {
      // overlap, add intersection
      intersection.add(mine.intersection(theirs));

      // Move the iterator of lower range [a..b], but not
      // the upper range as it may contain elements that will collide
      // with the next iterator. So, if mine=[0..115] and
      // theirs=[115..200], then intersection is 115 and move mine
      // but not theirs as theirs may collide with the next range
      // in thisIter.
      // move both iterators to next ranges
      if (mine.startsAfterNonDisjoint(theirs)) {
        j++;
      } else if (theirs.startsAfterNonDisjoint(mine)) {
        i++;
      }
    }
  }

  return intersection;
}

// Compute $\chi^2$ value for an interval
float chisquared_interval(IntervalSet& interval)
{
	// Count instances of each class
	map<ClassType, int> classCount;

	IntervalSet::iterator it = interval.begin(),
		end = interval.end();	
	for ( ; it != end; ++it) {
		int index = *it;
        ClassType classType = g_data[index].second;
		classCount[classType]++;
	}

	// Keep track of summation
	float chichiri = 0.0f;

	set<ClassType>::iterator sit = g_classTypes.begin(),
		send = g_classTypes.end();

	for ( ; sit != send; ++sit) {
		ClassType classType = *sit;

		// Uses notation from section 4.2.1,
		// "Discretization: An Enabling Technique" by H. Liu et al.

		float A_ij = static_cast<float>(classCount[classType]);

		float R_i  = static_cast<float>(interval.size());
		float C_j  = static_cast<float>(g_globalClassCount[classType]);
		float N    = static_cast<float>(g_data.size());

		float E_ij = (R_i * C_j) / N;

		float top = (A_ij - E_ij);

		chichiri += top*top / E_ij;
	}

	cout << "{" << chichiri << "} + ";

	return chichiri;
}

int main()
{
    IntervalSet s;
    s.insert(Interval<size_t>(1,10));
    s.insert(Interval<size_t>(10,20));
    s.insert(Interval<size_t>(20,30));
    s.insert(Interval<size_t>(30,40));

    print_set(s);

    printf("Matching intervals:\n");
    std::pair<IntervalSet::const_iterator,IntervalSet::const_iterator> ir = intersecting_intervals(s, Interval<size_t>(10,25));
    for(IntervalSet::const_iterator i=ir.first; i!=ir.second; ++i)
    {
        printf("%i-%i ", (int)i->start, (int)i->end);
    }
    printf("\n");

    intervalset_merge(s, Interval<size_t>(25, 50));
    print_set(s);
}

    bool clip(Ray& ray_, IntervalSet &bounds) const {
        Diff_Geom ipoints[2];
        float t[2];

        int numintersection = get_clip_points(ray_, ipoints,  t);
        if (numintersection == 0)
            return false;
        else if (numintersection == 2) {
            Diff_Geom *in = new Diff_Geom(ipoints[0].pos(), ipoints[0].normal(), ipoints[0].dNdx(), ipoints[0].dNdy(), t[0], ipoints[0].dPdx(), ipoints[0].dPdy(), ipoints[0].dDdx(), ipoints[0].dDdy(), ipoints[0].u(), ipoints[0].v(), ipoints[0].dTdx(), ipoints[0].dTdy(), true);
            Diff_Geom *out = new Diff_Geom(ipoints[1].pos(), ipoints[1].normal(), ipoints[1].dNdx(), ipoints[1].dNdy(), t[1], ipoints[1].dPdx(), ipoints[1].dPdy(), ipoints[1].dDdx(), ipoints[1].dDdy(), ipoints[1].u(), ipoints[1].v(), ipoints[1].dTdx(), ipoints[1].dTdy(), false);
            bounds.add(in, out);
            return true;
        } else {
            Diff_Geom *in, *out;
            IntervalSet capset;
            // check with cap plane
            Plane cap(m_vertex + m_axis * m_height, m_axis);
            if (cap.clip(ray_, capset) ) {
                in = capset.bounds()[0].data;
                in->set_u(in->u()/(PSCALE*m_radius));
                in->set_v(in->v()/(PSCALE*m_radius));
                in->set_dTdx(in->dTdx()* (1.f/(PSCALE*m_radius)));
                in->set_dTdy(in->dTdy()* (1.f/(PSCALE*m_radius)));
                if (in->t()< t[0]) {
                //if (p_diff_geom_in_->t() < t[0]) {
                    out = new Diff_Geom(ipoints[0].pos(), ipoints[0].normal(), ipoints[0].dNdx(), ipoints[0].dNdy(), t[0], ipoints[0].dPdx(), ipoints[0].dPdy(), ipoints[0].dDdx(), ipoints[0].dDdy(), ipoints[0].u(), ipoints[0].v(), ipoints[0].dTdx(), ipoints[0].dTdy(), false);
                } else {
                    out = in;
                    out->set_in(false);
                    in = new Diff_Geom(ipoints[0].pos(), ipoints[0].normal(), ipoints[0].dNdx(), ipoints[0].dNdy(), t[0], ipoints[0].dPdx(), ipoints[0].dPdy(), ipoints[0].dDdx(), ipoints[0].dDdy(), ipoints[0].u(), ipoints[0].v(), ipoints[0].dTdx(), ipoints[0].dTdy(), true);
                }
                delete (capset.bounds()[1].data);
                bounds.add(in, out);
                return true;
            } else {
                return false;
            }
        }
    }

int main()
{
    using namespace ch10;
    using namespace std;

    IntervalSet intervals;
    intervals.insert(1);
    intervals.insert(100);
    intervals.insert(200);
    intervals.insert(300);
    intervals.insert(400);
    intervals.insert(500);
    intervals.insert(600);
    intervals.insert(700);
    intervals.insert(800);
    intervals.insert(900);
    intervals.insert(1000);

    Histogram hist(intervals);

    cout << "Uniform distribution:" << endl;
    UniformRandomGenerator gen1(1, 1000);
    for(int i = 0; i < 100; i++)
    {
        hist.add(gen1.draw());
    }
    hist.print_formatted(cout);

    cout << "Exponential distribution:" << endl;
    ExponentialRandomGenerator gen2(1,200);
    hist = Histogram(intervals);
    for(int i = 0; i < 100; i++)
    {
        int result = gen2.draw();
        if(result >= 1 && result < 1000) hist.add(result);
    }
    hist.print_formatted(cout);

    cout << "Normal distribution:" << endl;
    NormalRandomGenerator gen3(500, 200);
    hist = Histogram(intervals);
    for(int i = 0; i < 100; i++)
    {
        int result = gen3.draw();
        if(result >= 1 && result < 1000) hist.add(result);
    }
    hist.print_formatted(cout);
}