Performs a grid search over a space of parameters for a given
        objective

        Warnings
        --------
        ``gridsearch`` is lazy and will not remove useless combinations
        from the search space, eg.

        >>> n_splines=np.arange(5,10), fit_splines=[True, False]

        will result in 10 loops, of which 5 are equivalent because
        ``fit_splines = False``

        Also, it is not recommended to search over a grid that alternates
        between known scales and unknown scales, as the scores of the
        candidate models will not be comparable.

        Parameters
        ----------
        X : array-like
          input data of shape (n_samples, m_features)

        y : array-like
          label data of shape (n_samples,)

        weights : array-like shape (n_samples,), optional
            sample weights

        return_scores : boolean, optional
            whether to return the hyperpamaters and score for each element
            in the grid

        keep_best : boolean, optional
            whether to keep the best GAM as self.

        objective : {'auto', 'AIC', 'AICc', 'GCV', 'UBRE'}, optional
            Metric to optimize.
            If `auto`, then grid search will optimize `GCV` for models with unknown
            scale and `UBRE` for models with known scale.

        progress : bool, optional
            whether to display a progress bar

        **kwargs
            pairs of parameters and iterables of floats, or
            parameters and iterables of iterables of floats.

            If no parameter are specified, ``lam=np.logspace(-3, 3, 11)`` is used.
            This results in a 11 points, placed diagonally across lam space.

            If grid is iterable of iterables of floats,
            the outer iterable must have length ``m_features``.
            the cartesian product of the subgrids in the grid will be tested.

            If grid is a 2d numpy array,
            each row of the array will be tested.

            The method will make a grid of all the combinations of the parameters
            and fit a GAM to each combination.


        Returns
        -------
        if ``return_scores=True``:
            model_scores: dict containing each fitted model as keys and corresponding
            objective scores as values
        else:
            self: ie possibly the newly fitted model

        Examples
        --------
        For a model with 4 terms, and where we expect 4 lam values,
        our search space for lam must have 4 dimensions.

        We can search the space in 3 ways:

        1. via cartesian product by specifying the grid as a list.
        our grid search will consider ``11 ** 4`` points:

        >>> lam = np.logspace(-3, 3, 11)
        >>> lams = [lam] * 4
        >>> gam.gridsearch(X, y, lam=lams)

        2. directly by specifying the grid as a np.ndarray.
        This is useful for when the dimensionality of the search space
        is very large, and we would prefer to execute a randomized search:

        >>> lams = np.exp(np.random.random(50, 4) * 6 - 3)
        >>> gam.gridsearch(X, y, lam=lams)

        3. copying grids for parameters with multiple dimensions.
        if we specify a 1D np.ndarray for lam, we are implicitly testing the
        space where all points have the same value

        >>> gam.gridsearch(lam=np.logspace(-3, 3, 11))

        is equivalent to:

        >>> lam = np.logspace(-3, 3, 11)
        >>> lams = np.array([lam] * 4)
        >>> gam.gridsearch(X, y, lam=lams)

