The R Project SVN R

Rev

Details | Last modification | View Log | RSS feed

Rev Author Line No. Line
2 r 1 
/*
1160 maechler 2 
 *  R : A Computer Language for Statistical Data Analysis
2 r 3 
 *  Copyright (C) 1995, 1996  Robert Gentleman and Ross Ihaka
22903 ripley 4 
 *  Copyright (C) 1997-2003   Robert Gentleman, Ross Ihaka and the
4562 pd 5 
 *                            R Development Core Team
2 r 6 
 *
7 
 *  This program is free software; you can redistribute it and/or modify
8 
 *  it under the terms of the GNU General Public License as published by
9 
 *  the Free Software Foundation; either version 2 of the License, or
10 
 *  (at your option) any later version.
11 
 *
12 
 *  This program is distributed in the hope that it will be useful,
13 
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
14 
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15 
 *  GNU General Public License for more details.
16 
 *
17 
 *  You should have received a copy of the GNU General Public License
18 
 *  along with this program; if not, write to the Free Software
5458 ripley 19 
 *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
1895 ihaka 20 
 *
21 
 *
22 
 *  Model Formula Manipulation
23 
 *
24 
 *  Can you say ``recurse your brains out'';
25 
 *  I knew you could. -- Mr Ro(ss)gers
2 r 26 
 */
27 
 
5187 hornik 28 
#ifdef HAVE_CONFIG_H
7701 hornik 29 
#include <config.h>
5187 hornik 30 
#endif
31 
 
2 r 32 
#include "Defn.h"
33 
 
34 
#define WORDSIZE (8*sizeof(int))
35 
 
36 
static SEXP tildeSymbol = NULL;
37 
static SEXP plusSymbol  = NULL;
38 
static SEXP minusSymbol = NULL;
39 
static SEXP timesSymbol = NULL;
40 
static SEXP slashSymbol = NULL;
41 
static SEXP colonSymbol = NULL;
42 
static SEXP powerSymbol = NULL;
43 
static SEXP dotSymbol   = NULL;
44 
static SEXP parenSymbol = NULL;
45 
static SEXP inSymbol    = NULL;
12976 pd 46 
/* unused static SEXP identSymbol = NULL; */
2 r 47 
 
48 
 
49 
static int intercept;		/* intercept term in the model */
50 
static int parity;		/* +/- parity */
51 
static int response;		/* response term in the model */
52 
static int nvar;		/* Number of variables in the formula */
53 
static int nwords;		/* # of words (ints) to code a term */
54 
static int nterm;		/* # of model terms */
55 
static SEXP varlist;		/* variables in the model */
2754 pd 56 
SEXP framenames;		/* variables names for specified frame */
57 
/* NOTE: framenames can't be static because it must be protected from
58 
   garbage collection. */
22903 ripley 59 
static Rboolean haveDot;	/* does RHS of formula contain `.'? */
2 r 60 
 
92 ihaka 61 
static int isZeroOne(SEXP x)
62 
{
1839 ihaka 63 
    if (!isNumeric(x)) return 0;
64 
    return (asReal(x) == 0.0 || asReal(x) == 1.0);
92 ihaka 65 
}
66 
 
67 
static int isZero(SEXP x)
68 
{
1839 ihaka 69 
    if (!isNumeric(x)) return 0;
70 
    return asReal(x) == 0.0;
92 ihaka 71 
}
72 
 
2 r 73 
static int isOne(SEXP x)
74 
{
1839 ihaka 75 
    if (!isNumeric(x)) return 0;
76 
    return asReal(x) == 1.0;
2 r 77 
}
78 
 
79 
 
1839 ihaka 80 
/* MatchVar determines whether two ``variables'' are */
81 
/* identical.  Expressions are identical if they have */
82 
/* the same list structure and their atoms are identical. */
83 
/* This is just EQUAL from lisp. */
84 
 
2 r 85 
static int MatchVar(SEXP var1, SEXP var2)
86 
{
2819 pd 87 
    /* For expedience, and sanity... */
88 
    if ( var1 == var2 )
89 
    	return 1;
1839 ihaka 90 
    /* Handle Nulls */
91 
    if (isNull(var1) && isNull(var2))
92 
	return 1;
93 
    if (isNull(var1) || isNull(var2))
2 r 94 
	return 0;
1839 ihaka 95 
    /* Non-atomic objects - compare CARs & CDRs */
96 
    if ((isList(var1) || isLanguage(var1)) &&
97 
	(isList(var2) || isLanguage(var2)))
98 
	return MatchVar(CAR(var1), CAR(var2)) &&
99 
	       MatchVar(CDR(var1), CDR(var2));
100 
    /* Symbols */
101 
    if (isSymbol(var1) && isSymbol(var2))
102 
	return (var1 == var2);
103 
    /* Literal Numerics */
104 
    if (isNumeric(var1) && isNumeric(var2))
105 
	return (asReal(var1) == asReal(var2));
2819 pd 106 
    /* Literal Strings */
107 
    if (isString(var1) && isString(var2))
22906 ripley 108 
	return (strcmp(CHAR(STRING_ELT(var1, 0)),
109 
		       CHAR(STRING_ELT(var2, 0))) == 0);
1839 ihaka 110 
    /* Nothing else matches */
111 
    return 0;
2 r 112 
}
113 
 
114 
 
1839 ihaka 115 
/* InstallVar locates a ``variable'' in the model */
116 
/* variable list;  adding it to the list if not found. */
2 r 117 
 
118 
static int InstallVar(SEXP var)
119 
{
1839 ihaka 120 
    SEXP v;
13257 maechler 121 
    int indx;
1839 ihaka 122 
    /* Check that variable is legitimate */
123 
    if (!isSymbol(var) && !isLanguage(var) && !isZeroOne(var))
5731 ripley 124 
	error("invalid term in model formula");
1839 ihaka 125 
    /* Lookup/Install it */
13257 maechler 126 
    indx = 0;
1839 ihaka 127 
    for (v = varlist; CDR(v) != R_NilValue; v = CDR(v)) {
13257 maechler 128 
	indx++;
1839 ihaka 129 
	if (MatchVar(var, CADR(v)))
13257 maechler 130 
	    return indx;
1839 ihaka 131 
    }
10172 luke 132 
    SETCDR(v, CONS(var, R_NilValue));
13257 maechler 133 
    return indx + 1;
1839 ihaka 134 
}
2 r 135 
 
136 
 
1839 ihaka 137 
/* If there is a dotsxp being expanded then we need to see */
138 
/* whether any of the variables in the data frame match with */
139 
/* the variable on the lhs. If so they shouldn't be included */
140 
/* in the factors */
2 r 141 
 
7758 ripley 142 
static void CheckRHS(SEXP v)
2 r 143 
{
1839 ihaka 144 
    int i, j;
145 
    SEXP s, t;
146 
    while ((isList(v) || isLanguage(v)) && v != R_NilValue) {
147 
	CheckRHS(CAR(v));
148 
	v = CDR(v);
149 
    }
150 
    if (isSymbol(v)) {
22906 ripley 151 
	for (i = 0; i < length(framenames); i++) {
152 
	    s = install(CHAR(STRING_ELT(framenames, i)));
1839 ihaka 153 
	    if (v == s) {
154 
		t=allocVector(STRSXP, length(framenames)-1);
22906 ripley 155 
		for (j = 0; j < length(t); j++) {
156 
		    if (j < i)
10172 luke 157 
			SET_STRING_ELT(t, j, STRING_ELT(framenames, j));
1839 ihaka 158 
		    else
10172 luke 159 
			SET_STRING_ELT(t, j, STRING_ELT(framenames, j+1));
2 r 160 
		}
1839 ihaka 161 
		framenames=t;
162 
	    }
2 r 163 
	}
1839 ihaka 164 
    }
2 r 165 
}
166 
 
167 
 
1839 ihaka 168 
/* ExtractVars recursively extracts the variables */
169 
/* in a model formula.  It calls InstallVar to do */
170 
/* the installation.  The code takes care of unary */
171 
/* + and minus.  No checks are made of the other */
172 
/* ``binary'' operators.  Maybe there should be some. */
2 r 173 
 
174 
static void ExtractVars(SEXP formula, int checkonly)
175 
{
1839 ihaka 176 
    int len, i;
177 
    SEXP v;
24985 ripley 178 
 
1839 ihaka 179 
    if (isNull(formula) || isZeroOne(formula))
180 
	return;
25912 ripley 181 
    if (isSymbol(formula)) {
1839 ihaka 182 
	if (!checkonly) {
22903 ripley 183 
	    if (formula == dotSymbol && framenames != R_NilValue) {
184 
		haveDot = TRUE;
185 
		for (i = 0; i < length(framenames); i++) {
186 
		    v = install(CHAR(STRING_ELT(framenames, i)));
20316 ripley 187 
		    if (!MatchVar(v, CADR(varlist))) InstallVar(v);
2 r 188 
		}
22903 ripley 189 
	    } else
2 r 190 
		InstallVar(formula);
191 
	}
1839 ihaka 192 
	return;
193 
    }
194 
    if (isLanguage(formula)) {
195 
	len = length(formula);
196 
	if (CAR(formula) == tildeSymbol) {
197 
	    if (response)
5731 ripley 198 
		error("invalid model formula");
1839 ihaka 199 
	    if (isNull(CDDR(formula))) {
200 
		response = 0;
201 
		ExtractVars(CADR(formula), 0);
202 
	    }
203 
	    else {
204 
		response = 1;
205 
		InstallVar(CADR(formula));
206 
		ExtractVars(CADDR(formula), 0);
207 
	    }
208 
	    return;
209 
	}
210 
	if (CAR(formula) == plusSymbol) {
211 
	    if (length(formula) > 1)
212 
		ExtractVars(CADR(formula), checkonly);
213 
	    if (length(formula) > 2)
214 
		ExtractVars(CADDR(formula), checkonly);
215 
	    return;
216 
	}
217 
	if (CAR(formula) == colonSymbol) {
218 
	    ExtractVars(CADR(formula), checkonly);
219 
	    ExtractVars(CADDR(formula), checkonly);
220 
	    return;
221 
	}
222 
	if (CAR(formula) == powerSymbol) {
223 
	    if (!isNumeric(CADDR(formula)))
5731 ripley 224 
		error("invalid power in formula");
1839 ihaka 225 
	    ExtractVars(CADR(formula), checkonly);
226 
	    return;
227 
	}
228 
	if (CAR(formula) == timesSymbol) {
229 
	    ExtractVars(CADR(formula), checkonly);
230 
	    ExtractVars(CADDR(formula), checkonly);
231 
	    return;
232 
	}
233 
	if (CAR(formula) == inSymbol) {
234 
	    ExtractVars(CADR(formula), checkonly);
235 
	    ExtractVars(CADDR(formula), checkonly);
236 
	    return;
237 
	}
238 
	if (CAR(formula) == slashSymbol) {
239 
	    ExtractVars(CADR(formula), checkonly);
240 
	    ExtractVars(CADDR(formula), checkonly);
241 
	    return;
242 
	}
243 
	if (CAR(formula) == minusSymbol) {
244 
	    if (len == 2) {
245 
		ExtractVars(CADR(formula), 1);
246 
	    }
247 
	    else {
248 
		ExtractVars(CADR(formula), checkonly);
249 
		ExtractVars(CADDR(formula), 1);
250 
	    }
251 
	    return;
252 
	}
253 
	if (CAR(formula) == parenSymbol) {
254 
	    ExtractVars(CADR(formula), checkonly);
255 
	    return;
256 
	}
257 
	InstallVar(formula);
258 
	return;
259 
    }
24985 ripley 260 
    error("invalid model formula in ExtractVars");
2 r 261 
}
262 
 
263 
 
1839 ihaka 264 
/* AllocTerm allocates an integer array for */
265 
/* bit string representation of a model term */
2 r 266 
 
