1 /******************************************************************************/
2 /* File: HD.cpp */
3 /* Created by: Rainer Dyckerhoff, Pavlo Mozharovskyi */
4 /* First published: 19.06.2015 */
5 /* Last revised: 03.07.2015 */
6 /* */
7 /* Contains functions that compute the exact halfspace depth of a point z */
8 /* w.r.t. n data points x[1],...x[n]. */
9 /* */
10 /******************************************************************************/
11
12 #include "stdafx.h"
13
14 #define _USE_MATH_DEFINES
15 #include <math.h>
16 #include <algorithm>
17
18
19 using namespace std;
20
21 /* Definition of constants */
22 const double eps_HD1 = 1e-8;
23 const double eps_HD2 = 1e-8;
24 const double eps_HDx = 1e-8;
25 const double eps_Cmb1 = 1e-8;
26 //const double eps_pivot = 1e-8;
27
28
29 /****************************************************************************/
30 /* */
31 /* 'norm2' computes the Euclidean norm of a vector x in d-space. */
32 /* */
33 /****************************************************************************/
34
norm2(double * x,int d)35 double norm2(double* x, int d) {
36 double result = 0;
37 for (int i = 0; i < d; i++) result += x[i] * x[i];
38 return sqrt(result);
39 }
40
41 /****************************************************************************/
42 /* */
43 /* 'intHD1' computes the integer hafspace depth of 0 w.r.t n data points */
44 /* in R^1. */
45 /* */
46 /****************************************************************************/
47
intHD1(double ** x,int n)48 int intHD1(double** x, int n) {
49 int cnt1 = 0, cnt2 = 0;
50 for (int i = 0; i < n; i++, x++) {
51 if (**x < eps_HD1) cnt1++;
52 if (**x > -eps_HD1) cnt2++;
53 }
54 return min(cnt1, cnt2);
55 }
56
57 /****************************************************************************/
58 /* */
59 /* 'intHD2' computes the integer hafspace depth of 0 w.r.t n data points */
60 /* in R^2. */
61 /* */
62 /* This is an implemetation of the algorithm of */
63 /* Rousseeuw, P.J.and Ruts, I. (1996). Algorithm AS 307: bivariate */
64 /* location depth. Journal of the Royal Statistical Society. Series C: */
65 /* Applied Statistics 45, 516?526. */
66 /* */
67 /****************************************************************************/
68
intHD2(double ** x,int n)69 int intHD2(double** x, int n) {
70 double* alpha = new double[n];
71 int nt = 0; // how many zeros in in x ?
72 int nh = 0; // how many points in the halfspace ?
73 // compute angles alpha[i] and initilize array angle
74 for (int i = 0; i < n; i++) {
75 if (hypot(x[i][0], x[i][1]) <= eps_HD2)
76 nt++;
77 else {
78 alpha[i - nt] = atan2(x[i][1], x[i][0]); // alpha in (-pi,pi]
79 // correction for points like (-1, -1e-16)
80 if (alpha[i - nt] < -M_PI + eps_HD2) alpha[i - nt] = M_PI;
81 if (alpha[i - nt] <= eps_HD2) nh++;
82 }
83 }
84 int nn = n - nt;
85 // sort angles
86 sort(alpha, alpha + nn);
87 // compute halfspace depth
88 int result = nh;
89 if (result > 0) {
90 int j = nh;
91 for (int i = 0; i < nh; i++) {
92 while ((j <= nn - 1) && (alpha[j] - M_PI <= alpha[i] + eps_HD2))
93 j++;
94 if (j - i <= result) result = j - i - 1;
95 }
96 j = 0;
97 for (int i = nh; i < nn; i++) {
98 while ((j <= nh - 1) && (alpha[j] + M_PI <= alpha[i] + eps_HD2))
99 j++;
100 if (j - (i - nn) <= result) result = j - (i - nn) - 1;
101 }
102 }
103 delete[] alpha;
104 return result + nt;
105 }
106
107 /****************************************************************************/
108 /* */
109 /* 'nHD_Rec' computes the integer halfspace depth of 0 w.r.t n data points */
110 /* in R^d. */
111 /* */
112 /* 'nHD_Rec' implements the recursive algorithm (k = 1) as described in */
113 /* Section 3.3 of "Exact computation of the halfspace depth" by */
114 /* Rainer Dyckerhoff and Pavlo Mozharovskyi (arXiv:1411:6927) */
115 /* */
116 /****************************************************************************/
117
nHD_Rec(double ** xx,int n,int d)118 int nHD_Rec(double** xx, int n, int d) {
119 if (d == 1) return intHD1(xx, n);
120 if (d == 2) return intHD2(xx, n);
121
122 double* y = new double[d - 1];
123 double* z = new double[d];
124 double** x = new double*[n];
125 for (int k = 0; k < n; k++) x[k] = new double[d - 1];
126
127 int result = n;
128 for (int i = 0; i < n; i++) {
129
130 int kmax = d;
131 double xmax = 0;
132 for (int k = 0; k < d; k++)
133 if (fabs(xx[i][k]) > xmax) {
134 xmax = fabs(xx[i][k]);
135 kmax = k;
136 }
137 if (xmax > eps_HDx) {
138 int nNull = 0, nPos = 0, nNeg = 0, m = 0;
139 for (int k = 0; k < d; k++) z[k] = xx[i][k] / xx[i][kmax];
140 // project data points
141 for (int j = 0; j < n; j++){
142 double alpha = -xx[j][kmax];
143 for (int k = 0; k < kmax; k++)
144 y[k] = xx[j][k] + alpha * z[k];
145 for (int k = kmax; k < d - 1; k++)
146 y[k] = xx[j][k + 1] + alpha * z[k + 1];
147 if (norm2(y, d - 1) > eps_HDx) {
148 for (int k = 0; k < d - 1; k++) x[m][k] = y[k];
149 m++;
150 }
151 else {
152 // in this case alpha = -sign(x_i*x_j)
153 if (alpha > eps_HDx) nPos++;
154 else if (alpha < -eps_HDx) nNeg++;
155 else nNull++;
156 }
157 }
158 result = min(result, nHD_Rec(x, m, d - 1) +
159 nNull + min(nPos, nNeg));
160 if (result == 0) break;
161 }
162 }
163 for (int k = 0; k < n; k++) delete[] x[k];
164 delete[] x;
165 delete[] z;
166 delete[] y;
167 return result;
168 }
169
170 /****************************************************************************/
171 /* */
172 /* 'HD_Rec' computes the halfspace depth of a point z w.r.t n data */
173 /* points in R^d. */
174 /* */
175 /* 'HD_Rec' does some preprocessing of the data and then calls */
176 /* 'nHD_Rec'. */
177 /* */
178 /* See also the description in the header file. */
179 /* */
180 /* HD_Rec calls the following routines: */
181 /* norm2 */
182 /* intHD1 */
183 /* intHD2 */
184 /* nHD_Rec */
185 /* */
186 /****************************************************************************/
187
HD_Rec(double * z,double ** xx,int n,int d)188 double HD_Rec(double* z, double** xx, int n, int d) {
189 if (n <= 0) throw invalid_argument("n <= 0");
190 if (d <= 0) throw invalid_argument("d <= 0");
191 // preprocess data
192 // subtract z from all data points x[i]
193 int m = 0;
194 double** x = new double*[n];
195 bool create = true;
196 for (int i = 0; i < n; i++) {
197 if (create)
198 x[m] = new double[d];
199
200 for (int j = 0; j < d; j++) x[m][j] = xx[i][j] - z[j];
201 create = norm2(x[m], d) >= eps_HDx;
202 if (create) m++;
203 }
204 int result = nHD_Rec(x, m, d) + (n - m);
205 // deallocate array x
206 if (!create) m++;
207 for (int i = 0; i < m; i++) delete[] x[i];
208 delete[] x;
209 return result / (double)n;
210 }
211
212 /****************************************************************************/
213 /* 'getRank' computes the rank of the matrix x */
214 /* */
215 /* 'getRank' is used in preprocessing the data in HD_comb and HD_Comb2. */
216 /* 'getRank' detects if the data points are contained in a lower */
217 /* dimensional space, by computing the rank of the matrix formed by the */
218 /* data points x[i]. */
219 /* */
220 /****************************************************************************/
221
getRank(double ** x,int n,int d,int * piv)222 int getRank(double** x, int n, int d, int* piv) {
223 int imax;
224 int pc = 0, rank = 0;
225 double amax;
226 // copy x to A
227 double** A = new double*[d];
228 for (int i = 0; i < d; i++) {
229 A[i] = new double[n];
230 for (int j = 0; j < n; j++) A[i][j] = x[j][i];
231 }
232 rank = 0;
233 for (int k = 0; k < min(n, d); k++) {
234 // k-th elimination step
235 do {
236 imax = k;
237 amax = fabs(A[k][pc]);
238 // find maximum element in column
239 for (int i = k + 1; i < d; i++) {
240 if (fabs(A[i][pc]) > amax) {
241 amax = fabs(A[i][pc]);
242 imax = i;
243 }
244 }
245 if (amax < eps_pivot) pc++;
246 } while ((amax < eps_pivot) && (pc < n));
247 if (pc < n) {
248 rank++;
249 piv[k] = pc;
250 // exchange rows
251 if (imax != k) {
252 for (int j = pc; j < n; j++) {
253 double tmp = A[k][j];
254 A[k][j] = A[imax][j];
255 A[imax][j] = tmp;
256 }
257 }
258 // elimination
259 for (int i = k + 1; i < d; i++) {
260 double factor = A[i][pc] / A[k][pc];
261 for (int j = pc + 1; j < n; j++)
262 A[i][j] -= factor * A[k][j];
263 }
264 if (++pc >= n) break;
265 }
266 else break;
267 }
268 for (int i = 0; i < d; i++) delete[] A[i];
269 delete[] A;
270
271 return rank;
272 }
273
274 /****************************************************************************/
275 /* 'project' projects the data points on a lower dimensional subspace. */
276 /* */
277 /* 'project' is used in preprocessing the data in HD_comb and HD_Comb2. */
278 /* If the data points x[i] are contained in a subspace of dimension 'rank' */
279 /* (as detected by a call to getRank), the representation of the data */
280 /* points w.r.t. a basis of this subspace is computed. This gives a */
281 /* representation of the data points in the Euclidean space of dimension */
282 /* rank. */
283 /* */
284 /****************************************************************************/
285
project(double ** x,int n,int d,int rank,int indices[])286 void project(double** x, int n, int d, int rank, int indices[]) {
287 double** z = new double*[n];
288 for (int k = 0; k < n; k++) {
289 z[k] = new double[rank];
290 for (int i = 0; i < rank; i++) {
291 z[k][i] = 0;
292 for (int l = 0; l < d; l++)
293 z[k][i] += x[k][l] * x[indices[i]][l];
294 }
295 }
296 for (int k = 0; k < n; k++) {
297 delete[] x[k];
298 x[k] = z[k];
299 }
300 delete[] z;
301 }
302
303 /****************************************************************************/
304 /* */
305 /* 'getNormal' computes the normal vector to the d-1 vectors passed */
306 /* in A. */
307 /* */
308 /* If the rank of A is equal to d-1, then the function returns 'true' and */
309 /* the normal vector is passed to the calling routine in 'normal[]'. */
310 /* If the rank of A is less than d-1, then the function returns 'false' */
311 /* and the value of 'normal[]' is undefined. */
312 /* */
313 /****************************************************************************/
314
getNormal(double ** A,int d,double * normal)315 bool getNormal(double** A, int d, double* normal) {
316 int imax, jmax;
317 int* colp = new int[d];
318 double amax;
319 for (int k = 0; k < d - 1; k++) {
320 imax = k;
321 jmax = k;
322 amax = fabs(A[k][k]);
323 colp[k] = k;
324 // find maximum element in column
325 for (int i = k + 1; i < d - 1; i++) {
326 if (fabs(A[i][k]) > amax) {
327 amax = fabs(A[i][k]);
328 imax = i;
329 }
330 }
331 // maximum equal to zero => complete pivoting
332 if (amax < eps_pivot) {
333 for (int j = k + 1; j < d; j++) {
334 for (int i = k; i < d - 1; i++) {
335 if (fabs(A[i][j]) > amax) {
336 amax = fabs(A[i][j]);
337 imax = i;
338 jmax = j;
339 }
340 }
341 }
342 if (amax < eps_pivot) {
343 delete[] colp;
344 return false;
345 }
346 // exchange columns
347 for (int i = 0; i < d - 1; i++) {
348 double tmp = A[i][k];
349 A[i][k] = A[i][jmax];
350 A[i][jmax] = tmp;
351 }
352 colp[k] = jmax;
353 }
354 // exchange rows
355 if (imax != k) {
356 for (int j = k; j < d; j++) {
357 double tmp = A[k][j];
358 A[k][j] = A[imax][j];
359 A[imax][j] = tmp;
360 }
361 }
362 // elimination
363 for (int i = k + 1; i < d - 1; i++) {
364 double factor = A[i][k] / A[k][k];
365 for (int j = k + 1; j < d; j++) A[i][j] -= factor * A[k][j];
366 }
367 }
368 // back substitution
369 colp[d - 1] = d - 1;
370 normal[d - 1] = -1;
371 for (int k = d - 2; k >= 0; k--) {
372 normal[k] = A[k][d - 1] / A[k][k];
373 for (int i = k - 1; i >= 0; i--) A[i][d - 1] -= normal[k] * A[i][k];
374 }
375 // reverse column permutations
376 for (int k = d - 1; k >= 0; k--) {
377 if (colp[k] != k) {
378 double temp = normal[k];
379 normal[k] = normal[colp[k]];
380 normal[colp[k]] = temp;
381 }
382 }
383 delete[] colp;
384
385 return true;
386 }
387
388 int nHD_Comb(double** xx, int n, int d);
389
390 /****************************************************************************/
391 /* */
392 /* 'HD1proj' performs the following steps: */
393 /* */
394 /* 1) All data points x[i] are projected in the direction p, */
395 /* i.e., z[i] = p'x[i] is computed. */
396 /* 2) The univariate integer halfspace depth of 0 is computed w.r.t. all */
397 /* the z[i] that are not equal to 0. */
398 /* 3) If there are more than d-1 values z[i] that are equal to zero, */
399 /* the respective points are projected on the orthogonal complement */
400 /* of p. Then, the integer halfspace depth of 0 w.r.t. these */
401 /* projected points is computed. */
402 /* 4) The sum of the values from step 2) and 3) is returned. */
403 /* */
404 /****************************************************************************/
405
HD1proj(double ** x,int n,int d,double * p,int indices[])406 int HD1proj(double** x, int n, int d, double* p, int indices[]) {
407 int cnt0 = 0, cnt1 = 0, cnt2 = 0, HDproj = 0;
408 int* plane = new int[n];
409 for (int i = 0; i < n; i++) {
410 double sum = 0;
411 for (int j = 0; j < d; j++) sum += p[j] * x[i][j];
412 if (sum > eps_HD1) cnt1++;
413 else if (sum < -eps_HD1) cnt2++;
414 else plane[cnt0++] = i;
415 }
416 if (cnt0 > d - 1) {
417 // recursion
418 double** z = new double*[cnt0];
419 for (int i = 0; i
420 < cnt0; i++) {
421 z[i] = new double[d - 1];
422 for (int j = 0; j < d - 1; j++) {
423 z[i][j] = 0;
424 for (int k = 0; k < d; k++)
425 z[i][j] += x[indices[j]][k] * x[plane[i]][k];
426 }
427 }
428 HDproj = nHD_Comb(z, cnt0, d - 1);
429 for (int i = 0; i < cnt0; i++) delete[] z[i];
430 delete[] z;
431 }
432 delete[] plane;
433 return min(cnt1, cnt2) + HDproj;
434 }
435
436 /****************************************************************************/
437 /* */
438 /* 'nHD_Comb' computes the integer halfspace depth of 0 w.r.t n data points */
439 /* in R^d. */
440 /* */
441 /* 'nHD_Comb' implements the combinatorial algorithm (k = d-1) as described */
442 /* in Section 3.1 of "Exact computation of the halfspace depth" by */
443 /* Rainer Dyckerhoff and Pavlo Mozharovskyi (arXiv:1411:6927) */
444 /* */
445 /****************************************************************************/
446
nHD_Comb(double ** xx,int n,int d)447 int nHD_Comb(double** xx, int n, int d) {
448 if (d == 1) return intHD1(xx, n);
449 if (d == 2) return intHD2(xx, n);
450
451 int result = n + 1;
452 double** a = new double*[d - 1];
453 for (int i = 0; i < d - 1; i++) a[i] = new double[d];
454 double* p = new double[d];
455 int* indices = new int[d - 1];
456 indices[0] = -1;
457 int pos = 0;
458 while (pos >= 0) {
459 indices[pos]++;
460 for (pos++; pos < d - 1; pos++) indices[pos] = indices[pos - 1] + 1;
461 pos--;
462 do {
463 for (int i = 0; i < d - 1; i++)
464 for (int j = 0; j < d; j++) a[i][j] = xx[indices[i]][j];
465 if (getNormal(a, d, p))
466 result = min(result, HD1proj(xx, n, d, p, indices));
467 indices[pos]++;
468 } while (indices[pos] < n - d + pos + 2);
469 do pos--; while (pos >= 0 && indices[pos] >= n - d + pos + 1);
470 }
471 for (int i = 0; i < d - 1; i++) delete[] a[i];
472 delete[] a;
473 delete[] p;
474 delete[] indices;
475 return result;
476 }
477
478 /****************************************************************************/
479 /* */
480 /* 'HD_Comb' computes the halfspace depth of a point z w.r.t n data */
481 /* points in R^d. */
482 /* */
483 /* 'HD_Comb' does some preprocessing of the data and then calls */
484 /* 'nHD_Comb'. */
485 /* */
486 /* See also the description in the header file. */
487 /* */
488 /* HD_Comb calls the following routines: */
489 /* norm2 */
490 /* intHD1 */
491 /* intHD2 */
492 /* getRank */
493 /* project */
494 /* getNormal */
495 /* HD1proj */
496 /* nHD_Rec */
497 /* */
498 /****************************************************************************/
499
HD_Comb(double * z,double ** xx,int n,int d)500 double HD_Comb(double* z, double** xx, int n, int d) {
501 if (n <= 0) throw invalid_argument("n <= 0");
502 if (d <= 0) throw invalid_argument("d <= 0");
503 // preprocess data
504 // subtract z from all data points x[i]
505 // check whether the data points are concentrated on a lower
506 // dimensional space
507 int m = 0, rank;
508 int* indices = new int[d];
509 double** x = new double*[n];
510 for (int i = 0; i < n; i++) {
511 x[m] = new double[d];
512 for (int j = 0; j < d; j++) x[m][j] = xx[i][j] - z[j];
513 if (norm2(x[m], d) >= eps_HDx) m++; else delete[] x[m];
514 }
515 if (m == 0) return 1.0;
516
517 rank = getRank(x, m, d, indices);
518 if (rank < d) project(x, m, d, rank, indices);
519 int result = nHD_Comb(x, m, rank) + (n - m);
520 // deallocate array x
521 for (int i = 0; i < m; i++) delete[] x[i];
522 delete[] x;
523 delete[] indices;
524 return result / (double)n;
525 }
526
527 /****************************************************************************/
528 /* */
529 /* 'getBasisComplement' computes a basis of the orthogonal complement of */
530 /* the d-2 vectors passed in A. */
531 /* */
532 /* If the rank of A is equal to d-2, then the function returns 'true' and */
533 /* the two basis vectors are passed to the calling routine in 'basis[][]'. */
534 /* If the rank of A is less than d-2, then the function returns 'false' */
535 /* and the value of 'basis[]' is undefined. */
536 /* */
537 /****************************************************************************/
538
getBasisComplement(double ** A,int d,double ** basis)539 bool getBasisComplement(double** A, int d, double** basis) {
540 int imax, jmax;
541 int* colp = new int[d];
542 double amax;
543 for (int k = 0; k < d - 2; k++) {
544 imax = k;
545 jmax = k;
546 amax = fabs(A[k][k]);
547 colp[k] = k;
548 // find maximum element in column
549 for (int i = k + 1; i < d - 2; i++) {
550 if (fabs(A[i][k]) > amax) {
551 amax = fabs(A[i][k]);
552 imax = i;
553 }
554 }
555 // maximum equal to zero => complete pivoting
556 if (amax < eps_pivot) {
557 for (int j = k + 1; j < d; j++) {
558 for (int i = k; i < d - 2; i++) {
559 if (fabs(A[i][j]) > amax) {
560 amax = fabs(A[i][j]);
561 imax = i;
562 jmax = j;
563 }
564 }
565 }
566 if (amax < eps_pivot) {
567 delete[] colp;
568 return false;
569 }
570 // exchange columns
571 for (int i = 0; i < d - 2; i++) {
572 double tmp = A[i][k];
573 A[i][k] = A[i][jmax];
574 A[i][jmax] = tmp;
575 }
576 colp[k] = jmax;
577 }
578 // exchange rows
579 if (imax != k) {
580 for (int j = k; j < d; j++) {
581 double tmp = A[k][j];
582 A[k][j] = A[imax][j];
583 A[imax][j] = tmp;
584 }
585 }
586 // elimination
587 for (int i = k + 1; i < d - 2; i++) {
588 double factor = A[i][k] / A[k][k];
589 for (int j = k + 1; j < d; j++) A[i][j] -= factor * A[k][j];
590 }
591 }
592 // back substitution
593 colp[d - 2] = d - 2;
594 colp[d - 1] = d - 1;
595 basis[0][d - 2] = -1;
596 basis[0][d - 1] = 0;
597 basis[1][d - 2] = 0;
598 basis[1][d - 1] = -1;
599 for (int k = d - 3; k >= 0; k--) {
600 basis[0][k] = A[k][d - 2] / A[k][k];
601 basis[1][k] = A[k][d - 1] / A[k][k];
602 for (int i = k - 1; i >= 0; i--) {
603 A[i][d - 2] -= basis[0][k] * A[i][k];
604 A[i][d - 1] -= basis[1][k] * A[i][k];
605 }
606 }
607 // reverse column permutations
608 for (int k = d - 1; k >= 0; k--) {
609 if (colp[k] != k) {
610 for (int l = 0; l < 2; l++) {
611 double temp = basis[l][k];
612 basis[l][k] = basis[l][colp[k]];
613 basis[l][colp[k]] = temp;
614 }
615 }
616 }
617 delete[] colp;
618 return true;
619 }
620
621 int nHD_Comb2(double** xx, int n, int d);
622
623 /****************************************************************************/
624 /* */
625 /* 'HD2proj' performs the following steps: */
626 /* */
627 /* 1) All data points x[i] are projected on the space spanned by the */
628 /* two vectors passed in p, i.e., y[i,1] = p[1]'x[i] and */
629 /* y[i,2] = p[2]'x[i] are computed. */
630 /* 2) The bivariate integer halfspace depth of 0 is computed w.r.t. all */
631 /* the y[i] that are not equal to (0,0). */
632 /* 3) If there are more than d-2 values y[i] that are equal to (0,0), */
633 /* the respective points are projected on the orthogonal complement */
634 /* of p. Then, the integer halfspace depth of 0 w.r.t. these */
635 /* projected points is computed. */
636 /* 4) The sum of the values from step 2) and 3) is returned. */
637 /* */
638 /****************************************************************************/
639
HD2proj(double ** x,int n,int d,double ** p,int * indices)640 int HD2proj(double** x, int n, int d, double** p, int* indices) {
641 double** y = new double*[n];
642 for (int i = 0; i < n; i++) y[i] = new double[2];
643
644 int cnt0 = 0, cnt1 = 0, HDproj = 0;
645 int* plane = new int[n];
646 for (int i = 0; i < n; i++) {
647 y[cnt1][0] = y[cnt1][1] = 0;
648 for (int j = 0; j < d; j++)
649 for (int k = 0; k < 2; k++) y[cnt1][k] += p[k][j] * x[i][j];
650 if (norm2(y[cnt1], 2) > eps_HD2) cnt1++;
651 else plane[cnt0++] = i;
652
653 }
654 if (cnt0 > d - 2) {
655 double** z = new double*[cnt0];
656 for (int i = 0; i < cnt0; i++) {
657 z[i] = new double[d - 2];
658 for (int j = 0; j < d - 2; j++) {
659 z[i][j] = 0;
660 for (int k = 0; k < d; k++)
661 z[i][j] += x[indices[j]][k] * x[plane[i]][k];
662 }
663 }
664 HDproj = nHD_Comb2(z, cnt0, d - 2);
665 for (int i = 0; i < cnt0; i++) delete[] z[i];
666 delete[] z;
667 }
668 int result = intHD2(y, cnt1) + HDproj;
669 delete[] plane;
670 for (int i = 0; i < n; i++) delete[] y[i];
671 delete[] y;
672 return result;
673 }
674
675 /****************************************************************************/
676 /* */
677 /* 'nHD_Comb2' computes the integer halfspace depth of 0 w.r.t n data */
678 /* points in R^d. */
679 /* */
680 /* 'nHD_Comb2' implements the combinatorial algorithm (k = d-2) as */
681 /* described in Section 3.2 of "Exact computation of the halfspace depth" */
682 /* by Rainer Dyckerhoff and Pavlo Mozharovskyi (arXiv:1411:6927) */
683 /* */
684 /****************************************************************************/
685
nHD_Comb2(double ** xx,int n,int d)686 int nHD_Comb2(double** xx, int n, int d) {
687 if (d == 1) return intHD1(xx, n);
688 if (d == 2) return intHD2(xx, n);
689
690 int result = n + 1;
691 double** a = new double*[d - 2];
692 for (int i = 0; i < d - 2; i++) a[i] = new double[d];
693 double** p = new double*[2];
694 for (int i = 0; i < 2; i++) p[i] = new double[d];
695 int* indices = new int[d - 2];
696
697 indices[0] = -1;
698 int pos = 0;
699 while (pos >= 0) {
700 indices[pos]++;
701 for (pos++; pos < d - 2; pos++) indices[pos] = indices[pos - 1] + 1;
702 pos--;
703 do {
704 for (int i = 0; i < d - 2; i++)
705 for (int j = 0; j < d; j++) a[i][j] = xx[indices[i]][j];
706 if (getBasisComplement(a, d, p))
707 result = min(result, HD2proj(xx, n, d, p, indices));
708 indices[pos]++;
709 } while (indices[pos] < n - d + pos + 3);
710 do pos--; while (pos >= 0 && indices[pos] >= n - d + pos + 2);
711 }
712 for (int i = 0; i < d - 2; i++) delete[] a[i];
713 delete[] a;
714 for (int i = 0; i < 2; i++) delete[] p[i];
715 delete[] p;
716 delete[] indices;
717 return result;
718 }
719
720 /****************************************************************************/
721 /* */
722 /* 'HD_Comb2' computes the halfspace depth of a point z w.r.t. n data */
723 /* points in R^d. */
724 /* */
725 /* 'HD_Comb2' does some preprocessing of the data and then calls */
726 /* 'nHD_Comb2'. */
727 /* */
728 /* See also the description in the header file. */
729 /* */
730 /* HD_Comb2 calls the following routines: */
731 /* norm2 */
732 /* intHD1 */
733 /* intHD2 */
734 /* getRank */
735 /* project */
736 /* getBasisComplement */
737 /* HD2proj */
738 /* nHD_Rec */
739 /* */
740 /****************************************************************************/
741
HD_Comb2(double * z,double ** xx,int n,int d)742 double HD_Comb2(double* z, double** xx, int n, int d) {
743 if (n <= 0) throw invalid_argument("n <= 0");
744 if (d <= 0) throw invalid_argument("d <= 0");
745 int m = 0, rank;
746 int* indices = new int[d];
747 double** x = new double*[n];
748 // preprocess data
749 // subtract z from all data points x[i]
750 // check whether the data points are concentrated on a lower
751 // dimensional spcae
752 for (int i = 0; i < n; i++) {
753 x[m] = new double[d];
754 for (int j = 0; j < d; j++) x[m][j] = xx[i][j] - z[j];
755 if (norm2(x[m], d) >= eps_HDx) m++; else delete[] x[m];
756 }
757 if (m == 0) return 1.0;
758 rank = getRank(x, m, d, indices);
759 if (rank < d) project(x, m, d, rank, indices);
760
761 int result = nHD_Comb2(x, m, rank) + (n - m);
762 // deallocate array x
763 for (int i = 0; i < m; i++) delete[] x[i];
764 delete[] x;
765 delete[] indices;
766 return result / (double)n;
767 }
768
769
770