TODO

* Add missing data handling to the just-pass-in-a-matrix bit of the high-level API

* Add parallel array handling to build_design_matrices

* Add parallel array handling of some sort to high-level API...

* Refactor build so that there are two stages
- first stage takes a set of factor evaluators, and returns a set of
  evaluated columns
- second stage handles interactions and categorical coding and assembles
  these together into design matrices
use case: any model where you actually want to get categorical data
out (like multinomial regression with a factor on the LHS, or CART
with factors on the right-hand side)
** first stage should also handle other "parallel" data, like weights, which need to participate in the missingness calculations
** possibly also support a "subset=" argument at this stage
** and for parallel vectors and subset=, allow a string as a value, and if seen then evaluate it as python code in the same context as formula data (like R's subset=(MyCol > 10))
** And do NaN/mask/missing data handling at this stage
*** Imputation?
*** numpy.ma

* Better NaN/masks/missing data handling in transforms. I think the
current ones will just blow up if there are any NaNs. (The previous
entry is about handling the term "x" where x has NAs; this entry is
about handling "center(x)" where x has NAs.) R's solution to this is
that scale(x) simply unconditionally ignores NAs when computing the
mean, regardless of the overall setting of na.action. That seems
reasonable...

* Advocacy
Potential users?
- statsmodels
- PyMC has a (closed) ticket requesting such features:
    http://code.google.com/p/pymc/issues/detail?id=162
- nipy, though they have their own thing...
- sklearn, which has regression (and might find it useful otherwise!)

* Do something smarter with mismatched pandas indexes
Right now we're conservative -- if you do ~ x + y and x and y don't
have identical indexes, then that's an error. It's possible we should
do something cleverer, though. Perhaps we should merge them somehow
(pandas.concat(..., join="outer")?). (This of course would require
that *all* items have indexes, though; right now you can mix plain
ndarrays and pandas objects.)

* Improve EvalFactor's stateful transform handling to follow . lookups
right now it can only detect stateful transforms when they are called
directly like
  scale(x)
but not if referenced through some module like
  mylib.scale(x)
In general we don't even want to try handling every possible function
lookup syntax (see next item for a safety check for that), but we
should allow for people to distribute non-builtin stateful
transforms.

* As a safety check for non-stateful transforms, we should always
evaluate each formula on just the first row of data alone, and make
sure the result matches what we got when evaluating it vectorized
(i.e., confirm f(x[0]) == f(x)[0], where f is our transform. However,
this is kind of tricky given that x might be pulled out of the
environment, the 'data' dict might have arbitrary objects,
etc. Hmm. Maybe intercept variable lookups and just munge those? This
is easy to do if someone's passing in a structured array or dataframe
and pulling all their data from it, or even if they use a dict with
well-behaved columns. But the problem is when people do things like:

  In [1]: logx = np.log(data["x"])

  # refers to data["y"] and logx together
  In [2]: lm("y ~ logx", data)

* More contrast tools
- Some sort of symbolic tools for user-defined contrasts -- take the
  comparisons that people want to compute in terms of linear
  combinations of level names, convert that to a matrix and do the
  pinv dance? We have the linear_contrast code already, but that's for
  describing constraints in terms of the coefficients you have -- it
  seems like people want to be able to describe constraints in terms
  of... I'm not sure what. Group means? The coefficients they could
  have had if they'd fit some other model? (Presumably the
  all-full-rank-dummy-coding-all-the-time model.) If I can ever figure
  out what this is (it has something to do with "estimable contrasts")
  then I'll implement it.
- Short-hands for Type II, Type III, and "remove this term and
  everything marginal to it" contrast tests?
  Might need to figure out the trick that car::Anova uses to do
  efficient Type II tests with two contrast matrices.
- Understand how coding matters for Type-III ANOVA. The tabs I had
  open last time I was looking at this:
    http://goanna.cs.rmit.edu.au/~fscholer/anova.php
    http://www.mail-archive.com/r-help@stat.math.ethz.ch/msg69781.html
    https://stat.ethz.ch/pipermail/r-help/2007-October/143047.html
    http://www.uni-kiel.de/psychologie/dwoll/r/ssTypes.php

* A good way to support magic functions like mgcv's s().
statsmodels wants this for things like
  y ~ arima(2, 3)
  y ~ garch(1, 1)
the cheap trick way of doing it is:
  class ArimaModelType(object):
    __patsy_magic__ = True
    ...
  def arima(n, m):
    return ArimaModelType(n, m)
and then in the factor type sniffing code detect these things and
separate them out from "real" factors. 

* make sure that pickling works
- And make sure that if we allow it at all, then it's sustainable!
  i.e. we'll be able to guarantee that if people pickle a ModelDesc or
  Design or whatever now, then they'll be able to get it back later.

* Should EvalEnvironment.capture make a copy of the scope dictionaries?
- The effect would be to prevent later changes in the enclosing scope
  from affecting predictions. Of course, we probably don't want to
  make a *deep* copy of the scope, so there's still no guarantees --
  changes to mutable objects within that scope would still be
  visible. Perhaps we *could* get away with making a deep copy of all
  mutable objects that are accessed during the initial build,
  though... I think we'd need to special-case and ignore any READONLY
  ndarrays, as a safety valve for people who have a giant data-set
  they're referring to. of course, even a deep copy isn't enough --
  they could call an immutable function which accesses mutable state.
- Josef points out that in long-running REPLs people often need to del
  local variables to let memory be released, and if all their formulas
  are going and making shallow copies of the environment then this
  will be impossible. So making a shallow copy is probably out.
- The other approach would be to extend the state dependency checking
  that we already want to do (to catch undeclared stateful
  transforms), and have it not only check that building an isolated
  row of data gives the same result as building the full list, but
  also that re-building that same row later at prediction time gives
  the same result as it did in the first place.

* Export information on which terms are marginal to which other ones
Marginality only makes sense within a numeric-interaction "bucket", so
this has to be computed in patsy.build and exported as part of
DesignMatrixColumnInfo. Then it can be used for Type II tests.

* Some way to specify the default contrast

* Support for R's magic "." term
- The "y ~ everything else" form
- The "what I had in this other ModelDesc" form (e.g., "y ~ . - a"
  to drop the 'a' predictor from an old model)
  - This will require that the formula->ModelDesc have access to the
    data or previous formula...

* More stateful transforms:
- Splines
- 'cut': numeric->factor by quantile dichotimization
- Orthogonal polynomials
- 'code': takes a Categorical (or coerces to one), and optionally
  a contrast, and and does the standard contrast-coding. And
  possibly this should replace _CatFactorEvaluator...

* Support for building sparse model matrices directly. (This should
be pretty straightforward when it comes to exploiting the intrinsic
sparsity of categorical factors; numeric factors that evaluate to a
sparse matrix directly might be slightly more complicated.)

* Real testing/support for formula syntax extensions
The tricky part here is making sure we produce something useful.
Use cases:
- multinomial log-linear modelling
- see below
Prior art:
- R package "lmer" interpets formulas like
    y ~ x1 + x2 + (1 | foo) + (1 + x | bar)
- The R [[http://cran.r-project.org/web/packages/Formula/vignettes/Formula.pdf][Formula]] package, which has two features:
  - you can write multivariate responses, like y1 + y2 ~ ... (in stock
    R, this is interpreted as addition (!)).
  - you can write multiple "parts" on each side, separated
    by |. Basically these are treated as a list of design matrix
    specifications, and there are ways to pull out the first, second
    etc. on each side.
  - R package "plm" uses Formula to allow formulas like:
       y ~ x1 + x2
       y ~ x1 + x2 | x3
       y ~ x1 + x2 | . + x3
    where the second part specifies "instrumental variables". I can't
    tell if the second part has an implicit intercept.
  - R package "frontier" uses Formula in a similar way, allowing
    formulas like
       y ~ x1 + x2
       y ~ x1 + x2 | x3
    where the first form computes a "error components frontier" and the
    latter computes an "efficiency effects frontier" (where the part
    after the | are regresses "used to explain the efficiency levels (Z
    variables)"). The part after the bar does have an implicit
    intercept.
  - package AER uses this in its "ivreg" command, which seems similar
    to plm. An example makes clear that "y ~ . | x1 + x2" works, and
    presumably the "." means the same thing as it would in "y ~ ." for
    lm.
  - package betareg does "beta regression", and a formula like "y ~
    x1 | x2" states that "x1" should be used for the "mean submodel"
    and "x2" should be used for the "precision submodel". Its betatree
    function extends this further to "y ~ x1 | x2 | c1 + c2" where
    "c1", "c2" are "partitioning variables". AFAICT this means that it
    does basically does a CART-style tree division of the data based
    on c1, c2, and then fits beta regression models x1 | x2 on each
    subset.
  - package "fdaMixed" uses formulas like
      Y | id ~ fixed | random
    where Y is a response variable, id is "a factor separating the
    samples", and fixed and random are linear models for the fixed
    and random effects. The 'id' part seems to be used to match
    multiple samples from the same random effects group?
  - package "growcurves" allows "y ~ fixed | random". If there is
    no |, then there is a second argument (random.only) which is
    consulted to determine whether the sole RHS argument is fixed or
    random. (Maybe 'y ~ x1 + x2 + random(x3 + x3)' would be a better
    syntax?)
  - package "games" uses a syntax like "y ~ x1 + x2 | 0 | x3 |
    z". There is another version with 8 entries instead of 4.
  - package "metafor" does effect-size calculations using the syntax
    "outcome ~ group | study" where each entry has to be a 2-level
    factor. (And the 'weights' argument gives the actual numbers.)
  - package "mhurdle" seems to describe a kind of multi-step process
    via three-part formulas
      y ~ x1 | x2 | x3
    where "the first part describes the selection process if any, the
    second part the regression equation, and the third part the
    purchase infrequency process". You can fill in 0 if you want to
    assume that some process doesn't actually apply (or leave out the
    last one altogether).
  - package "mlogit" uses three-part RHS formulas to specify different
    parts of a multinomial logit model. "the first one contains the
    alternative specific variables with generic coefficient, i.e. a
    unique coefficient for all the alternatives; the second one
    contains the individual specific variables for which one
    coefficient is estimated for all the alternatives except one of
    them ; the third one contains the alternative specific variables
    with alternative specific coefficients...If a standard formula is
    writen, it is assumed that there are only alternative specific
    variables with generic coefficients."
    The second RHS termlist has an intercept by default; for the other
    two termlists any intercept is ignored in any case.
  - package "polywog" does some clever polynomial basis function
    fitting thing, and uses formulas like
      y ~ x1 + x2 | z1 + z2
    to mean basically the equivalent of
      y ~ x1*x2 + z1 + z2
    i.e., the first termlist gets a super-rich non-linear interaction
    between all its entries, and the second is just entered linearly.

* Currently we don't distinguish between ordered and unordered categorical data.
Should that change?

* how should redundancy elimination and explicit factor matrices interact?
Example: If you do 1 + C(a, mat):C(b, mat), then currently it will
expand that to 1 + C(a, mat) + C(a, mat):C(b, mat), which is going to
be weird. Probably we should notice that the .contrast attribute in
these cases does not give us the option of full- versus reduced-rank
coding, and in redundancy.py we should note that such factors cannot
be "expanded".

* Profiling/optimization. There are lots of places where I use lazy
quadratic algorithms (or even exponential, in the case of the
non-redundant coding stuff). Perhaps worse is the heavy
multiplication used unconditionally to load data into the model
matrix. I'm pretty sure that at least most of the quadratic stuff
doesn't matter because it's n^2 where n is something like the
number of factors in an interaction term (and who has hundreds of
factors interacting in one term?), but it wouldn't hurt to run some
profiles to check. I think really what I mean is just, run timeit
on a 10-variable interaction to make sure it isn't completely
annoying.

* Possible optimization: let a stateful transform's memorize_chunk
function raise Stateless to indicate that actually, ha-ha, it turns
out that it doesn't need to memorize anything after all (b/c the
relevant data turns out to be specified explicitly in *args,
**kwargs).
Actually, this would be really useful for things like splines,
which need to do expensive quantile estimation, but not if knots
are specified.
Another use case: C(something_that's_already_categorical,
contrast=...). Note that this can't be detected until we do the first
round of evaluation.
A better interface would be memorize_needed(self, *args, **kwargs).
I guess we could even have memorize_passes_needed, but eh...

* Wacky idea: make factors into an actual stateful transform (one
that takes a dict-like object and returns a matrix or Categorical)
This would require:
- adding memorize_passes support to stateful transforms
- moving the factor memorization state inside an object (so it
  wouldn't be factors that would be stateful transforms, factors would
  be factories for stateful transforms)