1 #ifndef NMF_UTILS_H // include header only once
2 #define NMF_UTILS_H
3
4 #include <R.h>
5 #include <Rdefines.h>
6 #include <R_ext/Error.h>
7
8 extern "C" {
9
10 /** Returns the pointer address of 'x' as a character string*/
11 SEXP ptr_address (SEXP x);
12
13 /** Clone an object 'x'*/
14 SEXP clone_object (SEXP x);
15
16 /** pmin in place with 'y' being a single numeric value*/
17 SEXP ptr_pmax (SEXP x, SEXP y, SEXP skip);
18
19 /** Apply inequality constraints in place. */
20 SEXP ptr_neq_constraints(SEXP x, SEXP constraints, SEXP ratio=R_NilValue, SEXP value=R_NilValue);
21
22 /** Minimum per column*/
23 SEXP colMin(SEXP x);
24
25 /** Maximum per row*/
26 SEXP colMax(SEXP x);
27
28 /** Test if an external pointer is null.
29 *
30 * Function taken from the package bigmemory (v4.2.11).
31 */
ptr_isnil(SEXP address)32 SEXP ptr_isnil(SEXP address)
33 {
34 void *ptr = R_ExternalPtrAddr(address);
35 SEXP ret = PROTECT(NEW_LOGICAL(1));
36 LOGICAL_DATA(ret)[0] = (ptr==NULL) ? (Rboolean)TRUE : Rboolean(FALSE);
37 UNPROTECT(1);
38 return(ret);
39 }
40
41 }
42
43 // define the exported versions (for SEXP)
ptr_address(SEXP x)44 SEXP ptr_address (SEXP x){
45
46 SEXP ans = R_NilValue;
47 char tmp[15];
48 PROTECT(ans = allocVector(STRSXP, 1));
49 sprintf(tmp, "%p", (void *) x);
50 SET_STRING_ELT(ans, 0, mkChar(tmp));
51 UNPROTECT(1);
52 return ans;
53 }
54
clone_object(SEXP x)55 SEXP clone_object (SEXP x){
56
57 return Rf_duplicate(x);
58
59 }
60
ptr_pmax(SEXP x,SEXP y,SEXP skip=R_NilValue)61 SEXP ptr_pmax(SEXP x, SEXP y, SEXP skip=R_NilValue){
62
63 int n = length(x);
64 double* p_x = ( isNull(x) ? NULL : NUMERIC_POINTER(x) );
65 double lim = isNull(y) ? -1.0 : *NUMERIC_POINTER(y);
66
67 // backup skipped values
68 int n_skip = length(skip);
69 int ncol = isNull(GET_DIM(x)) ? 1 : INTEGER(GET_DIM(x))[1];
70 int nrow = n / ncol;
71 double* old_value = NULL;
72 int* p_skip = NULL;
73
74 if( !isNull(skip) && n_skip > 0 ){
75 old_value = (double*) R_alloc(n_skip*ncol, sizeof(double));
76 p_skip = INTEGER_POINTER(skip);
77 for(int k=ncol-1; k>=0; --k){
78 for(int i=n_skip-1; i>=0; --i){
79 //Rprintf("skip %i x %i\n", i, k);
80 int is = p_skip[i]-1;
81 double val = p_x[k*nrow + is];
82 old_value[k*n_skip + i] = val;
83 }
84 }
85 }
86
87 // apply limit inf to all values
88 double* p_x2 = p_x + n-1;
89 for(int i=n-1; i>=0; --i){
90 if( *p_x2 < lim )
91 *p_x2 = lim;
92 --p_x2;
93 }
94 p_x2 = NULL;
95
96 // restore skipped values
97 if( !isNull(skip) && n_skip > 0 ){
98 for(int k=ncol-1; k>=0; --k){
99 for(int i=n_skip-1; i>=0; --i){
100 //Rprintf("restore %i x %i\n", i, k);
101 int is = p_skip[i]-1;
102 p_x[k*nrow + is] = old_value[k*n_skip + i];
103 }
104 }
105 }
106
107
108 // return modified x
109 return x;
110 }
111
112 /** Apply inequality constraints in place. */
ptr_neq_constraints(SEXP x,SEXP constraints,SEXP ratio,SEXP value)113 SEXP ptr_neq_constraints(SEXP x, SEXP constraints, SEXP ratio, SEXP value){
114
115 double* p_x = ( isNull(x) ? NULL : NUMERIC_POINTER(x) );
116 double d_ratio = isNull(ratio) ? 0 : *NUMERIC_POINTER(ratio);
117 double* p_value = ( isNull(value) ? NULL : NUMERIC_POINTER(value) );
118 double eps = sqrt(DOUBLE_EPS);
119
120 // get dimensions
121 int ncol = isNull(GET_DIM(x)) ? 1 : INTEGER(GET_DIM(x))[1];
122 int nrow = isNull(GET_DIM(x)) ? length(x) : INTEGER(GET_DIM(x))[0];
123 int nc = length(constraints);
124 if( nc != ncol )
125 error("There must be as many elements in list `constraints` as columns in `x`.");
126
127 // apply each set of constraints (from first to last)
128 double* _xj = p_x; // pointer to marked column
129 double* _x_last = p_x + (ncol - 1) * nrow; // pointer to last column
130 for(int j=0; j<nc; ++j){
131 SEXP c_j = VECTOR_ELT(constraints, j);
132 int n = length(c_j);
133 int* p_i = INTEGER_POINTER(c_j);
134
135 // apply the constraint on each row in the set
136 for(int k=n-1; k>=0; --k){
137 double lim = d_ratio != 0.0 ? _xj[p_i[k]-1] / d_ratio - eps : 0.0;
138 if( lim < 0 )
139 lim = 0;
140
141 // apply constraints on each column
142 // pointer to current row in last column
143 double* _xi = _x_last + p_i[k]-1;
144 for(int l=ncol-1; l>=0; --l){
145 //Rprintf("Before: xi=%f > lim=%f ? => ", lim, *_xi);
146 if( l != j && *_xi > lim ){ // constrain column to 'lim'
147 *_xi = lim;
148 }else if( l == j && p_value != NULL ){ // constrain column to 'value'
149 *_xi = *p_value;
150 }
151 //Rprintf("xi=%f\n", *_xi);
152 // move to previous column
153 _xi -= nrow;
154 }
155 _xi = NULL;
156 }
157 // move to next marked column
158 _xj += nrow;
159 }
160
161 // return modified x
162 return x;
163 }
164
165
colMin(T * x,int n,int p,T * res,const T & NA_value)166 template<class T> inline void colMin(T* x, int n, int p, T* res, const T& NA_value){
167
168 // do nothing if there is no data or fill with NAs
169 if( n <= 0 ){
170 if( p <= 0 )
171 return;
172 for(int j=p-1; j>=0; --j, ++res)
173 *res = NA_value;
174 }
175
176 for(int j=p-1; j>=0; --j, ++res){
177 *res = *(x++);
178 for(int i=n-2; i>=0; --i, ++x){
179 if( *res > *x )
180 *res = *x;
181 }
182 }
183
184 }
185
colMax(T * x,int n,int p,T * res,const T & NA_value)186 template<class T> inline void colMax(T* x, int n, int p, T* res, const T& NA_value){
187
188 // do nothing if there is no data or fill with NAs
189 if( n <= 0 ){
190 if( p <= 0 )
191 return;
192 for(int j=p-1; j>=0; --j, ++res)
193 *res = NA_value;
194 }
195
196 for(int j=p-1; j>=0; --j, ++res){
197 *res = *(x++);
198 for(int i=n-2; i>=0; --i, ++x){
199 if( *res < *x )
200 *res = *x;
201 }
202 }
203
204 }
205
206 /**
207 * Minimum per column
208 */
colMin(SEXP x)209 SEXP colMin(SEXP x){
210
211 SEXP ans, dims;
212
213 // check that the argument is a matrix
214 dims = GET_DIM(x);
215 if (dims == R_NilValue)
216 error("a matrix-like object is required as argument to 'colMin'");
217 // check that it is a numeric data
218 if (!isNumeric(x))
219 error("a numeric object is required as argument to 'colMin'");
220
221 // get the dimension of the input matrix
222 int n = INTEGER(dims)[0];
223 int p = INTEGER(dims)[1];
224
225 if( TYPEOF(x) == REALSXP ){
226 // allocate memory for the result (a vector of length the number of columns of x)
227 PROTECT(ans = allocVector(REALSXP, p));
228 colMin(NUMERIC_POINTER(x), n, p, NUMERIC_POINTER(ans), NA_REAL);
229 UNPROTECT(1);
230 }
231 else{
232 // allocate memory for the result (a vector of length the number of columns of x)
233 PROTECT(ans = allocVector(INTSXP, p));
234 colMin(INTEGER_POINTER(x), n, p, INTEGER_POINTER(ans), NA_INTEGER);
235 UNPROTECT(1);
236 }
237
238 return ans;
239 }
240
241 /**
242 * Maximum per column
243 */
colMax(SEXP x)244 SEXP colMax(SEXP x){
245
246 SEXP ans, dims;
247
248 // check that the argument is a matrix
249 dims = GET_DIM(x);
250 if (dims == R_NilValue)
251 error("a matrix-like object is required as argument to 'colMax'");
252 // check that it is a numeric data
253 if (!isNumeric(x))
254 error("a numeric object is required as argument to 'colMax'");
255
256 // get the dimension of the input matrix
257 int n = INTEGER(dims)[0];
258 int p = INTEGER(dims)[1];
259
260 if( TYPEOF(x) == REALSXP ){
261 // allocate memory for the result (a vector of length the number of columns of x)
262 PROTECT(ans = allocVector(REALSXP, p));
263 colMax(NUMERIC_POINTER(x), n, p, NUMERIC_POINTER(ans), NA_REAL);
264 UNPROTECT(1);
265 }
266 else{
267 // allocate memory for the result (a vector of length the number of columns of x)
268 PROTECT(ans = allocVector(INTSXP, p));
269 colMax(INTEGER_POINTER(x), n, p, INTEGER_POINTER(ans), NA_INTEGER);
270 UNPROTECT(1);
271 }
272
273 return ans;
274 }
275
276 #endif //END ifdef NMF_UTILS_H
277