Changed implementation of WeightedRips to store simplex values (max distance between simplices' vertices) as an invisible layer on top of each simplex object, so that the data() field of WeightedRips has been freed for use by the users again.
// Copright (C) 1999-2006, Bernd Gaertner
// $Revision: 1.3 $
// $Date: 2006/11/16 08:01:52 $
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program 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 General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program; if not, write to the Free Software
// Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA,
// or download the License terms from prep.ai.mit.edu/pub/gnu/COPYING-2.0.
//
// Contact:
// --------
// Bernd Gaertner
// Institute of Theoretical Computer Science
// ETH Zuerich
// CAB G32.2
// CH-8092 Zuerich, Switzerland
// http://www.inf.ethz.ch/personal/gaertner
#include <cstdlib>
#include <cassert>
#include <cmath>
#include <iostream>
#include <list>
// Functions
// =========
inline double mb_sqr (double r) {return r*r;}
// Class Declarations
// ==================
// smallest enclosing ball of a set of n points in dimension d
class Miniball;
// smallest ball with a set of n <= d+1 points *on the boundary*
class Miniball_b;
// point in dimension d
class Point;
// Class Definitions
// =================
// Miniball_b
// ----------
class Miniball_b {
private:
// data members
int d; // dimension
int m, s; // size and number of support points
double* q0;
double* z;
double* f;
double** v;
double** a;
double** c;
double* sqr_r;
double* current_c; // refers to some c[j]
double current_sqr_r;
// if you want the copy constructor and assignment operator, please
// properly implement them yourself. The default copy/assignment
// semantics fails since a Miniball_b object stores pointers to
// dynamically allocated memory
Miniball_b (const Miniball_b& mb);
Miniball_b& operator=(const Miniball_b& mb);
public:
Miniball_b(int dim)
: d(dim)
{
q0 = new double[d];
z = new double[d+1];
f = new double[d+1];
v = new double*[d+1];
for (int i=0; i<d+1; ++i) v[i] = new double[d];
a = new double*[d+1];
for (int i=0; i<d+1; ++i) a[i] = new double[d];
c = new double*[d+1];
for (int i=0; i<d+1; ++i) c[i] = new double[d];
sqr_r = new double[d+1];
reset();
}
~Miniball_b()
{
delete[] sqr_r;
for (int i=0; i<d+1; ++i) delete[] c[i];
delete[] c;
for (int i=0; i<d+1; ++i) delete[] a[i];
delete[] a;
for (int i=0; i<d+1; ++i) delete[] v[i];
delete[] v;
delete[] f;
delete[] z;
delete[] q0;
}
// access
const double* center() const;
double squared_radius() const;
int size() const;
int support_size() const;
double excess (const Point& p) const;
// modification
void reset(); // generates empty sphere with m=s=0
bool push (const Point& p);
void pop ();
// checking
double slack() const;
};
// Miniball
// --------
class Miniball {
public:
// types
typedef std::list<Point>::iterator It;
typedef std::list<Point>::const_iterator Cit;
private:
// data members
int d; // dimension
std::list<Point> L; // internal point set
Miniball_b B; // the current ball
It support_end; // past-the-end iterator of support set
// private methods
void mtf_mb (It k);
void pivot_mb (It k);
void move_to_front (It j);
double max_excess (It t, It i, It& pivot) const;
public:
// creates an empty ball
Miniball(int dim) : d(dim), B(dim) {}
// copies p to the internal point set
void check_in (const Point& p);
// builds the smallest enclosing ball of the internal point set
void build ();
// returns center of the ball (undefined if ball is empty)
Point center() const;
// returns squared_radius of the ball (-1 if ball is empty)
double squared_radius () const;
// returns size of internal point set
int nr_points () const;
// returns begin- and past-the-end iterators for internal point set
Cit points_begin () const;
Cit points_end () const;
// returns size of support point set; this set has the following properties:
// - there are at most d+1 support points,
// - all support points are on the boundary of the computed ball, and
// - the smallest enclosing ball of the support point set equals the
// smallest enclosing ball of the internal point set
int nr_support_points () const;
// returns begin- and past-the-end iterators for internal point set
Cit support_points_begin () const;
Cit support_points_end () const;
// assesses the quality of the computed ball. The return value is the
// maximum squared distance of any support point or point outside the
// ball to the boundary of the ball, divided by the squared radius of
// the ball. If everything went fine, this will be less than e-15 and
// says that the computed ball approximately contains all the internal
// points and has all the support points on the boundary.
//
// The slack parameter that is set by the method says something about
// whether the computed ball is really the *smallest* enclosing ball
// of the support points; if everything went fine, this value will be 0;
// a positive value may indicate that the ball is not smallest possible,
// with the deviation from optimality growing with the slack
double accuracy (double& slack) const;
// returns true if the accuracy is below the given tolerance and the
// slack is 0
bool is_valid (double tolerance = 1e-15) const;
};
// Point (inline)
// --------------
class Point {
private:
int d;
double* coord;
public:
// default
Point(int dim)
: d (dim), coord(new double[dim])
{}
~Point ()
{
delete[] coord;
}
// copy from Point
Point (const Point& p)
: d (p.dim()), coord(new double[p.dim()])
{
for (int i=0; i<d; ++i)
coord[i] = p.coord[i];
}
// copy from double*
Point (int dim, const double* p)
: d (dim), coord(new double[dim])
{
for (int i=0; i<d; ++i)
coord[i] = p[i];
}
// assignment
Point& operator = (const Point& p)
{
assert (d == p.dim());
if (this != &p)
for (int i=0; i<d; ++i)
coord[i] = p.coord[i];
return *this;
}
// coordinate access
double& operator [] (int i)
{
return coord[i];
}
const double& operator [] (int i) const
{
return coord[i];
}
const double* begin() const
{
return coord;
}
const double* end() const
{
return coord+d;
}
// dimension access
int dim() const
{
return d;
}
};
// Class Implementations
// =====================
// Miniball
// --------
void Miniball::check_in (const Point& p)
{
assert (d == p.dim());
L.push_back(p);
}
void Miniball::build ()
{
B.reset();
support_end = L.begin();
pivot_mb (L.end());
}
void Miniball::mtf_mb (It i)
{
support_end = L.begin();
if ((B.size())==d+1) return;
for (It k=L.begin(); k!=i;) {
It j=k++;
if (B.excess(*j) > 0) {
if (B.push(*j)) {
mtf_mb (j);
B.pop();
move_to_front(j);
}
}
}
}
void Miniball::move_to_front (It j)
{
if (support_end == j)
support_end++;
L.splice (L.begin(), L, j);
}
void Miniball::pivot_mb (It i)
{
It t = ++L.begin();
mtf_mb (t);
double max_e, old_sqr_r = -1;
do {
It pivot;
max_e = max_excess (t, i, pivot);
if (max_e > 0) {
t = support_end;
if (t==pivot) ++t;
old_sqr_r = B.squared_radius();
B.push (*pivot);
mtf_mb (support_end);
B.pop();
move_to_front (pivot);
}
} while ((max_e > 0) && (B.squared_radius() > old_sqr_r));
}
double Miniball::max_excess (It t, It i, It& pivot) const
{
const double *c = B.center(), sqr_r = B.squared_radius();
double e, max_e = 0;
for (It k=t; k!=i; ++k) {
const double *p = (*k).begin();
e = -sqr_r;
for (int j=0; j<d; ++j)
e += mb_sqr(p[j]-c[j]);
if (e > max_e) {
max_e = e;
pivot = k;
}
}
return max_e;
}
Point Miniball::center () const
{
return Point(d, B.center());
}
double Miniball::squared_radius () const
{
return B.squared_radius();
}
int Miniball::nr_points () const
{
return L.size();
}
Miniball::Cit Miniball::points_begin () const
{
return L.begin();
}
Miniball::Cit Miniball::points_end () const
{
return L.end();
}
int Miniball::nr_support_points () const
{
return B.support_size();
}
Miniball::Cit Miniball::support_points_begin () const
{
return L.begin();
}
Miniball::Cit Miniball::support_points_end () const
{
return support_end;
}
double Miniball::accuracy (double& slack) const
{
double e, max_e = 0;
int n_supp=0;
Cit i;
for (i=L.begin(); i!=support_end; ++i,++n_supp)
if ((e = std::abs (B.excess (*i))) > max_e)
max_e = e;
// you've found a non-numerical problem if the following ever fails
assert (n_supp == nr_support_points());
for (i=support_end; i!=L.end(); ++i)
if ((e = B.excess (*i)) > max_e)
max_e = e;
slack = B.slack();
return (max_e/squared_radius());
}
bool Miniball::is_valid (double tolerance) const
{
double slack;
return ( (accuracy (slack) < tolerance) && (slack == 0) );
}
// Miniball_b
// ----------
const double* Miniball_b::center () const
{
return current_c;
}
double Miniball_b::squared_radius() const
{
return current_sqr_r;
}
int Miniball_b::size() const
{
return m;
}
int Miniball_b::support_size() const
{
return s;
}
double Miniball_b::excess (const Point& p) const
{
double e = -current_sqr_r;
for (int k=0; k<d; ++k)
e += mb_sqr(p[k]-current_c[k]);
return e;
}
void Miniball_b::reset ()
{
m = s = 0;
// we misuse c[0] for the center of the empty sphere
for (int j=0; j<d; ++j)
c[0][j]=0;
current_c = c[0];
current_sqr_r = -1;
}
void Miniball_b::pop ()
{
--m;
}
bool Miniball_b::push (const Point& p)
{
int i, j;
double eps = 1e-32;
if (m==0) {
for (i=0; i<d; ++i)
q0[i] = p[i];
for (i=0; i<d; ++i)
c[0][i] = q0[i];
sqr_r[0] = 0;
} else {
// set v_m to Q_m
for (i=0; i<d; ++i)
v[m][i] = p[i]-q0[i];
// compute the a_{m,i}, i< m
for (i=1; i<m; ++i) {
a[m][i] = 0;
for (j=0; j<d; ++j)
a[m][i] += v[i][j] * v[m][j];
a[m][i]*=(2/z[i]);
}
// update v_m to Q_m-\bar{Q}_m
for (i=1; i<m; ++i) {
for (j=0; j<d; ++j)
v[m][j] -= a[m][i]*v[i][j];
}
// compute z_m
z[m]=0;
for (j=0; j<d; ++j)
z[m] += mb_sqr(v[m][j]);
z[m]*=2;
// reject push if z_m too small
if (z[m]<eps*current_sqr_r) {
return false;
}
// update c, sqr_r
double e = -sqr_r[m-1];
for (i=0; i<d; ++i)
e += mb_sqr(p[i]-c[m-1][i]);
f[m]=e/z[m];
for (i=0; i<d; ++i)
c[m][i] = c[m-1][i]+f[m]*v[m][i];
sqr_r[m] = sqr_r[m-1] + e*f[m]/2;
}
current_c = c[m];
current_sqr_r = sqr_r[m];
s = ++m;
return true;
}
double Miniball_b::slack () const
{
double l[d+1], min_l=0;
l[0] = 1;
for (int i=s-1; i>0; --i) {
l[i] = f[i];
for (int k=s-1; k>i; --k)
l[i]-=a[k][i]*l[k];
if (l[i] < min_l) min_l = l[i];
l[0] -= l[i];
}
if (l[0] < min_l) min_l = l[0];
return ( (min_l < 0) ? -min_l : 0);
}
// Point
// -----
// Output
std::ostream& operator << (std::ostream& os, const Point& p)
{
os << "(";
int d = p.dim();
for (int i=0; i<d-1; ++i)
os << p[i] << ", ";
os << p[d-1] << ")";
return os;
}