Commit | Line | Data |
---|---|---|
cb51adb8 BA |
1 | #include <R.h> |
2 | #include <Rdefines.h> | |
3 | #include "sources/EMGLLF.h" | |
4 | #include "sources/utils/io.h" | |
5 | ||
6 | SEXP EMGLLF( | |
7 | SEXP M_, | |
8 | SEXP NIix_, | |
9 | SEXP alpha_, | |
10 | SEXP h_, | |
11 | SEXP epsilon_, | |
12 | SEXP maxiter_, | |
13 | SEXP symmNeighbs_, | |
14 | SEXP trace_ | |
15 | ) { | |
16 | // get parameters | |
17 | double alpha = NUMERIC_VALUE(alpha_); | |
18 | double h = NUMERIC_VALUE(h_); | |
19 | double epsilon = NUMERIC_VALUE(epsilon_); | |
20 | int maxiter = INTEGER_VALUE(maxiter_); | |
21 | int symmNeighbs = LOGICAL_VALUE(symmNeighbs_); | |
22 | int trace = LOGICAL_VALUE(trace_); | |
23 | ||
24 | // extract infos from M and get associate pointer | |
25 | SEXP dim = getAttrib(M_, R_DimSymbol); | |
26 | int nrow = INTEGER(dim)[0]; | |
27 | int ncol = INTEGER(dim)[1]; | |
28 | // M is always given by columns: easier to process in rows | |
29 | double* pM = transpose(REAL(M_), nrow, ncol); | |
30 | ||
31 | // extract NIix list vectors in a jagged array | |
32 | int* lengthNIix = (int*)malloc(nrow*sizeof(int)); | |
33 | int** NIix = (int**)malloc(nrow*sizeof(int*)); | |
34 | for (int i=0; i<nrow; i++) | |
35 | { | |
36 | lengthNIix[i] = LENGTH(VECTOR_ELT(NIix_,i)); | |
37 | SEXP tmp; | |
38 | PROTECT(tmp = AS_INTEGER(VECTOR_ELT(NIix_,i))); | |
39 | NIix[i] = (int*)malloc(lengthNIix[i]*sizeof(int)); | |
40 | for (int j=0; j<lengthNIix[i]; j++) | |
41 | NIix[i][j] = INTEGER(tmp)[j]; | |
42 | UNPROTECT(1); | |
43 | // WARNING: R indices start at 1, | |
44 | // so we must lower every index right now to avoid future bug | |
45 | for (int j=0; j<lengthNIix[i]; j++) | |
46 | NIix[i][j]--; | |
47 | } | |
48 | ||
49 | // Main call to core algorithm | |
50 | Parameters params = getVarsWithConvexOptim_core( | |
51 | pM, lengthNIix, NIix, nrow, ncol, alpha, h, epsilon, maxiter, symmNeighbs, trace); | |
52 | ||
53 | // free neighborhoods parameters arrays | |
54 | free(lengthNIix); | |
55 | for (int i=0; i<nrow; i++) | |
56 | free(NIix[i]); | |
57 | free(NIix); | |
58 | ||
59 | // copy matrix F into pF for output to R (1D matrices) | |
60 | SEXP f; | |
61 | PROTECT(f = allocMatrix(REALSXP, nrow, ncol)); | |
62 | double* pF = REAL(f); | |
63 | for (int i=0; i<nrow; i++) | |
64 | { | |
65 | for (int j=0; j<ncol; j++) | |
66 | pF[i+nrow*j] = params.f[i][j]; | |
67 | } | |
68 | // copy theta into pTheta for output to R | |
69 | SEXP theta; | |
70 | PROTECT(theta = allocVector(REALSXP, nrow)); | |
71 | double* pTheta = REAL(theta); | |
72 | for (int i=0; i<nrow; i++) | |
73 | pTheta[i] = params.theta[i]; | |
74 | ||
75 | // free params.f and params.theta | |
76 | free(params.theta); | |
77 | for (int i=0; i<nrow; i++) | |
78 | free(params.f[i]); | |
79 | free(params.f); | |
80 | ||
81 | // build return list with f and theta | |
82 | SEXP listParams, listNames; | |
83 | PROTECT(listParams = allocVector(VECSXP, 2)); | |
84 | char* lnames[2] = {"f", "theta"}; //lists labels | |
85 | PROTECT(listNames = allocVector(STRSXP,2)); | |
86 | for (int i=0; i<2; i++) | |
87 | SET_STRING_ELT(listNames,i,mkChar(lnames[i])); | |
88 | setAttrib(listParams, R_NamesSymbol, listNames); | |
89 | SET_VECTOR_ELT(listParams, 0, f); | |
90 | SET_VECTOR_ELT(listParams, 1, theta); | |
91 | ||
92 | UNPROTECT(4); | |
93 | return listParams; | |
94 | ||
95 | ||
96 | ||
97 | ||
98 | ||
99 | ||
100 | ||
101 | ||
102 | // Get matrices dimensions | |
103 | const mwSize n = mxGetDimensions(prhs[8])[0]; | |
104 | const mwSize p = mxGetDimensions(prhs[0])[0]; | |
105 | const mwSize m = mxGetDimensions(prhs[0])[1]; | |
106 | const mwSize k = mxGetDimensions(prhs[0])[2]; | |
107 | ||
108 | //////////// | |
109 | // INPUTS // | |
110 | //////////// | |
111 | ||
112 | // phiInit | |
113 | const mwSize* dimPhiInit = mxGetDimensions(prhs[0]); | |
114 | Real* brPhiInit = matlabToBrArray_real(mxGetPr(prhs[0]), dimPhiInit, 3); | |
115 | ||
116 | // rhoInit | |
117 | const mwSize* dimRhoInit = mxGetDimensions(prhs[1]); | |
118 | Real* brRhoInit = matlabToBrArray_real(mxGetPr(prhs[1]), dimRhoInit, 3); | |
119 | ||
120 | // piInit | |
121 | Real* piInit = mxGetPr(prhs[2]); | |
122 | ||
123 | // gamInit | |
124 | const mwSize* dimGamInit = mxGetDimensions(prhs[3]); | |
125 | Real* brGamInit = matlabToBrArray_real(mxGetPr(prhs[3]), dimGamInit, 2); | |
126 | ||
127 | // min number of iterations | |
128 | Int mini = ((Int*)mxGetData(prhs[4]))[0]; | |
129 | ||
130 | // max number of iterations | |
131 | Int maxi = ((Int*)mxGetData(prhs[5]))[0]; | |
132 | ||
133 | // gamma | |
134 | Real gamma = mxGetScalar(prhs[6]); | |
135 | ||
136 | // lambda | |
137 | Real lambda = mxGetScalar(prhs[7]); | |
138 | ||
139 | // X | |
140 | const mwSize* dimX = mxGetDimensions(prhs[8]); | |
141 | Real* brX = matlabToBrArray_real(mxGetPr(prhs[8]), dimX, 2); | |
142 | ||
143 | // Y | |
144 | const mwSize* dimY = mxGetDimensions(prhs[9]); | |
145 | Real* brY = matlabToBrArray_real(mxGetPr(prhs[9]), dimY, 2); | |
146 | ||
147 | // tau | |
148 | Real tau = mxGetScalar(prhs[10]); | |
149 | ||
150 | ///////////// | |
151 | // OUTPUTS // | |
152 | ///////////// | |
153 | ||
154 | // phi | |
155 | const mwSize dimPhi[] = {dimPhiInit[0], dimPhiInit[1], dimPhiInit[2]}; | |
156 | plhs[0] = mxCreateNumericArray(3,dimPhi,mxDOUBLE_CLASS,mxREAL); | |
157 | Real* phi = mxGetPr(plhs[0]); | |
158 | ||
159 | // rho | |
160 | const mwSize dimRho[] = {dimRhoInit[0], dimRhoInit[1], dimRhoInit[2]}; | |
161 | plhs[1] = mxCreateNumericArray(3,dimRho,mxDOUBLE_CLASS,mxREAL); | |
162 | Real* rho = mxGetPr(plhs[1]); | |
163 | ||
164 | // pi | |
165 | plhs[2] = mxCreateNumericMatrix(k,1,mxDOUBLE_CLASS,mxREAL); | |
166 | Real* pi = mxGetPr(plhs[2]); | |
167 | ||
168 | // LLF | |
169 | plhs[3] = mxCreateNumericMatrix(maxi,1,mxDOUBLE_CLASS,mxREAL); | |
170 | Real* LLF = mxGetPr(plhs[3]); | |
171 | ||
172 | // S | |
173 | const mwSize dimS[] = {p, m, k}; | |
174 | plhs[4] = mxCreateNumericArray(3,dimS,mxDOUBLE_CLASS,mxREAL); | |
175 | Real* S = mxGetPr(plhs[4]); | |
176 | ||
177 | //////////////////// | |
178 | // Call to EMGLLF // | |
179 | //////////////////// | |
180 | ||
181 | EMGLLF(brPhiInit,brRhoInit,piInit,brGamInit,mini,maxi,gamma,lambda,brX,brY,tau, | |
182 | phi,rho,pi,LLF,S, | |
183 | n,p,m,k); | |
184 | ||
185 | free(brPhiInit); | |
186 | free(brRhoInit); | |
187 | free(brGamInit); | |
188 | free(brX); | |
189 | free(brY); | |
190 | ||
191 | ||
192 | ||
193 | ||
194 | ||
195 | ||
196 | ||
197 | } |