DMD.h

 
// Description: Computes the DMD algorithm on the given snapshot matrix. The
//              implemented dynamic mode decomposition algorithm is derived from
//              Tu et. al's paper "On Dynamic Mode Decomposition: Theory and
//              Applications": https://arxiv.org/abs/1312.0041
//              This algorithm also works in the case that the first sample does
//              not start from t = 0.0 by incorporating a time offset.
 
#ifndef included_DMD_h
#define included_DMD_h
 
#include "ParametricDMD.h"
#include <vector>
#include <complex>
 
namespace CAROM {
 
class Matrix;
class Vector;
class ComplexEigenPair;
 
struct DMDInternal
{
    Matrix* snapshots_in;
 
    Matrix* snapshots_out;
 
    Matrix* basis;
 
    Matrix* basis_right;
 
    Matrix* S_inv;
 
    ComplexEigenPair* eigenpair;
};
 
class DMD
{
public:
 
    DMD(int dim, double dt, bool alt_output_basis = false,
        std::shared_ptr<Vector> state_offset = nullptr);
 
    DMD(std::string base_file_name);
 
    virtual ~DMD() { };
 
    virtual void setOffset(std::shared_ptr<Vector> & offset_vector, int order);
 
    virtual void takeSample(double* u_in, double t);
 
    virtual void train(double energy_fraction, const Matrix* W0 = NULL,
                       double linearity_tol = 0.0);
 
    virtual void train(int k, const Matrix* W0 = NULL, double linearity_tol = 0.0);
 
    void projectInitialCondition(const Vector & init, double t_offset = -1.0);
 
    std::unique_ptr<Vector> predict(double t, int deg = 0);
 
    double getTimeOffset() const;
 
    int getNumSamples() const
    {
        return d_snapshots.size();
    }
 
    int getDimension() const
    {
        return d_k;
    }
 
    std::unique_ptr<const Matrix> getSnapshotMatrix();
 
    virtual void load(std::string base_file_name);
 
    void load(const char* base_file_name);
 
    virtual void save(std::string base_file_name);
 
    void save(const char* base_file_name);
 
    void summary(std::string base_file_name);
 
protected:
    friend void getParametricDMD<DMD>(std::unique_ptr<DMD>& parametric_dmd,
                                      const std::vector<Vector>& parameter_points,
                                      std::vector<DMD*>& dmds,
                                      const Vector & desired_point,
                                      std::string rbf,
                                      std::string interp_method,
                                      double closest_rbf_val,
                                      bool reorthogonalize_W);
 
    DMD(int dim, bool alt_output_basis = false,
        std::shared_ptr<Vector> state_offset = nullptr);
 
    DMD(std::vector<std::complex<double>> & eigs,
        std::shared_ptr<Matrix> & phi_real,
        std::shared_ptr<Matrix> & phi_imaginary, int k, double dt,
        double t_offset, std::shared_ptr<Vector> & state_offset);
 
    DMD();
 
    DMD(const DMD& other);
 
    DMD&
    operator = (
        const DMD& rhs);
 
    std::pair<std::shared_ptr<Matrix>,std::shared_ptr<Matrix>> phiMultEigs(
                double t, int deg = 0);
 
    void constructDMD(const Matrix & f_snapshots,
                      int rank,
                      int num_procs,
                      const Matrix* W0,
                      double linearity_tol);
 
    virtual std::pair<Matrix*, Matrix*> computeDMDSnapshotPair(
        const Matrix & snapshots);
 
    virtual void computePhi(DMDInternal dmd_internal_obj);
 
    virtual std::complex<double> computeEigExp(std::complex<double> eig, double t);
 
    virtual void addOffset(Vector & result, double t = 0.0, int deg = 0);
 
    std::unique_ptr<const Matrix>
    createSnapshotMatrix(const std::vector<std::shared_ptr<Vector>> & snapshots);
 
    int d_rank;
 
    int d_num_procs;
 
    int d_dim;
 
    double d_dt = -1.0;
 
    double d_t_offset;
 
    std::vector<std::shared_ptr<Vector>> d_snapshots;
 
    std::vector<double> d_sampled_times;
 
    std::shared_ptr<Vector> d_state_offset;
 
    bool d_trained;
 
    bool d_init_projected;
 
    int d_num_singular_vectors;
 
    std::vector<double> d_sv;
 
    double d_energy_fraction;
 
    int d_k;
 
    std::shared_ptr<Matrix> d_basis;
 
    bool d_alt_output_basis = false;
 
    std::shared_ptr<Matrix> d_A_tilde;
 
    std::shared_ptr<Matrix> d_phi_real;
 
    std::shared_ptr<Matrix> d_phi_imaginary;
 
    std::unique_ptr<Matrix> d_phi_real_squared_inverse;
 
    std::unique_ptr<Matrix> d_phi_imaginary_squared_inverse;
 
    std::shared_ptr<Vector> d_projected_init_real;
 
    std::shared_ptr<Vector> d_projected_init_imaginary;
 
    std::vector<std::complex<double>> d_eigs;
 
};
 
}
 
#endif