MorseSmaleComplex.cpp

/*  _______________________________________________________________________

    DAKOTA: Design Analysis Kit for Optimization and Terascale Applications
    Copyright (c) 2010, Sandia National Laboratories.
    This software is distributed under the GNU Lesser General Public License.
    For more information, see the README file in the top Dakota directory.
    _______________________________________________________________________ */

//- Class:        MorseSmaleComplex
//- Description:  Class implementation for MorseSmaleComplex
//- Owner:        Dan Maljovec, Mohamed Ebeida
//- Version: $Id$

#include "MorseSmaleComplex.hpp"
#include <map>
#include <vector>
// From the Dionysus package
#include "topology/persistence-diagram.h"
#define USING_GL
#ifdef USING_GL
  #include <GL/glut.h>
#endif

//using namespace std;

///
//Vertex
//
Vertex::Vertex(int n, double *p, double _val, int _id)//n值为维数1(dim-1);p指向p[0],p[1];_val值为p[1];_id为数据中第i个序号
{
	d = n;
	x = new double[n];
	for(int i=0; i < n; i++)
		x[i] = p[i];

	e = -1;
	PC_UF_min = PC_UF_max = UF_min = UF_max = ID = _id;
	persistence = 0;
	val = _val;
	classification = REGULAR;
}

double Vertex::GetXi(int i)
{
	return x[i];
}

int Vertex::GetIthNeighbor(int i, KNN_Edge * E[])
{
	int neighborID = E[this->e]->end;
	int nextKNN_EdgeID = E[this->e]->nextKNN_Edge;
	for(int j = 0; j < i; j++)
	{
		if(nextKNN_EdgeID != -1)
		{
			neighborID = E[nextKNN_EdgeID]->end;
			nextKNN_EdgeID = E[nextKNN_EdgeID]->nextKNN_Edge;
		}
		else
		{
			neighborID = -1;
		}
	}
	return neighborID;
}

void Vertex::CreateANNpoint(ANNpoint &p)
{
	if(p != NULL)
		delete [] p;
	p = new ANNcoord[d];
	for(int i = 0; i < d; i++)
		p[i] = this->x[i];
}

void Vertex::Union_Max(Vertex * v, Vertex * V[])
{
	Vertex *a = V[this->Find_Max(V)];
	Vertex *b = V[v->Find_Max(V)];
	if(a->val > b->val)
	{
		b->UF_max = this->ID;
		b->PC_UF_max = a->ID;
	}
	else if(a->val == b->val)
	{
		if(a->ID > b->ID)
		{
			b->UF_max = this->ID;
			b->PC_UF_max = a->ID;
		}
		else
		{
			a->UF_max = v->ID;
			a->PC_UF_max = b->ID;
		}
	}
	else
	{
		a->UF_max = v->ID;
		a->PC_UF_max = b->ID;
	}
}
int Vertex::Find_Max(Vertex * V[])
{
	if(this->UF_max == this->ID)
		return this->ID;

	this->PC_UF_max = V[this->PC_UF_max]->Find_Max(V);	
	return this->PC_UF_max;
}

double Vertex::SDistance(Vertex *v)
{
	double sDist = 0;
	for(int i = 0; i < d; i++)
		sDist += (x[i]-v->x[i])*(x[i]-v->x[i]);
	return sDist;
}

