def stepwise_selection(variables, train_model, score_model, direction=
    'both', verbose=True):
    """ Variable selection using forward and/or backward selection 
    
    Input:
        variables: complete list of variables to consider in model building
        train_model: function that returns a fitted model for a given set of variables
        score_model: function that returns the score of a model; better models have lower scores
        direction: use it to limit stepwise selection to either 'forward' or 'backward'
    
    Returns:
        (best_model, best_variables) 
    """
    FORWARD = 'forward'
    BACKWARD = 'backward'
    directions = [FORWARD, BACKWARD]
    if direction.lower() == FORWARD:
        directions = [FORWARD]
    if direction.lower() == BACKWARD:
        directions = [BACKWARD]
    best_variables = [] if 'forward' in directions else list(variables)
    best_model = train_model(best_variables)
    best_score = score_model(best_model, best_variables)
    if verbose:
        print('Variables: ' + ', '.join(variables))
        print('Start: score={:.2f}, constant'.format(best_score))
    while True:
        step = [(best_score, None, best_model, 'unchanged')]
        if FORWARD in directions:
            for variable in variables:
                if variable in best_variables:
                    continue
                step_var = list(best_variables)
                step_var.append(variable)
                step_model = train_model(step_var)
                step_score = score_model(step_model, step_var)
                step.append((step_score, variable, step_model, 'add'))
        if 'backward' in directions:
            for variable in best_variables:
                step_var = list(best_variables)
                step_var.remove(variable)
                step_model = train_model(step_var)
                step_score = score_model(step_model, step_var)
                step.append((step_score, variable, step_model, 'remove'))
        step.sort(key=lambda x: x[0])
        best_score, chosen_variable, best_model, direction = step[0]
        if verbose:
            print('Step: score={:.2f}, {} {}'.format(best_score, direction,
                chosen_variable))
        if chosen_variable is None:
            break
        if direction == 'add':
            best_variables.append(chosen_variable)
        else:
            best_variables.remove(chosen_variable)
    return best_model, best_variables
