Merge pull request #31 from CLOCTools/v0.8.2

v0.9.0

Author: MYousuf3

CMakeLists.txt

@@ -11,7 +11,7 @@ set(CMAKE_TOOLCHAIN_FILE ${CMAKE_CURRENT_SOURCE_DIR}/vcpkg/scripts/buildsystems/
 ############ PROJECT SPECIFICATIONS
 ################################################################################
 # set project name and version (MAJOR.minor.patch)
-project(ldsCtrlEst VERSION 0.8.1 LANGUAGES CXX C)
+project(ldsCtrlEst VERSION 0.9.0 LANGUAGES CXX C)
 include(CTest)
 # some name variants I will be using:
 set(CMAKE_PROJECT_NAME_CAP LDSCTRLEST)#all caps

include/ldsCtrlEst_h/lds_gaussian_sys.h

@@ -86,19 +86,19 @@ class System : public lds::System {
   }
   /// Set output noise covariance
   void set_R(const Matrix& R) {
-    Reassign(R_,R);
+    Reassign(R_, R);
     do_recurse_Ke_ = true;
   };
 
   /// Set estimator gain
   void set_Ke(const Matrix& Ke) {
-    Reassign(Ke_,Ke);
+    Reassign(Ke_, Ke);
     // if users have set Ke, they must not want to calculate it online.
     do_recurse_Ke_ = false;
   };
   /// Set disturbance estimator gain
   void set_Ke_m(const Matrix& Ke_m) {
-    Reassign(Ke_m_,Ke_m);
+    Reassign(Ke_m_, Ke_m);
     // if users have set Ke, they must not want to calculate it online.
     do_recurse_Ke_ = false;
   };
@@ -113,13 +113,16 @@ class System : public lds::System {
     y_ = cx_ + d_;
   };
 
+  /// System output function: stateless
+  Vector h_(Vector x) override { return C_ * x + d_; };
+
   /// Recursively update estimator gain
   void RecurseKe() override;
 
   // Gaussian-output-specific
-  Matrix R_;           ///< covariance of output noise
+  Matrix R_;              ///< covariance of output noise
   bool do_recurse_Ke_{};  ///< whether to recursively calculate estimator gain
-};                      // System
+};  // System
 }  // namespace gaussian
 }  // namespace lds
 

include/ldsCtrlEst_h/lds_poisson_sys.h

@@ -83,6 +83,9 @@ class System : public lds::System {
     diag_y_.diag() = y_;
   };
 
+  /// System output function: stateless
+  Vector h_(Vector x) override { return exp(C_ * x + d_); };
+
   /**
    * Recursively recalculate estimator gain (Ke).
    *
@@ -100,10 +103,10 @@ class System : public lds::System {
 
  private:
   // Poisson-output-specific
-  Matrix diag_y_;                 ///< diagonal matrix with elements y
+  Matrix diag_y_;  ///< diagonal matrix with elements y
   std::poisson_distribution<size_t>
       pd_;  ///< poisson distribution for simulating data
-};          // System
+};  // System
 }  // namespace poisson
 }  // namespace lds
 

include/ldsCtrlEst_h/lds_sys.h

@@ -30,6 +30,7 @@
 #define LDSCTRLEST_LDS_SYS_H
 
 #include "lds.h"
+#include "lds_uniform_mats.h"
 
 namespace lds {
 /// Linear Dynamical System Type
@@ -94,6 +95,13 @@ class System {
    */
   virtual void h() = 0;
 
+  /**
+   * @brief      system output function (stateless)
+   * @param      x_t  state at time t
+   * @return     predicted state at time t + 1
+   */
+  virtual Vector h_(Vector x) = 0;
+
   /// Get number of inputs
   size_t n_u() const { return n_u_; };
   /// Get number of states
