241a09f52bbc4a76701fe60407d894c20905015f
3 #include "sources/EMGLLF.h"
4 #include "sources/utils/io.h"
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_
);
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
);
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
++)
36 lengthNIix
[i
] = LENGTH(VECTOR_ELT(NIix_
,i
));
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
];
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
++)
49 // Main call to core algorithm
50 Parameters params
= getVarsWithConvexOptim_core(
51 pM
, lengthNIix
, NIix
, nrow
, ncol
, alpha
, h
, epsilon
, maxiter
, symmNeighbs
, trace
);
53 // free neighborhoods parameters arrays
55 for (int i
=0; i
<nrow
; i
++)
59 // copy matrix F into pF for output to R (1D matrices)
61 PROTECT(f
= allocMatrix(REALSXP
, nrow
, ncol
));
63 for (int i
=0; i
<nrow
; i
++)
65 for (int j
=0; j
<ncol
; j
++)
66 pF
[i
+nrow
*j
] = params
.f
[i
][j
];
68 // copy theta into pTheta for output to R
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
];
75 // free params.f and params.theta
77 for (int i
=0; i
<nrow
; i
++)
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
);
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];
113 const mwSize
* dimPhiInit
= mxGetDimensions(prhs
[0]);
114 Real
* brPhiInit
= matlabToBrArray_real(mxGetPr(prhs
[0]), dimPhiInit
, 3);
117 const mwSize
* dimRhoInit
= mxGetDimensions(prhs
[1]);
118 Real
* brRhoInit
= matlabToBrArray_real(mxGetPr(prhs
[1]), dimRhoInit
, 3);
121 Real
* piInit
= mxGetPr(prhs
[2]);
124 const mwSize
* dimGamInit
= mxGetDimensions(prhs
[3]);
125 Real
* brGamInit
= matlabToBrArray_real(mxGetPr(prhs
[3]), dimGamInit
, 2);
127 // min number of iterations
128 Int mini
= ((Int
*)mxGetData(prhs
[4]))[0];
130 // max number of iterations
131 Int maxi
= ((Int
*)mxGetData(prhs
[5]))[0];
134 Real gamma
= mxGetScalar(prhs
[6]);
137 Real lambda
= mxGetScalar(prhs
[7]);
140 const mwSize
* dimX
= mxGetDimensions(prhs
[8]);
141 Real
* brX
= matlabToBrArray_real(mxGetPr(prhs
[8]), dimX
, 2);
144 const mwSize
* dimY
= mxGetDimensions(prhs
[9]);
145 Real
* brY
= matlabToBrArray_real(mxGetPr(prhs
[9]), dimY
, 2);
148 Real tau
= mxGetScalar(prhs
[10]);
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]);
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]);
165 plhs
[2] = mxCreateNumericMatrix(k
,1,mxDOUBLE_CLASS
,mxREAL
);
166 Real
* pi
= mxGetPr(plhs
[2]);
169 plhs
[3] = mxCreateNumericMatrix(maxi
,1,mxDOUBLE_CLASS
,mxREAL
);
170 Real
* LLF
= mxGetPr(plhs
[3]);
173 const mwSize dimS
[] = {p
, m
, k
};
174 plhs
[4] = mxCreateNumericArray(3,dimS
,mxDOUBLE_CLASS
,mxREAL
);
175 Real
* S
= mxGetPr(plhs
[4]);
181 EMGLLF(brPhiInit
,brRhoInit
,piInit
,brGamInit
,mini
,maxi
,gamma
,lambda
,brX
,brY
,tau
,