Main Page | Directories | File List | File Members | Related Pages

# Utils.c

Go to the documentation of this file.

    00001
    00009 #include "party.h"
    00010
    00011
    00023 void C_kronecker (const double *A, const int m, const int n,
    00024                   const double *B, const int r, const int s,
    00025                   double *ans) {
    00026
    00027     int i, j, k, l, mr, js, ir;
    00028     double y;
    00029
    00030     mr = m * r;
    00031     for (i = 0; i < m; i++) {
    00032         ir = i * r;
    00033         for (j = 0; j < n; j++) {
    00034             js = j * s;
    00035             y = A[j*m + i];
    00036             for (k = 0; k < r; k++) {
    00037                 for (l = 0; l < s; l++) {
    00038                     ans[(js + l) * mr + ir + k] = y * B[l * r + k];
    00039                 }
    00040             }
    00041         }
    00042     }
    00043 }
    00044
    00045
    00052 SEXP R_kronecker (SEXP A, SEXP B) {
    00053
    00054     /*  The Kronecker product, a real (mr x ns) matrix */
    00055     SEXP ans;
    00056     int *adim, *bdim;
    00057
    00058     if (!isReal(A) || !isReal(B))
    00059         error("R_kronecker: A and B are not of type REALSXP");
    00060
    00061     if (isMatrix(A)) {
    00062         adim = INTEGER(getAttrib(A, R_DimSymbol));
    00063     } else {
    00064         /* assume row vectors */
    00065         adim = Calloc(2, int);
    00066         adim[0] = 1;
    00067         adim[1] = LENGTH(A);
    00068     }
    00069
    00070     if (isMatrix(B)) {
    00071         bdim = INTEGER(getAttrib(B, R_DimSymbol));
    00072     } else {
    00073         /* assume row vectors */
    00074         bdim = Calloc(2, int);
    00075         bdim[0] = 1;
    00076         bdim[1] = LENGTH(B);
    00077     }
    00078
    00079     PROTECT(ans = allocMatrix(REALSXP,
    00080                               adim[0] * bdim[0],
    00081                               adim[1] * bdim[1]));
    00082     C_kronecker(REAL(A), adim[0], adim[1],
    00083                 REAL(B), bdim[0], bdim[1], REAL(ans));
    00084     if (!isMatrix(A)) Free(adim);
    00085     if (!isMatrix(B)) Free(bdim);
    00086     UNPROTECT(1);
    00087     return(ans);
    00088 }
    00089
    00090
    00097 SEXP CR_svd (SEXP x, SEXP svdmem) {
    00098
    00099     int p, i;
    00100     double *du, *dv;
    00101
    00102     if (!isMatrix(x) || !isReal(x))
    00103         error("x is not a real matrix");
    00104
    00105     du = REAL(GET_SLOT(svdmem, PL2_uSym));
    00106     dv = REAL(GET_SLOT(svdmem, PL2_vSym));
    00107     p = INTEGER(GET_SLOT(svdmem, PL2_pSym))[0];
    00108     for (i = 0; i < p*p; i++) {
    00109         du[i] = 0.0;
    00110         dv[i] = 0.0;
    00111     }
    00112     SET_SLOT(svdmem, PL2_svdSym, La_svd(GET_SLOT(svdmem, PL2_jobuSym),
    00113         GET_SLOT(svdmem, PL2_jobvSym), x, GET_SLOT(svdmem, PL2_sSym),
    00114         GET_SLOT(svdmem, PL2_uSym), GET_SLOT(svdmem, PL2_vSym),
    00115         GET_SLOT(svdmem, PL2_methodSym)));
    00116     return(R_NilValue);
    00117 }
    00118
    00119
    00128 void C_MPinv (SEXP x, double tol, SEXP svdmem, SEXP ans) {
    00129
    00130     SEXP svdx, d, u, vt, dummy;
    00131     int i, j, p, k, *positive;
    00132     double *dd, *du, *dvt, *dMPinv;
    00133     double *drank;
    00134
    00135     drank = REAL(GET_SLOT(ans, PL2_rankSym));
    00136     dMPinv = REAL(GET_SLOT(ans, PL2_MPinvSym));
    00137
    00138     dummy = CR_svd(x, svdmem);
    00139     svdx = GET_SLOT(svdmem, PL2_svdSym);
    00140     d = VECTOR_ELT(svdx, 0);
    00141     dd = REAL(d);
    00142     u = VECTOR_ELT(svdx, 1);
    00143     du = REAL(u);
    00144     vt = VECTOR_ELT(svdx, 2);
    00145     dvt = REAL(vt);
    00146     p = LENGTH(d);
    00147
    00148     if (tol * dd[0] > tol) tol = tol * dd[0];
    00149
    00150     positive = Calloc(p, int);
    00151
    00152     drank[0] = 0.0;
    00153     for (i = 0; i < p; i++) {
    00154         if (dd[i] > tol) {
    00155             positive[i] = 1;
    00156             drank[0] += 1.0;
    00157         }
    00158     }
    00159
    00160     for (j = 0; j < p; j++) {
    00161         if (positive[j]) {
    00162             for (i = 0; i < p; i++)
    00163                 du[j * p + i] *= (1 / dd[j]);
    00164         }
    00165     }
    00166
    00167     for (i = 0; i < p; i++) {
    00168         for (j = 0; j < p; j++) {
    00169             dMPinv[j * p + i] = 0.0;
    00170             for (k = 0; k < p; k++) {
    00171                 if (positive[k])
    00172                     dMPinv[j * p + i] += dvt[i * p + k] * du[p * k + j];
    00173             }
    00174         }
    00175     }
    00176
    00177     Free(positive);
    00178 }
    00179
    00187 SEXP R_MPinv (SEXP x, SEXP tol, SEXP svdmem) {
    00188
    00189     SEXP ans;
    00190     int p;
    00191
    00192     if (!isMatrix(x) || !isReal(x))
    00193         error("R_MPinv: x is not a real matrix");
    00194
    00195     if (nrow(x) != ncol(x))
    00196         error("R_MPinv: x is not a square matrix");
    00197
    00198     if (!isReal(tol) || LENGTH(tol) != 1)
    00199         error("R_MPinv: tol is not a scalar real");
    00200
    00201     p = nrow(x);
    00202     if (p != INTEGER(GET_SLOT(svdmem, PL2_pSym))[0])
    00203         error("R_MPinv: dimensions don't match");
    00204
    00205     PROTECT(ans = NEW_OBJECT(MAKE_CLASS("LinStatExpectCovarMPinv")));
    00206     SET_SLOT(ans, PL2_MPinvSym, PROTECT(allocMatrix(REALSXP, p, p)));
    00207     SET_SLOT(ans, PL2_rankSym, PROTECT(allocVector(REALSXP, 1)));
    00208
    00209     C_MPinv(x, REAL(tol)[0], svdmem, ans);
    00210
    00211     UNPROTECT(3);
    00212     return(ans);
    00213 }
    00214
    00222 double C_max(const double *x, const int n) {
    00223    double tmp = 0.0;
    00224    int i;
    00225
    00226    for (i = 0; i < n; i++) {
    00227        if (x[i] > tmp) tmp = x[i];
    00228    }
    00229    return(tmp);
    00230 }
    00231
    00232
    00238 SEXP R_max(SEXP x) {
    00239
    00240     SEXP ans;
    00241     int n;
    00242
    00243     if (!isReal(x))
    00244         error("R_max: x is not of type REALSXP");
    00245     n = LENGTH(x);
    00246     PROTECT(ans = allocVector(REALSXP, 1));
    00247     REAL(ans)[0] = C_max(REAL(x), n);
    00248     UNPROTECT(1);
    00249     return(ans);
    00250 }
    00251
    00252
    00259 void C_abs(double *x, int n) {
    00260
    00261     int i;
    00262     for (i = 0; i < n; i++) x[i] = fabs(x[i]);
    00263 }
    00264
    00265
    00271 SEXP R_abs(SEXP x) {
    00272
    00273     SEXP ans;
    00274     int n;
    00275
    00276     if (!isReal(x))
    00277         error("R_max: x is not of type REALSXP");
    00278     n = LENGTH(x);
    00279     PROTECT(ans = duplicate(x));
    00280     C_abs(REAL(ans), n);
    00281     UNPROTECT(1);
    00282     return(ans);
    00283 }
    00284
    00285
    00297 void C_matprod(double *x, int nrx, int ncx,
    00298                double *y, int nry, int ncy, double *z)
    00299 {
    00300     char *transa = "N", *transb = "N";
    00301     double one = 1.0, zero = 0.0;
    00302     int i;
    00303
    00304     if (nrx > 0 && ncx > 0 && nry > 0 && ncy > 0) {
    00305         F77_CALL(dgemm)(transa, transb, &nrx, &ncy, &ncx, &one,
    00306                         x, &nrx, y, &nry, &zero, z, &nrx);
    00307     } else /* zero-extent operations should return zeroes */
    00308         for(i = 0; i < nrx*ncy; i++) z[i] = 0;
    00309 }
    00310
    00311
    00318 SEXP R_matprod(SEXP x, SEXP y) {
    00319
    00320     SEXP ans;
    00321
    00322     int nrx, ncx, nry, ncy;
    00323
    00324     nrx = nrow(x);
    00325     ncx = ncol(x);
    00326     nry = nrow(y);
    00327     ncy = ncol(y);
    00328
    00329     if (ncx != nry)
    00330         error("R_matprod: dimensions don't match");
    00331     PROTECT(ans = allocMatrix(REALSXP, nrx, ncy));
    00332     C_matprod(REAL(x), nrx, ncx, REAL(y), nry, ncy, REAL(ans));
    00333     UNPROTECT(1);
    00334     return(ans);
    00335 }
    00336
    00337
    00349 void C_matprodT(double *x, int nrx, int ncx,
    00350                 double *y, int nry, int ncy, double *z)
    00351 {
    00352     char *transa = "N", *transb = "T";
    00353     double one = 1.0, zero = 0.0;
    00354     int i;
    00355
    00356     if (nrx > 0 && ncx > 0 && nry > 0 && ncy > 0) {
    00357         F77_CALL(dgemm)(transa, transb, &nrx, &nry, &ncy, &one,
    00358                         x, &nrx, y, &nry, &zero, z, &nrx);
    00359     } else /* zero-extent operations should return zeroes */
    00360         for(i = 0; i < nrx*nry; i++) z[i] = 0;
    00361 }
    00362
    00363
    00370 SEXP R_matprodT(SEXP x, SEXP y) {
    00371
    00372     SEXP ans;
    00373     int nrx, ncx, nry, ncy;
    00374
    00375     nrx = nrow(x);
    00376     ncx = ncol(x);
    00377     nry = nrow(y);
    00378     ncy = ncol(y);
    00379
    00380     if (ncx != ncy)
    00381         error("R_matprod: dimensions don't match");
    00382     PROTECT(ans = allocMatrix(REALSXP, nrx, nry));
    00383     C_matprodT(REAL(x), nrx, ncx, REAL(y), nry, ncy, REAL(ans));
    00384     UNPROTECT(1);
    00385     return(ans);
    00386 }
    00387
    00388
    00397 void C_SampleNoReplace(int *x, int m, int k, int *ans) {
    00398
    00399     int i, j, n = m;
    00400
    00401     for (i = 0; i < m; i++)
    00402         x[i] = i;
    00403     for (i = 0; i < k; i++) {
    00404         j = n * unif_rand();
    00405         ans[i] = x[j];
    00406         x[j] = x[--n];
    00407     }
    00408 }
    00409
    00410
    00416 SEXP R_permute(SEXP m) {
    00417
    00418     SEXP x, ans;
    00419     int n;
    00420
    00421     n = INTEGER(m)[0];
    00422     PROTECT(x = allocVector(INTSXP, n));
    00423     PROTECT(ans = allocVector(INTSXP, n));
    00424     C_SampleNoReplace(INTEGER(x), n, n, INTEGER(ans));
    00425     UNPROTECT(2);
    00426     return(ans);
    00427 }
    00428
    00429
    00436 SEXP R_rsubset(SEXP m, SEXP k) {
    00437
    00438     SEXP x, ans;
    00439     int n, j;
    00440
    00441     n = INTEGER(m)[0];
    00442     j = INTEGER(k)[0];
    00443     PROTECT(x = allocVector(INTSXP, n));
    00444     PROTECT(ans = allocVector(INTSXP, j));
    00445     C_SampleNoReplace(INTEGER(x), n, j, INTEGER(ans));
    00446     UNPROTECT(2);
    00447     return(ans);
    00448 }
    00449
    00450
    00458 int i_in_set(int i, int *iset, int p) {
    00459
    00460     int j, is = 0;
    00461
    00462     if (p == 0) return(0);
    00463
    00464     for (j = 0; j < p; j++) {
    00465         if (iset[j] == i) {
    00466             is = 1;
    00467             break;
    00468         }
    00469     }
    00470     return(is);
    00471 }
    00472
    00473 int C_i_in_set(int i, SEXP set) {
    00474     if (LENGTH(set) > 0)
    00475         return(i_in_set(i, INTEGER(set), LENGTH(set)));
    00476     else
    00477         return(0);
    00478 }
    00479
    00480 int nrow(SEXP x) {
    00481     return(INTEGER(getAttrib(x, R_DimSymbol))[0]);
    00482 }
    00483
    00484 int ncol(SEXP x) {
    00485     return(INTEGER(getAttrib(x, R_DimSymbol))[1]);
    00486 }
    00487
    00488 /* compute index of variable with smallest p-value
    00489    (and largest test statistic in case two or more p-values coincide --
    00490     should not happen anymore since we use 1 - (1 - p)^k for Bonferroni adjustment)
    00491 */
    00492 int C_whichmax(double *pvalue, double *teststat, int ninputs) {
    00493
    00494     int ans = -1, j;
    00495     double tmppval = 0.0, tmptstat = 0.0;
    00496
    00497     /* <FIXME> can we switch to the log scale here? </FIXME> */
    00498
    00499     tmppval = 0.0;
    00500     tmptstat = 0.0;
    00501     for (j = 0; j < ninputs; j++) {
    00502         if (pvalue[j] > tmppval) {
    00503             ans = j;
    00504             tmppval = pvalue[j];
    00505             tmptstat = teststat[j];
    00506         } else {
    00507             if (pvalue[j] == tmppval && teststat[j] > tmptstat) {
    00508                 ans = j;
    00509                 tmppval = pvalue[j];
    00510                 tmptstat = teststat[j];
    00511             }
    00512         }
    00513     }
    00514     return(ans);
    00515 }
    00516
    00517 SEXP R_whichmax(SEXP x, SEXP y) {
    00518     SEXP ans;
    00519
    00520     if (LENGTH(x) != LENGTH(y)) error("different length");
    00521     PROTECT(ans = allocVector(INTSXP, 1));
    00522     INTEGER(ans)[0] = C_whichmax(REAL(x), REAL(y), LENGTH(x));
    00523     UNPROTECT(1);
    00524     return(ans);
    00525 }
    00526
    00527 SEXP R_listplus(SEXP a, SEXP b, SEXP which) {
    00528
    00529     int na, nb, i, j, *iwhich;
    00530     double *dae, *dbe;
    00531     SEXP ae, be;
    00532
    00533     na = LENGTH(a);
    00534     nb = LENGTH(b);
    00535     if (na != nb) error("a and b are of different length");
    00536
    00537     iwhich = LOGICAL(which);
    00538
    00539     for (i = 0; i < na; i++) {
    00540         if (iwhich[i]) continue;
    00541
    00542         ae = VECTOR_ELT(a, i);
    00543         be = VECTOR_ELT(b, i);
    00544
    00545         if (LENGTH(ae) != LENGTH(be))
    00546             error("elements %d are of different length", i);
    00547
    00548         if (!isReal(ae) || !isReal(be))
    00549             error("elements %d are not of type double", i);
    00550
    00551         dae = REAL(ae);
    00552         dbe = REAL(be);
    00553         for (j = 0; j < LENGTH(ae); j++)
    00554             dae[j] += dbe[j];
    00555     }
    00556     return(a);
    00557 }
    00558
    00559 SEXP R_modify_response(SEXP x, SEXP vf) {
    00560
    00561     double *src, *tar;
    00562     int i, n;
    00563
    00564     src = REAL(x);
    00565     n = LENGTH(x);
    00566
    00567     tar = REAL(get_transformation(vf, 1));
    00568     for (i = 0; i < n; i++)
    00569         tar[i] = src[i];
    00570
    00571     tar = REAL(get_jointtransf(vf));
    00572     for (i = 0; i < n; i++)
    00573         tar[i] = src[i];
    00574
    00575     tar = REAL(get_variable(vf, 1));
    00576     for (i = 0; i < n; i++)
    00577         tar[i] = src[i];
    00578
    00579     return(R_NilValue);
    00580 }