def gridsearch(self, X, y, weights=None, return_scores=False,
                   keep_best=True, objective='auto', progress=True,
                   **param_grids):
        """
        Performs a grid search over a space of parameters for a given
        objective

        Warnings
        --------
        ``gridsearch`` is lazy and will not remove useless combinations
        from the search space, eg.

        >>> n_splines=np.arange(5,10), fit_splines=[True, False]

        will result in 10 loops, of which 5 are equivalent because
        ``fit_splines = False``

        Also, it is not recommended to search over a grid that alternates
        between known scales and unknown scales, as the scores of the
        candidate models will not be comparable.

        Parameters
        ----------
        X : array-like
          input data of shape (n_samples, m_features)

        y : array-like
          label data of shape (n_samples,)

        weights : array-like shape (n_samples,), optional
            sample weights

        return_scores : boolean, optional
            whether to return the hyperpamaters and score for each element
            in the grid

        keep_best : boolean, optional
            whether to keep the best GAM as self.

        objective : {'auto', 'AIC', 'AICc', 'GCV', 'UBRE'}, optional
            Metric to optimize.
            If `auto`, then grid search will optimize `GCV` for models with unknown
            scale and `UBRE` for models with known scale.

        progress : bool, optional
            whether to display a progress bar

        **kwargs
            pairs of parameters and iterables of floats, or
            parameters and iterables of iterables of floats.

            If no parameter are specified, ``lam=np.logspace(-3, 3, 11)`` is used.
            This results in a 11 points, placed diagonally across lam space.

            If grid is iterable of iterables of floats,
            the outer iterable must have length ``m_features``.
            the cartesian product of the subgrids in the grid will be tested.

            If grid is a 2d numpy array,
            each row of the array will be tested.

            The method will make a grid of all the combinations of the parameters
            and fit a GAM to each combination.


        Returns
        -------
        if ``return_scores=True``:
            model_scores: dict containing each fitted model as keys and corresponding
            objective scores as values
        else:
            self: ie possibly the newly fitted model

        Examples
        --------
        For a model with 4 terms, and where we expect 4 lam values,
        our search space for lam must have 4 dimensions.

        We can search the space in 3 ways:

        1. via cartesian product by specifying the grid as a list.
        our grid search will consider ``11 ** 4`` points:

        >>> lam = np.logspace(-3, 3, 11)
        >>> lams = [lam] * 4
        >>> gam.gridsearch(X, y, lam=lams)

        2. directly by specifying the grid as a np.ndarray.
        This is useful for when the dimensionality of the search space
        is very large, and we would prefer to execute a randomized search:

        >>> lams = np.exp(np.random.random(50, 4) * 6 - 3)
        >>> gam.gridsearch(X, y, lam=lams)

        3. copying grids for parameters with multiple dimensions.
        if we specify a 1D np.ndarray for lam, we are implicitly testing the
        space where all points have the same value

        >>> gam.gridsearch(lam=np.logspace(-3, 3, 11))

        is equivalent to:

        >>> lam = np.logspace(-3, 3, 11)
        >>> lams = np.array([lam] * 4)
        >>> gam.gridsearch(X, y, lam=lams)
        """
        # check if model fitted
        if not self._is_fitted:
            self._validate_params()
            self._validate_data_dep_params(X)

        y = check_y(y, self.link, self.distribution, verbose=self.verbose)
        X = check_X(X, verbose=self.verbose)
        check_X_y(X, y)

        if weights is not None:
            weights = np.array(weights).astype('f').ravel()
            weights = check_array(weights, name='sample weights',
                                  ndim=1, verbose=self.verbose)
            check_lengths(y, weights)
        else:
            weights = np.ones_like(y).astype('float64')

        # validate objective
        if objective not in ['auto', 'GCV', 'UBRE', 'AIC', 'AICc']:
            raise ValueError("objective mut be in "\
                             "['auto', 'GCV', 'UBRE', 'AIC', 'AICc'], '\
                             'but found objective = {}".format(objective))

        # check objective
        if self.distribution._known_scale:
            if objective == 'GCV':
                raise ValueError('GCV should be used for models with'\
                                 'unknown scale')
            if objective == 'auto':
                objective = 'UBRE'

        else:
            if objective == 'UBRE':
                raise ValueError('UBRE should be used for models with '\
                                 'known scale')
            if objective == 'auto':
                objective = 'GCV'

        # if no params, then set up default gridsearch
        if not bool(param_grids):
            param_grids['lam'] = np.logspace(-3, 3, 11)

        # validate params
        admissible_params = list(self.get_params()) + self._plural
        params = []
        grids = []
        for param, grid in list(param_grids.items()):

            # check param exists
            if param not in (admissible_params):
                raise ValueError('unknown parameter: {}'.format(param))

            # check grid is iterable at all
            if not (isiterable(grid) and (len(grid) > 1)): \
                raise ValueError('{} grid must either be iterable of '
                                 'iterables, or an iterable of lengnth > 1, '\
                                 'but found {}'.format(param, grid))

            # prepare grid
            if any(isiterable(g) for g in grid):

                # get required parameter shape
                target_len = len(flatten(getattr(self, param)))

                # check if cartesian product needed
                cartesian = (not isinstance(grid, np.ndarray) or grid.ndim != 2)

                # build grid
                grid = [np.atleast_1d(g) for g in grid]

                # check chape
                msg = '{} grid should have {} columns, '\
                      'but found grid with {} columns'.format(param, target_len, len(grid))
                if cartesian:
                    if len(grid) != target_len:
                        raise ValueError(msg)
                    grid = combine(*grid)

                if not all([len(subgrid) == target_len for subgrid in grid]):
                    raise ValueError(msg)

            # save param name and grid
            params.append(param)
            grids.append(grid)

        # build a list of dicts of candidate model params
        param_grid_list = []
        for candidate in combine(*grids):
            param_grid_list.append(dict(zip(params,candidate)))

        # set up data collection
        best_model = None # keep the best model
        best_score = np.inf
        scores = []
        models = []

        # check if our model has been fitted already and store it
        if self._is_fitted:
            models.append(self)
            scores.append(self.statistics_[objective])

            # our model is currently the best
            best_model = models[-1]
            best_score = scores[-1]

        # make progressbar optional
        if progress:
            pbar = ProgressBar()
        else:
            pbar = lambda x: x

        # loop through candidate model params
        for param_grid in pbar(param_grid_list):
            try:
                # try fitting
                # define new model
                gam = deepcopy(self)
                gam.set_params(self.get_params())
                gam.set_params(**param_grid)

                # warm start with parameters from previous build
                if models:
                    coef = models[-1].coef_
                    gam.set_params(coef_=coef, force=True, verbose=False)
                gam.fit(X, y, weights)

            except ValueError as error:
                msg = str(error) + '\non model with params:\n' + str(param_grid)
                msg += '\nskipping...\n'
                if self.verbose:
                    warnings.warn(msg)
                continue

            # record results
            models.append(gam)
            scores.append(gam.statistics_[objective])

            # track best
            if scores[-1] < best_score:
                best_model = models[-1]
                best_score = scores[-1]

        # problems
        if len(models) == 0:
            msg = 'No models were fitted.'
            if self.verbose:
                warnings.warn(msg)
            return self

        # copy over the best
        if keep_best:
            self.set_params(deep=True,
                            force=True,
                            **best_model.get_params(deep=True))
        if return_scores:
            return OrderedDict(zip(models, scores))
        else:
            return self