1# -*- coding: utf-8 -*-
2"""
3Created on Thu May 15 16:36:05 2014
4
5Author: Josef Perktold
6License: BSD-3
7
8"""
9
10import numpy as np
11
12
13class LinearConstraints(object):
14    """Class to hold linear constraints information
15
16    Affine constraints are defined as ``R b = q` where `R` is the constraints
17    matrix and `q` are the constraints values and `b` are the parameters.
18
19    This is in analogy to patsy's LinearConstraints class but can be pickled.
20
21    Parameters
22    ----------
23    constraint_matrix : ndarray
24        R matrix, 2-dim with number of columns equal to the number of
25        parameters. Each row defines one constraint.
26    constraint_values : ndarray
27        1-dim array of constant values
28    variable_names : list of strings
29        parameter names, used only for display
30    kwds : keyword arguments
31        keywords are attached to the instance.
32
33    """
34
35    def __init__(self, constraint_matrix, constraint_values,
36                 variable_names, **kwds):
37
38        self.constraint_matrix = constraint_matrix
39        self.constraint_values = constraint_values
40        self.variable_names = variable_names
41
42        # alias for patsy compatibility
43        self.coefs = constraint_matrix
44        self.constants = constraint_values
45
46        self.__dict__.update(kwds)
47        self.tuple = (self.constraint_matrix, self.constraint_values)
48
49    def __iter__(self):
50        yield from self.tuple
51
52    def __getitem__(self, idx):
53        return self.tuple[idx]
54
55    def __str__(self):
56        def prod_string(v, name):
57            v = np.abs(v)
58            if v != 1:
59                ss = str(v) + " * " + name
60            else:
61                ss = name
62            return ss
63
64        constraints_strings = []
65        for r, q in zip(*self):
66            ss = []
67            for v, name in zip(r, self.variable_names):
68                if v != 0 and ss == []:
69                    ss += prod_string(v, name)
70                elif v > 0:
71                    ss += " + " + prod_string(v, name)
72                elif v < 0:
73                    ss += " - " + prod_string(np.abs(v), name)
74            ss += " = " + str(q.item())
75            constraints_strings.append(''.join(ss))
76
77        return '\n'.join(constraints_strings)
78
79    @classmethod
80    def from_patsy(cls, lc):
81        """class method to create instance from patsy instance
82
83        Parameters
84        ----------
85        lc : instance
86            instance of patsy LinearConstraint, or other instances that have
87            attributes ``lc.coefs, lc.constants, lc.variable_names``
88
89        Returns
90        -------
91        instance of this class
92
93        """
94        return cls(lc.coefs, lc.constants, lc.variable_names)
95
96
97class TransformRestriction(object):
98    """Transformation for linear constraints `R params = q`
99
100    Note, the transformation from the reduced to the full parameters is an
101    affine and not a linear transformation if q is not zero.
102
103
104    Parameters
105    ----------
106    R : array_like
107        Linear restriction matrix
108    q : arraylike or None
109        values of the linear restrictions
110
111
112    Notes
113    -----
114    The reduced parameters are not sorted with respect to constraints.
115
116    TODO: error checking, eg. inconsistent constraints, how?
117
118    Inconsistent constraints will raise an exception in the calculation of
119    the constant or offset. However, homogeneous constraints, where q=0, will
120    can have a solution where the relevant parameters are constraint to be
121    zero, as in the following example::
122
123        b1 + b2 = 0 and b1 + 2*b2 = 0, implies that b2 = 0.
124
125    The transformation applied from full to reduced parameter space does not
126    raise and exception if the constraint does not hold.
127    TODO: maybe change this, what's the behavior in this case?
128
129
130    The `reduce` transform is applied to the array of explanatory variables,
131    `exog`, when transforming a linear model to impose the constraints.
132    """
133
134    def __init__(self, R, q=None):
135
136        # The calculations are based on Stata manual for makecns
137        R = self.R = np.atleast_2d(R)
138        if q is not None:
139            q = self.q = np.asarray(q)
140
141        k_constr, k_vars = R.shape
142        self.k_constr, self.k_vars = k_constr, k_vars
143        self.k_unconstr = k_vars - k_constr
144
145        m = np.eye(k_vars) - R.T.dot(np.linalg.pinv(R).T)
146        evals, evecs = np.linalg.eigh(m)
147
148        # This normalizes the transformation so the larges element is 1.
149        # It makes it easier to interpret simple restrictions, e.g. b1 + b2 = 0
150        # TODO: make this work, there is something wrong, does not round-trip
151        #       need to adjust constant
152        #evecs_maxabs = np.max(np.abs(evecs), 0)
153        #evecs = evecs / evecs_maxabs
154
155        self.evals = evals
156        self.evecs = evecs # temporarily attach as attribute
157        L = self.L = evecs[:, :k_constr]
158        self.transf_mat = evecs[:, k_constr:]
159
160        if q is not None:
161            # use solve instead of inv
162            #self.constant = q.T.dot(np.linalg.inv(L.T.dot(R.T)).dot(L.T))
163            try:
164                self.constant = q.T.dot(np.linalg.solve(L.T.dot(R.T), L.T))
165            except np.linalg.linalg.LinAlgError as e:
166                raise ValueError('possibly inconsistent constraints. error '
167                                 'generated by\n%r' % (e, ))
168        else:
169            self.constant = 0
170
171    def expand(self, params_reduced):
172        """transform from the reduced to the full parameter space
173
174        Parameters
175        ----------
176        params_reduced : array_like
177            parameters in the transformed space
178
179        Returns
180        -------
181        params : array_like
182            parameters in the original space
183
184        Notes
185        -----
186        If the restriction is not homogeneous, i.e. q is not equal to zero,
187        then this is an affine transform.
188        """
189        params_reduced = np.asarray(params_reduced)
190        return self.transf_mat.dot(params_reduced.T).T + self.constant
191
192    def reduce(self, params):
193        """transform from the full to the reduced parameter space
194
195        Parameters
196        ----------
197        params : array_like
198            parameters or data in the original space
199
200        Returns
201        -------
202        params_reduced : array_like
203            parameters in the transformed space
204
205        This transform can be applied to the original parameters as well
206        as to the data. If params is 2-d, then each row is transformed.
207        """
208        params = np.asarray(params)
209        return params.dot(self.transf_mat)
210
211
212def transform_params_constraint(params, Sinv, R, q):
213    """find the parameters that statisfy linear constraint from unconstrained
214
215    The linear constraint R params = q is imposed.
216
217    Parameters
218    ----------
219    params : array_like
220        unconstrained parameters
221    Sinv : ndarray, 2d, symmetric
222        covariance matrix of the parameter estimate
223    R : ndarray, 2d
224        constraint matrix
225    q : ndarray, 1d
226        values of the constraint
227
228    Returns
229    -------
230    params_constraint : ndarray
231        parameters of the same length as params satisfying the constraint
232
233    Notes
234    -----
235    This is the exact formula for OLS and other linear models. It will be
236    a local approximation for nonlinear models.
237
238    TODO: Is Sinv always the covariance matrix?
239    In the linear case it can be (X'X)^{-1} or sigmahat^2 (X'X)^{-1}.
240
241    My guess is that this is the point in the subspace that satisfies
242    the constraint that has minimum Mahalanobis distance. Proof ?
243    """
244
245    rsr = R.dot(Sinv).dot(R.T)
246
247    reduction = Sinv.dot(R.T).dot(np.linalg.solve(rsr, R.dot(params) - q))
248    return params - reduction
249
250
251def fit_constrained(model, constraint_matrix, constraint_values,
252                    start_params=None, fit_kwds=None):
253    # note: self is model instance
254    """fit model subject to linear equality constraints
255
256    The constraints are of the form   `R params = q`
257    where R is the constraint_matrix and q is the vector of constraint_values.
258
259    The estimation creates a new model with transformed design matrix,
260    exog, and converts the results back to the original parameterization.
261
262
263    Parameters
264    ----------
265    model: model instance
266        An instance of a model, see limitations in Notes section
267    constraint_matrix : array_like, 2D
268        This is R in the linear equality constraint `R params = q`.
269        The number of columns needs to be the same as the number of columns
270        in exog.
271    constraint_values :
272        This is `q` in the linear equality constraint `R params = q`
273        If it is a tuple, then the constraint needs to be given by two
274        arrays (constraint_matrix, constraint_value), i.e. (R, q).
275        Otherwise, the constraints can be given as strings or list of
276        strings.
277        see t_test for details
278    start_params : None or array_like
279        starting values for the optimization. `start_params` needs to be
280        given in the original parameter space and are internally
281        transformed.
282    **fit_kwds : keyword arguments
283        fit_kwds are used in the optimization of the transformed model.
284
285    Returns
286    -------
287    params : ndarray ?
288        estimated parameters (in the original parameterization
289    cov_params : ndarray
290        covariance matrix of the parameter estimates. This is a reverse
291        transformation of the covariance matrix of the transformed model given
292        by `cov_params()`
293        Note: `fit_kwds` can affect the choice of covariance, e.g. by
294        specifying `cov_type`, which will be reflected in the returned
295        covariance.
296    res_constr : results instance
297        This is the results instance for the created transformed model.
298
299
300    Notes
301    -----
302    Limitations:
303
304    Models where the number of parameters is different from the number of
305    columns of exog are not yet supported.
306
307    Requires a model that implement an offset option.
308    """
309    self = model   # internal alias, used for methods
310    if fit_kwds is None:
311        fit_kwds = {}
312
313    R, q = constraint_matrix, constraint_values
314    endog, exog = self.endog, self.exog
315
316    transf = TransformRestriction(R, q)
317
318    exogp_st = transf.reduce(exog)
319
320    offset = exog.dot(transf.constant.squeeze())
321    if hasattr(self, 'offset'):
322        offset += self.offset
323
324    if start_params is not None:
325        start_params =  transf.reduce(start_params)
326
327    #need copy, because we do not want to change it, we do not need deepcopy
328    import copy
329    init_kwds = copy.copy(self._get_init_kwds())
330
331    # TODO: refactor to combine with above or offset_all
332    if 'offset' in init_kwds:
333        del init_kwds['offset']
334
335    # using offset as keywords is not supported in all modules
336    mod_constr = self.__class__(endog, exogp_st, offset=offset, **init_kwds)
337    res_constr = mod_constr.fit(start_params=start_params, **fit_kwds)
338    params_orig = transf.expand(res_constr.params).squeeze()
339    cov_params = transf.transf_mat.dot(res_constr.cov_params()).dot(transf.transf_mat.T)
340
341    return params_orig, cov_params, res_constr
342
343
344def fit_constrained_wrap(model, constraints, start_params=None, **fit_kwds):
345    """fit_constraint that returns a results instance
346
347    This is a development version for fit_constrained methods or
348    fit_constrained as standalone function.
349
350    It will not work correctly for all models because creating a new
351    results instance is not standardized for use outside the `fit` methods,
352    and might need adjustements for this.
353
354    This is the prototype for the fit_constrained method that has been added
355    to Poisson and GLM.
356    """
357
358    self = model  # alias for use as method
359
360    #constraints = (R, q)
361    # TODO: temporary trailing underscore to not overwrite the monkey
362    #       patched version
363    # TODO: decide whether to move the imports
364    from patsy import DesignInfo
365    # we need this import if we copy it to a different module
366    #from statsmodels.base._constraints import fit_constrained
367
368    # same pattern as in base.LikelihoodModel.t_test
369    lc = DesignInfo(self.exog_names).linear_constraint(constraints)
370    R, q = lc.coefs, lc.constants
371
372    # TODO: add start_params option, need access to tranformation
373    #       fit_constrained needs to do the transformation
374    params, cov, res_constr = fit_constrained(self, R, q,
375                                              start_params=start_params,
376                                              fit_kwds=fit_kwds)
377    #create dummy results Instance, TODO: wire up properly
378    res = self.fit(start_params=params, maxiter=0,
379                   warn_convergence=False)  # we get a wrapper back
380    res._results.params = params
381    res._results.cov_params_default = cov
382    cov_type = fit_kwds.get('cov_type', 'nonrobust')
383    if cov_type == 'nonrobust':
384        res._results.normalized_cov_params = cov / res_constr.scale
385    else:
386        res._results.normalized_cov_params = None
387
388    k_constr = len(q)
389    res._results.df_resid += k_constr
390    res._results.df_model -= k_constr
391    res._results.constraints = LinearConstraints.from_patsy(lc)
392    res._results.k_constr = k_constr
393    res._results.results_constrained = res_constr
394    return res
395