267 
static SEXP AllocTerm()
268 
{
1839 ihaka 269 
    int i;
270 
    SEXP term = allocVector(INTSXP, nwords);
271 
    for (i = 0; i < nwords; i++)
272 
	INTEGER(term)[i] = 0;
273 
    return term;
2 r 274 
}
275 
 
276 
 
1839 ihaka 277 
/* SetBit sets bit ``whichBit'' to value ``value'' */
278 
/* in the bit string representation of a term. */
2 r 279 
 
280 
static void SetBit(SEXP term, int whichBit, int value)
281 
{
1839 ihaka 282 
    int word, offset;
283 
    word = (whichBit - 1) / WORDSIZE;
284 
    offset = (WORDSIZE - whichBit) % WORDSIZE;
285 
    if (value)
286 
	((unsigned *) INTEGER(term))[word] |= ((unsigned) 1 << offset);
287 
    else
288 
	((unsigned *) INTEGER(term))[word] &= ~((unsigned) 1 << offset);
2 r 289 
}
290 
 
291 
 
1839 ihaka 292 
/* GetBit gets bit ``whichBit'' from the */
293 
/* bit string representation of a term. */
2 r 294 
 
295 
static int GetBit(SEXP term, int whichBit)
296 
{
1839 ihaka 297 
    unsigned int word, offset;
298 
    word = (whichBit - 1) / WORDSIZE;
299 
    offset = (WORDSIZE - whichBit) % WORDSIZE;
300 
    return ((((unsigned *) INTEGER(term))[word]) >> offset) & 1;
2 r 301 
}
302 
 
303 
 
1839 ihaka 304 
/* OrBits computes a new (bit string) term */
305 
/* which contains the logical OR of the bits */
306 
/* in ``term1'' and ``term2''. */
2 r 307 
 
308 
static SEXP OrBits(SEXP term1, SEXP term2)
309 
{
1839 ihaka 310 
    SEXP term;
311 
    int i;
312 
    term = AllocTerm();
313 
    for (i = 0; i < nwords; i++)
314 
	INTEGER(term)[i] = INTEGER(term1)[i] | INTEGER(term2)[i];
315 
    return term;
2 r 316 
}
317 
 
318 
 
1839 ihaka 319 
/* BitCount counts the number of ``on'' */
320 
/* bits in a term */
321 
 
2 r 322 
static int BitCount(SEXP term)
323 
{
1839 ihaka 324 
    int i, sum;
325 
    sum = 0;
326 
    for (i = 1; i <= nvar; i++)
327 
	sum += GetBit(term, i);
328 
    return sum;
2 r 329 
}
330 
 
331 
 
1839 ihaka 332 
/* TermZero tests whether a (bit string) term is zero */
2 r 333 
 
334 
static int TermZero(SEXP term)
335 
{
1839 ihaka 336 
    int i, val;
337 
    val = 1;
338 
    for (i = 0; i < nwords; i++)
339 
	val = val && (INTEGER(term)[i] == 0);
340 
    return val;
2 r 341 
}
342 
 
343 
 
1839 ihaka 344 
/* TermEqual tests two (bit string) terms for equality. */
2 r 345 
 
346 
static int TermEqual(SEXP term1, SEXP term2)
347 
{
1839 ihaka 348 
    int i, val;
349 
    val = 1;
350 
    for (i = 0; i < nwords; i++)
351 
	val = val && (INTEGER(term1)[i] == INTEGER(term2)[i]);
352 
    return val;
2 r 353 
}
354 
 
355 
 
1839 ihaka 356 
/* StripTerm strips the specified term from */
357 
/* the given list.  This mutates the list. */
2 r 358 
 
359 
static SEXP StripTerm(SEXP term, SEXP list)
360 
{
1839 ihaka 361 
    SEXP tail;
362 
    if (TermZero(term))
363 
	intercept = 0;
364 
    if (list == R_NilValue)
2 r 365 
	return list;
1839 ihaka 366 
    tail = StripTerm(term, CDR(list));
367 
    if (TermEqual(term, CAR(list)))
368 
	return tail;
10172 luke 369 
    SETCDR(list, tail);
1839 ihaka 370 
    return list;
2 r 371 
}
372 
 
373 
 
1839 ihaka 374 
/* TrimRepeats removes duplicates of (bit string) terms */
375 
/* in a model formula by repeated use of ``StripTerm''. */
376 
/* Also drops zero terms. */
2 r 377 
 
378 
static SEXP TrimRepeats(SEXP list)
379 
{
1839 ihaka 380 
    if (list == R_NilValue)
381 
	return R_NilValue;
382 
    if (TermZero(CAR(list)))
383 
	return TrimRepeats(CDR(list));
10172 luke 384 
    SETCDR(list, TrimRepeats(StripTerm(CAR(list), CDR(list))));
1839 ihaka 385 
    return list;
2 r 386 
}
1839 ihaka 387 
 
388 
 
454 maechler 389 
 
390 
/*==========================================================================*/
2 r 391 
 
1839 ihaka 392 
/* Model Formula Manipulation */
393 
/* These functions take a numerically coded */
394 
/* formula and fully expand it. */
2 r 395 
 
454 maechler 396 
static SEXP EncodeVars(SEXP);/* defined below */
2 r 397 
 
398 
 
1839 ihaka 399 
/* PlusTerms expands ``left'' and ``right'' and */
400 
/* concatenates their terms (removing duplicates). */
2 r 401 
 
402 
static SEXP PlusTerms(SEXP left, SEXP right)
403 
{
1839 ihaka 404 
    PROTECT(left = EncodeVars(left));
405 
    right = EncodeVars(right);
406 
    UNPROTECT(1);
407 
    return TrimRepeats(listAppend(left, right));
2 r 408 
}
409 
 
410 
 
1839 ihaka 411 
/* InteractTerms expands ``left'' and ``right'' */
412 
/* and forms a new list of terms containing the bitwise */
413 
/* OR of each term in ``left'' with each term in ``right''. */
2 r 414 
 
415 
static SEXP InteractTerms(SEXP left, SEXP right)
416 
{
1839 ihaka 417 
    SEXP term, l, r, t;
418 
    PROTECT(left = EncodeVars(left));
419 
    PROTECT(right = EncodeVars(right));
420 
    PROTECT(term = allocList(length(left) * length(right)));
421 
    t = term;
422 
    for (l = left; l != R_NilValue; l = CDR(l))
423 
	for (r = right; r != R_NilValue; r = CDR(r)) {
10172 luke 424 
	    SETCAR(t, OrBits(CAR(l), CAR(r)));
1839 ihaka 425 
	    t = CDR(t);
426 
	}
427 
    UNPROTECT(3);
428 
    return TrimRepeats(term);
2 r 429 
}
430 
 
431 
 
1839 ihaka 432 
/* CrossTerms expands ``left'' and ``right'' */
433 
/* and forms the ``cross'' of the list of terms.  */
434 
/* Duplicates are removed. */
2 r 435 
 
436 
static SEXP CrossTerms(SEXP left, SEXP right)
437 
{
1839 ihaka 438 
    SEXP term, l, r, t;
439 
    PROTECT(left = EncodeVars(left));
440 
    PROTECT(right = EncodeVars(right));
441 
    PROTECT(term = allocList(length(left) * length(right)));
442 
    t = term;
443 
    for (l = left; l != R_NilValue; l = CDR(l))
444 
	for (r = right; r != R_NilValue; r = CDR(r)) {
10172 luke 445 
	    SETCAR(t, OrBits(CAR(l), CAR(r)));
1839 ihaka 446 
	    t = CDR(t);
447 
	}
448 
    UNPROTECT(3);
449 
    listAppend(right, term);
450 
    listAppend(left, right);
451 
    return TrimRepeats(left);
2 r 452 
}
453 
 
454 
 
1839 ihaka 455 
/* PowerTerms expands the ``left'' form and then */
456 
/* raises it to the power specified by the right term. */
457 
/* Allocation here is wasteful, but so what ... */
2 r 458 
 
459 
static SEXP PowerTerms(SEXP left, SEXP right)
460 
{
1839 ihaka 461 
    SEXP term, l, r, t;
13607 maechler 462 
    int i, ip;
463 
    ip = asInteger(right);
464 
    if (ip==NA_INTEGER || ip <= 1)
5731 ripley 465 
	error("Invalid power in formula");
2460 hornik 466 
    term = R_NilValue;		/* -Wall */
1839 ihaka 467 
    PROTECT(left = EncodeVars(left));
468 
    right = left;
13607 maechler 469 
    for (i=1; i < ip; i++)  {
1839 ihaka 470 
	PROTECT(right);
471 
	PROTECT(term = allocList(length(left) * length(right)));
472 
	t = term;
473 
	for (l = left; l != R_NilValue; l = CDR(l))
474 
	    for (r = right; r != R_NilValue; r = CDR(r)) {
10172 luke 475 
		SETCAR(t, OrBits(CAR(l), CAR(r)));
1839 ihaka 476 
		t = CDR(t);
477 
	    }
478 
	UNPROTECT(2);
479 
	right = TrimRepeats(term);
480 
    }
481 
    UNPROTECT(1);
482 
    return term;
2 r 483 
}
484 
 
485 
 
1839 ihaka 486 
/* InTerms expands ``left'' and ``right'' and */
487 
/* forms the ``nest'' of the the left in the */
488 
/* interaction of the right */
2 r 489 
 
490 
static SEXP InTerms(SEXP left, SEXP right)
491 
{
1839 ihaka 492 
    SEXP term, t;
493 
    int i;
494 
    PROTECT(left = EncodeVars(left));
495 
    PROTECT(right = EncodeVars(right));
496 
    PROTECT(term = AllocTerm());
497 
    /* Bitwise or of all terms on right */
498 
    for (t = right; t != R_NilValue; t = CDR(t)) {
499 
	for (i = 0; i < nwords; i++)
500 
	    INTEGER(term)[i] = INTEGER(term)[i] | INTEGER(CAR(t))[i];
501 
    }
502 
    /* Now bitwise or with each term on the left */
503 
    for (t = left; t != R_NilValue; t = CDR(t))
504 
	for (i = 0; i < nwords; i++)
505 
	    INTEGER(CAR(t))[i] = INTEGER(term)[i] | INTEGER(CAR(t))[i];
506 
    UNPROTECT(3);
507 
    return TrimRepeats(left);
2 r 508 
}
509 
 
1839 ihaka 510 
/* NestTerms expands ``left'' and ``right'' */
511 
/* and forms the ``nest'' of the list of terms.  */
512 
/* Duplicates are removed. */
2 r 513 
 
514 
static SEXP NestTerms(SEXP left, SEXP right)
515 
{
1839 ihaka 516 
    SEXP term, t;
517 
    int i;
518 
    PROTECT(left = EncodeVars(left));
519 
    PROTECT(right = EncodeVars(right));
520 
    PROTECT(term = AllocTerm());
521 
    /* Bitwise or of all terms on left */
522 
    for (t = left; t != R_NilValue; t = CDR(t)) {
523 
	for (i = 0; i < nwords; i++)
524 
	    INTEGER(term)[i] = INTEGER(term)[i] | INTEGER(CAR(t))[i];
525 
    }
526 
    /* Now bitwise or with each term on the right */
527 
    for (t = right; t != R_NilValue; t = CDR(t))
528 
	for (i = 0; i < nwords; i++)
529 
	    INTEGER(CAR(t))[i] = INTEGER(term)[i] | INTEGER(CAR(t))[i];
530 
    UNPROTECT(3);
531 
    listAppend(left, right);
532 
    return TrimRepeats(left);
2 r 533 
}
534 
 
535 
 
1839 ihaka 536 
/* DeleteTerms expands ``left'' and ``right'' */
537 
/* and then removes any terms which appear in */
538 
/* ``right'' from ``left''. */
2 r 539 
 
540 
static SEXP DeleteTerms(SEXP left, SEXP right)
541 
{
1839 ihaka 542 
    SEXP t;
543 
    PROTECT(left = EncodeVars(left));
544 
    parity = 1-parity;
545 
    PROTECT(right = EncodeVars(right));
546 
    parity = 1-parity;
547 
    for (t = right; t != R_NilValue; t = CDR(t))
548 
	left = StripTerm(CAR(t), left);
549 
    UNPROTECT(2);
550 
    return left;
2 r 551 
}
552 
 
