def recommender_function(model, X, treatment_column):
    """
    Calculate the recommender function for a set of patients based on a 
    previously fitted model. Implementation corresponds to the one 
    proposed in [1] (Eq. 6).
    
    Parameters
    ----------
    model: 
        Model that will be used to compute the log-hazards.
        It needs to be fitted previously.
    X: pandas DataFrame
        Data. Rows correspond to instances (i.e, patients).
        Columns correspond to features.
    treatment_column: string
        Which column in X corresponds to the treatment.

    Currently, it only supports comparison of two treatments.

    Returns
    -------
    rec_ij: NumPy array.
        The recommender function of all patients. For each patient:
        If rec_ij is positive, it means that treatment i lead to a higher 
        risk of death than treatment j.
        If rec_ij is negative, it means that treatment j lead to a higher
        risk of death than treatment i.
        
    References
    ----------
    [1] Katzman, Jared L., et al. "DeepSurv: personalized treatment recommender system using a Cox proportional hazards deep neural network." BMC medical research methodology 18.1 (2018): 24.
    """
    
    # Validate number of treatments.
    treatments = X[treatment_column].unique()
    n_treatments = len(treatments)
    
    if n_treatments == 1:
        raise ValueError("It is not possible to give a treatment recommendation with only one treatment value.")
    elif n_treatments > 2:
        raise ValueError(f"{n_treatments} found. Currently, only two treatments are supported for comparison.")
    
    # Create DataFrames to be used for prediction.
    X_treatment0 = X.copy(deep=True)
    X_treatment0[treatment_column] = treatments[0]
    
    X_treatment1 = X.copy(deep=True)
    X_treatment1[treatment_column] = treatments[1]
    
    # Calculate the log-hazards.
    h_i = model.predict(X_treatment0)
    h_j = model.predict(X_treatment1)
    
    # Calculate the recommender function (Eq. 6)
    rec_ij = h_i - h_j
    
    return rec_ij