def generate_random_lifetimes(hazard_rates, timelines, size=1, censor=None):
    """
    Based on the hazard rates, compute random variables from the survival function
      hazard_rates: (n,t) array of hazard rates
      timelines: (t,) the observation times
      size: the number to return, per hardard rate
      censor: If True, adds uniform censoring between timelines.max() and  0
              If a postive number, censors all events above that value.
              If (n,) np.array >=0 , censor elementwise.


    Returns
    -------
      survival_times: (size,n) array of random variables.
      (optional) censorship: if censor is true, returns (size,n) array with bool True
         if the death was observed (not right-censored)
    """
    n = hazard_rates.shape[1]
    survival_times = np.empty((n, size))
    cumulative_hazards = cumulative_integral(hazard_rates.values, timelines).T

    for i in range(size):
        u = random.rand(n, 1)
        e = -np.log(u)
        v = (e - cumulative_hazards) < 0
        cross = v.argmax(1)
        survival_times[:, i] = timelines[cross]
        survival_times[cross == 0, i] = np.inf

    if censor is not None:
        if isinstance(censor, bool):
            T = timelines.max()
            rv = T * random.uniform(size=survival_times.shape)
        else:
            rv = censor

        observed = np.less_equal(survival_times, rv)
        survival_times = np.minimum(rv, survival_times)
        return survival_times.T, observed.T
    else:
        return survival_times