Check if correlation matrix is consistent

src/semeio/fmudesign/create_design.py

@@ -17,6 +17,37 @@
 from semeio.fmudesign.iman_conover import ImanConover
 
 
+def is_consistent_correlation_matrix(matrix):
+    """
+    Check if a matrix is a consistent correlation matrix.
+
+    A correlation matrix is consistent if it has:
+    1. All diagonal elements equal to 1
+    2. Is positive semidefinite (all eigenvalues ≥ 0)
+
+    Args:
+        matrix: numpy array representing the correlation matrix
+
+    Returns:
+        bool: True if matrix is a consistent correlation matrix, False otherwise
+    """
+    # Check if diagonal elements are 1
+    if not np.allclose(np.diagonal(matrix), 1):
+        return False
+
+    # Check positive semidefiniteness using eigenvalues
+    try:
+        eigenvals = np.linalg.eigvals(matrix)
+        # Matrix is positive semidefinite if all eigenvalues are non-negative
+        # Using small tolerance to account for numerical errors
+        if not np.all(eigenvals > -1e-8):
+            return False
+    except np.linalg.LinAlgError:
+        return False
+
+    return True
+
+
 def nearest_correlation_matrix(matrix, *, weights=None, eps=1e-6, verbose=False):
     """Returns the correlation matrix nearest to `matrix`, weighted elementwise
     by `weights`.
@@ -822,9 +853,29 @@ def generate(self, realnums, parameters, seedvalues, corrdict, rng):
                     # Make correlation matrix symmetric by copying lower triangular part
                     correlations = np.triu(correlations.T, k=1) + np.tril(correlations)
 
-                    correlations = nearest_correlation_matrix(
-                        correlations, weights=None, eps=1e-6, verbose=False
-                    )
+                    if not is_consistent_correlation_matrix(correlations):
+                        print("\nWarning: Correlation matrix is not consistent")
+                        print("Requirements:")
+                        print("  - Ones on the diagonal")
+                        print("  - Positive semi-definite matrix")
+                        print("\nInput correlation matrix:")
+                        with np.printoptions(
+                            precision=2,
+                            suppress=True,
+                            formatter={"float": lambda x: f"{x:.2f}"},
+                        ):
+                            print(correlations)
+                        correlations = nearest_correlation_matrix(
+                            correlations, weights=None, eps=1e-6, verbose=False
+                        )
+                        print("\nAdjusted to nearest consistent correlation matrix:")
+                        with np.printoptions(
+                            precision=2,
+                            suppress=True,
+                            formatter={"float": lambda x: f"{x:.2f}"},
+                        ):
+                            print(correlations)
+                        print()
 
                     sampler = qmc.LatinHypercube(
                         d=len(multivariate_parameters), seed=rng

tests/fmudesign/test_designmatrix.py

@@ -68,8 +68,44 @@ def test_endpoint(tmpdir, monkeypatch):
     shutil.copy(str(designfile), ".")
     shutil.copy(Path(designfile).parent / dependency, ".")
 
-    subprocess.run(["fmudesign", str(designfile)], check=True)
+    result = subprocess.run(
+        ["fmudesign", str(designfile)], check=True, capture_output=True, text=True
+    )
 
+    expected_output = """Added sensitivity : seed
+    Added sensitivity : faults
+    Added sensitivity : velmodel
+    Added sensitivity : contacts
+    Added sensitivity : multz
+    Added sensitivity : sens6
+
+    Warning: Correlation matrix is not consistent
+    Requirements:
+    - Ones on the diagonal
+    - Positive semi-definite matrix
+
+    Input correlation matrix:
+    [[1.00 0.90 0.00 0.00]
+    [0.90 1.00 0.90 0.00]
+    [0.00 0.90 1.00 0.00]
+    [0.00 0.00 0.00 1.00]]
+
+    Adjusted to nearest consistent correlation matrix:
+    [[1.00 0.74 0.11 0.00]
+    [0.74 1.00 0.74 0.00]
+    [0.11 0.74 1.00 0.00]
+    [0.00 0.00 0.00 1.00]]
+
+    Added sensitivity : sens7
+    Added sensitivity : sens8
+    Provided number of background values (11) is smaller than number of realisations for sensitivity ('sens7', 'p10_p90') and parameter PARAM13. Will be filled with default values.
+    Provided number of background values (11) is smaller than number of realisations for sensitivity ('sens7', 'p10_p90') and parameter PARAM14. Will be filled with default values.
+    Provided number of background values (11) is smaller than number of realisations for sensitivity ('sens7', 'p10_p90') and parameter PARAM15. Will be filled with default values.
+    Provided number of background values (11) is smaller than number of realisations for sensitivity ('sens7', 'p10_p90') and parameter PARAM16. Will be filled with default values.
+    A total of 91 realizations were generated
+    Designmatrix written to generateddesignmatrix.xlsx"""
+
+    assert result.stdout.split() == expected_output.split()
     assert Path("generateddesignmatrix.xlsx").exists  # Default output file
     valid_designmatrix(pd.read_excel("generateddesignmatrix.xlsx", engine="openpyxl"))