void Vertex::Union_Min(Vertex * v, Vertex * V[])
{
	Vertex *a = V[this->Find_Min(V)];
	Vertex *b = V[v->Find_Min(V)];
	if(a->val < b->val)
	{
		b->UF_min = this->ID;
		b->PC_UF_min = a->ID;
	}
	else if(a->val == b->val)
	{
		if(a->ID < b->ID)
		{
			b->UF_min = this->ID;
			b->PC_UF_min = a->ID;
		}
		else
		{
			a->UF_min = v->ID;
			a->PC_UF_min = b->ID;
		}
	}
	else
	{
		a->UF_min = v->ID;
		a->PC_UF_min = b->ID;
	}
}
int Vertex::Find_Min(Vertex * V[])
{
	if(this->UF_min == this->ID)
		return this->ID;

	this->PC_UF_min = V[this->PC_UF_min]->Find_Min(V);	
	return this->PC_UF_min;
}
void Vertex::ResetExtrema()
{
	UF_max = UF_min = PC_UF_min = PC_UF_max = ID;
	persistence = 0;
}
double Vertex::Value() { return val; }
Vertex::Vertex(Vertex &v)
{
	classification = v.classification;
	d = v.d;
	e = v.e;
	ID = v.ID;
	PC_UF_max = v.PC_UF_max;
	PC_UF_min = v.PC_UF_min;
	UF_max = v.UF_max;
	UF_min = v.UF_max;
	val = v.val;
	x = new double[d];
	for(int i = 0; i < d; i++)
		x[i] = v.x[i];
	persistence = 0;
}
Vertex::~Vertex() { delete [] x; }
///
//KNN_Edge
//
KNN_Edge::KNN_Edge(Vertex * _start, Vertex * _end, KNN_Edge * E[], int id)
{
	ID = id;
	E[id] = this;

	start = _start->ID;
	end = _end->ID;

	this->nextKNN_Edge = _start->e;
	_start->e = id;
}
///
//Crystal
//
MS_Crystal::MS_Crystal(Vertex ** _V, int *_vIds, int _count, double _d)
{
	V = _V;
	v = _vIds;
	count = _count;
	minV = v[count-2];
	maxV = v[count-1];
	minP = 0;
	maxP = 100;

	d = _d;
}

MS_Crystal::~MS_Crystal() { delete [] v; }

MS_Crystal::MS_Crystal(MS_Crystal *a, MS_Crystal *b, bool mergeMax, int persistence)
{
	minP = a->maxP = b->maxP = persistence;
	V = a->V;
	v = new int[a->count+b->count];
	int i = 0;
	for( ; i < a->count; i++)
		v[i] = a->v[i];
	for( ; i < a->count+b->count; i++)
		v[i] = b->v[i-a->count];

	count = a->count+b->count;
	d = a->d;

	if(mergeMax)
	{
		if(a->minV != b->minV)
			printf("We have a problem.");
		else
			minV = a->minV;

		if(V[a->maxV]->persistence > V[b->maxV]->persistence)
			maxV = a->maxV;
		else
			maxV = b->maxV;
	}
	else
	{
		if(a->maxV != b->maxV)
			printf("We have a problem.");
		else
			maxV = a->maxV;

		if(V[a->minV]->persistence > V[b->minV]->persistence)
			minV = a->minV;
		else
			minV = b->minV;
	}

	maxP = 100;
}
MS_Crystal::MS_Crystal(MS_Crystal *a, bool mergeMax, int nuExtrema, int persistence)
{
	minP = a->maxP = persistence;
	V = a->V;
	v = new int[a->count+1];
	int i = 0;
	for( ; i < a->count; i++)
		v[i] = a->v[i];
	v[i] = nuExtrema;

	count = a->count+1;
	d = a->d;

	if(mergeMax)
	{
			minV = a->minV;
			maxV = nuExtrema;
	}
	else
	{
			maxV = a->maxV;
			minV = nuExtrema;
	}

	maxP = 100;
}
int MS_Crystal::size() {	return count;	}
int MS_Crystal::GetVertexID(int i) {	return v[i];	}
bool MS_Crystal::Exists(double persistence)
	{ return (persistence >= minP && persistence < maxP); }

