9b3049561e6747ed85e47ca0c6ffaf013dc9c7c4
[synclust.git] / src / adapters / a.convexSolver.c
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 }