553 
 
1839 ihaka 554 
/* EncodeVars performs  model expansion and bit string encoding. */
555 
/* This is the real workhorse of model expansion. */
2 r 556 
 
557 
static SEXP EncodeVars(SEXP formula)
558 
{
22896 ripley 559 
    SEXP term;
560 
    int len;
2 r 561 
 
1839 ihaka 562 
    if (isNull(formula))
563 
	return R_NilValue;
2 r 564 
 
1839 ihaka 565 
    if (isOne(formula)) {
566 
	if (parity) intercept = 1;
567 
	else intercept = 0;
568 
	return R_NilValue;
569 
    }
570 
    else if (isZero(formula)) {
571 
	if (parity) intercept = 0;
572 
	else intercept = 1;
573 
	return R_NilValue;
574 
    }
25912 ripley 575 
    if (isSymbol(formula)) {
1839 ihaka 576 
	if (formula == dotSymbol && framenames != R_NilValue) {
22896 ripley 577 
	    /* prior to 1.7.0 this made term.labels in reverse order. */
578 
	    SEXP r = R_NilValue, v = R_NilValue; /* -Wall */
579 
	    int i, j; char *c;
580 
 
581 
	    if (!LENGTH(framenames)) return r;
582 
	    for (i = 0; i < LENGTH(framenames); i++) {
20316 ripley 583 
		/* change in 1.6.0 do not use duplicated names */
584 
		c = CHAR(STRING_ELT(framenames, i));
585 
		for(j = 0; j < i; j++)
586 
		    if(!strcmp(c, CHAR(STRING_ELT(framenames, j))))
22896 ripley 587 
			error("duplicated name `%s' in data frame using `.'",
588 
			      c);
1839 ihaka 589 
		term = AllocTerm();
22896 ripley 590 
		SetBit(term, InstallVar(install(c)), 1);
591 
		if(i == 0) PROTECT(v = r = cons(term, R_NilValue));
592 
		else {SETCDR(v, CONS(term, R_NilValue)); v = CDR(v);}
1839 ihaka 593 
	    }
22896 ripley 594 
	    UNPROTECT(1);
1839 ihaka 595 
	    return r;
2 r 596 
	}
1839 ihaka 597 
	else {
598 
	    term = AllocTerm();
599 
	    SetBit(term, InstallVar(formula), 1);
600 
	    return CONS(term, R_NilValue);
92 ihaka 601 
	}
1839 ihaka 602 
    }
603 
    if (isLanguage(formula)) {
604 
	len = length(formula);
605 
	if (CAR(formula) == tildeSymbol) {
606 
	    if (isNull(CDDR(formula)))
607 
		return EncodeVars(CADR(formula));
608 
	    else
609 
		return EncodeVars(CADDR(formula));
2 r 610 
	}
1839 ihaka 611 
	if (CAR(formula) == plusSymbol) {
612 
	    if (len == 2)
613 
		return EncodeVars(CADR(formula));
614 
	    else
615 
		return PlusTerms(CADR(formula), CADDR(formula));
2 r 616 
	}
1839 ihaka 617 
	if (CAR(formula) == colonSymbol) {
618 
	    return InteractTerms(CADR(formula), CADDR(formula));
619 
	}
620 
	if (CAR(formula) == timesSymbol) {
621 
	    return CrossTerms(CADR(formula), CADDR(formula));
622 
	}
623 
	if (CAR(formula) == inSymbol) {
624 
	    return InTerms(CADR(formula), CADDR(formula));
625 
	}
626 
	if (CAR(formula) == slashSymbol) {
627 
	    return NestTerms(CADR(formula), CADDR(formula));
628 
	}
629 
	if (CAR(formula) == powerSymbol) {
630 
	    return PowerTerms(CADR(formula), CADDR(formula));
631 
	}
632 
	if (CAR(formula) == minusSymbol) {
633 
	    if (len == 2)
634 
		return DeleteTerms(R_NilValue, CADR(formula));
635 
	    return DeleteTerms(CADR(formula), CADDR(formula));
636 
	}
637 
	if (CAR(formula) == parenSymbol) {
638 
	    return EncodeVars(CADR(formula));
639 
	}
640 
	term = AllocTerm();
641 
	SetBit(term, InstallVar(formula), 1);
642 
	return CONS(term, R_NilValue);
643 
    }
24985 ripley 644 
    error("invalid model formula in EncodeVars");
1839 ihaka 645 
    return R_NilValue;/*NOTREACHED*/
2 r 646 
}
647 
 
648 
 
1839 ihaka 649 
/* TermCode decides on the encoding of a model term. */
650 
/* Returns 1 if variable ``whichBit'' in ``thisTerm'' */
651 
/* is to be encoded by contrasts and 2 if it is to be */
652 
/* encoded by dummy variables.  This is decided using */
15501 ripley 653 
/* the heuristic of Chambers and Heiberger described */
1839 ihaka 654 
/* in Statistical Models in S, Page 38. */
2 r 655 
 
656 
static int TermCode(SEXP termlist, SEXP thisterm, int whichbit, SEXP term)
657 
{
1839 ihaka 658 
    SEXP t;
659 
    int allzero, i;
2 r 660 
 
1839 ihaka 661 
    for (i = 0; i < nwords; i++)
662 
	INTEGER(term)[i] = INTEGER(CAR(thisterm))[i];
2 r 663 
 
1839 ihaka 664 
    /* Eliminate factor ``whichbit'' */
2 r 665 
 
1839 ihaka 666 
    SetBit(term, whichbit, 0);
2 r 667 
 
1839 ihaka 668 
    /* Search preceding terms for a match */
669 
    /* Zero is a possibility - it is a special case */
2 r 670 
 
1839 ihaka 671 
    allzero = 1;
672 
    for (i = 0; i < nwords; i++) {
673 
	if (INTEGER(term)[i]) {
674 
	    allzero = 0;
675 
	    break;
676 
	}
677 
    }
678 
    if (allzero)
679 
	return 1;
680 
 
681 
    for (t = termlist; t != thisterm; t = CDR(t)) {
2 r 682 
	allzero = 1;
683 
	for (i = 0; i < nwords; i++) {
1839 ihaka 684 
	    if ((~(INTEGER(CAR(t))[i])) & INTEGER(term)[i])
685 
		allzero = 0;
2 r 686 
	}
687 
	if (allzero)
1839 ihaka 688 
	    return 1;
689 
    }
690 
    return 2;
2 r 691 
}
454 maechler 692 
 
2 r 693 
 
1839 ihaka 694 
/* Internal code for the ``terms'' function */
695 
/* The value is a formula with an assortment */
696 
/* of useful attributes. */
2 r 697 
 
454 maechler 698 
/* .Internal(terms.formula(x, new.specials, abb, data, keep.order)) */
1839 ihaka 699 
 
22903 ripley 700 
static SEXP ExpandDots(SEXP object, SEXP value);
701 
 
