realset.cpp

非常著名的曲线拟合程序

// This is -*- C++ -*-
// $Id: RealSet.cpp,v 1.13 1999/07/20 20:52:09 trow Exp $

/* 
 * RealSet.cpp
 *
 * Copyright (C) 1998, 1999 EMC Capital Management, Inc.
 *
 * Developed by Jon Trowbridge <[e-mail]>.
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Library General Public
 * License as published by the Free Software Foundation; either
 * version 2 of the License, or (at your option) any later version.
 *
 * This library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Library General Public License for more details.
 *
 * You should have received a copy of the GNU Library General Public
 * License along with this library; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
 * 02111-1307, USA.
 */

#include <config.h>
#include <string>
#include <stdlib.h>
#include <algorithm>
#include <functional>
#include <iostream>
#include <iomanip>
#include <zdouble.h>
#include "RealSet.h"

// some random number...
const DataSet::type_t RealSet::typecode = 0xeaeff;

RealSet::RealSet() : DataSet_Typed<double>(RealSet::typecode), epsilon_(1e-8)
{
  clear();
}

RealSet::RealSet(const RealSet& ds)
  : DataSet_Typed<double>(RealSet::typecode)
{
  copy(ds);
}

RealSet&
RealSet::operator=(const RealSet& ds)
{
  copy(ds);
  return *this;
}

const double*
RealSet::sorted_data() const
{
  if (size() == 0) return 0; 
  if (sdata_.size()) return sdata_.data();

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

  if (sortedness() == 0) {
    std::sort(s, s+N, less<double>());
  } else if (sortedness() == -1) {
    for(unsigned i=0; i<N/2; ++i) {
      double t = s[i];
      s[i] = s[N-1-i];
      s[N-1-i] = t;
    }
  }

  return (const double*)s;
}

unsigned
RealSet::greater_than(double x) const
{
  if (size() == 0 || max() <= x) return 0;
  if (min() > x) return size();

  const double* d = sorted_data();
  const unsigned N = size();
  unsigned i=0, j=N-1;

  for(;;) {
    unsigned mid=(i+j)/2;
    if (mid == N-1) --mid;
    if (d[mid] <= x && x < d[mid+1])
      return N-mid-1;
    else if (d[mid] > x) 
      j = mid;
    else
      i = mid;
  }
}

unsigned
RealSet::less_than(double x) const
{
  if (size() == 0 || min() >= x) return 0;
  if (max() < x) return size();

  const double* d = sorted_data();
  const unsigned N = size();
  unsigned i=0, j=N-1;

  for(;;) {
    unsigned mid=(i+j)/2;
    if (mid == N-1) --mid;
    if (d[mid] < x && x <= d[mid+1])
      return mid+1;
    else if (d[mid] >= x)
      j = mid;
    else
      i = mid;
  }
}

void
RealSet::reserve(unsigned R)
{
  data_.reserve(R);
}

void
RealSet::resize(unsigned N)
{
  if (size() > N) {
    remove(N, size());
  } else {
    while(size() < N)
      add(0.0);
  }
}

void
RealSet::add(double x)
{
  update_add(x);
  data_.add(x);
  if (sortedness_ > 0)
    sdata_ = data_;
  else
    sdata_.clear();
}

void
RealSet::add_at(int i, double x)
{
#ifdef REALSET_BOUNDS_CHECKING
  if (! in_bounds(i))
    throw Exception_IndexOutOfBounds(i, min_index(), max_index());
#endif

  update_add(x);
  data_.add_at(i - min_index(),x);
  recalc_sortedness(); // sortedness will be wrong from update_add().
  if (sortedness_ > 0)
    sdata_ = data_;
  else
    sdata_.clear();
}

void
RealSet::remove(int i)
{
#ifdef REALSET_BOUNDS_CHECKING
  if (! in_bounds(i))
    throw Exception_IndexOutOfBounds(i, min_index(), max_index());
#endif

  update_remove(data(i));
  data_.remove_at(i - min_index());

  // Remember, removing an element preserves sortedness.  The only possible
  // change is to bring an unsorted list into a sorted state.
  if (sortedness_ == 0)
    recalc_sortedness();

  if (sortedness_ > 0)
    sdata_ = data_;
  else
    sdata_.clear();
}

void
RealSet::add(double* xs, size_t count)
{
  data_.add(xs,count);
  recalc();
}
void
RealSet::add_at(int i, double* xs, size_t count)
{
#ifdef REALSET_BOUNDS_CHECKING
  if (! in_bounds(i))
    throw Exception_IndexOutOfBounds(i,min_index(), max_index());
#endif
  if (i == max_index())
    add(xs, count);
  else {
    data_.add_at(i - min_index(),xs,count);
    recalc();
  }
}

