realset.cpp

非常著名的曲线拟合程序

{
  double* d = data_.data();
  const unsigned N = size();
  for(unsigned i=0; i<N; ++i)
    d[i] = fabs(d[i]-y);
  recalc();
}

void
RealSet::reverse()
{
  if (sortedness_ == 2) return;

  double* d = data_.data();
  const unsigned N = size();

  for(unsigned i=0; i<N/2; ++i) {
    double t = d[i];
    d[i] = d[N-i-1];
    d[N-1-i] = t;
  }
  sortedness_ = -sortedness_;
}

void
RealSet::sort(int dir)
{
  // If we are already sorted properly, ignore
  if (sortedness_ == 2 || 
      (sortedness_ == 1 && dir >= 0) ||
      (sortedness_ == -1 && dir < 0))
    return;

  // If we have sorted data already, use it!
  if (dir >= 0 && sdata_.size()) {
    data_ = sdata_;
    sortedness_ = 1;
    return;
  }

  // If we are being asked to "opposite-sort" already-sorted data, do so
  // via reverse()
  if ((dir >= 0 && sortedness_ == -1) || (dir < 0 && sortedness_ == 1)) {
    reverse();
  } else if (dir >= 0) {
    // Otherwise there is no avoiding it, so just sort the data...
    double* d = data_.data();
    std::sort(d, d+size(), less<double>());
    sortedness_ = 1;
  } else {
    double* d = data_.data();
    std::sort(d, d+size(), greater<double>());
    sortedness_ = -1;
  }

  if (sortedness_ > 0)
    sdata_ = data_;
}

void
RealSet::scramble_explicit(Random& rnd)
{
  if (size() == 0 || sortedness_ == 2) return;

  rnd.scramble(data_.data(), size());

  recalc_sortedness();
}

void
RealSet::permute(const Permutation& p)
{
  if (size() != p.size())
    throw Exception("Permutation and data are not of the same size");
  
  double* d = data_.data();
  p.permute(d);

  recalc_sortedness();
}

void
RealSet::rank()
{
  pair<double,int>* rankdata = new pair<double,int>[size()];
  double* d = data_.data();
  const unsigned N = size();

  for(unsigned i1=0; i1<N; ++i1) {
    rankdata[i1].first = d[i1];
    rankdata[i1].second = i1;
  }

  std::sort(rankdata, rankdata+N, less< pair<double,int> >());

  // set ranks, handling ties by the usual convention
  unsigned count=1;
  for(unsigned i=0; i<N; /* increment down below */ ) {
    unsigned j=i+1;
    unsigned count2=0;
    while (j < N && are_equal(rankdata[i].first, rankdata[j].first)) {
      ++j;
      ++count2;
    }
    while (i != j)
      rankdata[i++].first = (2*count + count2)/2.0; // we increment i here
    count += 1+count2;
  }

  for(unsigned i=0; i<N; ++i)
    d[rankdata[i].second] = rankdata[i].first;

  recalc(); // This could be optimized...

  delete [] rankdata;
}

void
RealSet::standardize()
{
  const double m = mean();
  const double s = sdev();
  const unsigned N = size();
  double* d = data_.data();
  for(unsigned i=0; i<N; ++i)
    d[i] = (d[i]-m)/s;
  if (have_sorted_data()) {
    double* sd = sdata_.data();
    for(unsigned i=0; i<N; ++i)
      sd[i] = (sd[i]-m)/s;
  }
  recalc_saving_sorted();
}

Permutation
RealSet::sorting_permutation(int dir) const
{
  Permutation p;
  const unsigned N = size();

  if (sortedness_ == 2 ||
      (dir >= 0 && sortedness_ == 1) ||
      (dir < 0 && sortedness_ == -1)) {

    p.set_identity(N);

  } else if ((dir > 0 && sortedness_ == -1) ||
	     (dir < 0 && sortedness_ == 1)) {

    p.set_reverse(N);

  } else {

    pair<double,unsigned>* tmp = new pair<double,unsigned>[N];
    unsigned* tmp2 = new unsigned[N];
    const double* d = data_.data();
    for (unsigned i=0; i<N; ++i) {
      tmp[i].first = d[i];
      tmp[i].second = i;
    }

    if (dir >= 0) 
      std::sort(tmp, tmp+N, less< pair<double,unsigned> >());
    else
      std::sort(tmp, tmp+N, greater< pair<double,unsigned> >());

    for(unsigned i=0; i<N; ++i)
      tmp2[tmp[i].second] = i;

    p.initialize(size(), tmp2);
    delete [] tmp;
    delete [] tmp2;
  }

  return p;
}

void
RealSet::add_compressed(const string& s)
{
  const char* buffer = s.c_str();
  const char* eob = buffer + s.size();
  while (buffer < eob) 
    add(decode_double(buffer));
}

