From 96a57eade7e885d49dcc874cd706808413c4dbbf Mon Sep 17 00:00:00 2001 From: "Siu Wun \"Tony\" Cheung" Date: Sun, 12 Nov 2023 14:03:37 -0800 Subject: [PATCH] DMDc implementation (#238) * Add MATLAB version * Update MATLAB script * Initialize DMDc scripts * Add heat conduction DMDC script * astyle * Minor updates * Astyle * Minor fix * Updates * Minor fix * Minor updates * Astyle * Minor Fixes * Fix some bugs * Astyle * Add examples * Astyle * Default empty mesh * Add reference map for general meshes * Remove mesh file * Astyle * Minor changes * Astyle * Add comments * Replace reduced dimension by energy fraction in MATLAB * Add delete * Fix time dependence * Astyle * Fix comment * Astyle * Update example results --- CMakeLists.txt | 2 + examples/dmd/heat_conduction_dmdc.cpp | 763 ++++++++++++++++++ lib/CMakeLists.txt | 1 + lib/algo/DMD.cpp | 1 + lib/algo/DMD.h | 8 +- lib/algo/DMDc.cpp | 1053 +++++++++++++++++++++++++ lib/algo/DMDc.h | 408 ++++++++++ unit_tests/matlab/dmdc.m | 52 ++ 8 files changed, 2286 insertions(+), 2 deletions(-) create mode 100644 examples/dmd/heat_conduction_dmdc.cpp create mode 100644 lib/algo/DMDc.cpp create mode 100644 lib/algo/DMDc.h create mode 100644 unit_tests/matlab/dmdc.m diff --git a/CMakeLists.txt b/CMakeLists.txt index fe262b982..af9036360 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -140,6 +140,7 @@ if (ENABLE_EXAMPLES) dg_advection nonlinear_elasticity heat_conduction + heat_conduction_dmdc parametric_heat_conduction de_parametric_heat_conduction_greedy wave_equation @@ -166,6 +167,7 @@ if (ENABLE_EXAMPLES) dmd dmd dmd + dmd dmd) list(LENGTH examples len1) diff --git a/examples/dmd/heat_conduction_dmdc.cpp b/examples/dmd/heat_conduction_dmdc.cpp new file mode 100644 index 000000000..7d76568c2 --- /dev/null +++ b/examples/dmd/heat_conduction_dmdc.cpp @@ -0,0 +1,763 @@ +/****************************************************************************** + * + * Copyright (c) 2013-2023, Lawrence Livermore National Security, LLC + * and other libROM project developers. See the top-level COPYRIGHT + * file for details. + * + * SPDX-License-Identifier: (Apache-2.0 OR MIT) + * + *****************************************************************************/ + +// libROM MFEM Example: DMDc for Heat_Conduction (adapted from ex16p.cpp) +// +// Compile with: make heat_conduction_dmdc +// +// ================================================================================= +// +// Sample runs and results for DMDc: +// +// Command 1: +// mpirun -np 8 heat_conduction_dmdc -s 1 -a 0.0 -k 1.0 -visit +// +// Output 1: +// Relative error of DMDc temperature (u) at t_final: 0.5 is 0.0021705658 +// +// Command 2: +// mpirun -np 8 heat_conduction_dmdc -s 1 -a 0.5 -k 0.5 -o 4 -tf 0.7 -vs 1 -visit +// +// Output 2: +// Relative error of DMDc temperature (u) at t_final: 0.7 is 0.00099736216 +// +// ================================================================================= +// +// Description: This example solves a time dependent nonlinear heat equation +// problem of the form du/dt = C(u) + s, with a non-linear diffusion +// operator C(u) = \nabla \cdot (\kappa + \alpha u) \nabla u +// and time-varying external inlet and outlet source. +// The inlet and the outlet is located at (0,0) and (0.5,0.5) +// in the reference domain [-1,1]^d, where the shut down time and +// the amplitude of the sources are the control variables. +// +// The example demonstrates the use of nonlinear operators (the +// class ConductionOperator defining C(u)), as well as their +// implicit time integration. Note that implementing the method +// ConductionOperator::ImplicitSolve is the only requirement for +// high-order implicit (SDIRK) time integration. Optional saving +// with ADIOS2 (adios2.readthedocs.io) is also illustrated. + +#include "mfem.hpp" +#include "algo/DMDc.h" +#include "linalg/Vector.h" +#include +#include +#include + +using namespace std; +using namespace mfem; + +/** After spatial discretization, the conduction model can be written as: + * + * du/dt = M^{-1}(-Ku) + * + * where u is the vector representing the temperature, M is the mass matrix, + * and K is the diffusion operator with diffusivity depending on u: + * (\kappa + \alpha u). + * + * Class ConductionOperator represents the right-hand side of the above ODE. + */ +class ConductionOperator : public TimeDependentOperator +{ +protected: + ParFiniteElementSpace &fespace; + Array ess_tdof_list; // this list remains empty for pure Neumann b.c. + + ParBilinearForm *M; + ParBilinearForm *K; + ParLinearForm *B; + + HypreParMatrix Mmat; + HypreParMatrix Kmat; + HypreParMatrix *T; // T = M + dt K + double current_dt; + + CGSolver M_solver; // Krylov solver for inverting the mass matrix M + HypreSmoother M_prec; // Preconditioner for the mass matrix M + + CGSolver T_solver; // Implicit solver for T = M + dt K + HypreSmoother T_prec; // Preconditioner for the implicit solver + + FunctionCoefficient &src_func; // Source function coefficient + double alpha, kappa; // Nonlinear parameters + + mutable Vector z; // auxiliary vector + HypreParVector *b; // source vector + +public: + ConductionOperator(ParFiniteElementSpace &f, FunctionCoefficient &s, + double alpha, double kappa, const Vector &u); + + virtual void Mult(const Vector &u, Vector &du_dt) const; + /** Solve the Backward-Euler equation: k = f(u + dt*k, t), for the unknown k. + This is the only requirement for high-order SDIRK implicit integration.*/ + virtual void ImplicitSolve(const double dt, const Vector &u, Vector &k); + + /// Update the diffusion BilinearForm K using the given true-dof vector `u`. + void SetSourceTime(const double t); + + /// Update the diffusion BilinearForm K using the given true-dof vector `u`. + void SetParameters(const Vector &u); + + virtual ~ConductionOperator(); +}; + +double InitialTemperature(const Vector &x); +double TimeWindowFunction(const double t, const double t_begin, + const double t_end); +double Amplitude(const double t, const int index); +double SourceFunction(const Vector &x, double t); + +Vector bb_min, bb_max; // Mesh bounding box +double amp_in = 0.2; +double t_end_in = 0.1; +double amp_out = 0.1; +double t_end_out = 0.3; +double dt = 1.0e-2; + +int main(int argc, char *argv[]) +{ + // 1. Initialize MPI. + int num_procs, myid; + MPI_Init(&argc, &argv); + MPI_Comm_size(MPI_COMM_WORLD, &num_procs); + MPI_Comm_rank(MPI_COMM_WORLD, &myid); + + // 2. Parse command-line options. + const char *mesh_file = ""; + int ser_ref_levels = 2; + int par_ref_levels = 1; + int order = 2; + int ode_solver_type = 1; + double t_final = 0.5; + double alpha = 1.0e-2; + double kappa = 0.5; + double ef = 0.9999; + int rdim = -1; + bool visualization = true; + bool visit = false; + int vis_steps = 5; + bool adios2 = false; + + int precision = 8; + cout.precision(precision); + + OptionsParser args(argc, argv); + args.AddOption(&mesh_file, "-m", "--mesh", + "Mesh file to use."); + args.AddOption(&ser_ref_levels, "-rs", "--refine-serial", + "Number of times to refine the mesh uniformly in serial."); + args.AddOption(&par_ref_levels, "-rp", "--refine-parallel", + "Number of times to refine the mesh uniformly in parallel."); + args.AddOption(&order, "-o", "--order", + "Order (degree) of the finite elements."); + args.AddOption(&ode_solver_type, "-s", "--ode-solver", + "ODE solver: 1 - Backward Euler, 2 - SDIRK2, 3 - SDIRK3,\n\t" + "\t 11 - Forward Euler, 12 - RK2, 13 - RK3 SSP, 14 - RK4."); + args.AddOption(&t_final, "-tf", "--t-final", + "Final time; start time is 0."); + args.AddOption(&dt, "-dt", "--time-step", + "Time step."); + args.AddOption(&alpha, "-a", "--alpha", + "Alpha coefficient."); + args.AddOption(&kappa, "-k", "--kappa", + "Kappa coefficient offset."); + args.AddOption(&_in, "-amp-in", "--amplitude-in", + "Amplitude of inlet source at (0,0)."); + args.AddOption(&_out, "-amp-out", "--amplitude-out", + "Amplitude of outlet source at (0.5,0.5)."); + args.AddOption(&t_end_in, "-t-end-in", "--t-end-in", + "End time of inlet source at (0,0)."); + args.AddOption(&t_end_out, "-t-end-out", "--t-end-out", + "End time of outlet source at (0.5,0.5)."); + args.AddOption(&visualization, "-vis", "--visualization", "-no-vis", + "--no-visualization", + "Enable or disable GLVis visualization."); + args.AddOption(&visit, "-visit", "--visit-datafiles", "-no-visit", + "--no-visit-datafiles", + "Save data files for VisIt (visit.llnl.gov) visualization."); + args.AddOption(&vis_steps, "-vs", "--visualization-steps", + "Visualize every n-th timestep."); + args.AddOption(&adios2, "-adios2", "--adios2-streams", "-no-adios2", + "--no-adios2-streams", + "Save data using adios2 streams."); + args.AddOption(&ef, "-ef", "--energy_fraction", + "Energy fraction for DMDc."); + args.AddOption(&rdim, "-rdim", "--rdim", + "Reduced dimension for DMDc."); + args.Parse(); + if (!args.Good()) + { + args.PrintUsage(cout); + MPI_Finalize(); + return 1; + } + + if (myid == 0) + { + args.PrintOptions(cout); + } + + // 3. Read the serial mesh from the given mesh file on all processors. We can + // handle triangular, quadrilateral, tetrahedral and hexahedral meshes + // with the same code. + Mesh *mesh; + if (mesh_file == "") + { + mesh = new Mesh(Mesh::MakeCartesian2D(2, 2, Element::QUADRILATERAL)); + } + else + { + mesh = new Mesh(mesh_file, 1, 1); + } + int dim = mesh->Dimension(); + + mesh->GetBoundingBox(bb_min, bb_max, max(order, 1)); + + // 4. Define the ODE solver used for time integration. Several implicit + // singly diagonal implicit Runge-Kutta (SDIRK) methods, as well as + // explicit Runge-Kutta methods are available. + ODESolver *ode_solver; + bool explicit_solver; + switch (ode_solver_type) + { + // Implicit L-stable methods + case 1: + ode_solver = new BackwardEulerSolver; + explicit_solver = false; + break; + case 2: + ode_solver = new SDIRK23Solver(2); + explicit_solver = false; + break; + case 3: + ode_solver = new SDIRK33Solver; + explicit_solver = false; + break; + // Explicit methods + case 11: + ode_solver = new ForwardEulerSolver; + explicit_solver = true; + break; + case 12: + ode_solver = new RK2Solver(0.5); + explicit_solver = true; + break; // midpoint method + case 13: + ode_solver = new RK3SSPSolver; + explicit_solver = true; + break; + case 14: + ode_solver = new RK4Solver; + explicit_solver = true; + break; + case 15: + ode_solver = new GeneralizedAlphaSolver(0.5); + explicit_solver = true; + break; + // Implicit A-stable methods (not L-stable) + case 22: + ode_solver = new ImplicitMidpointSolver; + explicit_solver = false; + break; + case 23: + ode_solver = new SDIRK23Solver; + explicit_solver = false; + break; + case 24: + ode_solver = new SDIRK34Solver; + explicit_solver = false; + break; + default: + cout << "Unknown ODE solver type: " << ode_solver_type << '\n'; + delete mesh; + return 3; + } + + // 5. Refine the mesh in serial to increase the resolution. In this example + // we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is + // a command-line parameter. + for (int lev = 0; lev < ser_ref_levels; lev++) + { + mesh->UniformRefinement(); + } + + // 6. Define a parallel mesh by a partitioning of the serial mesh. Refine + // this mesh further in parallel to increase the resolution. Once the + // parallel mesh is defined, the serial mesh can be deleted. + ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh); + delete mesh; + for (int lev = 0; lev < par_ref_levels; lev++) + { + pmesh->UniformRefinement(); + } + + // 7. Define the vector finite element space representing the current and the + // initial temperature, u_ref. + H1_FECollection fe_coll(order, dim); + ParFiniteElementSpace fespace(pmesh, &fe_coll); + + int fe_size = fespace.GlobalTrueVSize(); + if (myid == 0) + { + cout << "Number of temperature unknowns: " << fe_size << endl; + } + + ParGridFunction u_gf(&fespace); + + // 8. Set the initial conditions for u. All boundaries are considered + // natural. + FunctionCoefficient u_0(InitialTemperature); + u_gf.ProjectCoefficient(u_0); + Vector u; + u_gf.GetTrueDofs(u); + + // 9. Initialize the conduction operator and the VisIt visualization. + FunctionCoefficient s(SourceFunction); + ConductionOperator oper(fespace, s, alpha, kappa, u); + + u_gf.SetFromTrueDofs(u); + { + ostringstream mesh_name, sol_name; + mesh_name << "heat_conduction_dmdc-mesh." << setfill('0') << setw(6) << myid; + sol_name << "heat_conduction_dmdc-init." << setfill('0') << setw(6) << myid; + ofstream omesh(mesh_name.str().c_str()); + omesh.precision(precision); + pmesh->Print(omesh); + ofstream osol(sol_name.str().c_str()); + osol.precision(precision); + u_gf.Save(osol); + } + + VisItDataCollection visit_dc("Heat_Conduction", pmesh); + visit_dc.RegisterField("temperature", &u_gf); + if (visit) + { + visit_dc.SetCycle(0); + visit_dc.SetTime(0.0); + visit_dc.Save(); + } + + // Optionally output a BP (binary pack) file using ADIOS2. This can be + // visualized with the ParaView VTX reader. +#ifdef MFEM_USE_ADIOS2 + ADIOS2DataCollection* adios2_dc = NULL; + if (adios2) + { + std::string postfix(mesh_file); + postfix.erase(0, std::string("../data/").size() ); + postfix += "_o" + std::to_string(order); + postfix += "_solver" + std::to_string(ode_solver_type); + const std::string collection_name = "heat_conduction_dmdc-p-" + postfix + ".bp"; + + adios2_dc = new ADIOS2DataCollection(MPI_COMM_WORLD, collection_name, pmesh); + adios2_dc->SetParameter("SubStreams", std::to_string(num_procs/2) ); + adios2_dc->RegisterField("temperature", &u_gf); + adios2_dc->SetCycle(0); + adios2_dc->SetTime(0.0); + adios2_dc->Save(); + } +#endif + + socketstream sout; + if (visualization) + { + char vishost[] = "localhost"; + int visport = 19916; + sout.open(vishost, visport); + sout << "parallel " << num_procs << " " << myid << endl; + int good = sout.good(), all_good; + MPI_Allreduce(&good, &all_good, 1, MPI_INT, MPI_MIN, pmesh->GetComm()); + if (!all_good) + { + sout.close(); + visualization = false; + if (myid == 0) + { + cout << "Unable to connect to GLVis server at " + << vishost << ':' << visport << endl; + cout << "GLVis visualization disabled.\n"; + } + } + else + { + sout.precision(precision); + sout << "solution\n" << *pmesh << u_gf; + sout << flush; + } + } + + StopWatch fom_timer, dmd_training_timer, dmd_prediction_timer; + + fom_timer.Start(); + + // 10. Perform time-integration (looping over the time iterations, ti, with a + // time-step dt). + ode_solver->Init(oper); + double t = 0.0; + vector ts; + double f[2]; + fom_timer.Stop(); + + dmd_training_timer.Start(); + + // 11. Create DMDc object and take initial sample. + u_gf.SetFromTrueDofs(u); + f[0] = Amplitude(t, 0); + f[1] = Amplitude(t, 1); + + CAROM::DMDc* dmd_u; + dmd_u = new CAROM::DMDc(u.Size(), 2, dt); + dmd_u->takeSample(u.GetData(), t, f, false); + ts.push_back(t); + + dmd_training_timer.Stop(); + + bool last_step = false; + for (int ti = 1; !last_step; ti++) + { + fom_timer.Start(); + + if (t + dt >= t_final - dt/2) + { + last_step = true; + } + + // Assuming Euler method is used + if (explicit_solver) + { + oper.SetSourceTime(t); + } + else + { + oper.SetSourceTime(t + dt); + } + ode_solver->Step(u, t, dt); + + fom_timer.Stop(); + + dmd_training_timer.Start(); + + u_gf.SetFromTrueDofs(u); + f[0] = Amplitude(t, 0); + f[1] = Amplitude(t, 1); + dmd_u->takeSample(u.GetData(), t, f, last_step); + ts.push_back(t); + + dmd_training_timer.Stop(); + + if (last_step || (ti % vis_steps) == 0) + { + if (myid == 0) + { + cout << "step " << ti << ", t = " << t << endl; + } + + u_gf.SetFromTrueDofs(u); + if (visualization) + { + sout << "parallel " << num_procs << " " << myid << "\n"; + sout << "solution\n" << *pmesh << u_gf << flush; + } + + if (visit) + { + visit_dc.SetCycle(ti); + visit_dc.SetTime(t); + visit_dc.Save(); + } + +#ifdef MFEM_USE_ADIOS2 + if (adios2) + { + adios2_dc->SetCycle(ti); + adios2_dc->SetTime(t); + adios2_dc->Save(); + } +#endif + } + oper.SetParameters(u); + } + +#ifdef MFEM_USE_ADIOS2 + if (adios2) + { + delete adios2_dc; + } +#endif + + // 12. Save the final solution in parallel. This output can be viewed later + // using GLVis: "glvis -np -m heat_conduction_dmdc-mesh -g heat_conduction_dmdc-final". + { + ostringstream sol_name; + sol_name << "heat_conduction_dmdc-final." << setfill('0') << setw(6) << myid; + ofstream osol(sol_name.str().c_str()); + osol.precision(precision); + u_gf.Save(osol); + } + + // 13. Calculate the DMDc modes. + if (myid == 0 && rdim != -1 && ef != -1) + { + std::cout << "Both rdim and ef are set. ef will be ignored." << std::endl; + } + + dmd_training_timer.Start(); + + if (rdim != -1) + { + if (myid == 0) + { + std::cout << "Creating DMDc with rdim: " << rdim << std::endl; + } + dmd_u->train(rdim); + } + else if (ef != -1) + { + if (myid == 0) + { + std::cout << "Creating DMDc with energy fraction: " << ef << std::endl; + } + dmd_u->train(ef); + } + + dmd_training_timer.Stop(); + + Vector true_solution_u(u.Size()); + true_solution_u = u.GetData(); + + dmd_prediction_timer.Start(); + + // 14. Predict the state at t_final using DMDc. + if (myid == 0) + { + std::cout << "Predicting temperature using DMDc" << std::endl; + } + + CAROM::Vector* result_u = dmd_u->predict(ts[0]); + Vector initial_dmd_solution_u(result_u->getData(), result_u->dim()); + u_gf.SetFromTrueDofs(initial_dmd_solution_u); + + VisItDataCollection dmd_visit_dc("DMDc_Heat_Conduction", pmesh); + dmd_visit_dc.RegisterField("temperature", &u_gf); + if (visit) + { + dmd_visit_dc.SetCycle(0); + dmd_visit_dc.SetTime(0.0); + dmd_visit_dc.Save(); + } + + delete result_u; + + if (visit) + { + for (int i = 1; i < ts.size(); i++) + { + if (i == ts.size() - 1 || (i % vis_steps) == 0) + { + result_u = dmd_u->predict(ts[i]); + Vector dmd_solution_u(result_u->getData(), result_u->dim()); + u_gf.SetFromTrueDofs(dmd_solution_u); + + dmd_visit_dc.SetCycle(i); + dmd_visit_dc.SetTime(ts[i]); + dmd_visit_dc.Save(); + + delete result_u; + } + } + } + + dmd_prediction_timer.Stop(); + + result_u = dmd_u->predict(t_final); + + // 15. Calculate the relative error between the DMDc final solution and the true solution. + Vector dmd_solution_u(result_u->getData(), result_u->dim()); + Vector diff_u(true_solution_u.Size()); + subtract(dmd_solution_u, true_solution_u, diff_u); + + double tot_diff_norm_u = sqrt(InnerProduct(MPI_COMM_WORLD, diff_u, diff_u)); + double tot_true_solution_u_norm = sqrt(InnerProduct(MPI_COMM_WORLD, + true_solution_u, true_solution_u)); + + if (myid == 0) + { + + std::cout << "Relative error of DMDc temperature (u) at t_final: " << t_final << + " is " << tot_diff_norm_u / tot_true_solution_u_norm << std::endl; + printf("Elapsed time for solving FOM: %e second\n", fom_timer.RealTime()); + printf("Elapsed time for training DMDc: %e second\n", + dmd_training_timer.RealTime()); + printf("Elapsed time for predicting DMDc: %e second\n", + dmd_prediction_timer.RealTime()); + } + + // 16. Free the used memory. + delete ode_solver; + delete pmesh; + delete result_u; + + MPI_Finalize(); + + return 0; +} + +ConductionOperator::ConductionOperator(ParFiniteElementSpace &f, + FunctionCoefficient &s, + double al, double kap, const Vector &u) + : TimeDependentOperator(f.GetTrueVSize(), 0.0), fespace(f), src_func(s), + M(NULL), K(NULL), T(NULL), current_dt(0.0), + M_solver(f.GetComm()), T_solver(f.GetComm()), z(height) +{ + const double rel_tol = 1e-8; + + M = new ParBilinearForm(&fespace); + M->AddDomainIntegrator(new MassIntegrator()); + M->Assemble(0); // keep sparsity pattern of M and K the same + M->FormSystemMatrix(ess_tdof_list, Mmat); + + M_solver.iterative_mode = false; + M_solver.SetRelTol(rel_tol); + M_solver.SetAbsTol(0.0); + M_solver.SetMaxIter(100); + M_solver.SetPrintLevel(0); + M_prec.SetType(HypreSmoother::Jacobi); + M_solver.SetPreconditioner(M_prec); + M_solver.SetOperator(Mmat); + + alpha = al; + kappa = kap; + + T_solver.iterative_mode = false; + T_solver.SetRelTol(rel_tol); + T_solver.SetAbsTol(0.0); + T_solver.SetMaxIter(100); + T_solver.SetPrintLevel(0); + T_solver.SetPreconditioner(T_prec); + + SetParameters(u); + + B = new ParLinearForm(&fespace); + B->AddDomainIntegrator(new DomainLFIntegrator(src_func)); +} + +void ConductionOperator::Mult(const Vector &u, Vector &du_dt) const +{ + // Compute: + // du_dt = M^{-1}*-K(u) + // for du_dt + Kmat.Mult(u, z); + z.Neg(); // z = -z + z += *b; + M_solver.Mult(z, du_dt); +} + +void ConductionOperator::ImplicitSolve(const double dt, + const Vector &u, Vector &du_dt) +{ + // Solve the equation: + // du_dt = M^{-1}*[-K(u + dt*du_dt)] + // for du_dt + if (!T) + { + T = Add(1.0, Mmat, dt, Kmat); + current_dt = dt; + T_solver.SetOperator(*T); + } + MFEM_VERIFY(dt == current_dt, ""); // SDIRK methods use the same dt + Kmat.Mult(u, z); + z.Neg(); + z += *b; + T_solver.Mult(z, du_dt); +} + +void ConductionOperator::SetSourceTime(double t) +{ + src_func.SetTime(t); + B->Assemble(); + b = B->ParallelAssemble(); +} + +void ConductionOperator::SetParameters(const Vector &u) +{ + ParGridFunction u_alpha_gf(&fespace); + u_alpha_gf.SetFromTrueDofs(u); + for (int i = 0; i < u_alpha_gf.Size(); i++) + { + u_alpha_gf(i) = kappa + alpha*u_alpha_gf(i); + } + + delete K; + K = new ParBilinearForm(&fespace); + + GridFunctionCoefficient u_coeff(&u_alpha_gf); + + K->AddDomainIntegrator(new DiffusionIntegrator(u_coeff)); + K->Assemble(0); // keep sparsity pattern of M and K the same + K->FormSystemMatrix(ess_tdof_list, Kmat); + delete T; + T = NULL; // re-compute T on the next ImplicitSolve +} + +ConductionOperator::~ConductionOperator() +{ + delete T; + delete M; + delete K; + delete B; +} + +double InitialTemperature(const Vector &x) +{ + if (x.Norml2() < 0.5) + { + return 2.0; + } + else + { + return 1.0; + } +} + +double TimeWindowFunction(const double t, const double t_begin, + const double t_end) +{ + return 0.5 * (tanh((t - t_begin) / (5*dt)) - tanh((t - t_end) / (5*dt))); +} + +double Amplitude(const double t, const int index) +{ + if (index == 0) + { + return amp_in * TimeWindowFunction(t, 0.0, t_end_in); + } + else + { + return amp_out * TimeWindowFunction(t, 0.0, t_end_out); + } +} + +double SourceFunction(const Vector &x, double t) +{ + // map to the reference [-1,1] domain + Vector X(x.Size()), Y(x.Size()); + for (int i = 0; i < x.Size(); i++) + { + double center = (bb_min[i] + bb_max[i]) * 0.5; + X(i) = 2 * (x(i) - center) / (bb_max[i] - bb_min[i]); + Y(i) = X(i) - 0.5; + } + + double r1 = X.Norml2() / 0.01; + double r2 = Y.Norml2() / 0.01; + return Amplitude(t,0) * exp(-0.5*r1*r1) - Amplitude(t,1) * exp(-0.5*r2*r2); +} diff --git a/lib/CMakeLists.txt b/lib/CMakeLists.txt index a523dd67d..e65cca1f0 100644 --- a/lib/CMakeLists.txt +++ b/lib/CMakeLists.txt @@ -35,6 +35,7 @@ set(module_list linalg/svd/SVD linalg/svd/StaticSVD algo/DMD + algo/DMDc algo/AdaptiveDMD algo/NonuniformDMD algo/DifferentialEvolution diff --git a/lib/algo/DMD.cpp b/lib/algo/DMD.cpp index cb9f36b8a..4d84d15df 100755 --- a/lib/algo/DMD.cpp +++ b/lib/algo/DMD.cpp @@ -195,6 +195,7 @@ void DMD::takeSample(double* u_in, double t) } d_snapshots.push_back(sample); + Vector* sampled_time = new Vector(&t, 1, false); d_sampled_times.push_back(sampled_time); } diff --git a/lib/algo/DMD.h b/lib/algo/DMD.h index 5231c704b..09ac38e1d 100644 --- a/lib/algo/DMD.h +++ b/lib/algo/DMD.h @@ -114,19 +114,23 @@ class DMD virtual void takeSample(double* u_in, double t); /** + * @brief Train the DMD model with energy fraction criterion. + * * @param[in] energy_fraction The energy fraction to keep after doing SVD. * @param[in] W0 The initial basis to prepend to W. * @param[in] linearity_tol The tolerance for determining whether a column - of W is linearly independent with W0. + * of W is linearly independent with W0. */ virtual void train(double energy_fraction, const Matrix* W0 = NULL, double linearity_tol = 0.0); /** + * @brief Train the DMD model with specified reduced dimension. + * * @param[in] k The number of modes to keep after doing SVD. * @param[in] W0 The initial basis to prepend to W. * @param[in] linearity_tol The tolerance for determining whether a column - of W is linearly independent with W0. + * of W is linearly independent with W0. */ virtual void train(int k, const Matrix* W0 = NULL, double linearity_tol = 0.0); diff --git a/lib/algo/DMDc.cpp b/lib/algo/DMDc.cpp new file mode 100644 index 000000000..95fd97a55 --- /dev/null +++ b/lib/algo/DMDc.cpp @@ -0,0 +1,1053 @@ +/****************************************************************************** + * + * Copyright (c) 2013-2023, Lawrence Livermore National Security, LLC + * and other libROM project developers. See the top-level COPYRIGHT + * file for details. + * + * SPDX-License-Identifier: (Apache-2.0 OR MIT) + * + *****************************************************************************/ + +// Description: Implementation of the DMDc algorithm. + +#include "DMD.h" +#include "DMDc.h" + +#include "linalg/Matrix.h" +#include "linalg/Vector.h" +#include "linalg/scalapack_wrapper.h" +#include "utils/Utilities.h" +#include "utils/CSVDatabase.h" +#include "utils/HDFDatabase.h" +#include "mpi.h" + +#include + +/* Use C++11 built-in shared pointers if available; else fallback to Boost. */ +#if __cplusplus >= 201103L +#include +#else +#include +#endif + +/* Use automatically detected Fortran name-mangling scheme */ +#define zgetrf CAROM_FC_GLOBAL(zgetrf, ZGETRF) +#define zgetri CAROM_FC_GLOBAL(zgetri, ZGETRI) + +extern "C" { + // LU decomposition of a general matrix. + void zgetrf(int*, int*, double*, int*, int*, int*); + + // Generate inverse of a matrix given its LU decomposition. + void zgetri(int*, double*, int*, int*, double*, int*, int*); +} + +namespace CAROM { + +DMDc::DMDc(int dim, int dim_c, Vector* state_offset) +{ + CAROM_VERIFY(dim > 0); + CAROM_VERIFY(dim_c > 0); + + // Get the rank of this process, and the number of processors. + int mpi_init; + MPI_Initialized(&mpi_init); + if (mpi_init == 0) { + MPI_Init(nullptr, nullptr); + } + + MPI_Comm_rank(MPI_COMM_WORLD, &d_rank); + MPI_Comm_size(MPI_COMM_WORLD, &d_num_procs); + d_dim = dim; + d_dim_c = dim_c; + d_trained = false; + d_init_projected = false; + setOffset(state_offset); +} + +DMDc::DMDc(int dim, int dim_c, double dt, Vector* state_offset) +{ + CAROM_VERIFY(dim > 0); + CAROM_VERIFY(dim_c > 0); + CAROM_VERIFY(dt > 0.0); + + // Get the rank of this process, and the number of processors. + int mpi_init; + MPI_Initialized(&mpi_init); + if (mpi_init == 0) { + MPI_Init(nullptr, nullptr); + } + + MPI_Comm_rank(MPI_COMM_WORLD, &d_rank); + MPI_Comm_size(MPI_COMM_WORLD, &d_num_procs); + d_dim = dim; + d_dim_c = dim_c; + d_dt = dt; + d_trained = false; + d_init_projected = false; + setOffset(state_offset); +} + +DMDc::DMDc(std::string base_file_name) +{ + // Get the rank of this process, and the number of processors. + int mpi_init; + MPI_Initialized(&mpi_init); + if (mpi_init == 0) { + MPI_Init(nullptr, nullptr); + } + + MPI_Comm_rank(MPI_COMM_WORLD, &d_rank); + MPI_Comm_size(MPI_COMM_WORLD, &d_num_procs); + d_trained = true; + d_init_projected = true; + + load(base_file_name); +} + +DMDc::DMDc(std::vector> eigs, Matrix* phi_real, + Matrix* phi_imaginary, Matrix* B_tilde, int k, + double dt, double t_offset, Vector* state_offset) +{ + // Get the rank of this process, and the number of processors. + int mpi_init; + MPI_Initialized(&mpi_init); + if (mpi_init == 0) { + MPI_Init(nullptr, nullptr); + } + + MPI_Comm_rank(MPI_COMM_WORLD, &d_rank); + MPI_Comm_size(MPI_COMM_WORLD, &d_num_procs); + d_trained = true; + d_init_projected = false; + + d_eigs = eigs; + d_phi_real = phi_real; + d_phi_imaginary = phi_imaginary; + d_B_tilde = B_tilde; + d_k = k; + d_dt = dt; + d_t_offset = t_offset; + setOffset(state_offset); +} + +DMDc::~DMDc() +{ + for (auto snapshot : d_snapshots) + { + delete snapshot; + } + for (auto sampled_time : d_sampled_times) + { + delete sampled_time; + } + delete d_state_offset; + delete d_basis; + delete d_A_tilde; + delete d_B_tilde; + delete d_phi_real; + delete d_phi_imaginary; + delete d_phi_real_squared_inverse; + delete d_phi_imaginary_squared_inverse; + delete d_projected_init_real; + delete d_projected_init_imaginary; +} + +void DMDc::setOffset(Vector* offset_vector) +{ + d_state_offset = offset_vector; +} + +void DMDc::takeSample(double* u_in, double t, double* f_in, bool last_step) +{ + CAROM_VERIFY(u_in != 0); + CAROM_VERIFY(t >= 0.0); + Vector* sample = new Vector(u_in, d_dim, true); + + double orig_t = t; + if (d_snapshots.empty()) + { + d_t_offset = t; + t = 0.0; + } + else + { + t -= d_t_offset; + } + + // Erase any snapshots taken at the same or later time + while (!d_sampled_times.empty() && d_sampled_times.back()->item(0) >= t) + { + if (d_rank == 0) std::cout << "Removing existing snapshot at time: " << + d_t_offset + d_sampled_times.back()->item(0) << std::endl; + Vector* last_snapshot = d_snapshots.back(); + delete last_snapshot; + d_snapshots.pop_back(); + d_controls.pop_back(); + d_sampled_times.pop_back(); + } + + if (d_snapshots.empty()) + { + d_t_offset = orig_t; + t = 0.0; + } + else + { + CAROM_VERIFY(d_sampled_times.back()->item(0) < t); + } + d_snapshots.push_back(sample); + + if (!last_step) + { + Vector* control = new Vector(f_in, d_dim_c, false); + d_controls.push_back(control); + } + + Vector* sampled_time = new Vector(&t, 1, false); + d_sampled_times.push_back(sampled_time); +} + +void DMDc::train(double energy_fraction, const Matrix* B) +{ + const Matrix* f_snapshots = getSnapshotMatrix(); + const Matrix* f_controls = createSnapshotMatrix(d_controls); + CAROM_VERIFY(f_snapshots->numColumns() > 1); + CAROM_VERIFY(f_controls->numColumns() == f_snapshots->numColumns() - 1); + CAROM_VERIFY(energy_fraction > 0 && energy_fraction <= 1); + d_energy_fraction = energy_fraction; + constructDMDc(f_snapshots, f_controls, d_rank, d_num_procs, B); + + delete f_snapshots; +} + +void DMDc::train(int k, const Matrix* B) +{ + const Matrix* f_snapshots = getSnapshotMatrix(); + const Matrix* f_controls = createSnapshotMatrix(d_controls); + CAROM_VERIFY(f_snapshots->numColumns() > 1); + CAROM_VERIFY(f_controls->numColumns() == f_snapshots->numColumns() - 1); + CAROM_VERIFY(k > 0 && k <= f_snapshots->numColumns() - 1); + d_energy_fraction = -1.0; + d_k = k; + constructDMDc(f_snapshots, f_controls, d_rank, d_num_procs, B); + + delete f_snapshots; +} + +std::pair +DMDc::computeDMDcSnapshotPair(const Matrix* snapshots, const Matrix* controls, + const Matrix* B) +{ + CAROM_VERIFY(snapshots->numColumns() > 1); + CAROM_VERIFY(controls->numColumns() == snapshots->numColumns()-1); + + // TODO: Making two copies of the snapshot matrix has a lot of overhead. + // We need to figure out a way to do submatrix multiplication and to + // reimplement this algorithm using one snapshot matrix. + int input_control_dim = (B == NULL) ? controls->numRows() : 0; + Matrix* f_snapshots_in = new Matrix(snapshots->numRows() + input_control_dim, + snapshots->numColumns() - 1, snapshots->distributed()); + Matrix* f_snapshots_out = new Matrix(snapshots->numRows(), + snapshots->numColumns() - 1, snapshots->distributed()); + + // Break up snapshots into snapshots_in and snapshots_out + // snapshots_in = all columns of snapshots except last + // snapshots_out = all columns of snapshots except first + for (int i = 0; i < snapshots->numRows(); i++) + { + for (int j = 0; j < snapshots->numColumns() - 1; j++) + { + f_snapshots_in->item(i, j) = snapshots->item(i, j); + f_snapshots_out->item(i, j) = snapshots->item(i, j + 1); + if (d_state_offset) + { + f_snapshots_in->item(i, j) -= d_state_offset->item(i); + f_snapshots_out->item(i, j) -= d_state_offset->item(i); + } + } + } + + for (int i = 0; i < controls->numRows(); i++) + { + if (B == NULL) + { + for (int j = 0; j < snapshots->numColumns() - 1; j++) + { + f_snapshots_in->item(snapshots->numRows() + i, j) = controls->item(i, j); + } + } + else + { + Matrix* Bf = B->mult(controls); + *f_snapshots_out -= *Bf; + delete Bf; + } + } + + return std::pair(f_snapshots_in, f_snapshots_out); +} + +void +DMDc::constructDMDc(const Matrix* f_snapshots, + const Matrix* f_controls, + int d_rank, + int d_num_procs, + const Matrix* B) +{ + std::pair f_snapshot_pair = computeDMDcSnapshotPair( + f_snapshots, f_controls, B); + Matrix* f_snapshots_in = f_snapshot_pair.first; + Matrix* f_snapshots_out = f_snapshot_pair.second; + + int *row_offset = new int[d_num_procs + 1]; + row_offset[d_num_procs] = f_snapshots_in->numDistributedRows(); + row_offset[d_rank] = f_snapshots_in->numRows(); + + CAROM_VERIFY(MPI_Allgather(MPI_IN_PLACE, + 1, + MPI_INT, + row_offset, + 1, + MPI_INT, + MPI_COMM_WORLD) == MPI_SUCCESS); + for (int i = d_num_procs - 1; i >= 0; i--) { + row_offset[i] = row_offset[i + 1] - row_offset[i]; + } + + CAROM_VERIFY(row_offset[0] == 0); + + int d_blocksize = row_offset[d_num_procs] / d_num_procs; + if (row_offset[d_num_procs] % d_num_procs != 0) d_blocksize += 1; + + SLPK_Matrix svd_input; + + // Calculate svd of snapshots_in + initialize_matrix(&svd_input, f_snapshots_in->numColumns(), + f_snapshots_in->numDistributedRows(), + 1, d_num_procs, d_blocksize, d_blocksize); + + for (int rank = 0; rank < d_num_procs; ++rank) + { + scatter_block(&svd_input, 1, row_offset[rank] + 1, + f_snapshots_in->getData(), + f_snapshots_in->numColumns(), + row_offset[rank + 1] - row_offset[rank], rank); + } + + std::unique_ptr d_factorizer_in(new SVDManager); + + // This block does the actual ScaLAPACK call to do the factorization. + svd_init(d_factorizer_in.get(), &svd_input); + d_factorizer_in->dov = 1; + factorize(d_factorizer_in.get()); + free_matrix_data(&svd_input); + + // Compute how many basis vectors we will actually use. + d_num_singular_vectors = std::min(f_snapshots_in->numColumns(), + f_snapshots_in->numDistributedRows()); + for (int i = 0; i < d_num_singular_vectors; i++) + { + d_sv.push_back(d_factorizer_in->S[i]); + } + + int d_k_in; + if (d_energy_fraction != -1.0) + { + d_k_in = d_num_singular_vectors; + if (d_energy_fraction < 1.0) + { + double total_energy = 0.0; + for (int i = 0; i < d_num_singular_vectors; i++) + { + total_energy += d_factorizer_in->S[i]; + } + double current_energy = 0.0; + for (int i = 0; i < d_num_singular_vectors; i++) + { + current_energy += d_factorizer_in->S[i]; + if (current_energy / total_energy >= d_energy_fraction) + { + d_k_in = i + 1; + break; + } + } + } + } + + if (d_rank == 0) std::cout << "Using " << d_k_in << " basis vectors out of " << + d_num_singular_vectors << " for input." << std::endl; + + // Allocate the appropriate matrices and gather their elements. + Matrix* d_basis_in = new Matrix(f_snapshots_in->numRows(), d_k_in, + f_snapshots_in->distributed()); + Matrix* d_S_inv = new Matrix(d_k_in, d_k_in, false); + Matrix* d_basis_right = new Matrix(f_snapshots_in->numColumns(), d_k_in, false); + + for (int d_rank = 0; d_rank < d_num_procs; ++d_rank) { + // V is computed in the transposed order so no reordering necessary. + gather_block(&d_basis_in->item(0, 0), d_factorizer_in->V, + 1, row_offset[static_cast(d_rank)]+1, + d_k_in, row_offset[static_cast(d_rank) + 1] - + row_offset[static_cast(d_rank)], + d_rank); + + // gather_transposed_block does the same as gather_block, but transposes + // it; here, it is used to go from column-major to row-major order. + gather_transposed_block(&d_basis_right->item(0, 0), d_factorizer_in->U, 1, 1, + f_snapshots_in->numColumns(), d_k_in, d_rank); + } + + // Get inverse of singular values by multiplying by reciprocal. + for (int i = 0; i < d_k_in; ++i) + { + d_S_inv->item(i, i) = 1 / d_factorizer_in->S[static_cast(i)]; + } + + if (B == NULL) + { + // SVD on outputs + row_offset[d_num_procs] = f_snapshots_out->numDistributedRows(); + row_offset[d_rank] = f_snapshots_out->numRows(); + + CAROM_VERIFY(MPI_Allgather(MPI_IN_PLACE, + 1, + MPI_INT, + row_offset, + 1, + MPI_INT, + MPI_COMM_WORLD) == MPI_SUCCESS); + for (int i = d_num_procs - 1; i >= 0; i--) { + row_offset[i] = row_offset[i + 1] - row_offset[i]; + } + + CAROM_VERIFY(row_offset[0] == 0); + + int d_blocksize = row_offset[d_num_procs] / d_num_procs; + if (row_offset[d_num_procs] % d_num_procs != 0) d_blocksize += 1; + + SLPK_Matrix svd_output; + + // Calculate svd of snapshots_out + initialize_matrix(&svd_output, f_snapshots_out->numColumns(), + f_snapshots_out->numDistributedRows(), + 1, d_num_procs, d_blocksize, d_blocksize); + + for (int rank = 0; rank < d_num_procs; ++rank) + { + scatter_block(&svd_output, 1, row_offset[rank] + 1, + f_snapshots_out->getData(), + f_snapshots_out->numColumns(), + row_offset[rank + 1] - row_offset[rank], rank); + } + + std::unique_ptr d_factorizer_out(new SVDManager); + + // This block does the actual ScaLAPACK call to do the factorization. + svd_init(d_factorizer_out.get(), &svd_output); + d_factorizer_out->dov = 1; + factorize(d_factorizer_out.get()); + free_matrix_data(&svd_output); + + // Compute how many basis vectors we will actually use. + d_num_singular_vectors = std::min(f_snapshots_out->numColumns(), + f_snapshots_out->numDistributedRows()); + for (int i = 0; i < d_num_singular_vectors; i++) + { + d_sv.push_back(d_factorizer_out->S[i]); + } + + if (d_energy_fraction != -1.0) + { + d_k = d_num_singular_vectors; + if (d_energy_fraction < 1.0) + { + double total_energy = 0.0; + for (int i = 0; i < d_num_singular_vectors; i++) + { + total_energy += d_factorizer_out->S[i]; + } + double current_energy = 0.0; + for (int i = 0; i < d_num_singular_vectors; i++) + { + current_energy += d_factorizer_out->S[i]; + if (current_energy / total_energy >= d_energy_fraction) + { + d_k = i + 1; + break; + } + } + } + } + + if (d_rank == 0) std::cout << "Using " << d_k << " basis vectors out of " << + d_num_singular_vectors << " for output." << std::endl; + + // Allocate the appropriate matrices and gather their elements. + d_basis = new Matrix(f_snapshots_out->numRows(), d_k, + f_snapshots_out->distributed()); + for (int d_rank = 0; d_rank < d_num_procs; ++d_rank) { + // V is computed in the transposed order so no reordering necessary. + gather_block(&d_basis->item(0, 0), d_factorizer_out->V, + 1, row_offset[static_cast(d_rank)]+1, + d_k, row_offset[static_cast(d_rank) + 1] - + row_offset[static_cast(d_rank)], + d_rank); + } + } + else + { + d_basis = d_basis_in; + d_k = d_k_in; + } + + delete[] row_offset; + + // Calculate A_tilde and B_tilde + Matrix* d_basis_mult_f_snapshots_out = d_basis->transposeMult(f_snapshots_out); + Matrix* d_basis_mult_f_snapshots_out_mult_d_basis_right = + d_basis_mult_f_snapshots_out->mult(d_basis_right); + Matrix* d_A_tilde_orig = d_basis_mult_f_snapshots_out_mult_d_basis_right->mult( + d_S_inv); + + if (B == NULL) + { + Matrix* d_basis_in_state = new Matrix(f_snapshots->numRows(), + d_k_in, f_snapshots->distributed()); + Matrix* d_basis_in_control_transpose = new Matrix(d_k_in, f_controls->numRows(), + false); + for (int j = 0; j < d_k_in; j++) + { + for (int i = 0; i < f_snapshots->numRows(); i++) + { + d_basis_in_state->item(i, j) = d_basis_in->item(i, j); + } + for (int i = 0; i < f_controls->numRows(); i++) + { + d_basis_in_control_transpose->item(j, + i) = d_basis_in->item(f_snapshots->numRows() + i, + j); + } + } + Matrix* d_basis_state_rot = d_basis_in_state->transposeMult(d_basis); + d_A_tilde = d_A_tilde_orig->mult(d_basis_state_rot); + d_B_tilde = d_A_tilde_orig->mult(d_basis_in_control_transpose); + delete d_basis_in_state; + delete d_basis_in_control_transpose; + delete d_basis_state_rot; + } + else + { + d_A_tilde = d_A_tilde_orig; + d_B_tilde = d_basis->transposeMult(B); + } + + // Calculate the right eigenvalues/eigenvectors of A_tilde + ComplexEigenPair eigenpair = NonSymmetricRightEigenSolve(d_A_tilde); + d_eigs = eigenpair.eigs; + + //struct DMDInternal dmd_internal = {f_snapshots_in, f_snapshots_out, d_basis, d_basis_right, d_S_inv, &eigenpair}; + //computePhi(dmd_internal); + + d_phi_real = d_basis->mult(eigenpair.ev_real); + d_phi_imaginary = d_basis->mult(eigenpair.ev_imaginary); + + Vector* init = new Vector(d_basis->numRows(), true); + for (int i = 0; i < init->dim(); i++) + { + init->item(i) = f_snapshots_in->item(i, 0); + } + + // Calculate the projection initial_condition onto column space of d_basis. + project(init, f_controls); + + d_trained = true; + + delete d_basis_right; + delete d_S_inv; + delete d_basis_mult_f_snapshots_out; + delete d_basis_mult_f_snapshots_out_mult_d_basis_right; + delete d_A_tilde_orig; + delete f_snapshots_in; + delete f_snapshots_out; + delete eigenpair.ev_real; + delete eigenpair.ev_imaginary; + delete init; + + release_context(&svd_input); +} + +void +DMDc::project(const Vector* init, const Matrix* controls, double t_offset) +{ + Matrix* d_phi_real_squared = d_phi_real->transposeMult(d_phi_real); + Matrix* d_phi_real_squared_2 = d_phi_imaginary->transposeMult(d_phi_imaginary); + *d_phi_real_squared += *d_phi_real_squared_2; + + Matrix* d_phi_imaginary_squared = d_phi_real->transposeMult(d_phi_imaginary); + Matrix* d_phi_imaginary_squared_2 = d_phi_imaginary->transposeMult(d_phi_real); + *d_phi_imaginary_squared -= *d_phi_imaginary_squared_2; + + delete d_phi_real_squared_2; + delete d_phi_imaginary_squared_2; + + const int dprs_row = d_phi_real_squared->numRows(); + const int dprs_col = d_phi_real_squared->numColumns(); + double* inverse_input = new double[dprs_row * dprs_col * 2]; + for (int i = 0; i < d_phi_real_squared->numRows(); i++) + { + int k = 0; + for (int j = 0; j < d_phi_real_squared->numColumns() * 2; j++) + { + if (j % 2 == 0) + { + inverse_input[d_phi_real_squared->numColumns() * 2 * i + j] = + d_phi_real_squared->item(i, k); + } + else + { + inverse_input[d_phi_imaginary_squared->numColumns() * 2 * i + j] = + d_phi_imaginary_squared->item(i, k); + k++; + } + } + } + + // Call lapack routines to do the inversion. + // Set up some stuff the lapack routines need. + int info; + int mtx_size = d_phi_real_squared->numColumns(); + int lwork = mtx_size*mtx_size*std::max(10,d_num_procs); + int* ipiv = new int [mtx_size]; + double* work = new double [lwork]; + + // Now call lapack to do the inversion. + zgetrf(&mtx_size, &mtx_size, inverse_input, &mtx_size, ipiv, &info); + zgetri(&mtx_size, inverse_input, &mtx_size, ipiv, work, &lwork, &info); + + d_phi_real_squared_inverse = d_phi_real_squared; + d_phi_imaginary_squared_inverse = d_phi_imaginary_squared; + + for (int i = 0; i < d_phi_real_squared_inverse->numRows(); i++) + { + int k = 0; + for (int j = 0; j < d_phi_real_squared_inverse->numColumns() * 2; j++) + { + if (j % 2 == 0) + { + d_phi_real_squared_inverse->item(i, k) = + inverse_input[d_phi_real_squared_inverse->numColumns() * 2 * i + j]; + } + else + { + d_phi_imaginary_squared_inverse->item(i, k) = + inverse_input[d_phi_imaginary_squared_inverse->numColumns() * 2 * i + j]; + k++; + } + } + } + + // Initial condition + Vector* init_real = d_phi_real->transposeMult(init); + Vector* init_imaginary = d_phi_imaginary->transposeMult(init); + + Vector* d_projected_init_real_1 = d_phi_real_squared_inverse->mult(init_real); + Vector* d_projected_init_real_2 = d_phi_imaginary_squared_inverse->mult( + init_imaginary); + d_projected_init_real = d_projected_init_real_1->plus(d_projected_init_real_2); + + Vector* d_projected_init_imaginary_1 = d_phi_real_squared_inverse->mult( + init_imaginary); + Vector* d_projected_init_imaginary_2 = d_phi_imaginary_squared_inverse->mult( + init_real); + d_projected_init_imaginary = d_projected_init_imaginary_2->minus( + d_projected_init_imaginary_1); + + delete init_real; + delete init_imaginary; + delete d_projected_init_real_1; + delete d_projected_init_real_2; + delete d_projected_init_imaginary_1; + delete d_projected_init_imaginary_2; + + // Controls + Matrix* B_tilde_f = d_B_tilde->mult(controls); + Matrix* UBf = d_basis->mult(B_tilde_f); + Matrix* controls_real = d_phi_real->transposeMult(UBf); + Matrix* controls_imaginary = d_phi_imaginary->transposeMult(UBf); + + d_projected_controls_real = d_phi_real_squared_inverse->mult(controls_real); + Matrix* d_projected_controls_real_2 = d_phi_imaginary_squared_inverse->mult( + controls_imaginary); + *d_projected_controls_real += *d_projected_controls_real_2; + + d_projected_controls_imaginary = d_phi_imaginary_squared_inverse->mult( + controls_real); + Matrix* d_projected_controls_imaginary_2 = d_phi_real_squared_inverse->mult( + controls_imaginary); + *d_projected_controls_imaginary -= *d_projected_controls_imaginary_2; + + delete controls_real; + delete controls_imaginary; + delete d_projected_controls_real_2; + delete d_projected_controls_imaginary_2; + + delete [] inverse_input; + delete [] ipiv; + delete [] work; + + if (t_offset >= 0.0) + { + std::cout << "t_offset is updated from " << d_t_offset + << " to " << t_offset << std::endl; + d_t_offset = t_offset; + } + d_init_projected = true; +} + +Vector* +DMDc::predict(double t) +{ + CAROM_VERIFY(d_trained); + CAROM_VERIFY(d_init_projected); + CAROM_VERIFY(t >= 0.0); + + t -= d_t_offset; + + int n = round(t / d_dt); + //int n = min(round(t / d_dt), d_projected_controls_real->numColumns()); + + std::pair d_phi_pair = phiMultEigs(t); + Matrix* d_phi_mult_eigs_real = d_phi_pair.first; + Matrix* d_phi_mult_eigs_imaginary = d_phi_pair.second; + + Vector* d_predicted_state_real_1 = d_phi_mult_eigs_real->mult( + d_projected_init_real); + Vector* d_predicted_state_real_2 = d_phi_mult_eigs_imaginary->mult( + d_projected_init_imaginary); + Vector* d_predicted_state_real = d_predicted_state_real_1->minus( + d_predicted_state_real_2); + addOffset(d_predicted_state_real); + + delete d_phi_mult_eigs_real; + delete d_phi_mult_eigs_imaginary; + delete d_predicted_state_real_1; + delete d_predicted_state_real_2; + + Vector* f_control_real = new Vector(d_basis->numRows(), false); + Vector* f_control_imaginary = new Vector(d_basis->numRows(), false); + for (int k = 0; k < n; k++) + { + t -= d_dt; + std::pair d_phi_pair = phiMultEigs(t); + d_phi_mult_eigs_real = d_phi_pair.first; + d_phi_mult_eigs_imaginary = d_phi_pair.second; + + d_projected_controls_real->getColumn(k, *f_control_real); + d_projected_controls_imaginary->getColumn(k, *f_control_imaginary); + d_predicted_state_real_1 = d_phi_mult_eigs_real->mult( + f_control_real); + d_predicted_state_real_2 = d_phi_mult_eigs_imaginary->mult( + f_control_imaginary); + *d_predicted_state_real += *d_predicted_state_real_1; + *d_predicted_state_real -= *d_predicted_state_real_2; + + delete d_phi_mult_eigs_real; + delete d_phi_mult_eigs_imaginary; + delete d_predicted_state_real_1; + delete d_predicted_state_real_2; + } + + delete f_control_real; + delete f_control_imaginary; + return d_predicted_state_real; +} + +void +DMDc::addOffset(Vector*& result) +{ + if (d_state_offset) + { + *result += *d_state_offset; + } +} + +std::complex +DMDc::computeEigExp(std::complex eig, double t) +{ + return std::pow(eig, t / d_dt); +} + +std::pair +DMDc::phiMultEigs(double t) +{ + Matrix* d_eigs_exp_real = new Matrix(d_k, d_k, false); + Matrix* d_eigs_exp_imaginary = new Matrix(d_k, d_k, false); + + for (int i = 0; i < d_k; i++) + { + std::complex eig_exp = computeEigExp(d_eigs[i], t); + d_eigs_exp_real->item(i, i) = std::real(eig_exp); + d_eigs_exp_imaginary->item(i, i) = std::imag(eig_exp); + } + + Matrix* d_phi_mult_eigs_real = d_phi_real->mult(d_eigs_exp_real); + Matrix* d_phi_mult_eigs_real_2 = d_phi_imaginary->mult(d_eigs_exp_imaginary); + *d_phi_mult_eigs_real -= *d_phi_mult_eigs_real_2; + Matrix* d_phi_mult_eigs_imaginary = d_phi_real->mult(d_eigs_exp_imaginary); + Matrix* d_phi_mult_eigs_imaginary_2 = d_phi_imaginary->mult(d_eigs_exp_real); + *d_phi_mult_eigs_imaginary += *d_phi_mult_eigs_imaginary_2; + + delete d_eigs_exp_real; + delete d_eigs_exp_imaginary; + delete d_phi_mult_eigs_real_2; + delete d_phi_mult_eigs_imaginary_2; + + return std::pair(d_phi_mult_eigs_real, + d_phi_mult_eigs_imaginary); +} + +double +DMDc::getTimeOffset() const +{ + return d_t_offset; +} + +const Matrix* +DMDc::getSnapshotMatrix() +{ + return createSnapshotMatrix(d_snapshots); +} + +const Matrix* +DMDc::createSnapshotMatrix(std::vector snapshots) +{ + CAROM_VERIFY(snapshots.size() > 0); + CAROM_VERIFY(snapshots[0]->dim() > 0); + for (int i = 0 ; i < snapshots.size() - 1; i++) + { + CAROM_VERIFY(snapshots[i]->dim() == snapshots[i + 1]->dim()); + CAROM_VERIFY(snapshots[i]->distributed() == snapshots[i + 1]->distributed()); + } + + Matrix* snapshot_mat = new Matrix(snapshots[0]->dim(), snapshots.size(), + snapshots[0]->distributed()); + + for (int i = 0; i < snapshots[0]->dim(); i++) + { + for (int j = 0; j < snapshots.size(); j++) + { + snapshot_mat->item(i, j) = snapshots[j]->item(i); + } + } + + return snapshot_mat; +} + +void +DMDc::load(std::string base_file_name) +{ + CAROM_ASSERT(!base_file_name.empty()); + + char tmp[100]; + std::string full_file_name = base_file_name; + HDFDatabase database; + database.open(full_file_name, "r"); + + sprintf(tmp, "dt"); + database.getDouble(tmp, d_dt); + + sprintf(tmp, "t_offset"); + database.getDouble(tmp, d_t_offset); + + sprintf(tmp, "k"); + database.getInteger(tmp, d_k); + + sprintf(tmp, "num_eigs"); + int num_eigs; + database.getInteger(tmp, num_eigs); + + std::vector eigs_real; + std::vector eigs_imag; + + sprintf(tmp, "eigs_real"); + eigs_real.resize(num_eigs); + database.getDoubleArray(tmp, &eigs_real[0], num_eigs); + + sprintf(tmp, "eigs_imag"); + eigs_imag.resize(num_eigs); + database.getDoubleArray(tmp, &eigs_imag[0], num_eigs); + + for (int i = 0; i < num_eigs; i++) + { + d_eigs.push_back(std::complex(eigs_real[i], eigs_imag[i])); + } + database.close(); + + full_file_name = base_file_name + "_basis"; + d_basis = new Matrix(); + d_basis->read(full_file_name); + + full_file_name = base_file_name + "_A_tilde"; + d_A_tilde = new Matrix(); + d_A_tilde->read(full_file_name); + + full_file_name = base_file_name + "_B_tilde"; + d_B_tilde = new Matrix(); + d_B_tilde->read(full_file_name); + + full_file_name = base_file_name + "_phi_real"; + d_phi_real = new Matrix(); + d_phi_real->read(full_file_name); + + full_file_name = base_file_name + "_phi_imaginary"; + d_phi_imaginary = new Matrix(); + d_phi_imaginary->read(full_file_name); + + full_file_name = base_file_name + "_phi_real_squared_inverse"; + d_phi_real_squared_inverse = new Matrix(); + d_phi_real_squared_inverse->read(full_file_name); + + full_file_name = base_file_name + "_phi_imaginary_squared_inverse"; + d_phi_imaginary_squared_inverse = new Matrix(); + d_phi_imaginary_squared_inverse->read(full_file_name); + + full_file_name = base_file_name + "_projected_init_real"; + d_projected_init_real = new Vector(); + d_projected_init_real->read(full_file_name); + + full_file_name = base_file_name + "_projected_init_imaginary"; + d_projected_init_imaginary = new Vector(); + d_projected_init_imaginary->read(full_file_name); + + full_file_name = base_file_name + "_state_offset"; + if (Utilities::file_exist(full_file_name + ".000000")) + { + d_state_offset = new Vector(); + d_state_offset->read(full_file_name); + } + + d_init_projected = true; + d_trained = true; + + MPI_Barrier(MPI_COMM_WORLD); +} + +void +DMDc::load(const char* base_file_name) +{ + load(std::string(base_file_name)); +} + +void +DMDc::save(std::string base_file_name) +{ + CAROM_ASSERT(!base_file_name.empty()); + CAROM_VERIFY(d_trained); + + if (d_rank == 0) + { + char tmp[100]; + std::string full_file_name = base_file_name; + HDFDatabase database; + database.create(full_file_name); + + sprintf(tmp, "dt"); + database.putDouble(tmp, d_dt); + + sprintf(tmp, "t_offset"); + database.putDouble(tmp, d_t_offset); + + sprintf(tmp, "k"); + database.putInteger(tmp, d_k); + + sprintf(tmp, "num_eigs"); + database.putInteger(tmp, d_eigs.size()); + + std::vector eigs_real; + std::vector eigs_imag; + + for (int i = 0; i < d_eigs.size(); i++) + { + eigs_real.push_back(d_eigs[i].real()); + eigs_imag.push_back(d_eigs[i].imag()); + } + + sprintf(tmp, "eigs_real"); + database.putDoubleArray(tmp, &eigs_real[0], eigs_real.size()); + + sprintf(tmp, "eigs_imag"); + database.putDoubleArray(tmp, &eigs_imag[0], eigs_imag.size()); + database.close(); + } + + std::string full_file_name; + + if (d_basis != NULL) + { + full_file_name = base_file_name + "_basis"; + d_basis->write(full_file_name); + } + + if (d_A_tilde != NULL) + { + full_file_name = base_file_name + "_A_tilde"; + d_A_tilde->write(full_file_name); + } + + if (d_B_tilde != NULL) + { + full_file_name = base_file_name + "_B_tilde"; + d_B_tilde->write(full_file_name); + } + + full_file_name = base_file_name + "_phi_real"; + d_phi_real->write(full_file_name); + + full_file_name = base_file_name + "_phi_imaginary"; + d_phi_imaginary->write(full_file_name); + + full_file_name = base_file_name + "_phi_real_squared_inverse"; + d_phi_real_squared_inverse->write(full_file_name); + + full_file_name = base_file_name + "_phi_imaginary_squared_inverse"; + d_phi_imaginary_squared_inverse->write(full_file_name); + + full_file_name = base_file_name + "_projected_init_real"; + d_projected_init_real->write(full_file_name); + + full_file_name = base_file_name + "_projected_init_imaginary"; + d_projected_init_imaginary->write(full_file_name); + + if (d_state_offset != NULL) + { + full_file_name = base_file_name + "_state_offset"; + d_state_offset->write(full_file_name); + } + + MPI_Barrier(MPI_COMM_WORLD); +} + +void +DMDc::save(const char* base_file_name) +{ + save(std::string(base_file_name)); +} + +void +DMDc::summary(std::string base_file_name) +{ + if (d_rank == 0) + { + CSVDatabase* csv_db(new CSVDatabase); + + csv_db->putDoubleVector(base_file_name + "_singular_value.csv", d_sv, + d_num_singular_vectors); + csv_db->putComplexVector(base_file_name + "_eigenvalue.csv", d_eigs, + d_eigs.size()); + + delete csv_db; + } +} + +} diff --git a/lib/algo/DMDc.h b/lib/algo/DMDc.h new file mode 100644 index 000000000..4f0700265 --- /dev/null +++ b/lib/algo/DMDc.h @@ -0,0 +1,408 @@ +/****************************************************************************** + * + * Copyright (c) 2013-2023, Lawrence Livermore National Security, LLC + * and other libROM project developers. See the top-level COPYRIGHT + * file for details. + * + * SPDX-License-Identifier: (Apache-2.0 OR MIT) + * + *****************************************************************************/ + +// Description: Computes the DMDc algorithm on the given snapshot matrix and control matrix. The +// implemented dynamic mode decomposition with control algorithm is derived from +// Proctor et. al's paper "Dynamic mode decomposition with control": +// https://arxiv.org/abs/1409.6358 + +#ifndef included_DMDc_h +#define included_DMDc_h + +#include +#include + +namespace CAROM { + +class Matrix; +class Vector; +class ComplexEigenPair; + +/** + * Class DMDc implements the DMDc algorithm on a given snapshot matrix. + */ +class DMDc +{ +public: + + /** + * @brief Constructor. Basic DMDc with uniform time step size. + * + * @param[in] dim The full-order state dimension. + * @param[in] dim_c The control dimension. + * @param[in] dt The dt between samples. + * @param[in] state_offset The state offset. + */ + DMDc(int dim, int dim_c, double dt, Vector* state_offset = NULL); + + /** + * @brief Constructor. DMDc from saved models. + * + * @param[in] base_file_name The base part of the filename of the + * database to load when restarting from a save. + */ + DMDc(std::string base_file_name); + + /** + * @brief Destroy the DMDc object + */ + virtual ~DMDc(); + + /** + * @brief Set the state offset. + */ + virtual void setOffset(Vector* offset_vector); + + /** + * @brief Sample the new state, u_in. Any samples in d_snapshots + * taken at the same or later time will be erased. + * + * @pre u_in != 0 + * @pre t >= 0.0 + * + * @param[in] u_in The new state. + * @param[in] t The time of the newly sampled state. + * @param[in] f_in The control. + * @param[in] last_step Whether it is the last step. + */ + virtual void takeSample(double* u_in, double t, double* f_in, + bool last_step = false); + + /** + * @brief Train the DMDc model with energy fraction criterion. + * The control matrix B may be available and used in training. + * It is yet to be implemented and requires + * consideration in sparse matrix multiplication in general. + * (See Section III.B. in https://arxiv.org/pdf/1409.6358.pdf) + * In default, the control matrix is unknown, with B = NULL. + * + * @param[in] energy_fraction The energy fraction to keep after doing SVD. + * @param[in] B (Optional) The control matrix B. + */ + virtual void train(double energy_fraction, const Matrix* B = NULL); + + /** + * @brief Train the DMDc model with specified reduced dimension. + * The control matrix B may be available and used in training. + * It is yet to be implemented and requires + * consideration in sparse matrix multiplication in general. + * (See Section III.B. in https://arxiv.org/pdf/1409.6358.pdf) + * In default, the control matrix is unknown, with B = NULL. + * + * @param[in] k The number of modes to keep after doing SVD. + * @param[in] B (Optional) The control matrix B. + */ + virtual void train(int k, const Matrix* B = NULL); + + /** + * @brief Project U using d_phi, where U is the initial condition and the controls. + * Calculate pinv(phi) x U, or more precisely, + * (phi* x phi)^{-1} x phi* x U, where phi* is the conjugate transpose. + * + * @param[in] init The initial condition. + * @param[in] controls The controls. + * @param[in] t_offset The initial time offset. + */ + void project(const Vector* init, const Matrix* controls, + double t_offset = -1.0); + + /** + * @brief Predict state given a time. Uses the projected initial condition of the + * training dataset (the first column). + * + * @param[in] t The time of the output state + */ + Vector* predict(double t); + + /** + * @brief Get the time offset contained within d_t_offset. + */ + double getTimeOffset() const; + + /** + * @brief Returns the number of samples taken. + * + * @return The number of samples taken. + */ + int getNumSamples() const + { + return d_snapshots.size(); + } + + int getDimension() const + { + return d_k; + } + + /** + * @brief Get the snapshot matrix contained within d_snapshots. + */ + const Matrix* getSnapshotMatrix(); + + /** + * @brief Load the object state from a file. + * + * @param[in] base_file_name The base part of the filename to load the + * database from. + */ + virtual void load(std::string base_file_name); + + /** + * @brief Load the object state from a file. + * + * @param[in] base_file_name The base part of the filename to load the + * database from. + */ + void load(const char* base_file_name); + + /** + * @brief Save the object state to a file. + * + * @param[in] base_file_name The base part of the filename to save the + * database to. + */ + virtual void save(std::string base_file_name); + + /** + * @brief Save the object state to a file. + * + * @param[in] base_file_name The base part of the filename to save the + * database to. + */ + void save(const char* base_file_name); + + /** + * @brief Output the DMDc record in CSV files. + */ + void summary(std::string base_file_name); + +protected: + /** + * @brief Constructor. Variant of DMDc with non-uniform time step size. + * + * @param[in] dim The full-order state dimension. + * @param[in] dim_c The control dimension. + * @param[in] state_offset The state offset. + */ + DMDc(int dim, int dim_c, Vector* state_offset = NULL); + + /** + * @brief Constructor. Specified from DMDc components. + * + * @param[in] eigs d_eigs + * @param[in] phi_real d_phi_real + * @param[in] phi_imaginary d_phi_imaginary + * @param[in] B_tilde d_B_tilde + * @param[in] k d_k + * @param[in] dt d_dt + * @param[in] t_offset d_t_offset + * @param[in] state_offset d_state_offset + */ + DMDc(std::vector> eigs, Matrix* phi_real, + Matrix* phi_imaginary, Matrix* B_tilde, int k, + double dt, double t_offset, Vector* state_offset); + + /** + * @brief Unimplemented default constructor. + */ + DMDc(); + + /** + * @brief Unimplemented copy constructor. + */ + DMDc(const DMDc& other); + + /** + * @brief Unimplemented assignment operator. + */ + DMDc& + operator = ( + const DMDc& rhs); + + /** + * @brief Internal function to multiply d_phi with the eigenvalues. + */ + std::pair phiMultEigs(double t); + + /** + * @brief Construct the DMDc object. + */ + void constructDMDc(const Matrix* f_snapshots, + const Matrix* f_controls, + int rank, + int num_procs, + const Matrix* B); + + /** + * @brief Returns a pair of pointers to the minus and plus snapshot matrices + */ + virtual std::pair computeDMDcSnapshotPair( + const Matrix* snapshots, const Matrix* controls, const Matrix* B); + + /** + * @brief Compute the appropriate exponential function when predicting the solution. + */ + virtual std::complex computeEigExp(std::complex eig, double t); + + /** + * @brief Add the state offset when predicting the solution. + */ + virtual void addOffset(Vector*& result); + + /** + * @brief Get the snapshot matrix contained within d_snapshots. + */ + const Matrix* createSnapshotMatrix(std::vector snapshots); + + /** + * @brief The rank of the process this object belongs to. + */ + int d_rank; + + /** + * @brief The number of processors being run on. + */ + int d_num_procs; + + /** + * @brief The total dimension of the sample vector. + */ + int d_dim; + + /** + * @brief The total dimension of the control vector. + */ + int d_dim_c; + + /** + * @brief The time step size between samples. + */ + double d_dt = -1.0; + + /** + * @brief The time offset of the first sample. + */ + double d_t_offset; + + /** + * @brief std::vector holding the snapshots. + */ + std::vector d_snapshots; + + /** + * @brief std::vector holding the controls. + */ + std::vector d_controls; + + /** + * @brief The stored times of each sample. + */ + std::vector d_sampled_times; + + /** + * @brief State offset in snapshot. + */ + Vector* d_state_offset = NULL; + + /** + * @brief Whether the DMDc has been trained or not. + */ + bool d_trained; + + /** + * @brief Whether the initial condition has been projected. + */ + bool d_init_projected; + + /** + * @brief The maximum number of singular vectors. + */ + int d_num_singular_vectors; + + /** + * @brief std::vector holding the signular values. + */ + std::vector d_sv; + + /** + * @brief The energy fraction used to obtain the DMDc modes. + */ + double d_energy_fraction; + + /** + * @brief The number of columns used after obtaining the SVD decomposition. + */ + int d_k; + + /** + * @brief The left singular vector basis. + */ + Matrix* d_basis = NULL; + + /** + * @brief A_tilde + */ + Matrix* d_A_tilde = NULL; + + /** + * @brief B_tilde + */ + Matrix* d_B_tilde = NULL; + + /** + * @brief The real part of d_phi. + */ + Matrix* d_phi_real = NULL; + + /** + * @brief The imaginary part of d_phi. + */ + Matrix* d_phi_imaginary = NULL; + + /** + * @brief The real part of d_phi_squared_inverse. + */ + Matrix* d_phi_real_squared_inverse = NULL; + + /** + * @brief The imaginary part of d_phi_squared_inverse. + */ + Matrix* d_phi_imaginary_squared_inverse = NULL; + + /** + * @brief The real part of the projected initial condition. + */ + Vector* d_projected_init_real = NULL; + + /** + * @brief The imaginary part of the projected initial condition. + */ + Vector* d_projected_init_imaginary = NULL; + + /** + * @brief The real part of the projected controls. + */ + Matrix* d_projected_controls_real = NULL; + + /** + * @brief The imaginary part of the projected controls. + */ + Matrix* d_projected_controls_imaginary = NULL; + + /** + * @brief A vector holding the complex eigenvalues of the eigenmodes. + */ + std::vector> d_eigs; + +}; + +} + +#endif diff --git a/unit_tests/matlab/dmdc.m b/unit_tests/matlab/dmdc.m new file mode 100644 index 000000000..dd4fb513e --- /dev/null +++ b/unit_tests/matlab/dmdc.m @@ -0,0 +1,52 @@ +function dmdc(X, Y, ef, t, dt, varargin) +% Input: X: n by m snapshot matrix +% Y: n by (m-1) control snapshot matrix +% ef: energy fraction +% t: The time to predict +% dt: delta time +% B: (optional) known control matrix + + X_in = X(:,1:end-1); + X_out = X(:,2:end); + b0 = X(:,1); + + if nargin == 5 + X_in = [X_in; Y]; + elseif nargin == 6 + B = varargin{1}; + X_out = X_out - B*Y; + else + error('Invalid input'); + end + + [U, S, V] = svd(X_in, 'econ'); + r = sum(cumsum(diag(S))/sum(diag(S)) < ef) + + U = U(:,1:r); + S = S(1:r,1:r); + V = V(:,1:r); + + m = size(X,1); + U1 = U(1:m,:); + if nargin == 5 + U2 = U(m+1:end,:); + [U_out, S_out, ~] = svd(X_out, 'econ'); + r_out = sum(cumsum(diag(S_out))/sum(diag(S_out)) < ef) + U_out = U_out(:,1:r_out); + Atilde = U_out'*X_out*V*inv(S)*U1'*U_out; + Btilde = U_out'*X_out*V*inv(S)*U2'; + [W,ev] = eig(Atilde); + %Phi = X_out*V*inv(S)*U1'*U_out*W; + else + U_out = U1; + Atilde = U_out'*X_out*V*inv(S); + Btilde = U_out'*B; + [W,ev] = eig(Atilde); + %Phi = X_out*V*inv(S)*W; + end + + n = size(X_out,2); + pred = U_out*W*sum((diag(ev).^(n:-1:0)).*(W\[U_out'*b0, Btilde*Y]),2); + + norm(real(pred-X(:,end)))/norm(real(X(:,end))) +end