We present the most exhaustive and accurate inferences of the internal solar rotation rate and its evolution during solar cycle 23. A full solar cycle of MDI observations have been analyzed using our state of the art fitting methodology. Time series of various lengths have been fitted, from a single 4608-day long epoch (64 times 72 day or 12.6 yr) down to 64 separate segments for the "traditional" 72-day long epochs. We used time series of spherical harmonic coefficients computed by the MDI group but using an improved spatial decomposition. This decomposition now includes our best estimate of the image plate scale and of the MDI instrumental image distortion. The leakage matrix used for the fitting includes the distortion of the eigenfunctions by the solar differential rotation, and the undistorted leakage matrix was itself carefully reviewed and independently recomputed. Rotation inversions were carried out for all the available mode sets, fitted for that epoch and all available segments, including the MDI and GONG "pipe-line" sets. The improved inversions we used is an iterative methodology based on a least-squares regularization. It also implement a model grid optimization derived from the actual information in the input set. This optimized model grid is itself irregular, namely with a variable number of latitudes at different depths. We not only present the most accurate mean rotation rate, but also how its derivation may still be affected by uncertainties in the mode fitting (in particular the leakage matrix). We also focus on the change of the rotation rate with activity levels and how well these changes are significantly assessed at higher latitudes as well as deeper in the solar interior, down to the base of the convection zone.