examples/cech-complex/Miniball_dynamic_d.h
 author Dmitriy Morozov Thu, 27 Dec 2007 00:00:00 -0500 changeset 43 0f970f454094 parent 18 a4d403087fbc permissions -rw-r--r--
Added initial Avida code (computes 0-dimensional persistence of the hamming distance function)

//    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,
//
//    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;
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;
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
{
}

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();
}

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;
}

{
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 for the center of the empty sphere
for (int j=0; j<d; ++j)
c[j]=0;
current_c = c;
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[i] = q0[i];
sqr_r = 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 = 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 -= l[i];
}
if (l < min_l) min_l = l;
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;
}