Commit | Line | Data |
---|---|---|
15d1825d BA |
1 | #include <R.h> |
2 | #include <Rdefines.h> | |
3 | #include "sources/convexSolver.h" | |
4 | #include "sources/utils/algebra.h" | |
5 | ||
6 | // compute estimated ("repaired", "smoothed"...) variations from rows of M | |
7 | // NOTE: geographic coordinates dropped here, since they are unused | |
8 | SEXP getVarsWithConvexOptim( | |
9 | SEXP M_, | |
10 | SEXP NIix_, | |
11 | SEXP alpha_, | |
12 | SEXP h_, | |
13 | SEXP epsilon_, | |
14 | SEXP maxiter_, | |
15 | SEXP symmNeighbs_, | |
16 | SEXP trace_ | |
17 | ) { | |
18 | // get parameters | |
19 | double alpha = NUMERIC_VALUE(alpha_); | |
20 | double h = NUMERIC_VALUE(h_); | |
21 | double epsilon = NUMERIC_VALUE(epsilon_); | |
22 | int maxiter = INTEGER_VALUE(maxiter_); | |
23 | int symmNeighbs = LOGICAL_VALUE(symmNeighbs_); | |
24 | int trace = LOGICAL_VALUE(trace_); | |
25 | ||
26 | // extract infos from M and get associate pointer | |
27 | SEXP dim = getAttrib(M_, R_DimSymbol); | |
28 | int nrow = INTEGER(dim)[0]; | |
29 | int ncol = INTEGER(dim)[1]; | |
30 | // M is always given by columns: easier to process in rows | |
31 | double* pM = transpose(REAL(M_), nrow, ncol); | |
32 | ||
33 | // extract NIix list vectors in a jagged array | |
34 | int* lengthNIix = (int*)malloc(nrow*sizeof(int)); | |
35 | int** NIix = (int**)malloc(nrow*sizeof(int*)); | |
36 | for (int i=0; i<nrow; i++) | |
37 | { | |
38 | lengthNIix[i] = LENGTH(VECTOR_ELT(NIix_,i)); | |
39 | SEXP tmp; | |
40 | PROTECT(tmp = AS_INTEGER(VECTOR_ELT(NIix_,i))); | |
41 | NIix[i] = (int*)malloc(lengthNIix[i]*sizeof(int)); | |
42 | for (int j=0; j<lengthNIix[i]; j++) | |
43 | NIix[i][j] = INTEGER(tmp)[j]; | |
44 | UNPROTECT(1); | |
45 | // WARNING: R indices start at 1, | |
46 | // so we must lower every index right now to avoid future bug | |
47 | for (int j=0; j<lengthNIix[i]; j++) | |
48 | NIix[i][j]--; | |
49 | } | |
50 | ||
51 | // Main call to core algorithm | |
52 | Parameters params = getVarsWithConvexOptim_core( | |
53 | pM, lengthNIix, NIix, nrow, ncol, alpha, h, epsilon, maxiter, symmNeighbs, trace); | |
54 | ||
55 | // free neighborhoods parameters arrays | |
56 | free(lengthNIix); | |
57 | for (int i=0; i<nrow; i++) | |
58 | free(NIix[i]); | |
59 | free(NIix); | |
60 | ||
61 | // copy matrix F into pF for output to R (1D matrices) | |
62 | SEXP f; | |
63 | PROTECT(f = allocMatrix(REALSXP, nrow, ncol)); | |
64 | double* pF = REAL(f); | |
65 | for (int i=0; i<nrow; i++) | |
66 | { | |
67 | for (int j=0; j<ncol; j++) | |
68 | pF[i+nrow*j] = params.f[i][j]; | |
69 | } | |
70 | // copy theta into pTheta for output to R | |
71 | SEXP theta; | |
72 | PROTECT(theta = allocVector(REALSXP, nrow)); | |
73 | double* pTheta = REAL(theta); | |
74 | for (int i=0; i<nrow; i++) | |
75 | pTheta[i] = params.theta[i]; | |
76 | ||
77 | // free params.f and params.theta | |
78 | free(params.theta); | |
79 | for (int i=0; i<nrow; i++) | |
80 | free(params.f[i]); | |
81 | free(params.f); | |
82 | ||
83 | // build return list with f and theta | |
84 | SEXP listParams, listNames; | |
85 | PROTECT(listParams = allocVector(VECSXP, 2)); | |
86 | char* lnames[2] = {"f", "theta"}; //lists labels | |
87 | PROTECT(listNames = allocVector(STRSXP,2)); | |
88 | for (int i=0; i<2; i++) | |
89 | SET_STRING_ELT(listNames,i,mkChar(lnames[i])); | |
90 | setAttrib(listParams, R_NamesSymbol, listNames); | |
91 | SET_VECTOR_ELT(listParams, 0, f); | |
92 | SET_VECTOR_ELT(listParams, 1, theta); | |
93 | ||
94 | UNPROTECT(4); | |
95 | return listParams; | |
96 | } |