2 r 702 
SEXP do_termsform(SEXP call, SEXP op, SEXP args, SEXP rho)
703 
{
1839 ihaka 704 
    SEXP a, ans, v, pattern, formula, varnames, term, termlabs;
22903 ripley 705 
    SEXP specials, t, data, rhs;
1839 ihaka 706 
    int i, j, k, l, n, keepOrder;
2 r 707 
 
1839 ihaka 708 
    checkArity(op, args);
2 r 709 
 
1839 ihaka 710 
    /* Always fetch these values rather than trying */
711 
    /* to remember them between calls.  The overhead */
712 
    /* is minimal and we don't have to worry about */
713 
    /* intervening dump/restore problems. */
2 r 714 
 
1839 ihaka 715 
    tildeSymbol = install("~");
716 
    plusSymbol  = install("+");
717 
    minusSymbol = install("-");
718 
    timesSymbol = install("*");
719 
    slashSymbol = install("/");
720 
    colonSymbol = install(":");
721 
    powerSymbol = install("^");
722 
    dotSymbol   = install(".");
723 
    parenSymbol = install("(");
724 
    inSymbol = install("%in%");
12976 pd 725 
    /* identSymbol = install("I"); */
2 r 726 
 
1839 ihaka 727 
    /* Do we have a model formula? */
728 
    /* Check for unary or binary ~ */
2 r 729 
 
1839 ihaka 730 
    if (!isLanguage(CAR(args)) ||
731 
	CAR(CAR(args)) != tildeSymbol ||
732 
	(length(CAR(args)) != 2 && length(CAR(args)) != 3))
5731 ripley 733 
	error("argument is not a valid model");
2 r 734 
 
22903 ripley 735 
    haveDot = FALSE;
736 
 
1839 ihaka 737 
    PROTECT(ans = duplicate(CAR(args)));
2 r 738 
 
22903 ripley 739 
    /* The formula will be returned, modified if haveDot becomes TRUE */
2 r 740 
 
1839 ihaka 741 
    specials = CADR(args);
742 
    a = CDDR(args);
91 ihaka 743 
 
1839 ihaka 744 
    /* We use data to get the value to */
745 
    /* substitute for "." in formulae */
91 ihaka 746 
 
1839 ihaka 747 
    data = CAR(a);
748 
    a = CDR(a);
749 
    if (isNull(data) || isEnvironment(data))
750 
	framenames = R_NilValue;
751 
    else if (isFrame(data))
752 
	framenames = getAttrib(data, R_NamesSymbol);
753 
    else
5731 ripley 754 
	errorcall(call,"data argument is of the wrong type");
22896 ripley 755 
 
1839 ihaka 756 
    if (framenames != R_NilValue)
757 
	if (length(CAR(args))== 3)
758 
	    CheckRHS(CADR(CAR(args)));
2 r 759 
 
1839 ihaka 760 
    /* Preserve term order? */
91 ihaka 761 
 
1839 ihaka 762 
    keepOrder = asLogical(CAR(a));
763 
    if (keepOrder == NA_LOGICAL)
764 
	keepOrder = 0;
2 r 765 
 
10172 luke 766 
    if (specials == R_NilValue) {
22906 ripley 767 
	a = allocList(8);
10172 luke 768 
	SET_ATTRIB(ans, a);
769 
    }
770 
    else {
22906 ripley 771 
	a = allocList(9);
10172 luke 772 
	SET_ATTRIB(ans, a);
773 
    }
2 r 774 
 
1839 ihaka 775 
    /* Step 1: Determine the ``variables'' in the model */
776 
    /* Here we create an expression of the form */
777 
    /* list(...).  You can evaluate it to get */
778 
    /* the model variables or use substitute and then */
779 
    /* pull the result apart to get the variable names. */
2 r 780 
 
1839 ihaka 781 
    intercept = 1;
782 
    parity = 1;
783 
    response = 0;
784 
    PROTECT(varlist = lcons(install("list"), R_NilValue));
785 
    ExtractVars(CAR(args), 1);
786 
    UNPROTECT(1);
10172 luke 787 
    SETCAR(a, varlist);
788 
    SET_TAG(a, install("variables"));
1839 ihaka 789 
    a = CDR(a);
2 r 790 
 
1839 ihaka 791 
    nvar = length(varlist) - 1;
792 
    nwords = (nvar - 1) / WORDSIZE + 1;
22906 ripley 793 
 
1839 ihaka 794 
    /* Step 2: Recode the model terms in binary form */
795 
    /* and at the same time, expand the model formula. */
2 r 796 
 
1839 ihaka 797 
    /* FIXME: this includes specials in the model */
798 
    /* There perhaps needs to be a an extra pass */
799 
    /* through the model to delete any terms which */
800 
    /* contain specials.  Actually, specials should */
801 
    /* only enter additively so this should also be */
802 
    /* checked and abort forced if not. */
90 ihaka 803 
 
22906 ripley 804 
    /* BDR 2002-01-29: S does include specials, so code may rely on this */
805 
 
1839 ihaka 806 
    /* FIXME: this is also the point where nesting */
807 
    /* needs to be taken care of. */
90 ihaka 808 
 
1839 ihaka 809 
    PROTECT(formula = EncodeVars(CAR(args)));
2 r 810 
 
8505 pd 811 
    nvar = length(varlist) - 1; /* need to recompute, in case
812 
                                   EncodeVars stretched it */
813 
 
22907 ripley 814 
    /* Step 2a: Compute variable names */
815 
 
816 
    PROTECT(varnames = allocVector(STRSXP, nvar));
25911 ripley 817 
    for (v = CDR(varlist), i = 0; v != R_NilValue; v = CDR(v))
818 
	SET_STRING_ELT(varnames, i++, STRING_ELT(deparse1line(CAR(v), 0), 0));
22907 ripley 819 
 
22906 ripley 820 
    /* Step 2b: Remove any offset(s) */
821 
 
822 
    for (l = response, k = 0; l < nvar; l++)
823 
	if (!strncmp(CHAR(STRING_ELT(varnames, l)), "offset(", 7)) k++;
824 
    SETCAR(a, v = allocVector(INTSXP, k));
825 
    if (k > 0) {
26891 ripley 826 
	for (l = response, k = 0; l < nvar; l++)
827 
	    if (!strncmp(CHAR(STRING_ELT(varnames, l)), "offset(", 7))
828 
		INTEGER(v)[k++] = l + 1;
22906 ripley 829 
	call = formula; /* call is to be the previous value */
26872 ripley 830 
	for (l = response, k = 0; l < length(formula)+response; l++) {
831 
	    SEXP thisterm;
832 
	    Rboolean have_offset = FALSE;
833 
	    thisterm = (l > response) ? CDR(call): call;
834 
	    for (i = 1; i <= nvar; i++) {
835 
		if (GetBit(CAR(thisterm), i) &&
836 
		    !strncmp(CHAR(STRING_ELT(varnames, i-1)), "offset(", 7)) {
837 
		    have_offset = TRUE;
838 
		    break;
839 
		}
840 
	    }
841 
	    if (have_offset) {
22906 ripley 842 
		if (l == response) call = formula = CDR(formula);
843 
		else SETCDR(call, CDR(CDR(call)));
844 
	    } else if (l > response) call = CDR(call);
26872 ripley 845 
	}
22906 ripley 846 
	SET_TAG(a, install("offset"));
847 
	a = CDR(a);
848 
    }
849 
    nterm = length(formula);
850 
 
22896 ripley 851 
    /* Step 3: Reorder the model terms by BitCount, otherwise
852 
       preserving their order. */
2 r 853 
 
22896 ripley 854 
    if (!keepOrder) {
22903 ripley 855 
	SEXP sCounts;
856 
	int *counts, bitmax = 0;
857 
 
22896 ripley 858 
	PROTECT(pattern = allocVector(VECSXP, nterm));
22903 ripley 859 
	PROTECT(sCounts = allocVector(INTSXP, nterm));
860 
	counts = INTEGER(sCounts);
22896 ripley 861 
	for (call = formula, n = 0; call != R_NilValue; call = CDR(call)) {
862 
	    SET_VECTOR_ELT(pattern, n, CAR(call));
863 
	    counts[n++] = BitCount(CAR(call));
864 
	}
865 
	for (n = 0; n < nterm; n++)
866 
	    if(counts[n] > bitmax) bitmax = counts[n];
867 
	call = formula;
868 
	for (i = 0; i <= bitmax; i++) /* can order 0 occur? */
869 
	    for (n = 0; n < nterm; n++)
870 
		if (counts[n] == i) {
871 
		    SETCAR(call, VECTOR_ELT(pattern, n));
872 
		    SETLEVELS(CAR(call), i);
873 
		    call = CDR(call);
874 
		}
22903 ripley 875 
	UNPROTECT(2);
1839 ihaka 876 
    }
22896 ripley 877 
 
1839 ihaka 878 
 
879 
    /* Step 4: Compute the factor pattern for the model. */
880 
    /* 0 - the variable does not appear in this term. */
881 
    /* 1 - code the variable by contrasts in this term. */
882 
    /* 2 - code the variable by indicators in this term. */
883 
 
884 
    if (nterm > 0) {
10172 luke 885 
	SETCAR(a, pattern = allocMatrix(INTSXP, nvar, nterm));
886 
	SET_TAG(a, install("factors"));
1839 ihaka 887 
	a = CDR(a);
888 
	for (i = 0; i < nterm * nvar; i++)
889 
	    INTEGER(pattern)[i] = 0;
890 
	PROTECT(term = AllocTerm());
2 r 891 
	n = 0;
892 
	for (call = formula; call != R_NilValue; call = CDR(call)) {
1839 ihaka 893 
	    for (i = 1; i <= nvar; i++) {
894 
		if (GetBit(CAR(call), i))
895 
		    INTEGER(pattern)[i-1+n*nvar] =
896 
			TermCode(formula, call, i, term);
897 
	    }
898 
	    n++;
2 r 899 
	}
900 
	UNPROTECT(1);
1839 ihaka 901 
    }
902 
    else {
10172 luke 903 
	SETCAR(a, pattern = allocVector(INTSXP,0));
904 
	SET_TAG(a, install("factors"));
1839 ihaka 905 
	a = CDR(a);
906 
    }
2 r 907 
 
22906 ripley 908 
    /* Step 5: Compute term labels */
2 r 909 
 
1839 ihaka 910 
    PROTECT(termlabs = allocVector(STRSXP, nterm));
911 
    n = 0;
912 
    for (call = formula; call != R_NilValue; call = CDR(call)) {
913 
	l = 0;
914 
	for (i = 1; i <= nvar; i++) {
915 
	    if (GetBit(CAR(call), i)) {
916 
		if (l > 0)
917 
		    l += 1;
10172 luke 918 
		l += strlen(CHAR(STRING_ELT(varnames, i - 1)));
1839 ihaka 919 
	    }
2 r 920 
	}
10172 luke 921 
	SET_STRING_ELT(termlabs, n, allocString(l));
922 
	CHAR(STRING_ELT(termlabs, n))[0] = '\0';
1839 ihaka 923 
	l = 0;
924 
	for (i = 1; i <= nvar; i++) {
925 
	    if (GetBit(CAR(call), i)) {
926 
		if (l > 0)
10172 luke 927 
		    strcat(CHAR(STRING_ELT(termlabs, n)), ":");
22906 ripley 928 
		strcat(CHAR(STRING_ELT(termlabs, n)),
929 
		       CHAR(STRING_ELT(varnames, i - 1)));
1839 ihaka 930 
		l++;
931 
	    }
2 r 932 
	}
1839 ihaka 933 
	n++;
934 
    }
935 
    PROTECT(v = allocVector(VECSXP, 2));
10172 luke 936 
    SET_VECTOR_ELT(v, 0, varnames);
937 
    SET_VECTOR_ELT(v, 1, termlabs);
1839 ihaka 938 
    if (nterm > 0)
939 
	setAttrib(pattern, R_DimNamesSymbol, v);
2 r 940 
 
10172 luke 941 
    SETCAR(a, termlabs);
942 
    SET_TAG(a, install("term.labels"));
1839 ihaka 943 
    a = CDR(a);
2 r 944 
 
1839 ihaka 945 
    /* If there are specials stick them in here */
946 
 
947 
    if (specials != R_NilValue) {
948 
	i = length(specials);
949 
	PROTECT(v = allocList(i));
950 
	for (j = 0, t = v; j < i; j++, t = CDR(t)) {
10172 luke 951 
	    SET_TAG(t, install(CHAR(STRING_ELT(specials, j))));
952 
	    n = strlen(CHAR(STRING_ELT(specials, j)));
953 
	    SETCAR(t, allocVector(INTSXP, 0));
1839 ihaka 954 
	    k = 0;
955 
	    for (l = 0; l < nvar; l++) {
10172 luke 956 
		if (!strncmp(CHAR(STRING_ELT(varnames, l)),
957 
			    CHAR(STRING_ELT(specials, j)), n))
958 
		    if (CHAR(STRING_ELT(varnames, l))[n] == '(')
1839 ihaka 959 
			k++;
960 
	    }
961 
	    if (k > 0) {
10172 luke 962 
		SETCAR(t, allocVector(INTSXP, k));
1839 ihaka 963 
		k = 0;
964 
		for (l = 0; l < nvar; l++) {
10172 luke 965 
		    if (!strncmp(CHAR(STRING_ELT(varnames, l)),
966 
				CHAR(STRING_ELT(specials, j)), n))
967 
			if (CHAR(STRING_ELT(varnames, l))[n] == '('){
1839 ihaka 968 
			    INTEGER(CAR(t))[k++] = l+1;
2 r 969 
			}
970 
		}
1839 ihaka 971 
	    }
10172 luke 972 
	    else SETCAR(t, R_NilValue);
2 r 973 
	}
10172 luke 974 
	SETCAR(a, v);
975 
	SET_TAG(a, install("specials"));
2 r 976 
	a = CDR(a);
1839 ihaka 977 
	UNPROTECT(1);
978 
    }
2 r 979 
 
1839 ihaka 980 
    UNPROTECT(3);	/* keep termlabs until here */
2 r 981 
 
22903 ripley 982 
    /* Step 6: Fix up the formula by substituting for dot, which should be
983 
       the framenames joined by + */
22906 ripley 984 
 
22903 ripley 985 
    if (haveDot && LENGTH(framenames)) {
27063 ripley 986 
	PROTECT_INDEX ind;
987 
	PROTECT_WITH_INDEX(rhs = install(CHAR(STRING_ELT(framenames, 0))),
988 
			   &ind);
22903 ripley 989 
	for (i = 1; i < LENGTH(framenames); i++) {
27063 ripley 990 
	    REPROTECT(rhs = lang3(plusSymbol, rhs,
991 
				  install(CHAR(STRING_ELT(framenames, i)))),
992 
		      ind);
22903 ripley 993 
	}
994 
	if (!isNull(CADDR(ans)))
995 
	    SETCADDR(ans, ExpandDots(CADDR(ans), rhs));
996 
	else
997 
	    SETCADR(ans, ExpandDots(CADR(ans), rhs));
998 
	UNPROTECT(1);
999 
    }
1000 
 
10172 luke 1001 
    SETCAR(a, allocVector(INTSXP, nterm));
1839 ihaka 1002 
    n = 0;
1003 
    for (call = formula; call != R_NilValue; call = CDR(call))
1004 
	INTEGER(CAR(a))[n++] = LEVELS(CAR(call));
10172 luke 1005 
    SET_TAG(a, install("order"));
1839 ihaka 1006 
    a = CDR(a);
2 r 1007 
 
10172 luke 1008 
    SETCAR(a, allocVector(INTSXP, 1));
1839 ihaka 1009 
    INTEGER(CAR(a))[0] = (intercept != 0);
10172 luke 1010 
    SET_TAG(a, install("intercept"));
1839 ihaka 1011 
    a = CDR(a);
2 r 1012 
 
10172 luke 1013 
    SETCAR(a, allocVector(INTSXP, 1));
1839 ihaka 1014 
    INTEGER(CAR(a))[0] = (response != 0);
10172 luke 1015 
    SET_TAG(a, install("response"));
1839 ihaka 1016 
    a = CDR(a);
2 r 1017 
 
10172 luke 1018 
    SETCAR(a, mkString("terms"));
1019 
    SET_TAG(a, install("class"));
22906 ripley 1020 
    SETCDR(a, R_NilValue);  /* truncate if necessary */
10172 luke 1021 
    SET_OBJECT(ans, 1);
2 r 1022 
 
1839 ihaka 1023 
    UNPROTECT(2);
1024 
    return ans;
2 r 1025 
}
1026 
 
1839 ihaka 1027 
/* Update a model formula by the replacement of "." templates. */
2 r 1028 
 
