Package structure

class Data:
    """``Data`` class.

    When the object is instanced, the compression and recompression indices
    are automatically calculated with the default inputs of the methods.
    See the documentation of each method for more information regarding the
    default inputs.

    Attributes
    ----------
    rawData : pandas DataFrame
        Data from the incremental loading oedometer test. It includes three
        series with the following and strict order: effective vertical stress,
        axial strain and void ratio. The first row must include the on-table
        void ratio of the specimen.
    sigmaV : float
        In situ effective vertical stress of the specimen.
    strainPercent : bool, optional
        Boolean to specify if the axial strain of the row data is in percent or
        not. The default is True.
    reloading : bool, optional
        Boolean to specify if the test included a reloading stage. The default
        is True.
    secondUnloading : bool, optional
        Boolean to specify if the test included two reloading stages. The
        default is True.
    use_this_UR_stage : int, optional
        Index to specify which unloading-reloading stage will be used to
        compute the recompression index when the test data includes more than
        one stage.  Index notation as Python's standard (first stage → 0). The
        default is 0.

    Examples
    --------
    >>> urlCSV = ''.join(['https://raw.githubusercontent.com/eamontoyaa/',
    >>>                   'data4testing/main/pysigmap/testData.csv'])
    >>> urlXLSX = ''.join(['https://raw.githubusercontent.com/eamontoyaa/',
    >>>                    'data4testing/main/pysigmap/testData.xlsx'])
    >>> data = Data(pd.read_csv(urlCSV), sigmaV=75)
    >>> # or
    >>> data = Data(pd.read_excel(urlXLSX), sigmaV=75)
    >>> data.plot()
    >>> data.idxCc, data.idxCr
    (0.23829647651319996, 0.04873212611656529)
    >>> data.compressionIdx(range2fitCc=(1000, 8000))
    >>> data.plot()
    >>> data.idxCc
    0.22754962258271252
    >>> data.recompressionIdx(opt=2)
    >>> data.plot()
    >>> data.idxCr
    0.049481704199255676
    >>> data.recompressionIdx(opt=3)
    >>> data.plot()
    >>> data.idxCr
    0.05398842524225007
    """

    def __init__(
        self,
        rawData,
        sigmaV,
        strainPercent=True,
        reloading=True,
        secondUnloading=True,
        use_this_UR_stage=0,
    ):
        """Initialize the class."""
        self.raw = rawData
        self.sigmaV = sigmaV
        self.strainPercent = strainPercent
        self.reloading = reloading
        self.secondUnloading = secondUnloading
        self.use_this_UR_stage = use_this_UR_stage
        self.preprocessing()
        self.getBreakIndices()
        self.clean()
        self.idxCc = None
        self.idxCr = None
        # self.compressionIdx()
        # self.recompressionIdx()
        return

    def preprocessing(self):
        """
        Rename the series names and create the on-table void ratio attribute.

        Returns
        -------
        None.
        """
        # Standardizing series names
        self.raw.columns = ["stress", "strain", "e"]
        # Removing percentage format to strain values
        if self.strainPercent:
            self.raw["strain"] = self.raw["strain"].divide(100)
        # On-table (initial) void ratio
        self.e_0 = self.raw["e"].iloc[0]
        return

    def getBreakIndices(self):
        """
        Find the break indices of the compressibility curve.

        The break indices are the following:

            - brkIdx1: Start of the first unloading stage
            - brkIdx2: End of the first unloading stage
            - brkIdx3: Point on the compression range after the first reloading
            - brkIdx4: Index of the last point on the compression range

        Returns
        -------
        None.
        """
        idx_unloading_init = [  # brkIdx1: Start of the first unloading
            i
            for i in self.raw.index[1:-1]
            if (
                self.raw["stress"][i] > self.raw["stress"][i - 1]
                and self.raw["stress"][i] > self.raw["stress"][i + 1]
            )
        ]
        if self.secondUnloading:
            idx_unloading_init.pop(-1)
        self.idx_unloading_init = idx_unloading_init
        brkIdx1 = idx_unloading_init[self.use_this_UR_stage]

        if self.reloading:  # Cases 2-3: Reloading and second unloading
            idx_reloading_init = [  # brkIdx2: End of the first unloading
                i
                for i in self.raw.index[1:-1]
                if (
                    self.raw["stress"][i] < self.raw["stress"][i - 1]
                    and self.raw["stress"][i] < self.raw["stress"][i + 1]
                )
            ]
            brkIdx2 = idx_reloading_init[self.use_this_UR_stage]
            self.idx_reloading_init = idx_reloading_init

            # brkIdx3: Points on the NCL after the each reloading
            idx_reloading_ncl = [
                self.raw.query(f"stress == stress[{i}]").index[1]
                for i in idx_unloading_init
            ]
            brkIdx3 = idx_reloading_ncl[self.use_this_UR_stage]
            self.idx_reloading_ncl = idx_reloading_ncl
        else:  # Case 1: No reloading - No second unloading
            brkIdx2 = self.raw.index[-1]  # last point of the unluading
            brkIdx3 = None
        # brkIdx4: index of the last point on the NCL
        brkIdx4 = self.raw.query("stress == stress.max()").index[0]

        self.brkIdx1 = brkIdx1
        self.brkIdx2 = brkIdx2
        self.brkIdx3 = brkIdx3
        self.brkIdx4 = brkIdx4
        return

    def clean(self):
        """
        Clean the raw data removing the unloading and reloading stages.

        Returns
        -------
        None.
        """
        if self.reloading:
            idx_init = [-1] + self.idx_reloading_ncl
            idx_end = self.idx_unloading_init + [self.brkIdx4]
            self.cleaned = pd.concat(
                [
                    self.raw[idx_init[i] + 1 : idx_end[i] + 1]
                    for i in range(len(idx_init))
                ]
            )
        else:
            self.cleaned = self.raw[: self.brkIdx1 + 1]
        self.cleaned.reset_index(drop=True, inplace=True)
        sigmaLog = np.log10(self.cleaned["stress"][1:])
        # cs = CubicSpline(x=sigmaLog, y=self.cleaned["e"][1:])
        # self.eSigmaV = float(cs(np.log10(self.sigmaV)))
        interpolator = interp1d(x=sigmaLog, y=self.cleaned["e"][1:])
        self.eSigmaV = interpolator(np.log10(self.sigmaV))
        return

    def findStressIdx(self, stress2find, cleanedData=True):
        """
        Return the ceiling index of a specified stress.

        Parameters
        ----------
        stress2find : float
            Stress whose ceiling index is wanted.
        cleanedData : bool, optional
            Boolean to indicate if the index must be find in the row or cleaned
            data. The default is True.

        Returns
        -------
        idx : int
            The ceiling index of the stress input.
        """
        if stress2find == 0:
            return 1
        elif stress2find > self.raw["stress"].max():
            return None
        else:
            data4finding = self.cleaned if cleanedData else self.raw
            return data4finding.query(f"stress >= {stress2find}").index[0]

    def compressionIdx(self, range2fitCc=None):
        """
        Calculate the compression index (Cc).

        Parameters
        ----------
        range2fitCc : list, tuple or array (length=2), optional
            Initial and final pressures between which the first-order
            polynomial will be fit to the data on the compression range. If
            None, the Cc index will be calculated as the
            steepest slope of the cubic spline that passes through the data.
            The default is None.

        Returns
        -------
        None.
        """
        self.range2fitCc = range2fitCc
        if range2fitCc is None:  # Using a cubic spline
            sigmaLog = np.log10(self.cleaned["stress"][1:])
            cs = CubicSpline(x=sigmaLog, y=self.cleaned["e"][1:])
            sigmaCS = np.linspace(sigmaLog.iloc[0], sigmaLog.iloc[-1], 500)
            steepestSlopeIdx = np.argmin(cs(sigmaCS, 1))
            idxCc = cs(sigmaCS, 1)[steepestSlopeIdx]
            idxCcInt = (
                cs(sigmaCS[steepestSlopeIdx])
                - idxCc * sigmaCS[steepestSlopeIdx]
            )
            self.fitCc = False
        else:  # Using a linear fit of a first-order polynomial
            idxInitCc = self.findStressIdx(
                stress2find=range2fitCc[0], cleanedData=True
            )
            idxEndCc = self.findStressIdx(
                stress2find=range2fitCc[1], cleanedData=True
            )
            maskCc = np.full(len(self.cleaned), False)
            maskCc[idxInitCc:idxEndCc] = True
            self.maskCc = maskCc
            # -- Linear regresion
            sigmaCc = self.cleaned["stress"].iloc[maskCc]
            sigmaCclog = np.log10(sigmaCc)
            eCc = self.cleaned["e"].iloc[maskCc]
            idxCcInt, idxCc = polyfit(sigmaCclog, eCc, deg=1)
            r2Cc = r2_score(
                y_true=eCc, y_pred=polyval(sigmaCclog, [idxCcInt, idxCc])
            )
            self.r2Cc = r2Cc
            self.fitCc = True
        self.idxCc = abs(idxCc)
        self.idxCcInt = abs(idxCcInt)
        return

    def recompressionIdx(self, opt=1):
        """
        Calculate the recompression index (Cr).

        Parameters
        ----------
        opt : TYPE, optional
            Integer value to indicate which method will be used. Using only
            two points, the start and end of the unluading stage (opt=1); using
            all the points of the first unloading stage (opt=2); using the
            points of the first unloading and reloading stages (opt=3). The
            default is 1.

        Returns
        -------
        None.
        """
        maskCr = np.full(len(self.raw), False)
        if opt == 1:
            maskCr[self.brkIdx1] = True
            maskCr[self.brkIdx2] = True
        elif opt == 2:
            maskCr[self.brkIdx1 : self.brkIdx2 + 1] = True
        elif opt == 3:
            maskCr[self.brkIdx1 : self.brkIdx3 + 1] = True
        # -- Linear regresion
        sigmaCr = self.raw["stress"].iloc[maskCr]
        sigmaCrlog = np.log10(sigmaCr)
        eCr = self.raw["e"].iloc[maskCr]
        idxCrInt, idxCr = polyfit(sigmaCrlog, eCr, deg=1)
        r2Cr = r2_score(
            y_true=eCr, y_pred=polyval(sigmaCrlog, [idxCrInt, idxCr])
        )
        self.maskCr = maskCr
        self.r2Cr = r2Cr
        self.idxCr = abs(idxCr)
        self.idxCrInt = idxCrInt
        return

    def plot(self):
        """
        Plot the compressibility curve, Cc and Cr indices.

        Returns
        -------
        fig : matplotlib.figure.Figure
            Figure that includes the compressibility curve, Cc, Cr and the in
            situ effective vertical stress of the specimen.
        """
        # -- plotting
        fig = plt.figure(figsize=figsize)
        ax = fig.add_axes([0.08, 0.12, 0.55, 0.85])
        ax.plot(
            self.raw["stress"][1:],
            self.raw["e"][1:],
            ls=(0, (1, 1)),
            marker="o",
            lw=1.5,
            c="k",
            mfc="w",
            label="Experimental data",
        )
        ax.plot(
            self.sigmaV,
            self.eSigmaV,
            ls="",
            marker="|",
            c="r",
            ms=15,
            mfc="w",
            mew=1.5,
            label=str().join(
                ["$\sigma^\prime_\mathrm{v_0}=$ ", f"{self.sigmaV:.0f} kPa"]
            ),
        )
        # Compression index
        if self.idxCc is not None:
            x4Cc = np.linspace(
                self.cleaned["stress"].iloc[-2],
                self.cleaned["stress"].iloc[-1],
            )
            y4Cc = -self.idxCc * np.log10(x4Cc) + self.idxCcInt
            ax.axline(
                (x4Cc[0], y4Cc[0]),
                (x4Cc[-1], y4Cc[-1]),
                ls="-",
                lw=1.125,
                color=colors[1],
                label=str().join(["$C_\mathrm{c}=$", f"{self.idxCc:.3f}"]),
            )
            if self.fitCc:
                ax.plot(
                    self.cleaned["stress"].iloc[self.maskCc],
                    self.cleaned["e"].iloc[self.maskCc],
                    ls="",
                    marker="x",
                    color=colors[1],
                    label=f"Data for linear fit\n(R$^2={self.r2Cc:.3f}$)",
                )
        # Recompression index
        if self.idxCr is not None:
            x4Cr = np.linspace(
                self.raw["stress"].iloc[self.maskCr].min(),
                self.raw["stress"].iloc[self.maskCr].max(),
            )
            y4Cr = -self.idxCr * np.log10(x4Cr) + self.idxCrInt
            ax.plot(
                x4Cr,
                y4Cr,
                ls="-",
                lw=1.125,
                color=colors[2],
                label=str().join(["$C_\mathrm{r}=$", f"{self.idxCr:.3f}"]),
            )
            ax.plot(
                self.raw["stress"].iloc[self.maskCr],
                self.raw["e"].iloc[self.maskCr],
                ls="",
                marker="+",
                color=colors[2],
                label=f"Data for linear fit\n(R$^2={self.r2Cr:.3f}$)",
            )
        # other details
        ax.spines["top"].set_visible(False)
        ax.spines["right"].set_visible(False)
        ax.set(
            xscale="log",
            ylabel="Void ratio, $e$",
            xlabel=str().join(
                [
                    "Effective vertical stress, ",
                    "$\sigma^\prime_\mathrm{v}$ [kPa]",
                ]
            ),
        )
        ax.xaxis.set_major_formatter(mtick.ScalarFormatter())
        ax.yaxis.set_minor_locator(mtick.AutoMinorLocator())
        ax.grid(False)
        ax.legend(
            bbox_to_anchor=(1.125, 0.5),
            loc=6,
            title="$\\bf{Compressibility\ curve}$",
        )
        return fig

