feat: Introduce paper and figure generation script for PSID simulations, alongside minor test and documentation improvements.

Author: OmidS

source/PSID/tests/test_IPSID.py

@@ -158,6 +158,7 @@ def test_IPSID(self):
                 allZp2, allYp2, allXp2 = sId3_cp.predict(Y, U=U)
 
                 try:
+                    # Verify identification errors are within tolerance
                     np.testing.assert_array_less(
                         [err[p] for p in params],
                         5e-2,
@@ -179,6 +180,7 @@ def test_IPSID(self):
                         err_msg="Error too large for some params {}".format(params),
                     )
 
+                    # Verify decoding accuracy
                     np.testing.assert_allclose(allXp2, allXp, rtol=1e-3)
                     np.testing.assert_allclose(allZp2, allZp - ZMean, atol=1e-6)
                     np.testing.assert_allclose(allYp2, allYp - YMean, atol=1e-6)
@@ -263,6 +265,7 @@ def test_IPSID_enforce_stability(self):
                     else:
                         print("Ok")
                     print("Done")
+                    # Verify eigenvalues are within unit circle
                     np.testing.assert_array_less(eigsAbs, 1)
 
         if failInds:

source/PSID/tests/test_LSSM.py

@@ -81,6 +81,8 @@ def test_generateRealizationWithKF(self):
             plt.legend()
             """
             mse.append(np.mean((Y - Yp) ** 2))
+
+            # Verify stochastic and predictor realizations match
             np.testing.assert_almost_equal(Y, Yp, decimal=10)
             # np.testing.assert_equal(Y, Yp)
             # np.testing.assert_allclose(Y, Yp, rtol=1e-3)
@@ -105,6 +107,8 @@ def test_solveric(self):
         out = (np.array([[0.6479, 0.0861], [0.6025, -0.1736]]), True)
 
         P, has_solution = solveric(**inp)
+
+        # Verify solveric result against known solution
         np.testing.assert_equal(out[1], has_solution)
         np.testing.assert_allclose(out[0], P, rtol=1e-3)
 
@@ -128,6 +132,7 @@ def test_solve_discrete_are_iterative(self):
                     s.A.T, s.C.T, s.Q, s.R, s=s.S, return_log=True
                 )
 
+                # Verify iterative solver matches scipy implementation
                 np.testing.assert_almost_equal(Pp1, Pp2)
                 np.testing.assert_almost_equal(0, PpNormChanges[-1])
 
@@ -153,6 +158,7 @@ def test_LSSM_randomize_in_predictor_form(self):
                 etcBlock = np.copy(s.A_KC)
                 etcBlock[:n1, n1:] = 0
                 try:
+                    # Verify predictor form structure
                     np.testing.assert_almost_equal(trBlock, np.zeros_like(trBlock))
                     np.testing.assert_allclose(
                         (s.A_KC - etcBlock) / np.linalg.norm(etcBlock),
@@ -313,6 +319,7 @@ def test_kalmanFilterAndSmoothingForS0Nu0(self):
                     filtered_state_means_2, filtered_state_covariances_2 = kf.filter(Y)
                     smoothed_state_means_2, smoothed_state_covariances_2 = kf.smooth(Y)
 
+                    # Verify against pykalman
                     np.testing.assert_allclose(
                         filtered_state_means, filtered_state_means_2, rtol=1e-3
                     )
@@ -511,6 +518,7 @@ def test_forwardBackwardSmootherBeingTheSameAsRTS(self):
                 # """
 
                 try:
+                    # Verify Forward-Backward smoother matches RTS
                     np.testing.assert_allclose(allXp, allXpFB, rtol=1e-4)
                     np.testing.assert_allclose(allXs, allXsFB, rtol=1e-4)
                     np.testing.assert_allclose(allPs, allPsFB, rtol=1e-4)

source/PSID/tests/test_PSID.py

