The R Project SVN – R-packages – Diff – /trunk/SQLiteDF/src/sqlite_biglm.c


#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;
    const 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;
}

Rev 3821	Rev 4798
1	`#include "sqlite_dataframe.h"`	1	`#include "sqlite_dataframe.h"`
2	`#include <math.h>`	2	`#include <math.h>`
3		3
4	`#define TOL 1e-12`	4	`#define TOL 1e-12`
5	`#define ZERO 0.000000000000000000`	5	`#define ZERO 0.000000000000000000`
6	`#define EQUAL_ZERO(x) (abs((x)) < TOL)`	6	`#define EQUAL_ZERO(x) (abs((x)) < TOL)`
7		7
8	`void F77_NAME(includ)(int , int , double , double , double *,`	8	`void F77_NAME(includ)(int , int , double , double , double *,`
9	`double , double , double , double , int *);`	9	`double , double , double , double , int *);`
10		10
11	`/* EQUAL_ZERO gives error */`	11	`/* EQUAL_ZERO gives error */`
12	`static void __include(int np, int nrbar, double D, double rbar, double *thetab,`	12	`static void __include(int np, int nrbar, double D, double rbar, double *thetab,`
13	`double sserr, double xrow, double y, double weight) {`	13	`double sserr, double xrow, double y, double weight) {`
14	`int nextr = 0, i, k;`	14	`int nextr = 0, i, k;`
15	`long double xi, xk, di, wxi, dpi, cbar, sbar;`	15	`long double xi, xk, di, wxi, dpi, cbar, sbar;`
16	`for (i = 0; i < np; i++) {`	16	`for (i = 0; i < np; i++) {`
17	`if (weight == ZERO) return;`	17	`if (weight == ZERO) return;`
18	`xi = xrow[i];`	18	`xi = xrow[i];`
19	`if (xi == ZERO) { nextr += np - i; continue; }`	19	`if (xi == ZERO) { nextr += np - i; continue; }`
20	`di = D[i];`	20	`di = D[i];`
21	`wxi = weight * xi;`	21	`wxi = weight * xi;`
22	`dpi = di + wxi * wxi;`	22	`dpi = di + wxi * wxi;`
23	`cbar = di / dpi;`	23	`cbar = di / dpi;`
24	`sbar = wxi / dpi;`	24	`sbar = wxi / dpi;`
25	`weight *= cbar;`	25	`weight *= cbar;`
26	`D[i] = dpi;`	26	`D[i] = dpi;`
27		27
28	`if (i < np - 1) {`	28	`if (i < np - 1) {`
29	`for (k = i+1; k < np; k++) {`	29	`for (k = i+1; k < np; k++) {`
30	`xk = xrow[k];`	30	`xk = xrow[k];`
31	`xrow[k] = xk - xi * rbar[nextr];`	31	`xrow[k] = xk - xi * rbar[nextr];`
32	`rbar[nextr] = cbar * rbar[nextr] + sbar * xk;`	32	`rbar[nextr] = cbar * rbar[nextr] + sbar * xk;`
33	`nextr++;`	33	`nextr++;`
34	`}`	34	`}`
35	`xk = y;`	35	`xk = y;`
36	`y = xk - xi * thetab[i];`	36	`y = xk - xi * thetab[i];`
37	`thetab[i] = cbar * thetab[i] + sbar * xk;`	37	`thetab[i] = cbar * thetab[i] + sbar * xk;`
38	`}`	38	`}`
39		39
40	`sserr += weight y * y;`	40	`sserr += weight y * y;`
41	`}`	41	`}`
42	`}`	42	`}`
43		43
44		44
45	`SEXP sdf_do_biglm(SEXP sdfx, SEXP svecy, SEXP sdfx_dim, SEXP intercept) {`	45	`SEXP sdf_do_biglm(SEXP sdfx, SEXP svecy, SEXP sdfx_dim, SEXP intercept) {`
46	`SEXP sdflm, qr, D, rbar, thetab, tmp;`	46	`SEXP sdflm, qr, D, rbar, thetab, tmp;`
47	`int np, nrbar, nprotect = 0, intrcpt, i, ier;`	47	`int np, nrbar, nprotect = 0, intrcpt, i, ier;`
48	`double xrow, y, ptrD, ptrrbar, ptrthetab, sserr, weight;`	48	`double xrow, y, ptrD, ptrrbar, ptrthetab, sserr, weight;`
49	`sqlite3_stmt stmt1, stmt2;`	49	`sqlite3_stmt stmt1, stmt2;`
50	`char inamex, inamey, tblnamey, varnamey;`	50	`const char inamex, inamey, tblnamey, varnamey;`
51		51
52	`inamex = SDF_INAME(sdfx);`	52	`inamex = SDF_INAME(sdfx);`
53		53
54	`inamey = SDF_INAME(svecy);`	54	`inamey = SDF_INAME(svecy);`
55	`tblnamey = SVEC_TBLNAME(svecy);`	55	`tblnamey = SVEC_TBLNAME(svecy);`
56	`varnamey = SVEC_VARNAME(svecy);`	56	`varnamey = SVEC_VARNAME(svecy);`
57		57
58	`if (!USE_SDF1(inamex, TRUE, TRUE)) return R_NilValue;`	58	`if (!USE_SDF1(inamex, TRUE, TRUE)) return R_NilValue;`
59	`if (!USE_SDF1(inamey, TRUE, TRUE)) return R_NilValue;`	59	`if (!USE_SDF1(inamey, TRUE, TRUE)) return R_NilValue;`
60		60
61	`/* get number of columns for sdfx */`	61	`/* get number of columns for sdfx */`
62	`np = INTEGER(sdfx_dim)[1];`	62	`np = INTEGER(sdfx_dim)[1];`
63	`intrcpt = (LOGICAL(intercept)[0]) ? 1 : 0;`	63	`intrcpt = (LOGICAL(intercept)[0]) ? 1 : 0;`
64	`np += intrcpt;`	64	`np += intrcpt;`
65	`nrbar = (np*(np-1))/2;`	65	`nrbar = (np*(np-1))/2;`
66		66
67	`/* Rprintf("np: %d, nrbar %d\n", np, nrbar); */`	67	`/* Rprintf("np: %d, nrbar %d\n", np, nrbar); */`
68		68
69	`PROTECT(D = NEW_NUMERIC(np)); nprotect++;`	69	`PROTECT(D = NEW_NUMERIC(np)); nprotect++;`
70	`PROTECT(rbar = NEW_NUMERIC(nrbar)); nprotect++;`	70	`PROTECT(rbar = NEW_NUMERIC(nrbar)); nprotect++;`
71	`PROTECT(thetab = NEW_NUMERIC(np)); nprotect++;`	71	`PROTECT(thetab = NEW_NUMERIC(np)); nprotect++;`
72	`xrow = (double *)R_alloc(np, sizeof(double));`	72	`xrow = (double *)R_alloc(np, sizeof(double));`
73		73
74	`/* initialize bigqr values */`	74	`/* initialize bigqr values */`
75	`if (intrcpt) xrow[0] = 1.0;`	75	`if (intrcpt) xrow[0] = 1.0;`
76	`sserr = 0.0;`	76	`sserr = 0.0;`
77	`ptrD = REAL(D); ptrrbar = REAL(rbar); ptrthetab = REAL(thetab);`	77	`ptrD = REAL(D); ptrrbar = REAL(rbar); ptrthetab = REAL(thetab);`
78	`for (i = 0; i < np; i++) ptrD[i] = ptrrbar[i] = ptrthetab[i] = 0.0;`	78	`for (i = 0; i < np; i++) ptrD[i] = ptrrbar[i] = ptrthetab[i] = 0.0;`
79	`for ( ; i < nrbar; i++) ptrrbar[i] = 0.0;`	79	`for ( ; i < nrbar; i++) ptrrbar[i] = 0.0;`
80		80
81	`/* setup sqlite_stmt's */`	81	`/* setup sqlite_stmt's */`
82	`sprintf(g_sql_buf[0], "select * from [%s].sdf_data", inamex);`	82	`sprintf(g_sql_buf[0], "select * from [%s].sdf_data", inamex);`
83	`_sqlite_error(sqlite3_prepare(g_workspace, g_sql_buf[0], -1, &stmt1, NULL));`	83	`_sqlite_error(sqlite3_prepare(g_workspace, g_sql_buf[0], -1, &stmt1, NULL));`
84	`sprintf(g_sql_buf[0], "select [%s] from [%s].[%s]", varnamey, inamey, tblnamey);`	84	`sprintf(g_sql_buf[0], "select [%s] from [%s].[%s]", varnamey, inamey, tblnamey);`
85	`_sqlite_error(sqlite3_prepare(g_workspace, g_sql_buf[0], -1, &stmt2, NULL));`	85	`_sqlite_error(sqlite3_prepare(g_workspace, g_sql_buf[0], -1, &stmt2, NULL));`
86		86
87	`if (intrcpt) {`	87	`if (intrcpt) {`
88	`while (sqlite3_step(stmt1) == SQLITE_ROW) {`	88	`while (sqlite3_step(stmt1) == SQLITE_ROW) {`
89	`sqlite3_step(stmt2);`	89	`sqlite3_step(stmt2);`
90	`for (i = 1; i < np; i++)`	90	`for (i = 1; i < np; i++)`
91	`xrow[i] = sqlite3_column_double(stmt1, i);`	91	`xrow[i] = sqlite3_column_double(stmt1, i);`
92	`y = sqlite3_column_double(stmt2, 0);`	92	`y = sqlite3_column_double(stmt2, 0);`
93	`weight = 1.0;`	93	`weight = 1.0;`
94	`/__include(np, nrbar, ptrD, ptrrbar, ptrthetab, &sserr, xrow, y, 1.0); /`	94	`/__include(np, nrbar, ptrD, ptrrbar, ptrthetab, &sserr, xrow, y, 1.0); /`
95	`F77_CALL(includ)(&np, &nrbar, &weight, xrow, &y, ptrD, ptrrbar, ptrthetab, &sserr, &ier);`	95	`F77_CALL(includ)(&np, &nrbar, &weight, xrow, &y, ptrD, ptrrbar, ptrthetab, &sserr, &ier);`
96	`}`	96	`}`
97	`} else {`	97	`} else {`
98	`while (sqlite3_step(stmt1) == SQLITE_ROW) {`	98	`while (sqlite3_step(stmt1) == SQLITE_ROW) {`
99	`sqlite3_step(stmt2);`	99	`sqlite3_step(stmt2);`
100	`for (i = 0; i < np; i++)`	100	`for (i = 0; i < np; i++)`
101	`xrow[i] = sqlite3_column_double(stmt1, i+1);`	101	`xrow[i] = sqlite3_column_double(stmt1, i+1);`
102	`y = sqlite3_column_double(stmt2, 0);`	102	`y = sqlite3_column_double(stmt2, 0);`
103	`/__include(np, nrbar, ptrD, ptrrbar, ptrthetab, &sserr, xrow, y, 1.0); /`	103	`/__include(np, nrbar, ptrD, ptrrbar, ptrthetab, &sserr, xrow, y, 1.0); /`
104	`weight = 1;`	104	`weight = 1;`
105	`F77_CALL(includ)(&np, &nrbar, &weight, xrow, &y, ptrD, ptrrbar, ptrthetab, &sserr, &ier);`	105	`F77_CALL(includ)(&np, &nrbar, &weight, xrow, &y, ptrD, ptrrbar, ptrthetab, &sserr, &ier);`
106	`}`	106	`}`
107	`}`	107	`}`
108		108
109	`sqlite3_finalize(stmt1);`	109	`sqlite3_finalize(stmt1);`
110	`sqlite3_finalize(stmt2);`	110	`sqlite3_finalize(stmt2);`
111		111
112	`/* setup bigqr SEXP */`	112	`/* setup bigqr SEXP */`
113	`PROTECT(qr = NEW_LIST(6)); nprotect++;`	113	`PROTECT(qr = NEW_LIST(6)); nprotect++;`
114		114
115	`/* set names */`	115	`/* set names */`
116	`PROTECT(tmp = NEW_CHARACTER(6)); nprotect++;`	116	`PROTECT(tmp = NEW_CHARACTER(6)); nprotect++;`
117	`SET_STRING_ELT(tmp, 0, mkChar("D"));`	117	`SET_STRING_ELT(tmp, 0, mkChar("D"));`
118	`SET_STRING_ELT(tmp, 1, mkChar("rbar"));`	118	`SET_STRING_ELT(tmp, 1, mkChar("rbar"));`
119	`SET_STRING_ELT(tmp, 2, mkChar("thetab"));`	119	`SET_STRING_ELT(tmp, 2, mkChar("thetab"));`
120	`SET_STRING_ELT(tmp, 3, mkChar("ss"));`	120	`SET_STRING_ELT(tmp, 3, mkChar("ss"));`
121	`SET_STRING_ELT(tmp, 4, mkChar("checked"));`	121	`SET_STRING_ELT(tmp, 4, mkChar("checked"));`
122	`SET_STRING_ELT(tmp, 5, mkChar("tol"));`	122	`SET_STRING_ELT(tmp, 5, mkChar("tol"));`
123	`SET_NAMES(qr, tmp);`	123	`SET_NAMES(qr, tmp);`
124		124
125	`/* set values */`	125	`/* set values */`
126	`SET_VECTOR_ELT(qr, 0, D);`	126	`SET_VECTOR_ELT(qr, 0, D);`
127	`SET_VECTOR_ELT(qr, 1, rbar);`	127	`SET_VECTOR_ELT(qr, 1, rbar);`
128	`SET_VECTOR_ELT(qr, 2, thetab);`	128	`SET_VECTOR_ELT(qr, 2, thetab);`
129	`SET_VECTOR_ELT(qr, 3, ScalarReal(sserr));`	129	`SET_VECTOR_ELT(qr, 3, ScalarReal(sserr));`
130	`SET_VECTOR_ELT(qr, 4, ScalarLogical(FALSE));`	130	`SET_VECTOR_ELT(qr, 4, ScalarLogical(FALSE));`
131	`PROTECT(tmp = NEW_NUMERIC(np)); nprotect++;`	131	`PROTECT(tmp = NEW_NUMERIC(np)); nprotect++;`
132	`for (i = 0; i < np; i++) REAL(tmp)[i] = 0.0;`	132	`for (i = 0; i < np; i++) REAL(tmp)[i] = 0.0;`
133	`SET_VECTOR_ELT(qr, 5, tmp);`	133	`SET_VECTOR_ELT(qr, 5, tmp);`
134		134
135	`/* set class */`	135	`/* set class */`
136	`SET_CLASS(qr, mkString("bigqr"));`	136	`SET_CLASS(qr, mkString("bigqr"));`
137		137
138	`/* setup sdflm SEXP */`	138	`/* setup sdflm SEXP */`
139	`PROTECT(sdflm = NEW_LIST(3)); nprotect++;`	139	`PROTECT(sdflm = NEW_LIST(3)); nprotect++;`
140		140
141	`PROTECT(tmp = NEW_CHARACTER(3)); nprotect++;`	141	`PROTECT(tmp = NEW_CHARACTER(3)); nprotect++;`
142	`SET_STRING_ELT(tmp, 0, mkChar("qr"));`	142	`SET_STRING_ELT(tmp, 0, mkChar("qr"));`
143	`SET_STRING_ELT(tmp, 1, mkChar("X"));`	143	`SET_STRING_ELT(tmp, 1, mkChar("X"));`
144	`SET_STRING_ELT(tmp, 2, mkChar("intercept"));`	144	`SET_STRING_ELT(tmp, 2, mkChar("intercept"));`
145	`SET_NAMES(sdflm, tmp);`	145	`SET_NAMES(sdflm, tmp);`
146		146
147	`SET_VECTOR_ELT(sdflm, 0, qr);`	147	`SET_VECTOR_ELT(sdflm, 0, qr);`
148	`SET_VECTOR_ELT(sdflm, 1, duplicate(sdfx));`	148	`SET_VECTOR_ELT(sdflm, 1, duplicate(sdfx));`
149	`SET_VECTOR_ELT(sdflm, 2, ScalarLogical(intrcpt));`	149	`SET_VECTOR_ELT(sdflm, 2, ScalarLogical(intrcpt));`
150		150
151	`PROTECT(tmp = NEW_CHARACTER(2)); nprotect++;`	151	`PROTECT(tmp = NEW_CHARACTER(2)); nprotect++;`
152	`SET_STRING_ELT(tmp, 0, mkChar("sdflm"));`	152	`SET_STRING_ELT(tmp, 0, mkChar("sdflm"));`
153	`SET_STRING_ELT(tmp, 1, mkChar("biglm"));`	153	`SET_STRING_ELT(tmp, 1, mkChar("biglm"));`
154	`SET_CLASS(sdflm, tmp);`	154	`SET_CLASS(sdflm, tmp);`
155		155
156	`UNPROTECT(nprotect);`	156	`UNPROTECT(nprotect);`
157	`UNUSE_SDF2(inamex);`	157	`UNUSE_SDF2(inamex);`
158	`UNUSE_SDF2(inamey);`	158	`UNUSE_SDF2(inamey);`
159	`return sdflm;`	159	`return sdflm;`
160	`}`	160	`}`

The R Project SVN R-packages

(root)/trunk/SQLiteDF/src/sqlite_biglm.c – Rev 3821 → 4798