1029 
static SEXP ExpandDots(SEXP object, SEXP value)
1030 
{
1839 ihaka 1031 
    SEXP op;
2 r 1032 
 
1839 ihaka 1033 
    if (TYPEOF(object) == SYMSXP) {
1034 
	if (object == dotSymbol)
1035 
	    object = duplicate(value);
1036 
	return object;
1037 
    }
2 r 1038 
 
1839 ihaka 1039 
    if (TYPEOF(object) == LANGSXP) {
1040 
	if (TYPEOF(value) == LANGSXP) op = CAR(value);
1041 
	else op = NULL;
1042 
	PROTECT(object);
1043 
	if (CAR(object) == plusSymbol) {
1044 
	    if (length(object) == 2) {
10172 luke 1045 
		SETCADR(object, ExpandDots(CADR(object), value));
1839 ihaka 1046 
	    }
1047 
	    else if (length(object) == 3) {
10172 luke 1048 
		SETCADR(object, ExpandDots(CADR(object), value));
1049 
		SETCADDR(object, ExpandDots(CADDR(object), value));
1839 ihaka 1050 
	    }
1051 
	    else goto badformula;
2 r 1052 
	}
1839 ihaka 1053 
	else if (CAR(object) == minusSymbol) {
1054 
	    if (length(object) == 2) {
1055 
		if (CADR(object) == dotSymbol &&
1056 
		   (op == plusSymbol || op == minusSymbol))
10172 luke 1057 
		    SETCADR(object, lang2(parenSymbol,
1058 
					  ExpandDots(CADR(object), value)));
1839 ihaka 1059 
		else
10172 luke 1060 
		    SETCADR(object, ExpandDots(CADR(object), value));
1839 ihaka 1061 
	    }
1062 
	    else if (length(object) == 3) {
1063 
		if (CADR(object) == dotSymbol &&
1064 
		   (op == plusSymbol || op == minusSymbol))
10172 luke 1065 
		    SETCADR(object, lang2(parenSymbol,
1066 
					  ExpandDots(CADR(object), value)));
1839 ihaka 1067 
		else
10172 luke 1068 
		    SETCADR(object, ExpandDots(CADR(object), value));
1839 ihaka 1069 
		if (CADDR(object) == dotSymbol &&
1070 
		   (op == plusSymbol || op == minusSymbol))
10172 luke 1071 
		    SETCADDR(object, lang2(parenSymbol,
1072 
					   ExpandDots(CADDR(object), value)));
1839 ihaka 1073 
		else
10172 luke 1074 
		    SETCADDR(object, ExpandDots(CADDR(object), value));
1839 ihaka 1075 
	    }
1076 
	    else goto badformula;
1077 
	}
1078 
	else if (CAR(object) == timesSymbol || CAR(object) == slashSymbol) {
1079 
	    if (length(object) != 3)
1080 
		goto badformula;
1081 
	    if (CADR(object) == dotSymbol &&
1082 
	       (op == plusSymbol || op == minusSymbol))
10172 luke 1083 
		SETCADR(object, lang2(parenSymbol,
1084 
				      ExpandDots(CADR(object), value)));
1839 ihaka 1085 
	    else
10172 luke 1086 
		SETCADR(object, ExpandDots(CADR(object), value));
1839 ihaka 1087 
	    if (CADDR(object) == dotSymbol &&
1088 
	       (op == plusSymbol || op == minusSymbol))
10172 luke 1089 
		SETCADDR(object, lang2(parenSymbol,
1090 
				       ExpandDots(CADDR(object), value)));
1839 ihaka 1091 
	    else
10172 luke 1092 
		SETCADDR(object, ExpandDots(CADDR(object), value));
1839 ihaka 1093 
	}
1094 
	else if (CAR(object) == colonSymbol) {
1095 
	    if (length(object) != 3)
1096 
		goto badformula;
1097 
	    if (CADR(object) == dotSymbol &&
1098 
	       (op == plusSymbol || op == minusSymbol ||
1099 
		op == timesSymbol || op == slashSymbol))
10172 luke 1100 
		SETCADR(object, lang2(parenSymbol,
1101 
				      ExpandDots(CADR(object), value)));
1839 ihaka 1102 
	    else
10172 luke 1103 
		SETCADR(object, ExpandDots(CADR(object), value));
1839 ihaka 1104 
	    if (CADDR(object) == dotSymbol &&
1105 
	       (op == plusSymbol || op == minusSymbol))
10172 luke 1106 
		SETCADDR(object, lang2(parenSymbol,
1107 
				       ExpandDots(CADDR(object), value)));
1839 ihaka 1108 
	    else
10172 luke 1109 
		SETCADDR(object, ExpandDots(CADDR(object), value));
1839 ihaka 1110 
	}
1111 
	else if (CAR(object) == powerSymbol) {
1112 
	    if (length(object) != 3)
1113 
		goto badformula;
1114 
	    if (CADR(object) == dotSymbol &&
1115 
	       (op == plusSymbol || op == minusSymbol ||
1116 
		op == timesSymbol || op == slashSymbol ||
1117 
		op == colonSymbol))
10172 luke 1118 
		SETCADR(object, lang2(parenSymbol,
1119 
				      ExpandDots(CADR(object), value)));
1839 ihaka 1120 
	    else
10172 luke 1121 
		SETCADR(object, ExpandDots(CADR(object), value));
1839 ihaka 1122 
	    if (CADDR(object) == dotSymbol &&
1123 
	       (op == plusSymbol || op == minusSymbol))
10172 luke 1124 
		SETCADDR(object, lang2(parenSymbol,
1125 
				       ExpandDots(CADDR(object), value)));
1839 ihaka 1126 
	    else
10172 luke 1127 
		SETCADDR(object, ExpandDots(CADDR(object), value));
1839 ihaka 1128 
	}
1129 
	else {
1130 
	    op = object;
1131 
	    while(op != R_NilValue) {
10172 luke 1132 
		SETCAR(op, ExpandDots(CAR(op), value));
1839 ihaka 1133 
		op = CDR(op);
1134 
	    }
1135 
	}
1136 
	UNPROTECT(1);
1137 
	return object;
1138 
    }
1139 
    else return object;
2 r 1140 
 
1839 ihaka 1141 
 badformula:
5731 ripley 1142 
    error("invalid formula in update");
1839 ihaka 1143 
    return R_NilValue; /*NOTREACHED*/
2 r 1144 
}
1145 
 
1146 
SEXP do_updateform(SEXP call, SEXP op, SEXP args, SEXP rho)
1147 
{
1839 ihaka 1148 
    SEXP new, old, lhs, rhs;
2 r 1149 
 
1839 ihaka 1150 
    checkArity(op, args);
2 r 1151 
 
1839 ihaka 1152 
    /* Always fetch these values rather than trying */
1153 
    /* to remember them between calls.  The overhead */
1154 
    /* is minimal and we don't have to worry about */
1155 
    /* intervening dump/restore problems. */
2 r 1156 
 
1839 ihaka 1157 
    tildeSymbol = install("~");
1158 
    plusSymbol  = install("+");
1159 
    minusSymbol = install("-");
1160 
    timesSymbol = install("*");
1161 
    slashSymbol = install("/");
1162 
    colonSymbol = install(":");
1163 
    powerSymbol = install("^");
1164 
    dotSymbol   = install(".");
1165 
    parenSymbol = install("(");
1166 
    inSymbol = install("%in%");
12976 pd 1167 
    /* identSymbol = install("I"); */
2 r 1168 
 
1839 ihaka 1169 
    /* We must duplicate here because the */
1170 
    /* formulae may be part of the parse tree */
1171 
    /* and we don't want to modify it. */
2 r 1172 
 
1839 ihaka 1173 
    old = CAR(args);
10172 luke 1174 
    new = SETCADR(args, duplicate(CADR(args)));
2 r 1175 
 
1839 ihaka 1176 
    /* Check of new and old formulae. */
1177 
    if (TYPEOF(old) != LANGSXP ||
1178 
       (TYPEOF(new) != LANGSXP && CAR(old) != tildeSymbol) ||
1179 
       CAR(new) != tildeSymbol)
5731 ripley 1180 
	errorcall(call, "formula expected");
2 r 1181 
 
3279 pd 1182 
    if (length(old) == 3) {
1183 
	lhs = CADR(old);
1184 
	rhs = CADDR(old);
1185 
	/* We now check that new formula has a valid lhs.
1186 
	   If it doesn't, we add one and set it to the rhs of the old
1187 
	   formula. */
1188 
	if (length(new) == 2)
10172 luke 1189 
	    SETCDR(new, CONS(lhs, CDR(new)));
3279 pd 1190 
	/* Now we check the left and right sides of the new formula
1191 
	   and substitute the correct value for any "." templates.
1192 
	   We must parenthesize the rhs or we might upset arity and
1193 
	   precedence. */
1194 
	PROTECT(rhs);
10172 luke 1195 
	SETCADR(new, ExpandDots(CADR(new), lhs));
1196 
	SETCADDR(new, ExpandDots(CADDR(new), rhs));
3279 pd 1197 
	UNPROTECT(1);
1198 
    }
1199 
    else {
1200 
	/* The old formula had no lhs, so we only expand the rhs of the
1201 
	   new formula. */
1202 
	rhs = CADR(old);
1203 
	if (length(new) == 3)
10172 luke 1204 
	    SETCADDR(new, ExpandDots(CADDR(new), rhs));
3279 pd 1205 
	else
10172 luke 1206 
	    SETCADR(new, ExpandDots(CADR(new), rhs));
3279 pd 1207 
    }
2 r 1208 
 
1839 ihaka 1209 
    /* It might be overkill to zero the */
1210 
    /* the attribute list of the returned */
1211 
    /* value, but it can't hurt. */
2 r 1212 
 
10172 luke 1213 
    SET_ATTRIB(new, R_NilValue);
1839 ihaka 1214 
    return new;
2 r 1215 
}
454 maechler 1216 
 
2 r 1217 
 
1218 
/*
1219 
 *  model.frame
1220 
 *
1221 
 *  The argument "terms" contains the terms object generated from the
1222 
 *  model formula.  We first evaluate the "variables" attribute of
1223 
 *  "terms" in the "data" environment.  This gives us a list of basic
1224 
 *  variables to be in the model frame.  We do some basic sanity
1225 
 *  checks on these to ensure that resulting object make sense.
1226 
 *
1227 
 *  The argument "dots" gives additional things like "weights", "offsets"
1228 
 *  and "subset" which will also go into the model frame so that they can
1229 
 *  be treated in parallel.
1230 
 *
1231 
 *  Next we subset the data frame according to "subset" and finally apply
1232 
 *  "na.action" to get the final data frame.
1233 
 *
1234 
 *  Note that the "terms" argument is glued to the model frame as an
1235 
 *  attribute.  Code downstream appears to need this.
454 maechler 1236 
 *
2 r 1237 
 *  Q: Is this really needed, or can we get by with less info?
1238 
 */
1239 
 
20590 ripley 1240 
/* time to move more functionality back into compiled
9060 tlumley 1241 
   code (cycle of reincarnation) */
2 r 1242 
 
1858 ihaka 1243 
/* .Internal(model.frame(terms, rownames, variables, varnames, */
1244 
/*           dots, dotnames, subset, na.action)) */
454 maechler 1245 
 