@@ -120,6 +120,7 @@ def test_PSID(self):
                 allZp2, allYp2, allXp2 = sId3_cp.predict(Y, U=None)
 
                 try:
+                    # Verify identification errors are within tolerance
                     np.testing.assert_array_less(
                         [err[p] for p in params],
                         5e-2,
@@ -141,6 +142,7 @@ def test_PSID(self):
                         err_msg="Error too large for some params {}".format(params),
                     )
 
+                    # Verify decoding accuracy
                     np.testing.assert_allclose(allXp2, allXp, rtol=1e-3)
                     np.testing.assert_allclose(allZp2, allZp - ZMean, atol=1e-6)
                     np.testing.assert_allclose(allYp2, allYp - YMean, atol=1e-6)
@@ -192,6 +194,7 @@ def test_PSID_doest_change_inputs(self):
             UBU = copy.deepcopy(U)
 
             sId = SubspacePSID(Y, Z, nx=nx, n1=n1, i=horizon, time_first=True)
+            # Verify input data remains unchanged
             np.testing.assert_equal(Y, YBU)
             np.testing.assert_equal(Z, ZBU)
             np.testing.assert_equal(U, UBU)

source/PSID/tests/test_PSIDModel.py

@@ -130,6 +130,7 @@ def test_PSID(self):
                 err2.update(computeLSSMIdError(s, sId2.model))
 
                 try:
+                    # Verify identification errors are within tolerance
                     np.testing.assert_array_less(
                         [err[p] for p in params],
                         5e-2,
@@ -151,6 +152,7 @@ def test_PSID(self):
                         err_msg="Error too large for some params {}".format(params),
                     )
 
+                    # Verify decoding accuracy
                     np.testing.assert_allclose(allXp2, allXp, rtol=1e-3)
                     np.testing.assert_allclose(allZp2, allZp - ZMean, atol=1e-6)
                     np.testing.assert_allclose(allYp2, allYp - YMean, atol=1e-6)

source/PSID/tests/test_PrepModel.py

@@ -54,6 +54,7 @@ def test_preprocessing(self):
                         sm = PrepModel()
                         sm.fit(data, time_first=time_first, **args)
 
+                        # Verify input data remains unchanged
                         np.testing.assert_equal(data, dataCopy)
 
                         ddof = args["std_ddof"] if "std_ddof" in args else 1
@@ -66,16 +67,19 @@ def test_preprocessing(self):
                             newDataStd = np.std(newData, axis=1, ddof=ddof)
 
                         if args["zscore"] or args["remove_mean"]:
+                            # Verify mean is zero
                             np.testing.assert_almost_equal(
                                 newDataMean, np.zeros_like(newDataMean)
                             )
                         if args["zscore"]:
+                            # Verify std is one
                             np.testing.assert_almost_equal(
                                 newDataStd, np.ones_like(newDataStd)
                             )
 
                         recoveredData = sm.apply_inverse(newData, time_first=time_first)
 
+                        # Verify inverse transform recovers original data
                         np.testing.assert_almost_equal(recoveredData, dataCopy)
 
     def test_preprocessing_for_segmented_data(self):
@@ -133,17 +137,20 @@ def test_preprocessing_for_segmented_data(self):
                             newDataStd = np.std(newDataCat, axis=1, ddof=ddof)
 
                         if args["zscore"] or args["remove_mean"]:
+                            # Verify mean is zero
                             np.testing.assert_almost_equal(
                                 newDataMean, np.zeros_like(newDataMean)
                             )
                         if args["zscore"]:
+                            # Verify std is one
                             np.testing.assert_almost_equal(
                                 newDataStd, np.ones_like(newDataStd)
                             )
 
                         recoveredData = sm.apply_inverse(newData, time_first)
 
                         for t in range(len(dataCopy)):
+                            # Verify recovery of segmented data
                             np.testing.assert_almost_equal(
                                 recoveredData[t], dataCopy[t]
                             )