There are a number of important inequalities in both Harmonic Analysis and PDE which are of the form \[ \int_{\mathbb R^n}| Tf(x)| ^pw(x)\,dx\leq C\int_{\mathbb R^n}| Sf(x)| ^p w(x)\,dx,\tag \(*\) \] where typically \(T\) is an operator with some degree of singularity (e.g., certain singular integral operator) and \(S\) is an operator which is easier to handle (e.g., a maximal operator), and \(w\) is in some class of weights. As it is well known, the usual technique for proving such results is to establish a good-\(\lambda\) inequality between \(T\) and \(S\). Inspired by the extrapolation theory for \(A_p\) weights discovered by Rubio de Francia, D. Cruz-Uribe, J. M. Martell and C. Pérez [J. Funct. Anal. 213, No. 2, 412–439 (2004; Zbl 1052.42016)] presented another approach to derive inequalities like (\(\ast\)) without using the good-\(\lambda\) technique. Namely, assume that (\(\ast\)) holds for some fixed exponent \(p_0\in (0,\infty)\) and for all \(w\in A_\infty\) and all (reasonable) functions \(f\) for which the left-hand side is finite. Then, the authors of the paper mentioned above showed that there is a very general extrapolation principle that allows one to get the full range of exponents \(p\in(0,\infty)\). As a consequence of this general extrapolation principle, vector-valued inequalities were obtained in a very easy way without using the Banach-valued theory in the paper mentioned above. However, with the extrapolation results mentioned above, no useful estimate was obtained at the endpoint \(p=1\), since \(L^1(w)\) is usually not the suitable endpoint space for the operator \(S\). In the paper mentioned above, for a Calderón-Zygmund operator \(T\), the authors found that the suitable endpoint space is the Lorentz space \(L^{1,\infty}(w)\). However, this is not the case for the commutator, \([b, T]\), generated by \(T\) and \(b\in \text{BMO} (\mathbb R^n)\). The main purpose of the present paper is to provide a more general framework in which this kind of examples can also be treated. To be precise, in this paper, the authors present an extrapolation theory that allows the authors to obtain, from weighted \(L^p\) inequalities on pairs of operators \(S\) for \(p\) fixed and all \(A_\infty\) weights, estimates for the same pairs on very general rearrangement invariant quasi-Banach spaces with \(A_\infty\) weights and also modular inequalities with \(A_\infty\) weights. As above mentioned, vector-valued inequalities are obtained automatically, without the need of a Banach-valued theory. This provides a method to prove very fine estimates for a variety of operators which include singular and fractional integrals and their commutators. In particular, the authors obtain weighted, and vector-valued, extensions of the classical theorems of Boyd and Lorentz-Shimogaki. The key is to develop appropriate versions of Rubio de Francia’s algorithm.