1 /* _______________________________________________________________________
2
3 DAKOTA: Design Analysis Kit for Optimization and Terascale Applications
4 Copyright 2014-2020 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
5 This software is distributed under the GNU Lesser General Public License.
6 For more information, see the README file in the top Dakota directory.
7 _______________________________________________________________________ */
8
9 //- Class: MorseSmaleComplex
10 //- Description: Class implementation for MorseSmaleComplex
11 //- Owner: Dan Maljovec, Mohamed Ebeida
12 //- Version: $Id$
13
14 #include "MorseSmaleComplex.hpp"
15 #include <map>
16 #include <vector>
17 // From the Dionysus package
18 #include "topology/persistence-diagram.h"
19
20 #ifdef USING_GL
21 #include <GL/glut.h>
22 #endif
23
24 //using namespace std;
25
26 ///////////////////////////////////////////////////
27 //Vertex
28 //////////////////////////////////////////////////
Vertex(int n,double * p,double _val,int _id)29 Vertex::Vertex(int n, double *p, double _val, int _id)
30 {
31 d = n;
32 x = new double[n];
33 for(int i=0; i < n; i++)
34 x[i] = p[i];
35
36 e = -1;
37 PC_UF_min = PC_UF_max = UF_min = UF_max = ID = _id;
38 persistence = 0;
39 val = _val;
40 classification = REGULAR;
41 }
42
GetXi(int i)43 double Vertex::GetXi(int i)
44 {
45 return x[i];
46 }
47
GetIthNeighbor(int i,KNN_Edge * E[])48 int Vertex::GetIthNeighbor(int i, KNN_Edge * E[])
49 {
50 int neighborID = E[this->e]->end;
51 int nextKNN_EdgeID = E[this->e]->nextKNN_Edge;
52 for(int j = 0; j < i; j++)
53 {
54 if(nextKNN_EdgeID != -1)
55 {
56 neighborID = E[nextKNN_EdgeID]->end;
57 nextKNN_EdgeID = E[nextKNN_EdgeID]->nextKNN_Edge;
58 }
59 else
60 {
61 neighborID = -1;
62 }
63 }
64 return neighborID;
65 }
66
CreateANNpoint(ANNpoint & p)67 void Vertex::CreateANNpoint(ANNpoint &p)
68 {
69 if(p != NULL)
70 delete [] p;
71 p = new ANNcoord[d];
72 for(int i = 0; i < d; i++)
73 p[i] = this->x[i];
74 }
75
Union_Max(Vertex * v,Vertex * V[])76 void Vertex::Union_Max(Vertex * v, Vertex * V[])
77 {
78 Vertex *a = V[this->Find_Max(V)];
79 Vertex *b = V[v->Find_Max(V)];
80 if(a->val > b->val)
81 {
82 b->UF_max = this->ID;
83 b->PC_UF_max = a->ID;
84 }
85 else if(a->val == b->val)
86 {
87 if(a->ID > b->ID)
88 {
89 b->UF_max = this->ID;
90 b->PC_UF_max = a->ID;
91 }
92 else
93 {
94 a->UF_max = v->ID;
95 a->PC_UF_max = b->ID;
96 }
97 }
98 else
99 {
100 a->UF_max = v->ID;
101 a->PC_UF_max = b->ID;
102 }
103 }
Find_Max(Vertex * V[])104 int Vertex::Find_Max(Vertex * V[])
105 {
106 if(this->UF_max == this->ID)
107 return this->ID;
108
109 this->PC_UF_max = V[this->PC_UF_max]->Find_Max(V);
110 return this->PC_UF_max;
111 }
112
SDistance(Vertex * v)113 double Vertex::SDistance(Vertex *v)
114 {
115 double sDist = 0;
116 for(int i = 0; i < d; i++)
117 sDist += (x[i]-v->x[i])*(x[i]-v->x[i]);
118 return sDist;
119 }
120
Union_Min(Vertex * v,Vertex * V[])121 void Vertex::Union_Min(Vertex * v, Vertex * V[])
122 {
123 Vertex *a = V[this->Find_Min(V)];
124 Vertex *b = V[v->Find_Min(V)];
125 if(a->val < b->val)
126 {
127 b->UF_min = this->ID;
128 b->PC_UF_min = a->ID;
129 }
130 else if(a->val == b->val)
131 {
132 if(a->ID < b->ID)
133 {
134 b->UF_min = this->ID;
135 b->PC_UF_min = a->ID;
136 }
137 else
138 {
139 a->UF_min = v->ID;
140 a->PC_UF_min = b->ID;
141 }
142 }
143 else
144 {
145 a->UF_min = v->ID;
146 a->PC_UF_min = b->ID;
147 }
148 }
Find_Min(Vertex * V[])149 int Vertex::Find_Min(Vertex * V[])
150 {
151 if(this->UF_min == this->ID)
152 return this->ID;
153
154 this->PC_UF_min = V[this->PC_UF_min]->Find_Min(V);
155 return this->PC_UF_min;
156 }
ResetExtrema()157 void Vertex::ResetExtrema()
158 {
159 UF_max = UF_min = PC_UF_min = PC_UF_max = ID;
160 persistence = 0;
161 }
Value()162 double Vertex::Value() { return val; }
Vertex(Vertex & v)163 Vertex::Vertex(Vertex &v)
164 {
165 classification = v.classification;
166 d = v.d;
167 e = v.e;
168 ID = v.ID;
169 PC_UF_max = v.PC_UF_max;
170 PC_UF_min = v.PC_UF_min;
171 UF_max = v.UF_max;
172 UF_min = v.UF_max;
173 val = v.val;
174 x = new double[d];
175 for(int i = 0; i < d; i++)
176 x[i] = v.x[i];
177 persistence = 0;
178 }
~Vertex()179 Vertex::~Vertex() { delete [] x; }
180 ///////////////////////////////////////////////////
181 //KNN_Edge
182 //////////////////////////////////////////////////
KNN_Edge(Vertex * _start,Vertex * _end,KNN_Edge * E[],int id)183 KNN_Edge::KNN_Edge(Vertex * _start, Vertex * _end, KNN_Edge * E[], int id)
184 {
185 ID = id;
186 E[id] = this;
187
188 start = _start->ID;
189 end = _end->ID;
190
191 this->nextKNN_Edge = _start->e;
192 _start->e = id;
193 }
194 ///////////////////////////////////////////////////
195 //Crystal
196 //////////////////////////////////////////////////
MS_Crystal(Vertex ** _V,int * _vIds,int _count,double _d)197 MS_Crystal::MS_Crystal(Vertex ** _V, int *_vIds, int _count, double _d)
198 {
199 V = _V;
200 v = _vIds;
201 count = _count;
202 minV = v[count-2];
203 maxV = v[count-1];
204 minP = 0;
205 maxP = 100;
206
207 d = _d;
208 }
209
~MS_Crystal()210 MS_Crystal::~MS_Crystal() { delete [] v; }
211
MS_Crystal(MS_Crystal * a,MS_Crystal * b,bool mergeMax,int persistence)212 MS_Crystal::MS_Crystal(MS_Crystal *a, MS_Crystal *b, bool mergeMax, int persistence)
213 {
214 minP = a->maxP = b->maxP = persistence;
215 V = a->V;
216 v = new int[a->count+b->count];
217 int i = 0;
218 for( ; i < a->count; i++)
219 v[i] = a->v[i];
220 for( ; i < a->count+b->count; i++)
221 v[i] = b->v[i-a->count];
222
223 count = a->count+b->count;
224 d = a->d;
225
226 if(mergeMax)
227 {
228 if(a->minV != b->minV)
229 printf("We have a problem.");
230 else
231 minV = a->minV;
232
233 if(V[a->maxV]->persistence > V[b->maxV]->persistence)
234 maxV = a->maxV;
235 else
236 maxV = b->maxV;
237 }
238 else
239 {
240 if(a->maxV != b->maxV)
241 printf("We have a problem.");
242 else
243 maxV = a->maxV;
244
245 if(V[a->minV]->persistence > V[b->minV]->persistence)
246 minV = a->minV;
247 else
248 minV = b->minV;
249 }
250
251 maxP = 100;
252 }
MS_Crystal(MS_Crystal * a,bool mergeMax,int nuExtrema,int persistence)253 MS_Crystal::MS_Crystal(MS_Crystal *a, bool mergeMax, int nuExtrema, int persistence)
254 {
255 minP = a->maxP = persistence;
256 V = a->V;
257 v = new int[a->count+1];
258 int i = 0;
259 for( ; i < a->count; i++)
260 v[i] = a->v[i];
261 v[i] = nuExtrema;
262
263 count = a->count+1;
264 d = a->d;
265
266 if(mergeMax)
267 {
268 minV = a->minV;
269 maxV = nuExtrema;
270 }
271 else
272 {
273 maxV = a->maxV;
274 minV = nuExtrema;
275 }
276
277 maxP = 100;
278 }
size()279 int MS_Crystal::size() { return count; }
GetVertexID(int i)280 int MS_Crystal::GetVertexID(int i) { return v[i]; }
Exists(double persistence)281 bool MS_Crystal::Exists(double persistence)
282 { return (persistence >= minP && persistence < maxP); }
283
MaxV()284 int MS_Crystal::MaxV() { return maxV; }
MinV()285 int MS_Crystal::MinV() { return minV; }
286 ///////////////////////////////////////////////////
287 //MS Complex
288 //////////////////////////////////////////////////
swap(int i,int j,Saddle ** S)289 void swap(int i, int j, Saddle **S)
290 {
291 Saddle *temp = S[i];
292 S[i] = S[j];
293 S[j] = temp;
294 }
partition(Saddle ** S,int left,int right,int pivot)295 int partition(Saddle **S, int left, int right, int pivot)
296 {
297 swap(pivot, right, S);//pivot is now in right
298 int insertionIdx = left;
299 for(int i = left; i < right; i++)
300 {
301 if(S[i]->persistence < S[right]->persistence)
302 {
303 swap(i,insertionIdx, S);
304 insertionIdx++;
305 }
306 }
307 swap(insertionIdx, right, S);// put pivot in place
308 return insertionIdx;
309 }
quickSort(Saddle ** S,int left,int right)310 void quickSort(Saddle **S, int left, int right)
311 {
312 if(left < right)
313 {
314 int pivot = (right + left) / 2;
315 int nuPivot = partition(S, left, right, pivot);
316 quickSort(S,left, nuPivot-1);
317 quickSort(S, nuPivot+1, right);
318 }
319 }
SortSaddles(Saddle ** S,int n)320 void SortSaddles(Saddle **S, int n)
321 {
322 quickSort(S, 0, n-1);
323 }
MS_Complex(double * points,int dimension,int count,int _k,bool perturb)324 MS_Complex::MS_Complex(double *points, int dimension, int count, int _k, bool perturb)
325 {
326 perturbed = perturb;
327 d = dimension-1;
328 numKneighbors = _k;
329 szV = numV = count;
330 //szE = szV*numKneighbors;
331 //Edges should be bi-directional
332 szE = szV*szV;
333
334 V = new Vertex *[szV];
335 E = new KNN_Edge *[szE];
336
337 srand(8);
338 double *p = new double[dimension];
339 for(int i=0; i < count; i++)
340 {
341 double eps = (double)rand() / (double)RAND_MAX;
342 eps = eps *1e-6;
343 if(rand() > RAND_MAX / 2)
344 eps = -eps;
345
346 for(int j = 0; j < dimension; j++)
347 p[j] = points[i*dimension+j];
348 p[dimension-1] += (perturbed ? eps : 0);
349 V[i] = new Vertex(d, p, p[d], i);
350 }
351 delete [] p;
352
353 // std::cout << "numV=" << numV << std::endl;
354 // std::cout << "d=" << d << std::endl;
355 // std::cout << "numKneighbors=" << numKneighbors << std::endl;
356
357 KNN();
358 Compute();
359 maxDist = -1;
360 }
361
MS_Complex(MS_Complex & Complex,double * new_point)362 MS_Complex::MS_Complex(MS_Complex &Complex, double *new_point)
363 {
364 perturbed = Complex.perturbed;
365 d = Complex.d;
366 numKneighbors = Complex.numKneighbors;
367 szV = numV = Complex.numV + 1;
368 //szE = szV*numKneighbors;
369 szE = szV*szV;
370
371 V = new Vertex *[szV];
372 E = new KNN_Edge *[szE];
373
374 double *p = new double[d+1];
375 for(int i=0; i < Complex.numV; i++)
376 {
377 for(int j = 0; j < d; j++)
378 p[j] = Complex.V[i]->GetXi(j);
379 p[d] = Complex.V[i]->Value();
380 V[i] = new Vertex(d, p, p[d], i);
381 }
382
383 double eps = (double)rand() / (double)RAND_MAX;
384 eps = eps * 1e-6;
385 if(rand() > RAND_MAX / 2)
386 eps = -eps;
387
388 V[numV-1] = new Vertex(d, new_point, new_point[d] + (perturbed ? eps : 0), numV-1);
389 delete [] p;
390 KNN();
391 Compute();
392 maxDist = -1;
393 }
394
KNN()395 void MS_Complex::KNN()
396 {
397 int eIndex = 0;
398 int i,k;
399
400 ANNpointArray pa = new ANNpoint[numV];
401 for(i = 0; i < numV; i++)
402 {
403 pa[i] = NULL;
404 V[i]->CreateANNpoint(pa[i]);
405 }
406
407 SearchStructure = new ANNkd_tree(pa,numV,d);
408 ANNidxArray nn_idx;
409 ANNdistArray dists;
410 for(i = 0; i < numV; i++)
411 {
412 nn_idx = new ANNidx[numKneighbors+1];
413 dists = new ANNdist[numKneighbors+1];
414
415 SearchStructure->annkSearch(pa[i],numKneighbors+1,nn_idx,dists);
416
417 for(k=1;k<numKneighbors+1;k++)
418 {
419 if(!DoesEdgeExist(i, nn_idx[k], eIndex))
420 {
421 E[eIndex] = new KNN_Edge(V[i],V[nn_idx[k]],E, eIndex);
422 eIndex++;
423 }
424 if(!DoesEdgeExist(nn_idx[k],i, eIndex))
425 {
426 E[eIndex] = new KNN_Edge(V[nn_idx[k]],V[i],E, eIndex);
427 eIndex++;
428 }
429 }
430 delete [] nn_idx;
431 delete [] dists;
432 }
433
434 for(i = 0; i < numV; i++)
435 delete [] pa[i];
436
437 delete [] pa;
438 numE = eIndex;
439 }
440
Compute()441 void MS_Complex::Compute()
442 {
443 int i,j;
444 int maxCount=0;
445 int minCount=0;
446 double globalMax = V[0]->Value();
447 double globalMin = V[0]->Value();
448 for(i = 0; i < numV; i++)
449 {
450 Vertex *v = V[i];
451
452 if(globalMax < v->Value())
453 globalMax = v->Value();
454 if(globalMin > v->Value())
455 globalMin = v->Value();
456
457 //Vertex **neighbors = new Vertex *[numKneighbors];
458 std::vector<Vertex *> neighbors;
459 //for(int j = 0; j < numKneighbors; j++)
460 j = 0;
461 while(v->GetIthNeighbor(j,E) != -1)
462 {
463 //neighbors[j] = V[v->GetIthNeighbor(j,E)];
464 neighbors.push_back(V[v->GetIthNeighbor(j,E)]);
465 j++;
466 }
467
468 double maximum = v->Value();
469 double minimum = v->Value();
470 Vertex *steepestA = v;
471 Vertex *steepestD = v;
472
473 //for(int j = 0; j < numKneighbors; j++)
474 j = 0;
475 while(v->GetIthNeighbor(j,E) != -1)
476 {
477 Vertex *currentNeighbor = neighbors[j];
478
479 if(currentNeighbor->Value() > maximum ||
480 (currentNeighbor->Value() == maximum &&
481 currentNeighbor->ID > v->ID))
482 {
483 maximum = currentNeighbor->Value();
484 steepestA = currentNeighbor;
485 }
486
487 if(currentNeighbor->Value() < minimum ||
488 (currentNeighbor->Value() == minimum &&
489 currentNeighbor->ID < v->ID))
490 {
491 minimum = currentNeighbor->Value();
492 steepestD = currentNeighbor;
493 }
494 j++;
495 }
496 if(steepestA == v)
497 {
498 v->classification = LOCAL_MAX;
499 maxCount++;
500 }
501 else if(steepestD == v)
502 {
503 v->classification = LOCAL_MIN;
504 minCount++;
505 }
506
507 v->Union_Max(steepestA, V);
508 v->Union_Min(steepestD, V);
509
510 //delete [] neighbors;
511 }
512
513 //Go through each vertex, examine its neighbors, if maxima are different, create or update the entry in a map
514 // where the key is the two maxima (always use lower maxima as first in pair), if the value of this point is
515 // higher then update or create the entry with this point as the pseudo-saddle.
516 //When we are done looking at each point, read off the saddles.
517 std::map<std::pair<int, int>,int> maxSaddles = std::map<std::pair<int, int>,int>();
518 for(i = 0; i < numV; i++)
519 {
520 if(V[i]->classification == LOCAL_MAX || V[i]->classification == LOCAL_MIN)
521 continue;
522
523 //Vertex **neighbors = new Vertex *[numKneighbors];
524 std::vector<Vertex *> neighbors;
525 j = 0;
526
527 //for(int j = 0; j < numKneighbors; j++)
528 while(V[i]->GetIthNeighbor(j,E) != -1)
529 {
530 //neighbors[j] = V[V[i]->GetIthNeighbor(j,E)];
531 neighbors.push_back(V[V[i]->GetIthNeighbor(j,E)]);
532 j++;
533 }
534 //for(int Bindex = 0; Bindex < numKneighbors; Bindex++)
535 for(int Bindex = 0; Bindex < j; Bindex++)
536 {
537 int AmaxIndex = V[i]->Find_Max(V);
538 int BmaxIndex = neighbors[Bindex]->Find_Max(V);
539
540 if(AmaxIndex != BmaxIndex)
541 {
542 if(AmaxIndex > BmaxIndex)
543 {
544 int temp = AmaxIndex;
545 AmaxIndex = BmaxIndex;
546 BmaxIndex = temp;
547 }
548 std::pair<int,int> AB(AmaxIndex, BmaxIndex);
549 if(maxSaddles.find(AB) == maxSaddles.end() || V[maxSaddles[AB]]->Value() < V[i]->Value())
550 {
551 maxSaddles[AB] = i;
552 }
553 }
554 }
555 //delete [] neighbors;
556 }
557 //Saddles between Minima
558 std::map<std::pair<int, int>, int> minSaddles = std::map<std::pair<int, int>, int>();
559 for(i = 0; i < numV; i++)
560 {
561 if(V[i]->classification == LOCAL_MAX || V[i]->classification == LOCAL_MIN)
562 continue;
563
564 //Vertex **neighbors = new Vertex *[numKneighbors];
565 std::vector<Vertex *> neighbors;
566 j = 0;
567 //for(int j = 0; j < numKneighbors; j++)
568 while(V[i]->GetIthNeighbor(j,E) != -1)
569 {
570 //neighbors[j] = V[V[i]->GetIthNeighbor(j,E)];
571 neighbors.push_back(V[V[i]->GetIthNeighbor(j,E)]);
572 j++;
573 }
574
575 int Bindex = 0;
576 //for(; Bindex < numKneighbors; Bindex++)
577 for(; Bindex < j; Bindex++)
578 {
579 int AminIndex = V[i]->Find_Min(V);
580 int BminIndex = neighbors[Bindex]->Find_Min(V);
581
582 if(AminIndex != BminIndex)
583 {
584 if(AminIndex > BminIndex)
585 {
586 int temp = AminIndex;
587 AminIndex = BminIndex;
588 BminIndex = temp;
589 }
590 std::pair<int,int> AB(AminIndex, BminIndex);
591 if(minSaddles.find(AB) == minSaddles.end() || V[minSaddles[AB]]->Value() > V[i]->Value())
592 {
593 minSaddles[AB] = i;
594 }
595 }
596 }
597 //delete [] neighbors;
598 }
599
600 ///////////////////////
601 double eps = (globalMax-globalMin)*1e-7;
602 int numSaddles = maxSaddles.size() + minSaddles.size();
603 Saddle * *S = new Saddle*[numSaddles];
604
605 std::map<std::pair<int,int>, int>::iterator sIter = maxSaddles.begin();
606 //Create Saddles
607 int curIdx = 0;
608 for( ; sIter != maxSaddles.end(); sIter++)
609 {
610 S[curIdx] = new Saddle();
611 S[curIdx]->idx = sIter->second;
612
613 //Higher index breaks tie and gets to be more persistent
614 if(V[sIter->first.first]->Value() < V[sIter->first.second]->Value() ||
615 (V[sIter->first.first]->Value() == V[sIter->first.second]->Value() && sIter->first.first < sIter->first.second))
616 {
617 S[curIdx]->less = sIter->first.first;
618 S[curIdx]->more = sIter->first.second;
619 }
620 else
621 {
622 S[curIdx]->less = sIter->first.second;
623 S[curIdx]->more = sIter->first.first;
624 }
625
626 //We don't want zero persistence, so apply some small epsilon (epsilon
627 // should depend on the range of values for this complex not arbitrarily
628 // chosen as I have done here)
629 S[curIdx]->persistence = std::max(eps,fabs(V[sIter->second]->Value() - V[S[curIdx]->less]->Value()));
630
631 if(V[S[curIdx]->less]->persistence == 0 || V[S[curIdx]->less]->persistence > S[curIdx]->persistence)
632 V[S[curIdx]->less]->persistence = S[curIdx]->persistence;
633
634 curIdx++;
635 V[sIter->second]->classification = SADDLE;
636 }
637
638 sIter = minSaddles.begin();
639 for( ; sIter != minSaddles.end(); sIter++)
640 {
641 S[curIdx] = new Saddle();
642 S[curIdx]->idx = sIter->second;
643
644 //Higher index breaks tie and gets to be more persistent
645 if(V[sIter->first.first]->Value() > V[sIter->first.second]->Value() ||
646 (V[sIter->first.first]->Value() > V[sIter->first.second]->Value() && sIter->first.first < sIter->first.second))
647 {
648 S[curIdx]->less = sIter->first.first;
649 S[curIdx]->more = sIter->first.second;
650 }
651 else
652 {
653 S[curIdx]->less = sIter->first.second;
654 S[curIdx]->more = sIter->first.first;
655 }
656
657 S[curIdx]->persistence = std::max(eps,fabs(V[sIter->second]->Value() - V[S[curIdx]->less]->Value()));
658
659 if(V[S[curIdx]->less]->persistence == 0 || V[S[curIdx]->less]->persistence > S[curIdx]->persistence)
660 V[S[curIdx]->less]->persistence = S[curIdx]->persistence;
661
662 curIdx++;
663 V[sIter->second]->classification = SADDLE;
664 }
665
666 //Do this at the end, where we actually look at the global min and max values
667 // as opposed to just reporting the extrema without persistence, this is wrong
668 // as an extrema's persistence may get updated when we cancel saddles
669 // This will incorrectly set an extremum's persistence and report the wrong
670 // global extrema.
671 // for(int i = 0; i < numSaddles; i++)
672 // {
673 // if(V[S[i]->more]->persistence == 0)
674 // {
675 // printf("Range: %f\n", globalMax-globalMin);
676 // V[S[i]->more]->persistence = globalMax-globalMin;
677 // if(V[S[i]->more]->classification == LOCAL_MAX)
678 // {
679 // globalMaxIdx = S[i]->more;
680 // }
681 // else
682 // {
683 // globalMinIdx = S[i]->more;
684 // }
685 // }
686 // }
687
688 //Sort Saddles
689 SortSaddles(S, numSaddles);
690 Saddle *head = S[0];
691 for(i = 0; i < numSaddles - 1; i++)
692 {
693 S[i]->next = S[i+1];
694 }
695 if(numSaddles > 0)
696 S[numSaddles - 1]->next = NULL;
697 else
698 head = S[0] = NULL;
699
700 std::map< std::pair<int, int>, std::vector<Vertex *> > Crystals;
701 int cCount = 0;
702 for(i = 0; i < numV; i++)
703 {
704 if(V[i]->classification == LOCAL_MIN || V[i]->classification == LOCAL_MAX)
705 continue;
706
707 int minI = V[i]->Find_Min(V);
708 int maxI = V[i]->Find_Max(V);
709 std::pair<int,int> minMaxPair = std::pair<int,int>(minI,maxI);
710 if(Crystals.find(minMaxPair) == Crystals.end())
711 {
712 Crystals[minMaxPair] = std::vector<Vertex *>();
713 Crystals[minMaxPair].push_back(V[i]);
714 }
715 else
716 {
717 Crystals[minMaxPair].push_back(V[i]);
718 }
719 }
720
721 std::map<std::pair<int,int>, std::vector<Vertex *> >::iterator cIter = Crystals.begin();
722 numC = 10000;
723 std::vector<double> tempPersistencesList;
724 tempPersistencesList.push_back(0.0);
725 std::vector<MS_Crystal *> tempCrystalList;
726 int nextCID = 0;
727 while(cIter != Crystals.end())
728 {
729 int *vIds = new int[cIter->second.size()+2];
730 for(j = 0; j < cIter->second.size(); j++)
731 {
732 vIds[j] = cIter->second[j]->ID;
733 }
734 vIds[j++] = cIter->first.first;
735 vIds[j] = cIter->first.second;
736
737 //C[nextCID] = new MS_Crystal(V,vIds, cIter->second.size()+2,d);
738 tempCrystalList.push_back(new MS_Crystal(V,vIds, cIter->second.size()+2,d));
739
740 cIter++;
741 nextCID++;
742 }
743 //for( ; nextCID < numC; nextCID++)
744 // C[nextCID] = NULL;
745 //nextCID = Crystals.size();
746
747 ///PERSISTENCE Calculation
748 p_levels = numSaddles;
749 V_to_C = new int[numV*p_levels];//Stores the crystal id that a vertex belongs to at a given persistence level
750 //Assign every vertex to its original Crystal for every persistence level, update when new crystal emerges
751 for(i = 0; i < Crystals.size(); i++)
752 {
753 for(j = 0; j < tempCrystalList[i]->size(); j++)
754 {
755 int vIndex = tempCrystalList[i]->GetVertexID(j);
756 for(int p = 0; p < p_levels; p++)
757 V_to_C[vIndex*p_levels+p] = i;
758 }
759 }
760
761 int p_level = 0;
762 Saddle *current = head;
763 while(current != NULL)
764 {
765 p_level++;
766 tempPersistencesList.push_back(V[current->less]->persistence);
767 int idxToRemove = current->less;
768 int idxToReplace = current->more;
769 V[idxToRemove]->persistence =
770 V[current->idx]->persistence = std::max(eps,fabs(V[current->less]->Value() - V[current->idx]->Value()));
771 //printf("Combining %d and %d\n", current->idx, idxToRemove);
772
773 //Find the rest of the saddles that have this extremum and delete them, then reenter them into the
774 // queue with the new neighboring extremum
775 Saddle *findIter = current->next;
776 Saddle *trailingFind = current;
777 while(findIter != NULL)
778 {
779 if(findIter->less == idxToRemove)
780 {
781 if(idxToReplace != findIter->more)
782 {
783 Saddle *nuSaddle = new Saddle();
784 nuSaddle->idx = findIter->idx;
785 if(V[idxToReplace]->classification == LOCAL_MAX)
786 {
787 if(V[idxToReplace]->Value() > V[findIter->more]->Value() ||
788 (V[idxToReplace]->Value() == V[findIter->more]->Value() && idxToReplace > findIter->more))
789 {
790 nuSaddle->more = idxToReplace;
791 nuSaddle->less = findIter->more;
792 }
793 else
794 {
795 nuSaddle->less = idxToReplace;
796 nuSaddle->more = findIter->more;
797 }
798 }
799 else
800 {
801 if(V[idxToReplace]->Value() < V[findIter->more]->Value() ||
802 (V[idxToReplace]->Value() == V[findIter->more]->Value() && idxToReplace > findIter->more))
803 {
804 nuSaddle->more = idxToReplace;
805 nuSaddle->less = findIter->more;
806 }
807 else
808 {
809 nuSaddle->less = idxToReplace;
810 nuSaddle->more = findIter->more;
811 }
812 }
813 nuSaddle->persistence = std::max(eps,fabs(V[nuSaddle->less]->Value() - V[nuSaddle->idx]->Value()));
814 V[nuSaddle->less]->persistence = std::min(nuSaddle->persistence, V[nuSaddle->less]->persistence);
815
816 //Remove the dead saddle
817 trailingFind->next = findIter->next;
818 delete findIter;
819 findIter = trailingFind->next;
820
821 //Insert the next one (note we will always gain persistence, so we can start
822 // looking to place it after the one we deleted)
823 Saddle *insert = findIter;
824 Saddle *trailingInsert = trailingFind;
825 if(insert == NULL || nuSaddle->persistence < insert->persistence)
826 {
827 trailingInsert->next = nuSaddle;
828 nuSaddle->next = insert;
829 trailingFind = nuSaddle;
830 }
831 else
832 {
833 while(insert != NULL && nuSaddle->persistence > insert->persistence)
834 {
835 trailingInsert = insert;
836 insert = insert->next;
837 }
838 nuSaddle->next = insert;
839 trailingInsert->next = nuSaddle;
840 }
841 }
842 else
843 {
844 //Remove the dead saddle
845 trailingFind->next = findIter->next;
846 delete findIter;
847 findIter = trailingFind->next;
848 }
849 }
850 else if(findIter->more == idxToRemove)
851 {
852 if(idxToReplace != findIter->less)
853 {
854 Saddle *nuSaddle = new Saddle();
855 nuSaddle->idx = findIter->idx;
856 if(V[idxToReplace]->classification == LOCAL_MAX)
857 {
858 if(V[idxToReplace]->Value() > V[findIter->less]->Value() ||
859 (V[idxToReplace]->Value() == V[findIter->less]->Value() && idxToReplace > findIter->less))
860 {
861 nuSaddle->more = idxToReplace;
862 nuSaddle->less = findIter->less;
863 }
864 else
865 {
866 nuSaddle->more = idxToReplace;
867 nuSaddle->less = findIter->less;
868 }
869 }
870 else
871 {
872 if(V[idxToReplace]->Value() < V[findIter->less]->Value() ||
873 (V[idxToReplace]->Value() == V[findIter->less]->Value() && idxToReplace > findIter->less))
874 {
875 nuSaddle->more = idxToReplace;
876 nuSaddle->less = findIter->less;
877 }
878 else
879 {
880 nuSaddle->more = idxToReplace;
881 nuSaddle->less = findIter->less;
882 }
883 }
884 nuSaddle->persistence = std::max(eps,fabs(V[nuSaddle->less]->Value() - V[nuSaddle->idx]->Value()));
885 V[nuSaddle->less]->persistence = std::min(nuSaddle->persistence, V[nuSaddle->less]->persistence);
886
887 //Remove the dead saddle
888 trailingFind->next = findIter->next;
889 delete findIter;
890 findIter = trailingFind->next;
891
892 //Insert the next one (note we will always gain persistence, so we can start
893 // looking to place it after the one we deleted)
894 Saddle *insert = findIter;
895 Saddle *trailingInsert = trailingFind;
896 if(insert == NULL || nuSaddle->persistence < insert->persistence)
897 {
898 trailingInsert->next = nuSaddle;
899 nuSaddle->next = insert;
900 trailingFind = nuSaddle;
901 }
902 else
903 {
904 while(insert != NULL && nuSaddle->persistence > insert->persistence)
905 {
906 trailingInsert = insert;
907 insert = insert->next;
908 }
909 nuSaddle->next = insert;
910 trailingInsert->next = nuSaddle;
911 }
912 }
913 else
914 {
915 //Remove the dead saddle
916 trailingFind->next = findIter->next;
917 delete findIter;
918 findIter = trailingFind->next;
919 }
920 }
921 else
922 {
923 trailingFind = findIter;
924 findIter = findIter->next;
925 }
926 }
927 //Do the crystal additions here
928 int totalCrystals = tempCrystalList.size();
929 for(i = 0; i < totalCrystals; i++)
930 {
931 if (tempCrystalList[i] == NULL || !tempCrystalList[i]->Exists(p_level-1))
932 continue;
933
934 bool success = false;
935 if(V[idxToRemove]->classification == LOCAL_MAX && tempCrystalList[i]->MaxV() == idxToRemove)
936 {
937 for(j = 0; j < totalCrystals; j++)
938 {
939 if (tempCrystalList[j] == NULL || !tempCrystalList[j]->Exists(p_level-1) || i == j)
940 continue;
941
942 if(tempCrystalList[j]->MaxV() == idxToReplace && tempCrystalList[i]->MinV() == tempCrystalList[j]->MinV())
943 {
944 //C[nextCID] = new MS_Crystal(C[i], C[j], true, p_level);
945 MS_Crystal *nuCrystal = new MS_Crystal(tempCrystalList[i], tempCrystalList[j], true, p_level);
946 tempCrystalList.push_back(nuCrystal);
947 for(int k = 0; k < nuCrystal->size(); k++)
948 {
949 int vIndex = nuCrystal->GetVertexID(k);
950 for(int p = p_level; p < p_levels; p++)
951 V_to_C[vIndex*p_levels+p] = nextCID;
952 }
953 nextCID++;
954 success = true;
955 break;
956 }
957 }
958 if(!success)
959 {
960 //C[nextCID] = new MS_Crystal(C[i],true,idxToReplace,p_level);
961 MS_Crystal *nuCrystal = new MS_Crystal(tempCrystalList[i],true,idxToReplace,p_level);
962 tempCrystalList.push_back(nuCrystal);
963 for(int k = 0; k < nuCrystal->size(); k++)
964 {
965 int vIndex = nuCrystal->GetVertexID(k);
966 for(int p = p_level; p < p_levels; p++)
967 V_to_C[vIndex*p_levels+p] = nextCID;
968 }
969 nextCID++;
970 }
971 }
972 else if(V[idxToRemove]->classification == LOCAL_MIN && tempCrystalList[i]->MinV() == idxToRemove)
973 {
974 for(j = 0; j < totalCrystals; j++)
975 {
976 if (tempCrystalList[j] == NULL || !tempCrystalList[j]->Exists(p_level-1) || i == j)
977 continue;
978
979 if(tempCrystalList[j]->MinV() == idxToReplace && tempCrystalList[i]->MaxV() == tempCrystalList[j]->MaxV())
980 {
981 //C[nextCID] = new MS_Crystal(C[i], C[j], false, p_level);
982 MS_Crystal *nuCrystal = new MS_Crystal(tempCrystalList[i], tempCrystalList[j], false, p_level);
983 tempCrystalList.push_back(nuCrystal);
984 for(int k = 0; k < nuCrystal->size(); k++)
985 {
986 int vIndex = nuCrystal->GetVertexID(k);
987 for(int p = p_level; p < p_levels; p++)
988 V_to_C[vIndex*p_levels+p] = nextCID;
989 }
990 nextCID++;
991 success = true;
992 break;
993 }
994 }
995 if(!success)
996 {
997 //C[nextCID] = new MS_Crystal(C[i],false,idxToReplace,p_level);
998 MS_Crystal *nuCrystal = new MS_Crystal(tempCrystalList[i],false,idxToReplace,p_level);
999 tempCrystalList.push_back(nuCrystal);
1000 for(int k = 0; k < nuCrystal->size(); k++)
1001 {
1002 int vIndex = nuCrystal->GetVertexID(k);
1003 for(int p = p_level; p < p_levels; p++)
1004 V_to_C[vIndex*p_levels+p] = nextCID;
1005 }
1006 nextCID++;
1007 }
1008 }
1009 }
1010 Saddle *trash = current;
1011 current = current->next;
1012 delete trash;
1013 }
1014 numC= nextCID;
1015 szP = tempPersistencesList.size();
1016 persistences = new double[szP];
1017 C = new MS_Crystal *[tempCrystalList.size()];
1018 std::vector<MS_Crystal *>::iterator mscIter = tempCrystalList.begin();
1019 for(int i = 0; i < tempCrystalList.size(); i++)
1020 {
1021 C[i] = tempCrystalList[i];
1022 }
1023 for(int i = 0; i < tempPersistencesList.size(); i++)
1024 {
1025 persistences[i] = tempPersistencesList[i];
1026 }
1027
1028 globalMinIdx = 0;
1029 globalMaxIdx = 0;
1030 for(i = 0; i < numV; i++)
1031 {
1032 if(V[i]->classification == LOCAL_MIN && V[i]->Value() < V[globalMinIdx]->Value())
1033 globalMinIdx = i;
1034 else if(V[i]->classification == LOCAL_MAX && V[i]->Value() > V[globalMaxIdx]->Value())
1035 globalMaxIdx = i;
1036 }
1037
1038 for(i = 0; i < numV; i++)
1039 {
1040 if((V[i]->classification == LOCAL_MIN || V[i]->classification == LOCAL_MAX) && V[i]->persistence == 0 )
1041 {
1042 if(i == globalMinIdx || i == globalMaxIdx)
1043 V[i]->persistence = V[globalMaxIdx]->Value() - V[globalMinIdx]->Value();
1044 else
1045 {
1046 std::cerr << "ERROR: There exists an extrema with zero persistence!\n\tVertex "
1047 << i << "\n\tClassification = " << V[i]->classification
1048 << "\n\tValue = " << V[i]->Value() << std::endl;
1049 // << "\nPrinting points for reconstruction to badMSC.txt" << std::endl;
1050 // FILE *badMSC = fopen("badMSC.txt", "w");
1051 // for(int i = 0; i < numV; i++)
1052 // {
1053 // fprintf(badMSC, "%f\t%f\t%f\n", V[i]->GetXi(0), V[i]->GetXi(1), V[i]->Value());
1054 // }
1055 // fclose(badMSC);
1056 }
1057 }
1058 }
1059
1060 delete [] S;
1061 }
1062
CompareExtremaCount(MS_Complex & _C)1063 double MS_Complex::CompareExtremaCount(MS_Complex &_C)
1064 {
1065 int maxNumV = _C.numV;
1066 int minNumV = numV;
1067
1068 int countThisMin = 0;
1069 int countThatMin = 0;
1070
1071 int countThisMax = 0;
1072 int countThatMax = 0;
1073
1074 int countThisSaddle = 0;
1075 int countThatSaddle = 0;
1076
1077 int _ci = 0;
1078 for(int i = 0; i < maxNumV; i++, _ci++)
1079 {
1080 if(i < numV)
1081 {
1082 if(V[i]->classification == 0)
1083 countThisMin++;
1084 else if(V[i]->classification == 1)
1085 countThisMax++;
1086 else if(V[i]->classification == 2)
1087 countThisSaddle++;
1088 }
1089
1090 if(i < _C.numV)
1091 {
1092 if(_C.V[i]->classification == 0)
1093 countThatMin++;
1094 else if(_C.V[i]->classification == 1)
1095 countThatMax++;
1096 else if(_C.V[i]->classification == 2)
1097 countThatSaddle++;
1098 }
1099
1100 }
1101
1102 return abs(countThisMin - countThatMin + countThisMax - countThatMax + countThisSaddle - countThatSaddle);
1103 }
1104
CompareClassChange(MS_Complex & _C)1105 double MS_Complex::CompareClassChange(MS_Complex &_C)
1106 {
1107 int maxNumV = _C.numV;
1108 int minNumV = numV;
1109
1110 int countChanged = 0;
1111
1112 int _ci = 0;
1113 for(int i = 0; i < minNumV; i++, _ci++)
1114 {
1115 if(V[i]->classification != _C.V[_ci]->classification)
1116 countChanged++;
1117 }
1118
1119 return countChanged;
1120 }
1121
ComparePersistence(MS_Complex & _C)1122 double MS_Complex::ComparePersistence(MS_Complex &_C)
1123 {
1124 int maxNumV = _C.numV;
1125 int minNumV = numV;
1126
1127 int countChanged = 0;
1128
1129 int _ci = 0;
1130 double sumP = 0.0;
1131 for(int i = 0; i < minNumV; i++, _ci++)
1132 {
1133 sumP += fabs(V[i]->persistence - _C.V[_ci]->persistence);
1134 }
1135
1136 return (sumP)/(double)(minNumV);
1137 }
1138
ComparePersistenceNoSaddles(MS_Complex & _C)1139 double MS_Complex::ComparePersistenceNoSaddles(MS_Complex &_C)
1140 {
1141 int maxNumV = _C.numV;
1142 int minNumV = numV;
1143
1144 int countChanged = 0;
1145
1146 double sumP = 0.0;
1147 for(int i = 0; i < minNumV; i++)
1148 {
1149 if(V[i]->classification == SADDLE || _C.V[i]->classification == SADDLE)
1150 continue;
1151 sumP += fabs(V[i]->persistence - _C.V[i]->persistence);
1152 }
1153
1154 return (sumP)/(double)(minNumV);
1155 }
1156
1157
CompareBottleneck(MS_Complex & _C)1158 double MS_Complex::CompareBottleneck(MS_Complex &_C)
1159 {
1160 PersistenceDiagram<double> pDia(1);
1161 for(int i = 0; i < numV; i++)
1162 {
1163 if(V[i]->classification == LOCAL_MIN && i != globalMinIdx)
1164 {
1165 double birth = V[i]->Value();
1166 double death = V[i]->Value() + V[i]->persistence;
1167 PDPoint<double> nu(birth, death);
1168 //std::cout << "MIN: " << nu << std::endl;
1169 pDia.push_back(nu);
1170 }
1171 else if(V[i]->classification == LOCAL_MAX && i != globalMaxIdx)
1172 {
1173 double birth = V[i]->Value() - V[i]->persistence;
1174 double death = V[i]->Value();
1175 PDPoint<double> nu(birth, death);
1176 //std::cout << "MAX: " << nu << std::endl;
1177 pDia.push_back(nu);
1178 }
1179 }
1180 PDPoint<double> globalMinMax(V[globalMinIdx]->Value(), V[globalMaxIdx]->Value());
1181 pDia.push_back(globalMinMax);
1182 PersistenceDiagram<double> pDia2(1);
1183 for(int i = 0; i < _C.numV; i++)
1184 {
1185 if(_C.V[i]->classification == LOCAL_MIN && i != _C.globalMinIdx)
1186 {
1187 double birth = _C.V[i]->Value();
1188 double death = _C.V[i]->Value() + _C.V[i]->persistence;
1189 PDPoint<double> nu(birth, death);
1190 //std::cout << "MIN: " << nu << std::endl;
1191 pDia2.push_back(nu);
1192 }
1193 else if(_C.V[i]->classification == LOCAL_MAX && i != _C.globalMaxIdx)
1194 {
1195 double birth = _C.V[i]->Value() - _C.V[i]->persistence;
1196 double death = _C.V[i]->Value();
1197 PDPoint<double> nu(birth, death);
1198 //std::cout << "MAX: " << nu << std::endl;
1199 pDia2.push_back(nu);
1200 }
1201 }
1202 PDPoint<double> globalMinMax2(_C.V[_C.globalMinIdx]->Value(), _C.V[_C.globalMaxIdx]->Value());
1203 pDia.push_back(globalMinMax2);
1204 //std::cout << "SIZES: " << pDia.size() << " " << pDia2.size() << std::endl;
1205 double ret_value = bottleneck_distance(pDia, pDia2);
1206 return ret_value;
1207 }
1208
CompareLabels(MS_Complex & _C,double p)1209 double MS_Complex::CompareLabels(MS_Complex &_C, double p)
1210 {
1211 int diff_labels = 0;
1212 double testPersistence = p*persistences[szP-1];
1213 int p_level = 0;
1214
1215 while(persistences[p_level] < testPersistence && p_level < szP-1)
1216 {
1217 p_level++;
1218 }
1219
1220 //double testPersistence2 = p*_C.persistences[_C.szP-1];
1221 int p_level2 = 0;
1222 while(_C.persistences[p_level2] < testPersistence && p_level2 < _C.szP-1)
1223 {
1224 p_level2++;
1225 }
1226
1227 for(int i = 0; i < numV && i < _C.numV; i++)
1228 {
1229 int idx1 = V_to_C[i*p_levels+p_level];
1230 int idx2 = _C.V_to_C[i*_C.p_levels+p_level2];
1231
1232 if(idx1 != idx2)
1233 diff_labels++;
1234 //if(idx1 < numC && idx2 < _C.numC)
1235 // if(C[idx1]->MaxV() != _C.C[idx2]->MaxV() || C[idx1]->MinV() != _C.C[idx2]->MinV())
1236 // diff_labels++;
1237 }
1238
1239 return diff_labels;
1240 }
1241
CloseToExtrema(double * v,double p)1242 double MS_Complex::CloseToExtrema(double *v, double p)
1243 {
1244 bool first = true;
1245 double minsDistance;
1246 double testPersistence = p*persistences[szP-1];
1247 for(int i = 0; i < numV; i++)
1248 {
1249 if(V[i]->classification == REGULAR || V[i]->persistence < testPersistence)
1250 continue;
1251
1252 double sDistance = 0.0;
1253 for(int j = 0; j < d; j++)
1254 sDistance += (V[i]->GetXi(j)-v[j])*(V[i]->GetXi(j)-v[j]);
1255
1256 if(first || minsDistance > sDistance)
1257 {
1258 first = false;
1259 minsDistance = sDistance;
1260 }
1261 }
1262 if(minsDistance == 0) //obviously don't want to sample the same point twice
1263 return 0;
1264 return 1/(minsDistance);
1265 }
1266
FarFromExtrema(double * v,double p)1267 double MS_Complex::FarFromExtrema(double *v, double p)
1268 {
1269 bool first = true;
1270 double minsDistance;
1271 double testPersistence = p*persistences[szP-1];
1272 int p_level = 0;
1273
1274 for(int i = 0; i < numV; i++)
1275 {
1276 if(V[i]->classification == REGULAR || V[i]->classification == SADDLE || V[i]->persistence < testPersistence)
1277 continue;
1278
1279 double sDistance = 0.0;
1280 for(int j = 0; j < d; j++)
1281 sDistance += (V[i]->GetXi(j)-v[j])*(V[i]->GetXi(j)-v[j]);
1282
1283 if(first || minsDistance > sDistance)
1284 {
1285 first = false;
1286 minsDistance = sDistance;
1287 }
1288
1289 }
1290 return minsDistance;
1291 }
1292
FarFromCP(double * v,double p)1293 double MS_Complex::FarFromCP(double *v, double p)
1294 {
1295 bool first = true;
1296 double minsDistance;
1297 double testPersistence = p*persistences[szP-1];
1298 int p_level = 0;
1299
1300 for(int i = 0; i < numV; i++)
1301 {
1302 if(V[i]->classification == REGULAR || V[i]->persistence < testPersistence)
1303 continue;
1304
1305 double sDistance = 0.0;
1306 for(int j = 0; j < d; j++)
1307 sDistance += (V[i]->GetXi(j)-v[j])*(V[i]->GetXi(j)-v[j]);
1308
1309 if(first || minsDistance > sDistance)
1310 {
1311 first = false;
1312 minsDistance = sDistance;
1313 }
1314
1315 }
1316 return minsDistance;
1317 }
1318
FarFromSaddle(double * v,double p)1319 double MS_Complex::FarFromSaddle(double *v, double p)
1320 {
1321 bool first = true;
1322 double minsDistance;
1323 double testPersistence = p*persistences[szP-1];
1324 int p_level = 0;
1325
1326 for(int i = 0; i < numV; i++)
1327 {
1328 if(V[i]->classification == REGULAR || V[i]->classification == LOCAL_MIN || V[i]->classification == LOCAL_MAX || V[i]->persistence < testPersistence)
1329 continue;
1330
1331 double sDistance = 0.0;
1332 for(int j = 0; j < d; j++)
1333 sDistance += (V[i]->GetXi(j)-v[j])*(V[i]->GetXi(j)-v[j]);
1334
1335 if(first || minsDistance > sDistance)
1336 {
1337 first = false;
1338 minsDistance = sDistance;
1339 }
1340
1341 }
1342 return minsDistance;
1343 }
1344
~MS_Complex()1345 MS_Complex::~MS_Complex()
1346 {
1347 for(int i = 0; i < numV; i++)
1348 delete V[i];
1349 delete [] V;
1350 for(int i = 0; i < numE; i++)
1351 delete E[i];
1352 delete [] E;
1353 for(int i = 0; i < numC; i++)
1354 delete C[i];
1355 delete [] C;
1356
1357 delete [] V_to_C;
1358 delete [] persistences;
1359 delete SearchStructure;
1360 }
1361
Destroy()1362 void MS_Complex::Destroy()
1363 {
1364 for(int i = 0; i < numV; i++)
1365 delete V[i];
1366 delete [] V;
1367 for(int i = 0; i < numE; i++)
1368 delete E[i];
1369 delete [] E;
1370 for(int i = 0; i < numC; i++)
1371 delete C[i];
1372 delete [] C;
1373
1374 delete [] V_to_C;
1375 delete [] persistences;
1376 delete SearchStructure;
1377 }
1378
GetIthHighestPersistence(int i)1379 int MS_Complex::GetIthHighestPersistence(int i)
1380 {
1381 int numCPs = 0;
1382 for(int j = 0; j < numV; j++)
1383 {
1384 if(V[j]->classification != REGULAR)
1385 numCPs++;
1386 }
1387 if(numCPs <= i)
1388 {
1389 return -1;
1390 }
1391 Saddle **AllExtrema = new Saddle *[numCPs];
1392
1393 int currentIdx = 0;
1394 for(int j = 0; j < numV && currentIdx < numCPs; j++)
1395 {
1396 if(V[j]->classification != REGULAR)
1397 {
1398 AllExtrema[currentIdx] = new Saddle();
1399 AllExtrema[currentIdx]->idx = j;
1400 AllExtrema[currentIdx]->persistence = V[j]->persistence;
1401 AllExtrema[currentIdx]->more = -1;
1402 AllExtrema[currentIdx]->less = -1;
1403 AllExtrema[currentIdx]->next = NULL;
1404 currentIdx++;
1405 }
1406 }
1407
1408 SortSaddles(AllExtrema, numCPs);
1409
1410 int retValue;
1411 if((numCPs-1) >= i && i >= 0)
1412 retValue = AllExtrema[(numCPs-1)-i]->idx;
1413 else
1414 retValue = -1;
1415
1416 delete [] AllExtrema;
1417 return retValue;
1418 }
1419
GetIthHighestPersistenceSaddle(int i)1420 int MS_Complex::GetIthHighestPersistenceSaddle(int i)
1421 {
1422 int numS = 0;
1423
1424 for(int j = 0; j < numV; j++)
1425 if(V[j]->classification == SADDLE)
1426 numS++;
1427
1428 if(numS <= i)
1429 {
1430 return -1;
1431 }
1432
1433 Saddle **AllExtrema = new Saddle *[numS];
1434
1435 int currentIdx = 0;
1436 for(int j = 0; j < numV && currentIdx < numS; j++)
1437 {
1438 if(V[j]->classification != REGULAR)
1439 {
1440 AllExtrema[currentIdx] = new Saddle();
1441 AllExtrema[currentIdx]->idx = j;
1442 AllExtrema[currentIdx]->persistence = V[j]->persistence;
1443 AllExtrema[currentIdx]->more = -1;
1444 AllExtrema[currentIdx]->less = -1;
1445 AllExtrema[currentIdx]->next = NULL;
1446 currentIdx++;
1447 }
1448 }
1449
1450 SortSaddles(AllExtrema, numS);
1451 return AllExtrema[numS-i]->idx;
1452 }
1453
IsCritical(int i)1454 bool MS_Complex::IsCritical(int i)
1455 {
1456 int j = i;
1457 return j < numV && j >= 0 && V[j]->classification != REGULAR;
1458 }
1459
IsExtrema(int i)1460 bool MS_Complex::IsExtrema(int i)
1461 {
1462 int j = i;
1463 return j < numV && j >= 0 && (V[j]->classification == LOCAL_MIN || V[j]->classification == LOCAL_MAX);
1464 }
1465
IsSaddle(int i)1466 bool MS_Complex::IsSaddle(int i)
1467 {
1468 int j = i;
1469 return j < numV && j >= 0 && (V[j]->classification == SADDLE);
1470 }
1471
GetVertex(int i)1472 Vertex * MS_Complex::GetVertex(int i)
1473 {
1474 int j = i;
1475 return (j < numV && j >= 0) ? V[j] : NULL;
1476 }
1477
ScoreTOPOB(MS_Complex & C1,double * x)1478 double ScoreTOPOB(MS_Complex &C1, double *x)
1479 {
1480 MS_Complex C2 = MS_Complex(C1,x);
1481 return C1.CompareBottleneck(C2);
1482 }
1483
ScoreTOPOP(MS_Complex & C1,double * x)1484 double ScoreTOPOP(MS_Complex &C1, double *x)
1485 {
1486 MS_Complex C2 = MS_Complex(C1,x);
1487 return C1.ComparePersistenceNoSaddles(C2);
1488 }
1489
ScoreTOPOHP(int dimension,int knn,double * training,double * trainingY,int n_training,double * candidates,double * candidateY,int n_candidates)1490 std::vector<int> ScoreTOPOHP(int dimension, int knn,
1491 double *training, double *trainingY, int n_training,
1492 double *candidates, double *candidateY, int n_candidates)
1493 {
1494 double *p = new double[(dimension+1)*(n_training+n_candidates)];
1495 for(int i = 0; i < n_training; i++)
1496 {
1497 for(int j = 0; j < dimension; j++)
1498 {
1499 p[i*(dimension+1)+j] = training[i*dimension+j];
1500 }
1501 p[i*(dimension+1)+dimension] = trainingY[i];
1502 }
1503 for(int i = 0; i < n_candidates; i++)
1504 {
1505 int ioffset = i+n_training;
1506 for(int j = 0; j < dimension; j++)
1507 {
1508 p[ioffset*(dimension+1)+j] = candidates[i*dimension+j];
1509 }
1510 p[ioffset*(dimension+1)+dimension] = candidateY[i];
1511 }
1512 MS_Complex Complex(p, dimension+1, n_training+n_candidates, knn);
1513 // std::cout << "DAN AMSC:\n";
1514 // Complex.Print(std::cout);
1515
1516 std::vector<int> returnOrder;
1517
1518 int count = 0;
1519 for(int i = 0; i < n_candidates; i++)
1520 {
1521 int tempIdx = Complex.GetIthHighestPersistence(count++) - n_training;
1522
1523 //If we don't find any CPs, then just sort them arbitrarily
1524 // for simplicity, I am returning points of zero persistence (regular)
1525 // in order they were given.
1526 if(tempIdx+n_training == -1)
1527 break;
1528
1529 //Probably overshooting on the iterations, but should not hurt to do so,
1530 // we won't add it if it is invalid anyway
1531 while(tempIdx < 0 && count < Complex.numV)
1532 tempIdx = Complex.GetIthHighestPersistence(count++) - n_training;
1533
1534 //Conditional indicates we actually found something useful above
1535 // and didn't just run out of rope
1536 if(count < Complex.numV)
1537 returnOrder.push_back(tempIdx);
1538 }
1539
1540 for(int i = 0; i < n_candidates; i++)
1541 {
1542 bool found = false;
1543 for(int j = 0; j < returnOrder.size(); j++)
1544 {
1545 if(i == returnOrder[j])
1546 {
1547 found = true;
1548 break;
1549 }
1550 }
1551 if(!found)
1552 returnOrder.push_back(i);
1553 }
1554
1555 delete p;
1556 return returnOrder;
1557 }
1558
Print(std::ostream & out)1559 void MS_Complex::Print(std::ostream &out)
1560 {
1561 out << "Minima: " << std::endl << "\t";
1562 for(int i = 0; i < numV; i++)
1563 {
1564 if(V[i]->classification == LOCAL_MIN)
1565 {
1566 for(int j = 0; j < d; j++)
1567 out << V[i]->GetXi(j) << ' ';
1568 out << V[i]->Value() << " (p=" << V[i]->persistence << ")" << std::endl << "\t";
1569 }
1570 }
1571 out << std::endl;
1572 out << "Maxima: " << std::endl;
1573 for(int i = 0; i < numV; i++)
1574 {
1575 if(V[i]->classification == LOCAL_MAX)
1576 {
1577 for(int j = 0; j < d; j++)
1578 out << V[i]->GetXi(j) << ' ';
1579 out << V[i]->Value() << " (p=" << V[i]->persistence << ")" << std::endl << "\t";
1580 }
1581 }
1582 out << std::endl;
1583 out << "Saddles: " << std::endl;
1584 for(int i = 0; i < numV; i++)
1585 {
1586 if(V[i]->classification == SADDLE)
1587 {
1588 for(int j = 0; j < d; j++)
1589 out << V[i]->GetXi(j) << ' ';
1590 out << V[i]->Value() << " (p=" << V[i]->persistence << ")" << std::endl << "\t";
1591 }
1592 }
1593 out << std::endl;
1594 }
1595
CountExtrema(double p)1596 int MS_Complex::CountExtrema(double p)
1597 {
1598 double testPersistence = p*persistences[szP-1];
1599 int count = 0;
1600
1601 for(int i = 0; i < numV; i++)
1602 {
1603 if((V[i]->classification == LOCAL_MAX || V[i]->classification == LOCAL_MIN) && V[i]->persistence >= testPersistence)
1604 count++;
1605 }
1606 return count;
1607 }
CountMaxima(double p)1608 int MS_Complex::CountMaxima(double p)
1609 {
1610 double testPersistence = p*persistences[szP-1];
1611 int count = 0;
1612
1613 for(int i = 0; i < numV; i++)
1614 {
1615 if(V[i]->classification == LOCAL_MAX && V[i]->persistence >= testPersistence)
1616 count++;
1617 }
1618 return count;
1619 }
1620
CountMinima(double p)1621 int MS_Complex::CountMinima(double p)
1622 {
1623 double testPersistence = p*persistences[szP-1];
1624 int count = 0;
1625
1626 for(int i = 0; i < numV; i++)
1627 {
1628 if(V[i]->classification == LOCAL_MIN && V[i]->persistence >= testPersistence)
1629 count++;
1630 }
1631 return count;
1632 }
CountSaddles(double p)1633 int MS_Complex::CountSaddles(double p)
1634 {
1635 double testPersistence = p*persistences[szP-1];
1636 int count = 0;
1637
1638 for(int i = 0; i < numV; i++)
1639 {
1640 if(V[i]->classification == SADDLE && V[i]->persistence >= testPersistence)
1641 count++;
1642 }
1643 return count;
1644 }
1645
1646 #ifdef USING_GL
Draw(double gMin,double gMax,bool flatMode)1647 void MS_Complex::Draw(double gMin, double gMax,bool flatMode)
1648 {
1649 if(d != 2)
1650 return;
1651
1652 glLineWidth(1.0);
1653 glBegin(GL_LINES);
1654 for(int i = 0; i < numE; i++)
1655 {
1656 glColor3f(0.25,0.25,0.25);
1657 glVertex3f(V[E[i]->start]->GetXi(0),V[E[i]->start]->GetXi(1), flatMode ? 0 : ((V[E[i]->start]->Value()-gMin)/(gMax-gMin) - 1./2.));
1658 glVertex3f(V[E[i]->end]->GetXi(0), V[E[i]->end]->GetXi(1), flatMode ? 0 : ((V[E[i]->end]->Value()-gMin)/(gMax-gMin) - 1./2.));
1659 }
1660 glEnd();
1661
1662 for(int i = 0; i < numV; i++)
1663 {
1664 Vertex *curV = V[i];
1665 if(curV->classification == 0)
1666 glColor3f(0,0.5,1);
1667 else if(curV->classification == 1)
1668 glColor3f(1,0,0);
1669 else if(curV->classification == 2)
1670 glColor3f(0,1,0);
1671 else
1672 glColor3f(0,0,0);
1673
1674 glTranslatef(curV->GetXi(0),curV->GetXi(1), flatMode ? 0 : ((curV->Value()-gMin)/(gMax-gMin) - 1./2.));
1675 glutSolidSphere(curV->classification == 3 ? 0.005 : 0.01,16,16);
1676 glTranslatef(-curV->GetXi(0),-curV->GetXi(1), -(flatMode ? 0 : ((curV->Value()-gMin)/(gMax-gMin) - 1./2.)));
1677 }
1678 for(int i = 0; i < numV; i++)
1679 {
1680 Vertex *curV = V[i];
1681 glDisable(GL_LIGHTING);
1682 glEnable(GL_BLEND);
1683 glColor3f(0,0,0);
1684 glLineWidth(6.0);
1685 for(int k = 0; k < numKneighbors; k++)
1686 {
1687 int nextIdx = curV->GetIthNeighbor(k, E);
1688 Vertex *nextV = V[nextIdx];
1689 if(nextIdx == curV->NeighborMax() || nextIdx == curV->NeighborMin())
1690 {
1691 glBegin(GL_LINES);
1692 if(i == nextV->NeighborMax())
1693 glColor4f(1,0,0,1);
1694 else if(i == nextV->NeighborMin())
1695 glColor4f(0.8,1.,1,1);
1696 else
1697 glColor4f(0,0,0,0);
1698
1699 glVertex3f(curV->GetXi(0),curV->GetXi(1), flatMode ? 0 : ((curV->Value()-gMin)/(gMax-gMin) - 1./2.));
1700 if(nextIdx == curV->NeighborMax())
1701 glColor4f(1,0,0,1);
1702 else if(nextIdx == curV->NeighborMin())
1703 glColor4f(0.8,1.,1,1);
1704 else
1705 glColor4f(0,0,0,0);
1706
1707 glVertex3f(nextV->GetXi(0),nextV->GetXi(1), flatMode ? 0 : ((nextV->Value()-gMin)/(gMax-gMin) - 1./2.));
1708 glEnd();
1709 }
1710 }
1711 glEnable(GL_LIGHTING);
1712 }
1713 glLineWidth(1.0);
1714 glDisable(GL_BLEND);
1715 }
1716
DrawBurst()1717 void MS_Complex::DrawBurst()
1718 {
1719 double maxSdist = 0;
1720 if(maxDist < 0)
1721 {
1722 for(int vi = 0; vi < numV; vi++)
1723 {
1724 for(int vj = vi+1; vj < numV; vj++)
1725 {
1726 double sdist = V[vi]->SDistance(V[vj]);
1727 if(sdist > maxSdist)
1728 maxSdist = sdist;
1729 }
1730 }
1731 maxDist = sqrt(maxSdist);
1732 }
1733
1734 double factor = M_PI/(2.*maxDist);
1735
1736 double gMin = V[globalMinIdx]->Value();
1737 double gMax = V[globalMaxIdx]->Value();
1738
1739 std::vector<int> mins;
1740 std::vector<int> maxs;
1741 std::vector<int> saddles;
1742 for(int i = 0; i < numV; i++)
1743 {
1744 if(V[i]->classification == LOCAL_MIN)
1745 mins.push_back(i);
1746 if(V[i]->classification == LOCAL_MAX)
1747 maxs.push_back(i);
1748 if(V[i]->classification == SADDLE)
1749 saddles.push_back(i);
1750 }
1751
1752 double x = 0;
1753 double y = 0;
1754 for(int i = 0; i < mins.size(); i++)
1755 {
1756 int midx = mins[i];
1757
1758 x = (double)(i+1)/(double)(mins.size()+1);
1759 y = 0.5*((V[midx]->Value() - gMin)/(gMax-gMin)) - 0.5;
1760
1761 glColor3f(0,0,1);
1762 glTranslatef( x, 0.5, y);
1763 glutSolidSphere(0.005,16,16);
1764 glTranslatef(-x, -0.5, -y);
1765
1766 glColor3f(0,0,0);
1767 for(int j = 0; j < numV; j++)
1768 {
1769 if(V[j]->Find_Min(V) == midx && midx != j)
1770 {
1771 glBegin(GL_LINES);
1772 glVertex3f(x,0.5,y);
1773 double radius = 0.5*(V[j]->Value() - V[midx]->Value())/(gMax-gMin);
1774 double dist = sqrt(V[midx]->SDistance(V[j]));
1775 double theta = dist*factor + M_PI/2.;
1776 double xtemp = x + radius*cos(theta);
1777 double ytemp = y + radius*sin(theta);
1778 glVertex3f(xtemp, 0.5, ytemp);
1779 glEnd();
1780 glTranslatef(xtemp, 0.5, ytemp);
1781 glutSolidSphere(0.0025,16,16);
1782 glTranslatef(-xtemp, -0.5, -ytemp);
1783 }
1784 }
1785 }
1786 for(int i = 0; i < maxs.size(); i++)
1787 {
1788 int midx = maxs[i];
1789
1790 x = (double)(i+1)/(double)(maxs.size()+1);
1791 y = 0.5*((V[midx]->Value() - gMin)/(gMax-gMin));
1792
1793 glColor3f(1,0,0);
1794 glTranslatef( x, 0.5, y);
1795 glutSolidSphere(0.005,16,16);
1796 glTranslatef(-x, -0.5, -y);
1797
1798 glColor3f(0,0,0);
1799 for(int j = 0; j < numV; j++)
1800 {
1801 if(V[j]->Find_Max(V) == midx && midx != j)
1802 {
1803 glBegin(GL_LINES);
1804 glVertex3f(x,0.5,y);
1805 double radius = 0.5*(V[midx]->Value() - V[j]->Value())/(gMax-gMin);
1806 double dist = sqrt(V[midx]->SDistance(V[j]));
1807 double theta = dist*factor + M_PI/2.;
1808 double xtemp = x + radius*cos(theta);
1809 double ytemp = y - radius*sin(theta);
1810 glVertex3f(xtemp, 0.5, ytemp);
1811 glEnd();
1812 glTranslatef(xtemp, 0.5, ytemp);
1813 glutSolidSphere(0.0025,16,16);
1814 glTranslatef(-xtemp, -0.5, -ytemp);
1815 }
1816 }
1817 }
1818
1819 }
1820 #endif
1821