int MS_Crystal::MaxV() { return maxV; }
int MS_Crystal::MinV() { return minV; }
///
//MS Complex
//
void swap(int i, int j, Saddle **S)
{
	Saddle *temp = S[i];
	S[i] = S[j];
	S[j] = temp;
}
int partition(Saddle **S, int left, int right, int pivot)
{
	swap(pivot, right, S);//pivot is now in right
	int insertionIdx = left;
	for(int i = left; i < right; i++)
	{
		if(S[i]->persistence < S[right]->persistence)
		{
			swap(i,insertionIdx, S);
			insertionIdx++;
		}
	}
	swap(insertionIdx, right, S);// put pivot in place
	return insertionIdx;
}
void quickSort(Saddle **S, int left, int right)
{
	if(left < right)
	{
		int pivot = (right + left) / 2;
		int nuPivot = partition(S, left, right, pivot);
		quickSort(S,left, nuPivot-1);
		quickSort(S, nuPivot+1, right);
	}
}
void SortSaddles(Saddle **S, int n)
{
	quickSort(S, 0, n-1);
}
MS_Complex::MS_Complex(double  *points, int dimension, int count, int _k, bool perturb)
{
  perturbed = perturb;
	d = dimension-1;
	numKneighbors = _k;
	szV = numV = count;
	//szE = szV*numKneighbors;
  //Edges should be bi-directional
	szE = szV*szV;

	V = new Vertex *[szV];
	E = new KNN_Edge *[szE];

  srand(8);	//随机产生eps，值大小为1e-7
	double *p = new double[dimension];
	for(int i=0; i < count; i++)
	{
    double eps = (double)rand() / (double)RAND_MAX;
    eps = eps *1e-6;
    if(rand() > RAND_MAX / 2)
      eps = -eps;

	for(int j = 0; j < dimension; j++)
		p[j] = points[i*dimension+j];
    p[dimension-1] += (perturbed ? eps : 0);
		V[i] = new Vertex(d, p, p[d], i);
	}	
	delete [] p;

//  std::cout << "numV=" << numV << std::endl;
//  std::cout << "d=" << d << std::endl;
//  std::cout << "numKneighbors=" << numKneighbors << std::endl;

	KNN();
	Compute();
  maxDist = -1;
}

MS_Complex::MS_Complex(MS_Complex &Complex, double  *new_point)
{
  perturbed = Complex.perturbed;
  d = Complex.d;
  numKneighbors = Complex.numKneighbors;
	szV = numV = Complex.numV + 1;
	//szE = szV*numKneighbors;
  szE = szV*szV;

	V = new Vertex *[szV];
	E = new KNN_Edge *[szE];
	
	double *p = new double[d+1];
	for(int i=0; i < Complex.numV; i++)
	{
		for(int j = 0; j < d; j++)
			p[j] = Complex.V[i]->GetXi(j);
    p[d] = Complex.V[i]->Value();
		V[i] = new Vertex(d, p, p[d], i);
	}	

  double eps = (double)rand() / (double)RAND_MAX;
  eps = eps * 1e-6;
  if(rand() > RAND_MAX / 2)
    eps = -eps;

  V[numV-1] = new Vertex(d, new_point, new_point[d] + (perturbed ? eps : 0), numV-1);
	delete [] p;
	KNN();
	Compute();
  maxDist = -1;
}

void MS_Complex::KNN()
{
	int eIndex = 0;
	int i,k;

	ANNpointArray pa = new ANNpoint[numV];
	for(i = 0; i < numV; i++)
	{
		pa[i] = NULL;
		V[i]->CreateANNpoint(pa[i]);
	}

	SearchStructure = new ANNkd_tree(pa,numV,d);
	ANNidxArray nn_idx;
	ANNdistArray dists;
	for(i = 0; i < numV; i++)
	{		
		nn_idx = new ANNidx[numKneighbors+1];
		dists = new ANNdist[numKneighbors+1];

		SearchStructure->annkSearch(pa[i],numKneighbors+1,nn_idx,dists);

		for(k=1;k<numKneighbors+1;k++)
		{
      if(!DoesEdgeExist(i, nn_idx[k], eIndex))
      {
			  E[eIndex] = new KNN_Edge(V[i],V[nn_idx[k]],E, eIndex);
			  eIndex++;
      }
      if(!DoesEdgeExist(nn_idx[k],i, eIndex))
      {
			  E[eIndex] = new KNN_Edge(V[nn_idx[k]],V[i],E, eIndex);
			  eIndex++;
      }
		}
		delete [] nn_idx;
		delete [] dists;
	}

	for(i = 0; i < numV; i++)
		delete [] pa[i];

	delete [] pa;
	numE = eIndex;
}

void MS_Complex::Compute()
{
  int i,j;
	int maxCount=0;
	int minCount=0;
	double globalMax = V[0]->Value();
	double globalMin = V[0]->Value();
	for(i = 0; i < numV; i++)
	{
		Vertex *v = V[i];

		if(globalMax < v->Value())
			globalMax = v->Value();
		if(globalMin > v->Value())
			globalMin = v->Value();

		//Vertex **neighbors = new Vertex *[numKneighbors];
    std::vector<Vertex *> neighbors;
		//for(int j = 0; j < numKneighbors; j++)
    j = 0;
    while(v->GetIthNeighbor(j,E) != -1)
		{
			//neighbors[j] = V[v->GetIthNeighbor