@@ -149,7 +157,7 @@ class System {
   /// Set input matrix
   void set_B(const Matrix& B) { Reassign(B_, B); };
   /// Set process disturbance
-  void set_m(const Vector& m, bool do_force_assign=false) {
+  void set_m(const Vector& m, bool do_force_assign = false) {
     Reassign(m0_, m);
     if ((!do_adapt_m) || do_force_assign) {
       Reassign(m_, m);
@@ -180,6 +188,10 @@ class System {
   /// Reset system variables
   void Reset();
 
+  std::vector<UniformMatrixList<kMatFreeDim2>> nstep_pred_block(
+      UniformMatrixList<kMatFreeDim2> u, UniformMatrixList<kMatFreeDim2> z,
+      size_t n_pred = 1);
+
   /// Print system variables to stdout
   void Print();
 
@@ -220,7 +232,7 @@ class System {
 
   Matrix Ke_;    ///< estimator gain
   Matrix Ke_m_;  ///< estimator gain for process disturbance
-};               // System
+};  // System
 
 }  // namespace lds
 

include/ldsCtrlEst_h/lds_uniform_mats.h

@@ -38,13 +38,15 @@ class UniformMatrixList : public std::vector<Matrix> {
  private:
   // TODO(mfbolus): would rather *uncomment* the below for sake of conversion
   // using std::vector<Matrix>::vector;
+  // don't allow push_back to be used since it doesn't check dims
+  using std::vector<Matrix>::push_back;
+
+ public:
   using std::vector<Matrix>::operator=;
   using std::vector<Matrix>::operator[];
   using std::vector<Matrix>::begin;
   using std::vector<Matrix>::end;
   using std::vector<Matrix>::size;
-
- public:
   using std::vector<Matrix>::at;
   /**
    * @brief      Constructs a new UniformMatrixList.
@@ -139,6 +141,12 @@ class UniformMatrixList : public std::vector<Matrix> {
    * @return     reference to object
    */
   UniformMatrixList<D>& operator=(UniformMatrixList<D>&& that) noexcept;
+  /**
+   * @brief      appends a matrix to the list
+   *
+   * @param      mat   input matrix
+   */
+  void append(const Matrix& mat);
 
  private:
   void CheckDimensions(std::array<size_t, 2> dim);
@@ -167,9 +175,11 @@ inline void UniformMatrixList<D>::Swap(Matrix& that, size_t n) {
     return;
   }
   // if checks pass, perform swap
-  Matrix tmp = std::move((*this)[n]);
-  (*this)[n] = std::move(that);
-  that = std::move(tmp);
+  // not moving, since it causes memory issues.
+  // so this method isn't a memory-saver as designed for now
+  Matrix tmp = (*this)[n];
+  (*this)[n] = that;
+  that = tmp;
 
   if (D == kMatFreeDim1) {
     this->dim_[n][0] = (*this)[n].n_rows;
@@ -179,6 +189,14 @@ inline void UniformMatrixList<D>::Swap(Matrix& that, size_t n) {
   }
 }
 
+template <MatrixListFreeDim D>
+void UniformMatrixList<D>::append(const Matrix& mat) {
+  std::array<size_t, 2> dim({mat.n_rows, mat.n_cols});
+  CheckDimensions(dim);
+  std::vector<Matrix>::push_back(mat);
+  dim_.push_back(dim);
+}
+
 template <MatrixListFreeDim D>
 inline UniformMatrixList<D>& UniformMatrixList<D>::operator=(
     const UniformMatrixList<D>& that) {

python/MANIFEST.in

@@ -0,0 +1,6 @@
+include ldsctrlest/*.so
+include ldsctrlest/*.pyd
+include ldsctrlest/*.lib
+include ldsctrlest/*.pdb
+include ldsctrlest/*.dylib
+include ldsctrlest/*.dll

python/carma

@@ -78,7 +78,6 @@
 
 # compare fit to original without state noise
 sys_hat = glds.System(fit)
-sys_hat.Q = np.zeros_like(sys_hat.Q)
 y_hat, x_hat, _ = sys_hat.simulate_block(u)
 y_imp_hat = sys_hat.simulate_imp(n_samp_imp)
 
@@ -96,7 +95,6 @@
 axs[1].legend([l1[0], l2[0]], ["ground truth", "estimated"])
 axs[1].set(ylabel="Impulse Response (a.u.)", xlabel="Time (s)")
 fig.tight_layout()
-fig
 
 # %%
 # SSID plot var explained

python/examples/eg_glds_fit.py

@@ -53,12 +53,12 @@
 for trial in range(n_trials):
     u.append(np.zeros((n_u, n_samp)))
     for k in range(1, n_samp):
-        u_t = 0.975*u[trial][:,k-1] + 1e-1*np.random.normal(size=n_u)
+        u_t = 0.975 * u[trial][:, k - 1] + 1e-1 * np.random.normal(size=n_u)
         u[trial][:, k] = u_t
 
 y, x, z = sys.simulate_block(u)
 
-n_samp_imp = int(np.ceil(0.1/dt))
+n_samp_imp = int(np.ceil(0.1 / dt))
 y_imp = sys.simulate_imp(n_samp_imp)
 t_imp = np.arange(0, n_samp_imp * dt, dt)
 
@@ -77,7 +77,6 @@
 
 # compare fit to original without state noise
 sys_hat = plds.System(fit)
-sys_hat.Q = np.zeros_like(sys_hat.Q)
 y_hat, x_hat, _ = sys_hat.simulate_block(u)
 y_imp_hat = sys_hat.simulate_imp(n_samp_imp)
 
@@ -87,7 +86,7 @@
 
 fig, axs = plt.subplots(1, 2)
 axs[0].semilogy(sing_vals[:n_h], "-o", color=[0.5, 0.5, 0.5])
-axs[0].semilogy(sing_vals[:n_h], color='k', linewidth=2)
+axs[0].semilogy(sing_vals[:n_h], color="k", linewidth=2)
 axs[0].set(ylabel="Singular Values", xlabel="Singular Value Index")
 
 l1 = axs[1].plot(t_imp, y_imp[0].T, "-", c="k", linewidth=2)
@@ -112,13 +111,21 @@
 axs[0].plot(t, y[eg_trial][0, :] / dt, "k-")
 axs[0].plot(t, y_hat[eg_trial][0, :] / dt, "-", c="gray", linewidth=2)
 axs[0].legend(["measurement", "fit"])
-axs[0].set(ylabel="Output 1 (events/s)", xlabel="Time (s)", title=f"proportion var explained (training): {pve[0]:0.3f}")
+axs[0].set(
+    ylabel="Output 1 (events/s)",
+    xlabel="Time (s)",
+    title=f"proportion var explained (training): {pve[0]:0.3f}",
+)
 
 axs[1].plot(t, y[eg_trial][1, :] / dt, "k-")
 axs[1].plot(t, y_hat[eg_trial][1, :] / dt, "-", c="gray", linewidth=2)
-axs[1].set(ylabel="Output 2 (events/s)", xlabel="Time (s)", title=f"proportion var explained (training): {pve[1]:0.3f}")
+axs[1].set(
+    ylabel="Output 2 (events/s)",
+    xlabel="Time (s)",
+    title=f"proportion var explained (training): {pve[1]:0.3f}",
+)
 
-axs[2].plot(t, u[eg_trial].T, 'k')
+axs[2].plot(t, u[eg_trial].T, "k")
 axs[2].set(ylabel="Input (a.u.)", xlabel="Time (s)")
 
 fig.tight_layout()
@@ -131,18 +138,20 @@
 
 # %%
 # Refit by E-M
-calc_dynamics = True #calculate dynamics (A, B mats)
-calc_Q = True #calculate process noise cov (Q)
-calc_init = True #calculate initial conditions
-calc_output = True #calculate output (C)
-calc_measurement = True #calculate output noise (R)
+calc_dynamics = True  # calculate dynamics (A, B mats)
+calc_Q = True  # calculate process noise cov (Q)
+calc_init = True  # calculate initial conditions
+calc_output = True  # calculate output (C)
+calc_measurement = True  # calculate output noise (R)
 max_iter = 50
 tol = 1e-2
 
 em = plds.FitEM(fit, u_train, z_train)
 
 start = time.perf_counter()
-fit_em = em.Run(calc_dynamics, calc_Q, calc_init, calc_output, calc_measurement, max_iter, tol)
+fit_em = em.Run(
+    calc_dynamics, calc_Q, calc_init, calc_output, calc_measurement, max_iter, tol
+)
 stop = time.perf_counter()
 print(f"Finished EM fit in {(stop-start)*1000} ms.")
 
@@ -162,9 +171,9 @@
 axs[0].legend(["measurement", "EM re-estimated"])
 axs[0].set(ylabel="Output (events/s)")
 
-axs[1].plot(t, z[eg_trial][0, :], 'k')
+axs[1].plot(t, z[eg_trial][0, :], "k")
 
-axs[2].plot(t, u[eg_trial].T, 'k')
+axs[2].plot(t, u[eg_trial].T, "k")
 axs[2].set(ylabel="Input (a.u.)", xlabel="Time (s)")
 
 fig.tight_layout()
@@ -175,7 +184,7 @@
 fig, axs = plt.subplots(1, 2)
 
 l1 = axs[1].plot(t_imp, y_imp[0].T, "-", c="k", linewidth=2)
-l2 = axs[1].plot(t_imp, y_imp_hat_em[0].T, "-", c='gray', linewidth=2)
+l2 = axs[1].plot(t_imp, y_imp_hat_em[0].T, "-", c="gray", linewidth=2)
 axs[1].legend([l1[0], l2[0]], ["ground truth", "EM re-estimated"])
 axs[1].set(ylabel="Impulse Response (a.u.)", xlabel="Time (s)")
 fig.tight_layout()
@@ -192,13 +201,21 @@
 axs[0].plot(t, y[eg_trial][0, :] / dt, "k-")
 axs[0].plot(t, y_hat_em[eg_trial][0, :] / dt, "-", c="gray", linewidth=2)
 axs[0].legend(["measurement", "EM re-estimated"])
-axs[0].set(ylabel="Output 1 (a.u.)", xlabel="Time (s)", title=f"EM-refit proportion var explained (training): {pve_em[0]:0.3f}")
+axs[0].set(
+    ylabel="Output 1 (a.u.)",
+    xlabel="Time (s)",
+    title=f"EM-refit proportion var explained (training): {pve_em[0]:0.3f}",
+)
 
 axs[1].plot(t, y[eg_trial][1, :] / dt, "k-")
 axs[1].plot(t, y_hat_em[eg_trial][1, :] / dt, "-", c="gray", linewidth=2)
-axs[1].set(ylabel="Output 2 (a.u.)", xlabel="Time (s)", title=f"EM-refit proportion var explained (training): {pve_em[1]:0.3f}")
+axs[1].set(
+    ylabel="Output 2 (a.u.)",
+    xlabel="Time (s)",
+    title=f"EM-refit proportion var explained (training): {pve_em[1]:0.3f}",
+)
 
-axs[2].plot(t, u[eg_trial].T, 'k')
+axs[2].plot(t, u[eg_trial].T, "k")
 axs[2].set(ylabel="Input (a.u.)", xlabel="Time (s)")
 
 fig.tight_layout()