void
RealSet::write_xml(ostream& out, bool compress) const
{
  out << "<realset label=\"" << label() << "\" size=" << size();
  if (compress)
    out << " compressed";
  out << ">" << endl;

  const double* d = data_.data();
  const unsigned N = size();

  if (compress) {
    // This is bigger than necessary.
    char buffer[16];
    memset(buffer,0,16);

    for(unsigned i=0; i<N; ++i) {
      encode_double(d[i], buffer);
      out << buffer;
      if (i % 7 == 6)
	out << endl;
    }
    if (N % 7)
      out << endl;
  } else {
    int old_prec = out.precision();
    out.precision(17);
    for(unsigned i=0; i<N; ++i)
      out << d[i] << endl;
    out.precision(old_prec);
  }

  out << "</realset>" << endl;
}

void
RealSet::update_add(double x)
{
  // update statistics /////

  const unsigned N = size();

  if (N == 0 || x < min_)
    min_ = x;
  if (N == 0 || x > max_)
    max_ = x;

  double old_mean = mean_;
  mean_ += (x - mean_)/(N+1);
  if (N > 0) {
    sumsq_ += (x-mean_)*(x-old_mean); // cool, huh?

    // Update sortedness flag
    if (sortedness_) {
      const double* d = data_.data_const();
      double prev = d[N-1];
      if (sortedness_ == 2) { // constant
	if (prev < x)
	  sortedness_ = 1;
	else if (prev > x)
	  sortedness_ = -1;
      } else if (sortedness_ == 1) { // ascending
	if (prev > x)
	  sortedness_ = 0;
      } else if (sortedness_ == -1) { // descending
	if (prev < x)
	  sortedness_ = 0;
      }
    }
  }
}

void
RealSet::update_remove(double x)
{
  const unsigned N = size()-1;
  double old_mean = mean_;
  mean_ += (mean_ - x)/N;
  if (N > 0)
    sumsq_ -= (x-old_mean)*(x-mean_);

  if (N != 0) {
    const double* d = data_.data_const();
    const unsigned NN = size();
    if (min_ == x) {
      bool skip_x = true;
      min_ = d[0];
      if (x == d[0]) {
	min_ = d[1];
	skip_x = false;
      }
      for(unsigned i=1; i<NN; ++i) {
	if (d[i] == x && skip_x) {
	  skip_x = false;
	  continue;
	}
	if (min_ > d[i])
	  min_ = d[i];
      }
    }
    if (max_ == x) {
      bool skip_x = true;
      max_ = d[0];
      if (x == d[0]) {
	max_ = d[1];
	skip_x = false;
      }
      for(unsigned i=1; i<NN; ++i) {
	if (d[i] == x && skip_x) {
	  skip_x = false;
	  continue;
	}
	if (max_ < d[i])
	  max_ = d[i];
      }
    }
  }
  // Leaves sortedness in fubared state.
}

// There is some code duplication between these and add(), but calling a
// recalc() should be a lot more efficient than just repeatedly add()-ing
// the numbers if you've been fiddling directly with data_.

void
RealSet::recalc_sortedness()
{
  sortedness_ = 2;
  const double* d = data_.data_const();
  const unsigned N = data_.size();

  for(unsigned i=1; i<N && sortedness_; ++i) {
    double x = d[i];
    double prev = d[i-1];
    if (sortedness_ == 2) {
      if (prev < x)
	sortedness_ = 1;
      else if (prev > x)
	sortedness_ = -1;
    } else if (sortedness_ == 1) {
      if (prev > x)
	sortedness_ = 0;
    } else if (sortedness_ == -1) {
      if (prev < x)
	sortedness_ = 0;
    }
  }
}

void
RealSet::recalc()
{
  recalc_saving_sorted();
  sdata_.clear();
}

void
RealSet::recalc_saving_sorted()
{
  mean_ = sumsq_ = 0;
  sortedness_ = 2;
  if (size() == 0) return;

  const double* d = data_.data_const();
  const unsigned N = data_.size();

  min_ = max_ = d[0];
  double old_mean;
  for (unsigned i=0; i<N; ++i) {
    double x = d[i];
    if (min_ > x) min_ = x;
    if (max_ < x) max_ = x;
    old_mean = mean_;
    mean_ += (x - mean_)/(i+1);
    if (i>0) {
      sumsq_ += (x-mean_)*(x-old_mean);
      if (sortedness_) {
	double prev = d[i-1];
	if (sortedness_ == 2) {
	  if (prev < x)
	    sortedness_ = 1;
	  else if (prev > x)
	    sortedness_ = -1;
	} else if (sortedness_ == 1) {
	  if (prev > x)
	    sortedness_ = 0;
	} else if (sortedness_ == -1) {
	  if (prev < x)
	    sortedness_ = 0;
	}
      }
    }
  }

  sdata_.clear();
}

void
RealSet::copy(const RealSet& ds)
{
  if (this == &ds) return; // can't copy ourselves.

  clear();
  set_label(ds.label());
  data_ = ds.data_;
  sdata_ = ds.sdata_;
  min_ = ds.min_;
  max_ = ds.max_;
  mean_ = ds.mean_;
  sumsq_ = ds.sumsq_;
  sortedness_ = ds.sortedness_;
  epsilon_ = ds.epsilon_;
}



// $Id: RealSet.cpp,v 1.13 1999/07/20 20:52:09 trow Exp $