Upload relases v0.1.10

See the History page in the documentation for details

.gitignore

@@ -124,6 +124,5 @@ dmypy.json
 # Pyre type checker
 .pyre/
 
-illustrative_examples/tests_safet_mutapcija.py
 illustrative_examples/*
-!illustrative_examples/illustrative_examples.py
+.vscode/*

docs/conf.py

@@ -52,7 +52,7 @@
 project = "pySigmaP"
 copyright = "2020, E. A. Montoya-Araque, A. J. Aparicio-Ortube, D. G. Zapata-Medina, L. G. Arboleda-Monsalve and Universidad Nacional de Colombia"
 author = "E. A. Montoya-Araque\\\A. J. Aparicio-Ortube\\\D. G. Zapata-Medina\\\L. G. Arboleda-Monsalve"
-release = "0.1.9"
+release = "0.1.10"
 
 # -- Options for HTML output -------------------------------------------
 

docs/history.rst

@@ -40,7 +40,7 @@ History
 0.1.7 (2020-12-26)
 ------------------
 
-* Use mstools instead of scikit-learn for linear regression coefficient of determination.
+* Use mstools instead of scikit-learn for linear regression coefficient of determination.
 
 0.1.8 (2021-06-28)
 ------------------
@@ -51,5 +51,21 @@ History
 ------------------
 
 * Update documentation.
-* Include feature to the `Data` class to read data with multiple unloading-reloading stages.
+* Include feature to the ``Data`` class to read data with multiple unloading-reloading stages.
 * Remove dependency of external LaTeX installation.
+
+0.1.10 (2023-10-31)
+-------------------
+
+* Include the ``range2fitCS`` to the ``getSigmaP`` method of the ``Casagrande`` class to limit the stress range for the cubic spline in the calculation of the maximum curvature point.
+* The maximum curvature point is now calculated using a function within the ``casagrande`` module
+  that determines the maximum value of the curvature that is not at the ends, instead
+  of the ``find_peaks`` SciPy's function. However, if the value is not determined, the maximum
+  curvature point is calculated as the absolute maximum value of the curvature.
+* Improve some figures:
+
+    - The void ratio of the :math:`\sigma'_\mathrm{v}` value is now determined by
+      projecting :math:`\sigma'_\mathrm{v}` over a linear interpolation instead of a cubic spline.
+    - :math:`C_\mathrm{c}` and :math:`C_\mathrm{r}` are not calculated inmediatelly when loading the data.
+      Instead, they are calculated when theis respective methods are called.
+    - Other minor improvements.

pyproject.toml

@@ -14,12 +14,12 @@ readme = "README.rst"
 license = {file = "LICENSE"}
 classifiers = ["License :: OSI Approved :: MIT License"]
 dynamic = ["version", "description"]
-requires-python = ">=3.6"
+requires-python = ">=3.7"
 dependencies = [
-    "numpy >=1.19.1",
-    "matplotlib >=3.2.2",
-    "scipy >= 1.5.0",
-    "pandas >= 1.1.1",
+    "numpy >=1.23.5",
+    "matplotlib >=3.7.1",
+    "scipy >= 1.11.3",
+    "pandas >= 1.5.3",
     "mstools >= 0.1.0"
 ]
 

pysigmap/__init__.py

@@ -3,7 +3,7 @@
 """Application software for calculating the preconsolidation pressure from the incremental loading oedometer testing"""
 
 __all__ = ["data", "casagrande", "energy", "bilog", "pachecosilva", "boone"]
-__version__ = "0.1.9"
+__version__ = "0.1.10"
 
 # Global variables
 figsize = [10, 5.3]

pysigmap/bilog.py

@@ -23,7 +23,7 @@
 # -- Required modules
 import numpy as np
 from numpy.polynomial.polynomial import polyfit, polyval
-from scipy.interpolate import CubicSpline
+from scipy.interpolate import CubicSpline, interp1d
 import matplotlib as mpl
 import matplotlib.pyplot as plt
 import matplotlib.ticker as mtick
@@ -168,7 +168,9 @@ def transformY(y, opt=1):
         volLog = transformY(self.data.cleaned["vol"][1:], opt)
         cs = CubicSpline(x=sigmaLog, y=volLog)
         # Specific volume at sigma V
-        volSigmaV = cs(transformX(self.data.sigmaV, opt))
+        # volSigmaV = cs(transformX(self.data.sigmaV, opt))
+        interpolator = interp1d(x=sigmaLog, y=volLog)
+        volSigmaV = interpolator(transformX(self.data.sigmaV, opt))
 
         # -- Compression range (CR)
         self.maskCR = np.full(len(self.data.cleaned), False)

pysigmap/casagrande.py

@@ -75,7 +75,9 @@ def __init__(self, data):
         self.data = data
         return
 
-    def getSigmaP(self, mcp=None, range2fitFOP=None, loglog=True):
+    def getSigmaP(
+        self, mcp=None, range2fitFOP=None, range2fitCS=None, loglog=True
+    ):
         """
         Return the value of the preconsolidation pressure.
 
@@ -105,32 +107,47 @@ def getSigmaP(self, mcp=None, range2fitFOP=None, loglog=True):
 
         def transform(x, reverse=False):
             if reverse:  # Remove a logaritmic scale
-                if loglog:
-                    return 10**10**x
-                else:
-                    return 10**x
+                return 10**10**x if loglog else 10**x
             else:  # Set a logaritmic scale
-                if loglog:
-                    return np.log10(np.log10(x))
-                else:
-                    return np.log10(x)
-
-        sigmaLog = np.log10(self.data.cleaned["stress"][1:])
-        cs = CubicSpline(x=sigmaLog, y=self.data.cleaned["e"][1:])
+                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:  # Using a cubic spline
+        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]
+            # 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
 
-        else:  # Using a fourth order polynomial (FOP)
+        elif mcp is None:  # Using a fourth order polynomial (FOP)
             # -- Indices to fit the FOP
             idxInitFOP = self.data.findStressIdx(
                 stress2find=range2fitFOP[0], cleanedData=True
@@ -160,9 +177,18 @@ def transform(x, reverse=False):
             )
             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]
+            # 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
 
@@ -248,6 +274,24 @@ def transform(x, reverse=False):
                 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(
@@ -268,9 +312,6 @@ def transform(x, reverse=False):
                 )
                 allLayers += l4 + l5
             allLayers += l6
-        else:  # Cubic spline
-            # allLayers = l1 + l2
-            pass
         # Other plots
         l7 = ax1.plot(
             self.data.sigmaV,
@@ -349,6 +390,33 @@ def transform(x, reverse=False):
         return fig
 
 
+def find_max_not_at_ends(vector):
+    """
+    Find the index of the maximum value of a vector that is not at the ends.
+
+    Parameters
+    ----------
+    vector : array_like
+        Vector to find the maximum value.
+
+    Returns
+    -------
+    max_index : int
+        Index of the maximum value of the vector that is not at the ends.
+
+    """
+
+    while True:
+        # Exclude the first and last element
+        vector = vector[1:-1]
+        # Find the index of the maximum value
+        max_index = np.argmax(vector)
+        # If the maximum value is not at the ends of the new vector, break the loop
+        if max_index not in [0, len(vector) - 1]:
+            break
+    return max_index + 1
+
+
 # %%
 """
 2-Clause BSD License.

pysigmap/data.py

@@ -9,7 +9,7 @@
 # -- Required modules
 import numpy as np
 from numpy.polynomial.polynomial import polyfit, polyval
-from scipy.interpolate import CubicSpline
+from scipy.interpolate import CubicSpline, interp1d
 import matplotlib as mpl
 import matplotlib.pyplot as plt
 import matplotlib.ticker as mtick
@@ -109,8 +109,10 @@ def __init__(
         self.preprocessing()
         self.getBreakIndices()
         self.clean()
-        self.compressionIdx()
-        self.recompressionIdx()
+        self.idxCc = None
+        self.idxCr = None
+        # self.compressionIdx()
+        # self.recompressionIdx()
         return
 
     def preprocessing(self):
@@ -210,8 +212,10 @@ def clean(self):
             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)))
+        # 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):
@@ -368,49 +372,52 @@ def plot(self):
             ),
         )
         # Compression index
-        x4Cc = np.linspace(
-            self.cleaned["stress"].iloc[-4], self.cleaned["stress"].iloc[-1]
-        )
-        y4Cc = -self.idxCc * np.log10(x4Cc) + self.idxCcInt
-        ax.plot(
-            x4Cc,
-            y4Cc,
-            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",
+        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=f"Data for linear fit\n(R$^2={self.r2Cc:.3f}$)",
+                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
-        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}$)",
-        )
+        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)