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 }