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,
        Vector* state_offset = NULL);
 
    DMD(std::string base_file_name);
 
    virtual ~DMD();
 
    virtual void setOffset(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);
 
    Vector* predict(double t, int deg = 0);
 
    double getTimeOffset() const;
 
    int getNumSamples() const
    {
        return d_snapshots.size();
    }
 
    int getDimension() const
    {
        return d_k;
    }
 
    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>(DMD*& parametric_dmd,
                                      std::vector<Vector*>& parameter_points,
                                      std::vector<DMD*>& dmds,
                                      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, Vector* state_offset = NULL);
 
    DMD(std::vector<std::complex<double>> eigs, Matrix* phi_real,
        Matrix* phi_imaginary, int k,
        double dt, double t_offset, Vector* state_offset);
 
    DMD();
 
    DMD(const DMD& other);
 
    DMD&
    operator = (
        const DMD& rhs);
 
    std::pair<Matrix*, 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);
 
    const Matrix* createSnapshotMatrix(std::vector<Vector*> snapshots);
 
    int d_rank;
 
    int d_num_procs;
 
    int d_dim;
 
    double d_dt = -1.0;
 
    double d_t_offset;
 
    std::vector<Vector*> d_snapshots;
 
    std::vector<Vector*> d_sampled_times;
 
    Vector* d_state_offset = NULL;
 
    bool d_trained;
 
    bool d_init_projected;
 
    int d_num_singular_vectors;
 
    std::vector<double> d_sv;
 
    double d_energy_fraction;
 
    int d_k;
 
    Matrix* d_basis = NULL;
 
    bool d_alt_output_basis = false;
 
    Matrix* d_A_tilde = NULL;
 
    Matrix* d_phi_real = NULL;
 
    Matrix* d_phi_imaginary = NULL;
 
    Matrix* d_phi_real_squared_inverse = NULL;
 
    Matrix* d_phi_imaginary_squared_inverse = NULL;
 
    Vector* d_projected_init_real = NULL;
 
    Vector* d_projected_init_imaginary = NULL;
 
    std::vector<std::complex<double>> d_eigs;
 
};
 
}
 
#endif