void
RealSet::remove(int i, int j)
{
  if (i > j)
    throw Exception("Illegal range");

#ifdef REALSET_BOUNDS_CHECKING
  if (! in_bounds(j))
    throw Exception_IndexOutOfBounds(j,min_index(),max_index());
#endif

  data_.remove_at(i-min_index(), j-min_index());
  recalc();
}

void
RealSet::set(int i, double x)
{
#ifdef REALSET_BOUNDS_CHECKING
  if (! in_bounds(i))
    throw Exception_IndexOutOfBounds(i,min_index(),max_index());
#endif

  data_.data()[i - min_index()] = x;
  recalc();
}

void
RealSet::set(int i, const double* xs, size_t count)
{
#ifdef REALSET_BOUNDS_CHECKING
  if (! in_bounds(i))
    throw Exception_IndexOutOfBounds(i,min_index(),max_index());
  if (! in_bounds(i+count-1))
    throw Exception_IndexOutOfBounds(i+count-1,min_index(),max_index());
#endif
    
  double* array = data_.data();
  unsigned j = 0;
  while (j < count) {
    array[i+j - min_index()] = xs[j];
    ++j;
  }
  recalc();
}

void
RealSet::clear()
{
  set_offset(0);
  data_.clear();
  sdata_.clear();
  mean_ = sumsq_ = 0;
  sortedness_ = 2; // constant
}

void
RealSet::append(const RealSet& ds)
{
  // clearly, this could be optimized...

  // By saving N before we iterate, we should properly handle
  // the case where &ds == this.
  const double* d = ds.data();
  const unsigned N = ds.size();
  for(unsigned i=0; i<N; ++i)
    add(d[i]);
}

/*
void
RealSet::apply(double (*fn)(double))
{
  double* d = data_.data();
  for(size_t i=0; i<size(); ++i)
    d[i] = fn(d[i]);
  recalc();
}
*/

void
RealSet::linear_transform(double a, double b)
{
  double* d = data_.data();
  const unsigned N = size();

  for(unsigned i=0; i<N; ++i) 
    d[i] = a*d[i]+b;

  
  double m1 = a*min_+b, m2 = a*max_+b;
  min_ = m1<m2 ? m1 : m2;
  max_ = m1<m2 ? m2 : m1;

  mean_ = a*mean_+b;
  sumsq_ *= a*a;
  if (a<0 && sortedness_ != 2)
    sortedness_ = -sortedness_;

  if (sdata_.size()) {
    double* s = sdata_.data();
    // linearly transform our sdata buffer
    for(unsigned i=0; i<N; ++i)
      s[i] = a*s[i]+b;

    // if a<0, reverse the sdata to preserve sortedness
    if (a<0) {
      for(unsigned i=0; i<N/2; ++i) {
	double t=s[i];
	s[i]=s[N-1-i];
	s[N-1-i]=t;
	}
    }
  }
}

void
RealSet::log_transform()
{
  // We should throw an exception here...
  if (min() <= 0) {
    clear();
    return;
  }

  const unsigned N = size();
  double* d = data_.data();
  for (unsigned i=0; i<N; ++i)
    d[i] = log(d[i]);

  if (have_sorted_data()) {
    double* sd = sdata_.data();
    for (unsigned i=0; i<N; ++i)
      sd[i] = log(sd[i]);
  }

  recalc_saving_sorted();
}

void
RealSet::exp_transform()
{
  const unsigned N = size();
  double* d = data_.data();
  for (unsigned i=0; i<N; ++i)
    d[i] = exp(d[i]);

  if (have_sorted_data()) {
    double* sd = sdata_.data();
    for (unsigned i=0; i<N; ++i)
      sd[i] = exp(sd[i]);
  }
  
  recalc_saving_sorted();
}

void
RealSet::logit_transform()
{
  // We should throw an exception here...
  if (min() <= 0 || max() >= 1) {
    clear();
    return;
  }

  // This could be optimized to preserve sorted data.
  double* d = data_.data();
  const unsigned N = size();
  for(unsigned i=0; i<N; ++i)
    d[i] = log(d[i]/(1-d[i]));

  recalc();
}

void
RealSet::fabs_transform()
{
  // This could be optimized to preserve sorted data.

  // Check for the cases that require less work
  if (max() <= 0) {
    linear_transform(-1, 0);
  } else if (min() < 0) {
    const unsigned N = size();
    double* d = data_.data();
    for(unsigned i=0; i<N; ++i)
      d[i] = fabs(d[i]);
    recalc();
  }
}

void
RealSet::deviations_transform(double y)