The R Project SVN R-packages

#include "sqlite_dataframe.h"
#include <math.h>

#define TOL 1e-12
#define ZERO 0.000000000000000000
#define EQUAL_ZERO(x) (abs((x)) < TOL)

void F77_NAME(includ)(int *, int *, double *, double *, double *,
              double *, double *, double *, double *, int *);

/* EQUAL_ZERO gives error */
static void __include(int np, int nrbar, double *D, double *rbar, double *thetab, 
        double *sserr, double *xrow, double y, double weight) {
    int nextr = 0, i, k;
    long double xi, xk, di, wxi, dpi, cbar, sbar;
    for (i = 0; i < np; i++) {
        if (weight == ZERO) return;
        xi = xrow[i];
        if (xi == ZERO) { nextr += np - i; continue; }
        di = D[i];
        wxi = weight * xi;
        dpi = di + wxi * wxi;
        cbar = di / dpi;
        sbar = wxi / dpi;
        weight *= cbar;
        D[i] = dpi;

        if (i < np - 1) {
            for (k = i+1; k < np; k++) {
                xk = xrow[k];
                xrow[k] = xk - xi * rbar[nextr];
                rbar[nextr] = cbar * rbar[nextr] + sbar * xk;
                nextr++;
            }
            xk = y;
            y = xk - xi * thetab[i];
            thetab[i] = cbar * thetab[i] + sbar * xk;
        }

        *sserr += weight * y * y;
    }
}


SEXP sdf_do_biglm(SEXP sdfx, SEXP svecy, SEXP sdfx_dim, SEXP intercept) {
    SEXP sdflm, qr, D, rbar, thetab, tmp;
    int np, nrbar, nprotect = 0, intrcpt, i, ier;
    double *xrow, y, *ptrD, *ptrrbar, *ptrthetab, sserr, weight;
    sqlite3_stmt *stmt1, *stmt2;
    char *inamex, *inamey, *tblnamey, *varnamey;

    inamex = SDF_INAME(sdfx);

    inamey = SDF_INAME(svecy);
    tblnamey = SVEC_TBLNAME(svecy);
    varnamey = SVEC_VARNAME(svecy);

    if (!USE_SDF1(inamex, TRUE, TRUE)) return R_NilValue;
    if (!USE_SDF1(inamey, TRUE, TRUE)) return R_NilValue;

    /* get number of columns for sdfx */
    np = INTEGER(sdfx_dim)[1]; 
    intrcpt = (LOGICAL(intercept)[0]) ? 1 : 0;
    np += intrcpt;
    nrbar = (np*(np-1))/2;

    /* Rprintf("np: %d, nrbar %d\n", np, nrbar); */

    PROTECT(D = NEW_NUMERIC(np)); nprotect++;
    PROTECT(rbar = NEW_NUMERIC(nrbar)); nprotect++;
    PROTECT(thetab = NEW_NUMERIC(np)); nprotect++;
    xrow = (double *)R_alloc(np, sizeof(double));

    /* initialize bigqr values */
    if (intrcpt) xrow[0] = 1.0;
    sserr = 0.0;
    ptrD = REAL(D); ptrrbar = REAL(rbar); ptrthetab = REAL(thetab);
    for (i = 0; i < np; i++) ptrD[i] = ptrrbar[i] = ptrthetab[i] = 0.0;
    for ( ; i < nrbar; i++) ptrrbar[i] = 0.0;

    /* setup sqlite_stmt's */
    sprintf(g_sql_buf[0], "select * from [%s].sdf_data", inamex);
    _sqlite_error(sqlite3_prepare(g_workspace, g_sql_buf[0], -1, &stmt1, NULL));
    sprintf(g_sql_buf[0], "select [%s] from [%s].[%s]", varnamey, inamey, tblnamey);
    _sqlite_error(sqlite3_prepare(g_workspace, g_sql_buf[0], -1, &stmt2, NULL));

    if (intrcpt) {
        while (sqlite3_step(stmt1) == SQLITE_ROW) {
            sqlite3_step(stmt2);
            for (i = 1; i < np; i++) 
                xrow[i] = sqlite3_column_double(stmt1, i);
            y = sqlite3_column_double(stmt2, 0);
            weight = 1.0;
            /*__include(np, nrbar, ptrD, ptrrbar, ptrthetab, &sserr, xrow, y, 1.0); */
            F77_CALL(includ)(&np, &nrbar, &weight, xrow, &y, ptrD, ptrrbar, ptrthetab, &sserr, &ier);
        }
    } else { 
        while (sqlite3_step(stmt1) == SQLITE_ROW) {
            sqlite3_step(stmt2);
            for (i = 0; i < np; i++)
                xrow[i] = sqlite3_column_double(stmt1, i+1);
            y = sqlite3_column_double(stmt2, 0);
            /*__include(np, nrbar, ptrD, ptrrbar, ptrthetab, &sserr, xrow, y, 1.0); */
            weight = 1;
            F77_CALL(includ)(&np, &nrbar, &weight, xrow, &y, ptrD, ptrrbar, ptrthetab, &sserr, &ier);
        }
    }

    sqlite3_finalize(stmt1);
    sqlite3_finalize(stmt2);

    /* setup bigqr SEXP */
    PROTECT(qr = NEW_LIST(6)); nprotect++;

    /* set names */
    PROTECT(tmp = NEW_CHARACTER(6)); nprotect++;
    SET_STRING_ELT(tmp, 0, mkChar("D"));
    SET_STRING_ELT(tmp, 1, mkChar("rbar"));
    SET_STRING_ELT(tmp, 2, mkChar("thetab"));
    SET_STRING_ELT(tmp, 3, mkChar("ss"));
    SET_STRING_ELT(tmp, 4, mkChar("checked"));
    SET_STRING_ELT(tmp, 5, mkChar("tol"));
    SET_NAMES(qr, tmp);

    /* set values */
    SET_VECTOR_ELT(qr, 0, D);
    SET_VECTOR_ELT(qr, 1, rbar);
    SET_VECTOR_ELT(qr, 2, thetab);
    SET_VECTOR_ELT(qr, 3, ScalarReal(sserr));
    SET_VECTOR_ELT(qr, 4, ScalarLogical(FALSE));
    PROTECT(tmp = NEW_NUMERIC(np)); nprotect++;
    for (i = 0; i < np; i++) REAL(tmp)[i] = 0.0;
    SET_VECTOR_ELT(qr, 5, tmp);

    /* set class */
    SET_CLASS(qr, mkString("bigqr"));

    /* setup sdflm SEXP */
    PROTECT(sdflm = NEW_LIST(3)); nprotect++;
    
    PROTECT(tmp = NEW_CHARACTER(3)); nprotect++;
    SET_STRING_ELT(tmp, 0, mkChar("qr"));
    SET_STRING_ELT(tmp, 1, mkChar("X"));
    SET_STRING_ELT(tmp, 2, mkChar("intercept"));
    SET_NAMES(sdflm, tmp);

    SET_VECTOR_ELT(sdflm, 0, qr);
    SET_VECTOR_ELT(sdflm, 1, duplicate(sdfx));
    SET_VECTOR_ELT(sdflm, 2, ScalarLogical(intrcpt));

    PROTECT(tmp = NEW_CHARACTER(2)); nprotect++;
    SET_STRING_ELT(tmp, 0, mkChar("sdflm"));
    SET_STRING_ELT(tmp, 1, mkChar("biglm"));
    SET_CLASS(sdflm, tmp);

    UNPROTECT(nprotect);
    UNUSE_SDF2(inamex);
    UNUSE_SDF2(inamey);
    return sdflm;
}