2 r 1246 
SEXP do_modelframe(SEXP call, SEXP op, SEXP args, SEXP rho)
1247 
{
1858 ihaka 1248 
    SEXP terms, data, names, variables, varnames, dots, dotnames, na_action;
4790 hornik 1249 
    SEXP ans, row_names, subset, tmp;
1839 ihaka 1250 
    char buf[256];
9060 tlumley 1251 
    int i, j, nr, nc;
1252 
    int nvars, ndots, nactualdots;
2 r 1253 
 
1839 ihaka 1254 
    checkArity(op, args);
1255 
    terms = CAR(args); args = CDR(args);
1858 ihaka 1256 
    row_names = CAR(args); args = CDR(args);
1257 
    variables = CAR(args); args = CDR(args);
1258 
    varnames = CAR(args); args = CDR(args);
1839 ihaka 1259 
    dots = CAR(args); args = CDR(args);
1858 ihaka 1260 
    dotnames = CAR(args); args = CDR(args);
1839 ihaka 1261 
    subset = CAR(args); args = CDR(args);
1262 
    na_action = CAR(args); args = CDR(args);
2 r 1263 
 
1858 ihaka 1264 
    /* Argument Sanity Checks */
2 r 1265 
 
1858 ihaka 1266 
    if (!isNewList(variables))
5731 ripley 1267 
	errorcall(call, "invalid variables");
1858 ihaka 1268 
    if (!isString(varnames))
5731 ripley 1269 
	errorcall(call, "invalid variable names");
1858 ihaka 1270 
    if ((nvars = length(variables)) != length(varnames))
5731 ripley 1271 
	errorcall(call, "number of variables != number of variable names");
20590 ripley 1272 
 
1858 ihaka 1273 
    if (!isNewList(dots))
5731 ripley 1274 
	errorcall(call, "invalid extra variables");
1858 ihaka 1275 
    if ((ndots = length(dots)) != length(dotnames))
5731 ripley 1276 
	errorcall(call, "number of variables != number of variable names");
10538 rgentlem 1277 
    if ( ndots && !isString(dotnames))
1278 
	errorcall(call, "invalid extra variable names");
20590 ripley 1279 
 
15501 ripley 1280 
    /*  check for NULL extra arguments -- moved from interpreted code */
2 r 1281 
 
15501 ripley 1282 
    nactualdots = 0;
1283 
    for (i = 0; i < ndots; i++){
20590 ripley 1284 
	if (VECTOR_ELT(dots, i) != R_NilValue)
9060 tlumley 1285 
	    nactualdots++;
1286 
    }
1287 
 
1839 ihaka 1288 
    /* Assemble the base data frame. */
20590 ripley 1289 
 
9060 tlumley 1290 
    PROTECT(data = allocVector(VECSXP, nvars + nactualdots));
1291 
    PROTECT(names = allocVector(STRSXP, nvars + nactualdots));
2 r 1292 
 
1858 ihaka 1293 
    for (i = 0; i < nvars; i++) {
10172 luke 1294 
	SET_VECTOR_ELT(data, i, VECTOR_ELT(variables, i));
1295 
	SET_STRING_ELT(names, i, STRING_ELT(varnames, i));
1839 ihaka 1296 
    }
15501 ripley 1297 
    for (i = 0,j = 0; i < ndots; i++) {
1298 
	if (VECTOR_ELT(dots, i) == R_NilValue)
9060 tlumley 1299 
	    continue;
25630 ripley 1300 
	if(strlen(CHAR(STRING_ELT(dotnames, i))) + 3 > 256)
1301 
	    error("overlong names in %s", CHAR(STRING_ELT(dotnames, i)));
1302 
	sprintf(buf, "(%s)", CHAR(STRING_ELT(dotnames, i)));
10172 luke 1303 
	SET_VECTOR_ELT(data, nvars + j, VECTOR_ELT(dots, i));
1304 
	SET_STRING_ELT(names, nvars + j,  mkChar(buf));
9060 tlumley 1305 
	j++;
1839 ihaka 1306 
    }
1858 ihaka 1307 
    setAttrib(data, R_NamesSymbol, names);
1839 ihaka 1308 
    UNPROTECT(2);
1864 ihaka 1309 
 
1839 ihaka 1310 
    /* Sanity checks to ensure that the the answer can become */
1311 
    /* a data frame.  Be deeply suspicious here! */
2 r 1312 
 
1839 ihaka 1313 
    nc = length(data);
2460 hornik 1314 
    nr = 0;			/* -Wall */
3279 pd 1315 
    if (nc > 0) {
10172 luke 1316 
	nr = nrows(VECTOR_ELT(data, 0));
1839 ihaka 1317 
	for (i = 0; i < nc; i++) {
10172 luke 1318 
	    ans = VECTOR_ELT(data, i);
1839 ihaka 1319 
	    if (TYPEOF(ans) < LGLSXP ||
1320 
		TYPEOF(ans) > REALSXP)
5731 ripley 1321 
		errorcall(call, "invalid variable type");
1839 ihaka 1322 
	    if (nrows(ans) != nr)
5731 ripley 1323 
		errorcall(call, "variable lengths differ");
2 r 1324 
	}
24239 ripley 1325 
    } else nr = length(row_names);
1326 
 
1839 ihaka 1327 
    PROTECT(data);
1328 
    PROTECT(subset);
2 r 1329 
 
1839 ihaka 1330 
    /* Turn the data "list" into a "data.frame" */
1331 
    /* so that subsetting methods will work. */
1332 
    /* To do this we must attach "class"  and */
1333 
    /* "row.names" attributes */
2 r 1334 
 
1839 ihaka 1335 
    PROTECT(tmp = mkString("data.frame"));
1336 
    setAttrib(data, R_ClassSymbol, tmp);
1337 
    UNPROTECT(1);
1338 
    if (length(row_names) == nr) {
1339 
	setAttrib(data, R_RowNamesSymbol, row_names);
1340 
    }
1341 
    else {
1342 
	PROTECT(row_names = allocVector(STRSXP, nr));
1343 
	for (i=0; i<nr; i++) {
1344 
	    sprintf(buf, "%d", i+1);
10172 luke 1345 
	    SET_STRING_ELT(row_names, i, mkChar(buf));
2 r 1346 
	}
1839 ihaka 1347 
	setAttrib(data, R_RowNamesSymbol, row_names);
1348 
	UNPROTECT(1);
1349 
    }
2 r 1350 
 
1839 ihaka 1351 
    /* Do the subsetting, if required. */
3475 pd 1352 
    /* Need to save and restore 'most' attributes */
2 r 1353 
 
1839 ihaka 1354 
    if (subset != R_NilValue) {
20590 ripley 1355 
	PROTECT(tmp=install("[.data.frame"));
18405 hornik 1356 
	PROTECT(tmp=LCONS(tmp,list4(data,subset,R_MissingArg,mkFalse())));
4669 pd 1357 
	data = eval(tmp, rho);
1358 
	UNPROTECT(2);
1839 ihaka 1359 
    }
1858 ihaka 1360 
    UNPROTECT(2);
1839 ihaka 1361 
    PROTECT(data);
2 r 1362 
 
1839 ihaka 1363 
    /* finally, we run na.action on the data frame */
8202 pd 1364 
    /* usually, this will be na.omit */
2 r 1365 
 
1839 ihaka 1366 
    if (na_action != R_NilValue) {
14333 ripley 1367 
	/* some na.actions need this to distinguish responses from
1368 
	   explanatory variables */
1369 
	setAttrib(data, install("terms"), terms);
1839 ihaka 1370 
	if (isString(na_action) && length(na_action) > 0)
10172 luke 1371 
	    na_action = install(CHAR(STRING_ELT(na_action, 0)));
1839 ihaka 1372 
	PROTECT(na_action);
1373 
	PROTECT(tmp = lang2(na_action, data));
3475 pd 1374 
	PROTECT(ans = eval(tmp, rho));
3076 pd 1375 
	if (!isNewList(ans) || length(ans) != length(data))
5731 ripley 1376 
	    errorcall(call, "invalid result from na.action");
20590 ripley 1377 
	/* need to transfer _all but dim_ attributes, possibly lost
1378 
	   by subsetting in na.action.  */
3076 pd 1379 
	for ( i = length(ans) ; i-- ; )
10172 luke 1380 
	  	copyMostAttrib(VECTOR_ELT(data, i),VECTOR_ELT(ans, i));
8202 pd 1381 
 
3475 pd 1382 
	UNPROTECT(3);
1839 ihaka 1383 
    }
1384 
    else ans = data;
1385 
    UNPROTECT(1);
1386 
    PROTECT(ans);
1387 
 
26987 ripley 1388 
    /* Finally, tack on a terms attribute
1389 
       Now done at R level.
1390 
       setAttrib(ans, install("terms"), terms); */
1839 ihaka 1391 
    UNPROTECT(1);
1392 
    return ans;
2 r 1393 
}
1394 
 
1395 
	/* Internal code for the ~ operator */
1396 
	/* Just returns the unevaluated call */
1397 
	/* No longer needed??? */
1398 
 
1399 
SEXP do_tilde(SEXP call, SEXP op, SEXP args, SEXP rho)
1400 
{
10538 rgentlem 1401 
    if (isObject(call))
1402 
        return duplicate(call);
1403 
    else {
1404 
        SEXP class;
1405 
        PROTECT(call = duplicate(call));
1406 
        PROTECT(class = allocVector(STRSXP, 1));
1407 
        SET_STRING_ELT(class, 0, mkChar("formula"));
1408 
        setAttrib(call, R_ClassSymbol, class);
1409 
        setAttrib(call, R_DotEnvSymbol, rho);
1410 
        UNPROTECT(2);
1411 
        return call;
1412 
    }
2 r 1413 
}
1414 
 
1415 
 
1416 
	/* The code below is related to model expansion */
1417 
	/* and is ultimately called by do_modelmatrix. */
1418 
 
1839 ihaka 1419 
static void firstfactor(double *x, int nrx, int ncx,
1420 
			double *c, int nrc, int ncc, int *v)
2 r 1421 
{
1839 ihaka 1422 
    double *cj, *xj;
1423 
    int i, j;
2 r 1424 
 
1839 ihaka 1425 
    for (j = 0; j < ncc; j++) {
1426 
	xj = &x[j*nrx];
1427 
	cj = &c[j*nrc];
1428 
	for (i = 0; i < nrx; i++)
10834 ripley 1429 
	    if(v[i] == NA_INTEGER) xj[i] = NA_REAL;
1430 
	    else xj[i] = cj[v[i]-1];
1839 ihaka 1431 
    }
2 r 1432 
}
1433 
 
1839 ihaka 1434 
static void addfactor(double *x, int nrx, int ncx,
1435 
		      double *c, int nrc, int ncc, int *v)
2 r 1436 
{
1839 ihaka 1437 
    int i, j, k;
1438 
    double *ck, *xj, *yj;
2 r 1439 
 
1839 ihaka 1440 
    for (k = ncc - 1; k >= 0; k--) {
1441 
	for (j = 0; j < ncx; j++) {
1442 
	    xj = &x[j*nrx];
1443 
	    yj = &x[(k*ncx+j)*nrx];
1444 
	    ck = &c[k*nrc];
1445 
	    for (i = 0; i < nrx; i++)
10834 ripley 1446 
	    if(v[i] == NA_INTEGER) yj[i] = NA_REAL;
1447 
	    else yj[i] = ck[v[i]-1] * xj[i];
2 r 1448 
	}
1839 ihaka 1449 
    }
2 r 1450 
}
1451 
 
1452 
static void firstvar(double *x, int nrx, int ncx, double *c, int nrc, int ncc)
1453 
{
1839 ihaka 1454 
    double *cj, *xj;
1455 
    int i, j;
2 r 1456 
 
1839 ihaka 1457 
    for (j = 0; j < ncc; j++) {
1458 
	xj = &x[j*nrx];
1459 
	cj = &c[j*nrc];
1460 
	for (i = 0; i < nrx; i++)
1461 
	    xj[i] = cj[i];
1462 
    }
2 r 1463 
}
1464 
 
1465 
static void addvar(double *x, int nrx, int ncx, double *c, int nrc, int ncc)
1466 
{
1839 ihaka 1467 
    int i, j, k;
1468 
    double *ck, *xj, *yj;
2 r 1469 
 
1839 ihaka 1470 
    for (k = ncc - 1; k >= 0; k--) {
1471 
	for (j = 0; j < ncx; j++) {
1472 
	    xj = &x[j*nrx];
1473 
	    yj = &x[(k*ncx+j)*nrx];
1474 
	    ck = &c[k*nrc];
1475 
	    for (i = 0; i < nrx; i++)
1476 
		yj[i] = ck[i] * xj[i];
2 r 1477 
	}
1839 ihaka 1478 
    }
2 r 1479 
}
1480 
 
