1 /*
2 * R : A Computer Language for Statistical Data Analysis
3
4 * Copyright (C) 1999-2014 The R Core Team
5 *
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation; either version 2 of the License, or
9 * (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program; if not, a copy is available at
18 * https://www.R-project.org/Licenses/.
19 */
20
21 #include <R_ext/Boolean.h>
22 #include <Rinternals.h>
23 #include "statsR.h"
24
cutree(SEXP merge,SEXP which)25 SEXP cutree(SEXP merge, SEXP which)
26 {
27 /* Return grouping vector from cutting a (binary) (cluster) tree
28 * into which[j] groups.
29 * merge = (n-1) x 2 matrix, described in help(hclust) */
30 SEXP ans;
31 int n, k, l, nclust, m1, m2, j, mm = 0;
32 Rboolean found_j, *sing;
33 int *m_nr, *z, *i_merge, *i_which, *i_ans;
34
35 PROTECT(merge = coerceVector(merge, INTSXP)); i_merge = INTEGER(merge);
36 PROTECT(which = coerceVector(which, INTSXP)); i_which = INTEGER(which);
37
38 n = nrows(merge)+1;
39 /* using 1-based indices ==> "--" */
40 sing = (Rboolean *) R_alloc(n, sizeof(Rboolean)); sing--;
41 m_nr = (int *) R_alloc(n, sizeof(int)); m_nr--;
42 z = (int *) R_alloc(n, sizeof(int)); z--;
43 PROTECT(ans = allocMatrix(INTSXP, n, LENGTH(which)));
44 i_ans = INTEGER(ans);
45
46 for(k = 1; k <= n; k++) {
47 sing[k] = TRUE;/* is k-th obs. still alone in cluster ? */
48 m_nr[k] = 0;/* containing last merge-step number of k-th obs. */
49 }
50
51 for(k = 1; k <= n-1; k++) {
52 /* k-th merge, from n-k+1 to n-k atoms: (m1,m2) = merge[ k , ] */
53 m1 = i_merge[k-1];
54 m2 = i_merge[n-1+k-1];
55
56 if(m1 < 0 && m2 < 0) {/* merging atoms [-m1] and [-m2] */
57 m_nr[-m1] = m_nr[-m2] = k;
58 sing[-m1] = sing[-m2] = FALSE;
59 }
60 else if(m1 < 0 || m2 < 0) {/* the other >= 0 */
61 if(m1 < 0) { j = -m1; m1 = m2; } else j = -m2;
62 /* merging atom j & cluster m1 */
63 for(l = 1; l <= n; l++)
64 if (m_nr[l] == m1)
65 m_nr[l] = k;
66 m_nr[j] = k;
67 sing[j] = FALSE;
68 }
69 else { /* both m1, m2 >= 0 */
70 for(l=1; l <= n; l++) {
71 if( m_nr[l] == m1 || m_nr[l] == m2 )
72 m_nr[l] = k;
73 }
74 }
75
76 /* does this k-th merge belong to a desired group size which[j] ?
77 * if yes, find j (maybe multiple ones): */
78 found_j = FALSE;
79 for(j = 0; j < LENGTH(which); j++) {
80 if(i_which[j] == n - k) {
81 if(!found_j) { /* first match (and usually only one) */
82 found_j = TRUE;
83 for(l = 1; l <= n; l++)
84 z[l] = 0;
85 nclust = 0;
86 mm = j*n; /*may want to copy this column of ans[] */
87 for(l = 1, m1 = mm; l <= n; l++, m1++) {
88 if(sing[l])
89 i_ans[m1] = ++nclust;
90 else {
91 if (z[m_nr[l]] == 0)
92 z[m_nr[l]] = ++nclust;
93 i_ans[m1] = z[m_nr[l]];
94 }
95 }
96 }
97 else { /* found_j: another which[j] == n-k : copy column */
98 for(l = 1, m1 = j*n, m2 = mm; l <= n; l++, m1++, m2++)
99 i_ans[m1] = i_ans[m2];
100 }
101 } /* if ( match ) */
102 } /* for(j .. which[j] ) */
103 } /* for(k ..) {merge} */
104
105 /* Dealing with trivial case which[] = n : */
106 for(j = 0; j < LENGTH(which); j++)
107 if(i_which[j] == n)
108 for(l = 1, m1 = j*n; l <= n; l++, m1++)
109 i_ans[m1] = l;
110
111 UNPROTECT(3);
112 return(ans);
113 }
114