class Casagrande:
    """``Casagrande`` class.

    When the object is instanced, the method ``getSigmaP()`` calculates the
    preconsolidation pressure by the method proposed by Casagrande (1963)
    based on the compression index obtained with the ``Data`` class and the
    parameters of the method. See the method documentation for more
    information.

    Attributes
    ----------
    data : Object instanced from the ``Data`` class.
        Contains the data structure from the oedometer test. See the class
        documentation for more information.

    Examples
    --------
    >>> urlCSV = ''.join(['https://raw.githubusercontent.com/eamontoyaa/',
    >>>                   'data4testing/main/pysigmap/testData.csv'])
    >>> data = Data(pd.read_csv(urlCSV), sigmaV=75)
    >>> method = Casagrande(data)
    >>> method.getSigmaP(range2fitFOP=None, loglog=True)
    >>> method.sigmaP, method.ocr
    (925.6233444375316, 12.34164459250042)
    >>> method.getSigmaP(range2fitFOP=[20, 5000], loglog=True)
    >>> method.sigmaP, method.ocr
    (651.5675456910274, 8.687567275880365)
    >>> method.getSigmaP(mcp=200)
    >>> method.sigmaP, method.ocr
    (498.6050168481866, 6.6480668913091545)
    """

    def __init__(self, data):
        """Initialize the class."""
        self.data = data
        return

    def getSigmaP(
        self, mcp=None, range2fitFOP=None, range2fitCS=None, loglog=True
    ):
        """
        Return the value of the preconsolidation pressure.

        Parameters
        ----------
        mcp : float, optional
            Variable to manually specify the maximun curvature point. If a
            value is input, the other parameters are not taken into account.
            The default is None.
        range2fitFOP : list, tuple or array (length=2), optional
            Initial and final pressures between which the fourth-order
            polynomial (FOP) will be fitted to the compressibility curve to
            calculate the maximum curvature point (MCP). If None, the MCP will
            be obtained using a cubic spline that passes through the data. The
            default is None.
        loglog : bool, optional
            Boolean to specify if the vertical effective stress will be
            transformed applying a logarithm twice to fit the FOT. The default
            is True.

        Returns
        -------
        fig : matplotlib figure
            Figure with the development of the method and the results.

        """

        def transform(x, reverse=False):
            if reverse:  # Remove a logaritmic scale
                return 10**10**x if loglog else 10**x
            else:  # Set a logaritmic scale
                return np.log10(np.log10(x)) if loglog else np.log10(x)

        # -- Indices to fit the Cubis Spline (CS)
        if range2fitCS is None:
            idxInitCS, idxEndCS = 1, len(self.data.cleaned["stress"]) + 1
        else:
            idxInitCS = self.data.findStressIdx(
                stress2find=range2fitCS[0], cleanedData=True
            )
            idxEndCS = self.data.findStressIdx(
                stress2find=range2fitCS[1], cleanedData=True
            )
        sigmaLog = np.log10(self.data.cleaned["stress"][idxInitCS:idxEndCS])
        cs = CubicSpline(
            x=sigmaLog, y=self.data.cleaned["e"][idxInitCS:idxEndCS]
        )
        sigmaCS = np.linspace(sigmaLog.iloc[0], sigmaLog.iloc[-1], 100)
        if range2fitFOP is None and mcp is None:  # Using a cubic spline
            x4FOP = 10**sigmaCS
            # -- Curvature function k(x) = f''(x)/(1+(f'(x))²)³/²
            curvature = abs(cs(sigmaCS, 2)) / (1 + cs(sigmaCS, 1) ** 2) ** (
                3 / 2
            )
            # maxCurvIdx = find_peaks(
            #     curvature, distance=self.data.cleaned["stress"].max()
            # )[0][0]
            try:
                maxCurvIdx = find_max_not_at_ends(curvature)
            except Exception:
                print(
                    "Maximun curvature not found matematicaly.",
                    "Choosing the absolute maximum value.",
                    sep="\n",
                )
                maxCurvIdx = np.argmax(curvature)
            self.sigmaMC = 10 ** sigmaCS[maxCurvIdx]  # Max. Curvature point
            self.eMC = cs(sigmaCS[maxCurvIdx])  # Void ratio at MC

        elif mcp is None:  # Using a fourth order polynomial (FOP)
            # -- Indices to fit the FOP
            idxInitFOP = self.data.findStressIdx(
                stress2find=range2fitFOP[0], cleanedData=True
            )
            idxEndFOP = self.data.findStressIdx(
                stress2find=range2fitFOP[1], cleanedData=True
            )

            # -- fittig a polynomial to data without unloads
            sigmaFOP = self.data.cleaned["stress"][idxInitFOP:idxEndFOP]
            sigmaFOPlog = transform(sigmaFOP)
            eFOP = self.data.cleaned["e"][idxInitFOP:idxEndFOP]
            p0, p1, p2, p3, p4 = polyfit(sigmaFOPlog, eFOP, deg=4)
            r2FOP = r2_score(
                y_true=eFOP, y_pred=polyval(sigmaFOPlog, [p0, p1, p2, p3, p4])
            )
            x4FOP = np.linspace(sigmaFOP.iloc[0], sigmaFOP.iloc[-1], 1000)
            x4FOPlog = transform(x4FOP)
            y4FOP = polyval(x4FOPlog, [p0, p1, p2, p3, p4])

            # -- Curvature function k(x) = f''(x)/(1+(f'(X))²)³/²
            firstDer = (
                p1
                + 2 * p2 * x4FOPlog
                + 3 * p3 * x4FOPlog**2
                + 4 * p4 * x4FOPlog**3
            )
            secondDer = 2 * p2 + 6 * p3 * x4FOPlog + 12 * p4 * x4FOPlog**2
            curvature = abs(secondDer) / (1 + firstDer**2) ** 1.5
            # maxCurvIdx = find_peaks(
            #     curvature, distance=self.data.cleaned["stress"].max()
            # )[0][0]
            try:
                maxCurvIdx = find_max_not_at_ends(curvature)
            except Exception:
                print(
                    "Maximun curvature not found matematicaly.",
                    "Choosing the absolute maximum value.",
                    sep="\n",
                )
                maxCurvIdx = np.argmax(curvature)
            self.sigmaMC = transform(x4FOPlog[maxCurvIdx], True)  # Max. Curv.
            self.eMC = y4FOP[maxCurvIdx]  # Void ratio at max. curvature

        if mcp is not None:
            self.sigmaMC = mcp  # Max. Curvature point
            self.eMC = cs(np.log10(self.sigmaMC))  # Void ratio at MC

        # -- Bisector line
        slopeMP = cs(np.log10(self.sigmaMC), nu=1)  # Slope at MC
        y1, x1 = self.eMC, np.log10(self.sigmaMC)
        x2 = np.log10(
            np.linspace(self.sigmaMC, self.data.cleaned["stress"].iloc[-1], 50)
        )
        y2 = slopeMP * (x2 - x1) + y1
        slopeBisect = np.tan(
            0.5 * np.arctan(slopeMP)
        )  # slope of bisector line
        y2bis = slopeBisect * (x2 - x1) + y1

        # -- Preconsolidation pressure
        self.sigmaP = 10 ** (
            (y1 - slopeBisect * x1 - self.data.idxCcInt)
            / (-self.data.idxCc - slopeBisect)
        )
        self.eSigmaP = slopeBisect * (np.log10(self.sigmaP) - x1) + y1
        self.ocr = self.sigmaP / self.data.sigmaV

        # -- plotting
        fig = plt.figure(figsize=figsize)
        ax1 = fig.add_axes([0.08, 0.12, 0.55, 0.85])
        l1 = ax1.plot(
            self.data.raw["stress"][1:],
            self.data.raw["e"][1:],
            ls=(0, (1, 1)),
            marker="o",
            lw=1.5,
            c="k",
            mfc="w",
            label="Compressibility curve",
        )
        # Lines of the Casagrande's method
        ax1.plot(
            [self.sigmaMC, self.data.cleaned["stress"].iloc[-1]],
            [self.eMC, self.eMC],
            ls="--",
            lw=1.125,
            c="k",
        )  # hztl line
        ax1.plot(10**x2, y2, ls="--", lw=1.125, color="k")  # tangent line
        ax1.plot(
            10**x2, y2bis, ls="--", lw=1.125, color="k"
        )  # bisector line
        # Compression index (Cc)
        x4Cc = np.linspace(self.sigmaP, self.data.cleaned["stress"].iloc[-1])
        y4Cc = -self.data.idxCc * np.log10(x4Cc) + self.data.idxCcInt
        l2 = ax1.plot(
            x4Cc,
            y4Cc,
            ls="-",
            lw=1.5,
            color=colors[1],
            label=str().join(["$C_\mathrm{c}=$", f"{self.data.idxCc:.3f}"]),
        )
        allLayers = l1 + l2
        if self.data.fitCc:
            l3 = ax1.plot(
                self.data.cleaned["stress"].iloc[self.data.maskCc],
                self.data.cleaned["e"].iloc[self.data.maskCc],
                ls="",
                marker="x",
                color=colors[1],
                label=f"Data for linear fit\n(R$^2={self.data.r2Cc:.3f}$)",
            )
            allLayers.insert(2, l3[0])
        if mcp is None:  # Curvature
            ax2 = ax1.twinx()  # second y axis for curvature function
            l6 = ax2.plot(
                x4FOP,
                curvature,
                ls="--",
                c=colors[0],
                lw=1.125,
                mfc="w",
                label="Curvature",
            )
            if range2fitCS is not None:  # Cubic spline
                l4 = ax1.plot(
                    10**sigmaCS,
                    cs(sigmaCS),
                    ls="--",
                    lw=1.125,
                    color=colors[2],
                    label="Cubic spline",
                )
                l5 = ax1.plot(
                    self.data.cleaned["stress"][idxInitCS:idxEndCS],
                    self.data.cleaned["e"][idxInitCS:idxEndCS],
                    ls="",
                    marker="+",
                    c=colors[2],
                    label="Data for cubic spline",
                )
                allLayers += l4 + l5
            # if range2fitFOP is None:  # Curvature
            if range2fitFOP is not None:  # FOP fit
                l4 = ax1.plot(
                    x4FOP,
                    y4FOP,
                    ls="--",
                    lw=1.125,
                    color=colors[2],
                    label="$4^\mathrm{th}$-order polynomial",
                )
                l5 = ax1.plot(
                    sigmaFOP,
                    eFOP,
                    ls="",
                    marker="+",
                    c=colors[2],
                    label=f"Data for linear fit\n(R$^2={r2FOP:.3f}$)",
                )
                allLayers += l4 + l5
            allLayers += l6
        # Other plots
        l7 = ax1.plot(
            self.data.sigmaV,
            self.data.eSigmaV,
            ls="",
            marker="|",
            c="r",
            ms=15,
            mfc="w",
            mew=1.5,
            label=str().join(
                [
                    "$\sigma^\prime_\mathrm{v_0}=$ ",
                    f"{self.data.sigmaV:.0f} kPa",
                ]
            ),
        )
        l8 = ax1.plot(
            self.sigmaMC,
            self.eMC,
            ls="",
            marker="^",
            c=colors[0],
            mfc="w",
            mew=1.5,
            ms=7,
            label="Max. curvature point",
        )
        l9 = ax1.plot(
            self.sigmaP,
            self.eSigmaP,
            ls="",
            marker="o",
            mfc="w",
            c=colors[0],
            ms=7,
            mew=1.5,
            label=str().join(
                [
                    "$\sigma^\prime_\mathrm{p}=$ ",
                    f"{self.sigmaP:.0f} kPa\n",
                    f"OCR= {self.ocr:.1f}",
                ]
            ),
        )
        allLayers += l7 + l8 + l9
        # Legend
        labs = [layer.get_label() for layer in allLayers]
        ax1.legend(
            allLayers,
            labs,
            bbox_to_anchor=(1.125, 0.5),
            loc=6,
            title="$\\bf{Casagrande\ method}$",
        )
        # Other details
        ax1.spines["top"].set_visible(False)
        ax1.spines["right"].set_visible(False)
        ax1.set(
            xscale="log",
            ylabel="Void ratio, $e$",
            xlabel=str().join(
                [
                    "Effective vertical stress, ",
                    "$\sigma^\prime_\mathrm{v}$ [kPa]",
                ]
            ),
        )
        ax1.xaxis.set_major_formatter(mtick.ScalarFormatter())
        ax1.yaxis.set_minor_locator(mtick.AutoMinorLocator())
        ax1.grid(False)
        if mcp is None:  # Curvature
            ax2.yaxis.set_minor_locator(mtick.AutoMinorLocator())
            ax2.set(ylabel="Curvature, $k$")
            ax2.spines["top"].set_visible(False)
        return fig