4662 pd 1481 
#define BUFSIZE 4096
2 r 1482 
 
1483 
static char *AppendString(char *buf, char *str)
1484 
{
1839 ihaka 1485 
    while (*str)
1486 
	*buf++ = *str++;
1487 
    *buf = '\0';
1488 
    return buf;
2 r 1489 
}
1490 
 
1491 
static char *AppendInteger(char *buf, int i)
1492 
{
1839 ihaka 1493 
    sprintf(buf, "%d", i);
1494 
    while(*buf) buf++;
1495 
    return buf;
2 r 1496 
}
1497 
 
1858 ihaka 1498 
static SEXP ColumnNames(SEXP x)
1499 
{
1500 
    SEXP dn = getAttrib(x, R_DimNamesSymbol);
1501 
    if (dn == R_NilValue)
1502 
	return R_NilValue;
1503 
    else
10172 luke 1504 
	return VECTOR_ELT(dn, 1);
1858 ihaka 1505 
}
1506 
 
2 r 1507 
SEXP do_modelmatrix(SEXP call, SEXP op, SEXP args, SEXP rho)
1508 
{
12976 pd 1509 
    SEXP expr, factors, terms, vars, vnames, assign;
1858 ihaka 1510 
    SEXP xnames, tnames, rnames;
1839 ihaka 1511 
    SEXP count, contrast, contr1, contr2, nlevs, ordered, columns, x;
1512 
    SEXP variable, var_i;
13607 maechler 1513 
    int fik, first, i, j, k, kk, ll, n, nc, nterms, nVar;
1514 
    int intrcept, jstart, jnext, risponse, indx, rhs_response;
21144 tlumley 1515 
    char buf[BUFSIZE]="\0", *bufp, *addp;
1516 
 
2 r 1517 
 
1839 ihaka 1518 
    checkArity(op, args);
2 r 1519 
 
1839 ihaka 1520 
    /* Get the "terms" structure and extract */
1521 
    /* the intercept and response attributes. */
2 r 1522 
 
1839 ihaka 1523 
    terms = CAR(args);
2 r 1524 
 
13607 maechler 1525 
    intrcept = asLogical(getAttrib(terms, install("intercept")));
1526 
    if (intrcept == NA_INTEGER)
1527 
	intrcept = 0;
2 r 1528 
 
13607 maechler 1529 
    risponse = asLogical(getAttrib(terms, install("response")));
1530 
    if (risponse == NA_INTEGER)
1531 
	risponse = 0;
2 r 1532 
 
1839 ihaka 1533 
    /* Get the factor pattern matrix.  We duplicate this because */
1534 
    /* we may want to alter it if we are in the no-intercept case. */
13607 maechler 1535 
    /* Note: the values of "nVar" and "nterms" are the REAL number of */
1839 ihaka 1536 
    /* variables in the model data frame and the number of model terms. */
2 r 1537 
 
13607 maechler 1538 
    nVar = nterms = 0;		/* -Wall */
1839 ihaka 1539 
    PROTECT(factors = duplicate(getAttrib(terms, install("factors"))));
1540 
    if (length(factors) == 0) {
25000 ripley 1541 
	/* if (intrcept == 0)
1542 
	   errorcall(call, "illegal model (zero parameters).");*/
13607 maechler 1543 
	nVar = 1;
1839 ihaka 1544 
	nterms = 0;
1545 
    }
1546 
    else if (isInteger(factors) && isMatrix(factors)) {
13607 maechler 1547 
	nVar = nrows(factors);
1839 ihaka 1548 
	nterms = ncols(factors);
1549 
    }
5731 ripley 1550 
    else errorcall(call, "invalid terms argument");
2 r 1551 
 
1839 ihaka 1552 
    /* Get the variable names from the factor matrix */
2 r 1553 
 
1839 ihaka 1554 
    vnames = getAttrib(factors, R_DimNamesSymbol);
2460 hornik 1555 
    if (length(factors) > 0) {
1556 
	if (length(vnames) < 1 ||
13607 maechler 1557 
	    (nVar - intrcept > 0 && !isString(VECTOR_ELT(vnames, 0))))
5731 ripley 1558 
	    errorcall(call, "invalid terms argument");
10172 luke 1559 
	vnames = VECTOR_ELT(vnames, 0);
2460 hornik 1560 
    }
2 r 1561 
 
1839 ihaka 1562 
    /* Get the variables from the model frame.  First perform */
1563 
    /* elementary sanity checks.  Notes:  1) We need at least */
1564 
    /* one variable (lhs or rhs) to compute the number of cases. */
1565 
    /* 2) We don't type-check the response. */
2 r 1566 
 
1839 ihaka 1567 
    vars = CADR(args);
13607 maechler 1568 
    if (!isNewList(vars) || length(vars) < nVar)
5731 ripley 1569 
	errorcall(call, "invalid model frame");
1839 ihaka 1570 
    if (length(vars) == 0)
5731 ripley 1571 
	errorcall(call, "don't know how many cases");
10172 luke 1572 
    n = nrows(VECTOR_ELT(vars, 0));
1846 ihaka 1573 
    rnames = getAttrib(vars, R_RowNamesSymbol);
2 r 1574 
 
1839 ihaka 1575 
    /* This section of the code checks the types of the variables */
1576 
    /* in the model frame.  Note that it should really only check */
1577 
    /* the variables if they appear in a term in the model. */
2 r 1578 
 
13607 maechler 1579 
    PROTECT(variable = allocVector(VECSXP, nVar));
1580 
    PROTECT(nlevs = allocVector(INTSXP, nVar));
1581 
    PROTECT(ordered = allocVector(LGLSXP, nVar));
1582 
    PROTECT(columns = allocVector(INTSXP, nVar));
91 ihaka 1583 
 
13607 maechler 1584 
    for (i = 0; i < nVar; i++) {
10172 luke 1585 
	var_i = SET_VECTOR_ELT(variable, i, VECTOR_ELT(vars, i));
1839 ihaka 1586 
	if (nrows(var_i) != n)
5731 ripley 1587 
	    errorcall(call, "variable lengths differ");
13607 maechler 1588 
	/*if (i == risponse - 1) {
1839 ihaka 1589 
	    LOGICAL(ordered)[0] = 0;
1590 
	    INTEGER(nlevs)[0] = 0;
1591 
	    INTEGER(columns)[0] = 0;
2 r 1592 
	}
7162 ripley 1593 
	else */
1594 
	if (isOrdered(var_i)) {
1839 ihaka 1595 
	    LOGICAL(ordered)[i] = 1;
27024 ripley 1596 
	    if((INTEGER(nlevs)[i] = nlevels(var_i)) < 1)
1597 
		errorcall(call, "variable %d has no levels", i+1);
27033 ripley 1598 
	    /* will get updated later when contrasts are set */
1839 ihaka 1599 
	    INTEGER(columns)[i] = ncols(var_i);
1600 
	}
1601 
	else if (isUnordered(var_i)) {
1602 
	    LOGICAL(ordered)[i] = 0;
27024 ripley 1603 
	    if((INTEGER(nlevs)[i] = nlevels(var_i)) < 1)
1604 
		errorcall(call, "variable %d has no levels", i+1);
27033 ripley 1605 
	    /* will get updated later when contrasts are set */
1839 ihaka 1606 
	    INTEGER(columns)[i] = ncols(var_i);
1607 
	}
15501 ripley 1608 
	else if (isLogical(var_i)) {
27033 ripley 1609 
	    /* currently this cannot happen as R code turns
1610 
	       logical into factor when setting contrasts */
15501 ripley 1611 
	    LOGICAL(ordered)[i] = 0;
1612 
	    INTEGER(nlevs)[i] = 2;
1613 
	    INTEGER(columns)[i] = ncols(var_i);
1614 
	}
1839 ihaka 1615 
	else if (isNumeric(var_i)) {
10172 luke 1616 
	    SET_VECTOR_ELT(variable, i, coerceVector(var_i, REALSXP));
1617 
	    var_i = VECTOR_ELT(variable, i);
1839 ihaka 1618 
	    LOGICAL(ordered)[i] = 0;
1619 
	    INTEGER(nlevs)[i] = 0;
1620 
	    INTEGER(columns)[i] = ncols(var_i);
1621 
	}
1895 ihaka 1622 
	else
5731 ripley 1623 
	    errorcall(call, "invalid variable type");
1839 ihaka 1624 
    }
2 r 1625 
 
1839 ihaka 1626 
    /* If there is no intercept we look through the factor pattern */
1627 
    /* matrix and adjust the code for the first factor found so that */
1628 
    /* it will be coded by dummy variables rather than contrasts. */
2 r 1629 
 
13607 maechler 1630 
    if (!intrcept) {
1839 ihaka 1631 
	for (j = 0; j < nterms; j++) {
13607 maechler 1632 
	    for (i = risponse; i < nVar; i++) {
1839 ihaka 1633 
		if (INTEGER(nlevs)[i] > 1
20590 ripley 1634 
		    && INTEGER(factors)[i + j * nVar] > 0) {
13607 maechler 1635 
		    INTEGER(factors)[i + j * nVar] = 2;
1839 ihaka 1636 
		    goto alldone;
2 r 1637 
		}
1839 ihaka 1638 
	    }
2 r 1639 
	}
1839 ihaka 1640 
    }
1641 
 alldone:
1642 
    ;
2 r 1643 
 
1858 ihaka 1644 
    /* Compute the required contrast or dummy variable matrices. */
1645 
    /* We set up a symbolic expression to evaluate these, substituting */
1646 
    /* the required arguments at call time.  The calls have the following */
1647 
    /* form: (contrast.type nlevs contrasts) */
2 r 1648 
 
13607 maechler 1649 
    PROTECT(contr1 = allocVector(VECSXP, nVar));
1650 
    PROTECT(contr2 = allocVector(VECSXP, nVar));
397 ihaka 1651 
 
1839 ihaka 1652 
    PROTECT(expr = allocList(3));
10172 luke 1653 
    SET_TYPEOF(expr, LANGSXP);
1654 
    SETCAR(expr, install("contrasts"));
1655 
    SETCADDR(expr, allocVector(LGLSXP, 1));
91 ihaka 1656 
 
1839 ihaka 1657 
    /* FIXME: We need to allow a third argument to this function */
1658 
    /* which allows us to specify contrasts directly.  That argument */
1659 
    /* would be used here in exactly the same way as the below. */
1660 
    /* I.e. we would search the list of constrast specs before */
1661 
    /* we try the evaluation below. */
397 ihaka 1662 
 
13607 maechler 1663 
    for (i = 0; i < nVar; i++) {
1839 ihaka 1664 
	if (INTEGER(nlevs)[i]) {
1665 
	    k = 0;
1666 
	    for (j = 0; j < nterms; j++) {
13607 maechler 1667 
		if (INTEGER(factors)[i + j * nVar] == 1)
1839 ihaka 1668 
		    k |= 1;
13607 maechler 1669 
		else if (INTEGER(factors)[i + j * nVar] == 2)
1839 ihaka 1670 
		    k |= 2;
1671 
	    }
10172 luke 1672 
	    SETCADR(expr, VECTOR_ELT(variable, i));
1839 ihaka 1673 
	    if (k & 1) {
1674 
		LOGICAL(CADDR(expr))[0] = 1;
10172 luke 1675 
		SET_VECTOR_ELT(contr1, i, eval(expr, rho));
1839 ihaka 1676 
	    }
1677 
	    if (k & 2) {
1678 
		LOGICAL(CADDR(expr))[0] = 0;
10172 luke 1679 
		SET_VECTOR_ELT(contr2, i, eval(expr, rho));
1839 ihaka 1680 
	    }
2 r 1681 
	}
1839 ihaka 1682 
    }
2 r 1683 
 
7162 ripley 1684 
    /* By convention, an rhs term identical to the response generates nothing
1685 
       in the model matrix (but interactions involving the response do). */
1686 
 
1687 
    rhs_response = -1;
13607 maechler 1688 
    if (risponse > 0) /* there is a response specified */
7162 ripley 1689 
	for (j = 0; j < nterms; j++)
13607 maechler 1690 
	    if (INTEGER(factors)[risponse - 1 + j * nVar]) {
1691 
		for (i = 0, k = 0; i < nVar; i++)
1692 
		    k += INTEGER(factors)[i + j * nVar] > 0;
7162 ripley 1693 
		if (k == 1) {
1694 
		    rhs_response = j;
1695 
		    break;
1696 
		}
1697 
	    }
1698 
 
20590 ripley 1699 
 
1858 ihaka 1700 
    /* We now have everything needed to build the design matrix. */
1701 
    /* The first step is to compute the matrix size and to allocate it. */
1702 
    /* Note that "count" holds a count of how many columns there are */
1703 
    /* for each term in the model and "nc" gives the total column count. */
2 r 1704 
 
1839 ihaka 1705 
    PROTECT(count = allocVector(INTSXP, nterms));
13607 maechler 1706 
    if (intrcept)
1858 ihaka 1707 
	nc = 1;
1708 
    else
1709 
	nc = 0;
1839 ihaka 1710 
    for (j = 0; j < nterms; j++) {
8464 tlumley 1711 
	if (j == rhs_response) {
25614 ripley 1712 
	    warning("the response appeared on the rhs and was dropped");
1713 
	    INTEGER(count)[j] = 0;  /* need this initialised */
8464 tlumley 1714 
	    continue;
1715 
	}
1839 ihaka 1716 
	k = 1;
13607 maechler 1717 
	for (i = 0; i < nVar; i++) {
1718 
	    if (INTEGER(factors)[i + j * nVar]) {
1839 ihaka 1719 
		if (INTEGER(nlevs)[i]) {
13607 maechler 1720 
		    switch(INTEGER(factors)[i + j * nVar]) {
1839 ihaka 1721 
		    case 1:
10172 luke 1722 
			k *= ncols(VECTOR_ELT(contr1, i));
1839 ihaka 1723 
			break;
1724 
		    case 2:
10172 luke 1725 
			k *= ncols(VECTOR_ELT(contr2, i));
1839 ihaka 1726 
			break;
1727 
		    }
2 r 1728 
		}
1839 ihaka 1729 
		else k *= INTEGER(columns)[i];
1730 
	    }
2 r 1731 
	}
1839 ihaka 1732 
	INTEGER(count)[j] = k;
1733 
	nc = nc + k;
1734 
    }
2 r 1735 
 
1858 ihaka 1736 
    /* Record which columns of the design matrix are associated */
1737 
    /* with which model terms. */
2 r 1738 
 
1839 ihaka 1739 
    PROTECT(assign = allocVector(INTSXP, nc));
1740 
    k = 0;
13607 maechler 1741 
    if (intrcept) INTEGER(assign)[k++] = 0;
1839 ihaka 1742 
    for (j = 0; j < nterms; j++)
1743 
	for (i = 0; i < INTEGER(count)[j]; i++)
1744 
	    INTEGER(assign)[k++] = j+1;
2 r 1745 
 
454 maechler 1746 
 
1839 ihaka 1747 
    /* Create column labels for the matrix columns. */
2 r 1748 
 
1839 ihaka 1749 
    PROTECT(xnames = allocVector(STRSXP, nc));
2 r 1750 
 
1858 ihaka 1751 
    /* Here we loop over the terms in the model and, within each */
1752 
    /* term, loop over the corresponding columns of the design */
1753 
    /* matrix, assembling the names. */
1754 
 
1755 
    /* FIXME : The body within these two loops should be embedded */
1756 
    /* in its own function. */
1757 
 
1839 ihaka 1758 
    k = 0;
13607 maechler 1759 
    if (intrcept)
10172 luke 1760 
	SET_STRING_ELT(xnames, k++, mkChar("(Intercept)"));
2 r 1761 
 
1839 ihaka 1762 
    for (j = 0; j < nterms; j++) {
7162 ripley 1763 
	if (j == rhs_response) continue;
1839 ihaka 1764 
	for (kk = 0; kk < INTEGER(count)[j]; kk++) {
1765 
	    first = 1;
13257 maechler 1766 
	    indx = kk;
1839 ihaka 1767 
	    bufp = &buf[0];
13607 maechler 1768 
	    for (i = 0; i < nVar; i++) {
1769 
		ll = INTEGER(factors)[i + j * nVar];
2460 hornik 1770 
		if (ll) {
20590 ripley 1771 
		    var_i = VECTOR_ELT(variable, i);
1858 ihaka 1772 
		    if (!first)
11025 ripley 1773 
			bufp = AppendString(bufp, ":");
1839 ihaka 1774 
		    first = 0;
15567 ripley 1775 
		    if (isFactor(var_i) || isLogical(var_i)) {
1839 ihaka 1776 
			if (ll == 1) {
10172 luke 1777 
			    x = ColumnNames(VECTOR_ELT(contr1, i));
1778 
			    ll = ncols(VECTOR_ELT(contr1, i));
2 r 1779 
			}
1839 ihaka 1780 
			else {
10172 luke 1781 
			    x = ColumnNames(VECTOR_ELT(contr2, i));
1782 
			    ll = ncols(VECTOR_ELT(contr2, i));
1839 ihaka 1783 
			}
10172 luke 1784 
			addp = CHAR(STRING_ELT(vnames, i));
7162 ripley 1785 
			if(strlen(buf) + strlen(addp) < BUFSIZE)
1786 
			    bufp = AppendString(bufp, addp);
20590 ripley 1787 
			else
7162 ripley 1788 
			    warningcall(call, "term names will be truncated");
1789 
			if (x == R_NilValue) {
1790 
			    if(strlen(buf) + 10 < BUFSIZE)
13257 maechler 1791 
				bufp = AppendInteger(bufp, indx % ll + 1);
20590 ripley 1792 
			    else
7162 ripley 1793 
				warningcall(call, "term names will be truncated");
1794 
			} else {
13257 maechler 1795 
			    addp = CHAR(STRING_ELT(x, indx % ll));
7162 ripley 1796 
			    if(strlen(buf) + strlen(addp) < BUFSIZE)
1797 
				bufp = AppendString(bufp, addp);
20590 ripley 1798 
			    else
7162 ripley 1799 
				warningcall(call, "term names will be truncated");
1800 
			}
1839 ihaka 1801 
		    }
1802 
		    else {
1858 ihaka 1803 
			x = ColumnNames(var_i);
1839 ihaka 1804 
			ll = ncols(var_i);
10172 luke 1805 
			addp = CHAR(STRING_ELT(vnames, i));
7162 ripley 1806 
			if(strlen(buf) + strlen(addp) < BUFSIZE)
1807 
			    bufp = AppendString(bufp, addp);
20590 ripley 1808 
			else
7162 ripley 1809 
			    warningcall(call, "term names will be truncated");
1839 ihaka 1810 
			if (ll > 1) {
7162 ripley 1811 
			    if (x == R_NilValue) {
1812 
				if(strlen(buf) + 10 < BUFSIZE)
13257 maechler 1813 
				    bufp = AppendInteger(bufp, indx % ll + 1);
20590 ripley 1814 
				else
1815 
				    warningcall(call, "term names will be truncated");
7162 ripley 1816 
			    } else {
13257 maechler 1817 
				addp = CHAR(STRING_ELT(x, indx % ll));
7162 ripley 1818 
				if(strlen(buf) + strlen(addp) < BUFSIZE)
1819 
				    bufp = AppendString(bufp, addp);
20590 ripley 1820 
				else
7162 ripley 1821 
				    warningcall(call, "term names will be truncated");
1822 
			    }
1839 ihaka 1823 
			}
1824 
		    }
13257 maechler 1825 
		    indx /= ll;
2 r 1826 
		}
1839 ihaka 1827 
	    }
10172 luke 1828 
	    SET_STRING_ELT(xnames, k++, mkChar(buf));
2 r 1829 
	}
1839 ihaka 1830 
    }
2 r 1831 
 
1839 ihaka 1832 
    /* Allocate and compute the design matrix. */
2 r 1833 
 
1839 ihaka 1834 
    PROTECT(x = allocMatrix(REALSXP, n, nc));
2 r 1835 
 
1839 ihaka 1836 
    /* a) Begin with a column of 1s for the intercept. */
2 r 1837 
 
13607 maechler 1838 
    if ((jnext = jstart = intrcept) != 0) {
1839 ihaka 1839 
	for (i = 0; i < n; i++) {
1840 
	    REAL(x)[i] = 1.0;
2 r 1841 
	}
1839 ihaka 1842 
    }
454 maechler 1843 
 
1839 ihaka 1844 
    /* b) Now loop over the model terms */
2 r 1845 
 
2460 hornik 1846 
    contrast = R_NilValue;	/* -Wall */
1839 ihaka 1847 
    for (k = 0; k < nterms; k++) {
7162 ripley 1848 
	if (k == rhs_response) continue;
13607 maechler 1849 
	for (i = 0; i < nVar; i++) {
2762 pd 1850 
	    if (INTEGER(columns)[i] == 0)
1851 
		continue;
10172 luke 1852 
	    var_i = VECTOR_ELT(variable, i);
13607 maechler 1853 
	    fik = INTEGER(factors)[i + k * nVar];
1839 ihaka 1854 
	    if (fik) {
1855 
		switch(fik) {
1856 
		case 1:
10172 luke 1857 
		    contrast = VECTOR_ELT(contr1, i);
1839 ihaka 1858 
		    break;
1859 
		case 2:
10172 luke 1860 
		    contrast = VECTOR_ELT(contr2, i);
1839 ihaka 1861 
		    break;
2 r 1862 
		}
1839 ihaka 1863 
		if (jnext == jstart) {
1864 
		    if (INTEGER(nlevs)[i] > 0) {
15567 ripley 1865 
			int adj = isLogical(var_i)?1:0;
1858 ihaka 1866 
			firstfactor(&REAL(x)[jstart * n], n, jnext - jstart,
1839 ihaka 1867 
				    REAL(contrast), nrows(contrast),
15567 ripley 1868 
				    ncols(contrast), INTEGER(var_i)+adj);
1839 ihaka 1869 
			jnext = jnext + ncols(contrast);
1870 
		    }
1871 
		    else {
1858 ihaka 1872 
			firstvar(&REAL(x)[jstart * n], n, jnext - jstart,
1839 ihaka 1873 
				 REAL(var_i), n, ncols(var_i));
1874 
			jnext = jnext + ncols(var_i);
1875 
		    }
1876 
		}
1877 
		else {
1878 
		    if (INTEGER(nlevs)[i] > 0) {
15567 ripley 1879 
			int adj = isLogical(var_i)?1:0;
1858 ihaka 1880 
			addfactor(&REAL(x)[jstart * n], n, jnext - jstart,
1839 ihaka 1881 
				  REAL(contrast), nrows(contrast),
15567 ripley 1882 
				  ncols(contrast), INTEGER(var_i)+adj);
1839 ihaka 1883 
			jnext = jnext + (jnext - jstart)*(ncols(contrast) - 1);
1884 
		    }
1885 
		    else {
1858 ihaka 1886 
			addvar(&REAL(x)[jstart * n], n, jnext - jstart,
1839 ihaka 1887 
			       REAL(var_i), n, ncols(var_i));
1858 ihaka 1888 
			jnext = jnext + (jnext - jstart) * (ncols(var_i) - 1);
1839 ihaka 1889 
		    }
1890 
		}
1891 
	    }
2 r 1892 
	}
1839 ihaka 1893 
	jstart = jnext;
1894 
    }
1895 
    PROTECT(tnames = allocVector(VECSXP, 2));
10172 luke 1896 
    SET_VECTOR_ELT(tnames, 0, rnames);
1897 
    SET_VECTOR_ELT(tnames, 1, xnames);
1839 ihaka 1898 
    setAttrib(x, R_DimNamesSymbol, tnames);
1899 
    setAttrib(x, install("assign"), assign);
1900 
    UNPROTECT(13);
1901 
    return x;
2 r 1902 
}