.. _program_listing_file_include_Suzuki-Trotter-Evolver_UnitaryEvolver.hpp:

Program Listing for File UnitaryEvolver.hpp
===========================================

|exhale_lsh| :ref:`Return to documentation for file <file_include_Suzuki-Trotter-Evolver_UnitaryEvolver.hpp>` (``include/Suzuki-Trotter-Evolver/UnitaryEvolver.hpp``)

.. |exhale_lsh| unicode:: U+021B0 .. UPWARDS ARROW WITH TIP LEFTWARDS

.. code-block:: cpp

   
   #include <complex>
   #include <cmath>
   #include <vector>
   #include <Eigen/Eigenvalues>
   #include <Eigen/SparseCore>
   using std::complex;
   using std::vector;
   using Eigen::Dynamic;
   using std::pow;
   using std::conj;
   
   namespace Suzuki_Trotter_Evolver {
   
   template<int n = Dynamic, int m = Dynamic> using DMatrix =
       Eigen::Matrix<complex<double>, n, m>;
   
   typedef Eigen::SparseMatrix<complex<double>> SMatrix;
   
   template<int n_ctrl = Dynamic,
            int dim = Dynamic,
            typename Matrix = DMatrix<dim, dim>>
   struct UnitaryEvolver {
       private: const complex<double> minus_i = complex<double>(0.0, -1.0);
       public:
   
       static const int dim_x_n_ctrl = (dim == Dynamic || n_ctrl == Dynamic)
                          ? Dynamic
                          : dim * n_ctrl;
       size_t length;
   
       Eigen::Array<complex<double>, dim, 1> d0;
   
       vector<Eigen::Array<complex<double>, dim, 1>> ds;
   
       Matrix u0;
   
       Matrix u0_inverse;
   
       vector<Matrix> us;
   
       vector<Matrix> us_individual;
   
       vector<Matrix> us_inverse_individual;
   
       vector<Matrix> hs;
   
       Matrix u0_inverse_u_last;
       
       template<typename T = Matrix>
       UnitaryEvolver(std::enable_if_t<std::is_same<T, DMatrix<dim, dim>>::value,
                      DMatrix<dim, dim>> drift_hamiltonian,
                      DMatrix<dim_x_n_ctrl, dim> control_hamiltonians
                     ) {
           typedef Eigen::SelfAdjointEigenSolver<DMatrix<dim, dim>> EigSolver;
           EigSolver eigs0(drift_hamiltonian);
           d0 = minus_i*eigs0.eigenvalues().array();
           u0 = eigs0.eigenvectors();
           u0_inverse = u0.adjoint();
   
           Eigen::Index l = drift_hamiltonian.rows();
           length = control_hamiltonians.rows()/control_hamiltonians.cols();
           
           if (length == 0) {
               u0_inverse_u_last = DMatrix<dim, dim>::Identity(l, l);
               us.push_back(u0);
           } else {
               hs.push_back((DMatrix<dim, dim>)(control_hamiltonians.block(0, 0, l, l)));
               EigSolver eigs(hs.back());
               ds.push_back((Eigen::Array<complex<double>, dim, 1>)(
                   minus_i*eigs.eigenvalues().array()));
               DMatrix<dim, dim> u = eigs.eigenvectors();
               us_individual.push_back(u);
               us_inverse_individual.push_back(u.adjoint());
               us.push_back(u.adjoint() * u0);
               us.push_back(u); // append twice
               for (size_t i = 1; i < length; i++) {
                   hs.push_back((DMatrix<dim, dim>)(control_hamiltonians.block(i*l, 0, l, l)));
                   eigs = EigSolver(hs.back());
                   ds.push_back((Eigen::Array<complex<double>, dim, 1>)(
                       minus_i*eigs.eigenvalues().array()));
                   u = eigs.eigenvectors();
                   us_individual.push_back(u);
                   us_inverse_individual.push_back(u.adjoint());
                   us.back() = u.adjoint() * us.back();
                   us.push_back(u);
               }
               u0_inverse_u_last.noalias() = u0_inverse * u;
           }
       };
   
       // Sparse initialisation
       template<typename T = Matrix>
       UnitaryEvolver(std::enable_if_t<std::is_same<T, SMatrix>::value,
                      DMatrix<dim, dim>> drift_hamiltonian,
                      DMatrix<dim_x_n_ctrl, dim> control_hamiltonians) {
           typedef Eigen::SelfAdjointEigenSolver<DMatrix<dim, dim>> EigSolver;
           EigSolver eigs0(drift_hamiltonian);
           d0 = minus_i*eigs0.eigenvalues().array();
           u0 = eigs0.eigenvectors().sparseView();
           u0_inverse = u0.adjoint();
   
           Eigen::Index l = drift_hamiltonian.rows();
           length = control_hamiltonians.rows()/control_hamiltonians.cols();
           if (length == 0) {
               u0_inverse_u_last = DMatrix<dim, dim>::Identity(l, l).sparseView();
               us.push_back(u0);
           } else {
               hs.push_back(control_hamiltonians.block(0, 0, l, l).sparseView());
               EigSolver eigs((DMatrix<dim, dim>)(hs.back()));
               ds.push_back((Eigen::Array<complex<double>, dim, 1>)(
                                              minus_i*eigs.eigenvalues().array()));
               Matrix u = eigs.eigenvectors().sparseView();
               us_individual.push_back(u);
               us_inverse_individual.push_back(u.adjoint());
               us.push_back((u.adjoint() * u0).pruned());
               us.push_back(u); // append twice
               for (size_t i = 1; i < length; i++) {
                   hs.push_back(control_hamiltonians.block(i*l, 0, l, l).sparseView());
                   eigs = EigSolver((DMatrix<dim, dim>)(hs.back()));
                   ds.push_back((Eigen::Array<complex<double>, dim, 1>)(
                                              minus_i*eigs.eigenvalues().array()));
                   u = eigs.eigenvectors().sparseView();
                   us_individual.push_back(u);
                   us_inverse_individual.push_back(u.adjoint());
                   us.back() = (u.adjoint() * us.back()).pruned();
                   us.push_back(u);
               }
               u0_inverse_u_last = u0_inverse * u;
           }
       };
   
       UnitaryEvolver(size_t l,
                      Eigen::Array<complex<double>, dim, 1> d0,
                      vector<Eigen::Array<complex<double>, dim, 1>> ds,
                      Matrix u0,
                      Matrix u0_inverse,
                      vector<Matrix> us,
                      vector<Matrix> us_individual,
                      vector<Matrix> us_inverse_individual,
                      vector<Matrix> hs,
                      Matrix u0_inverse_u_last)
           : length{l},
             d0{d0},
             ds{ds},
             u0{u0},
             u0_inverse{u0_inverse},
             us{us},
             us_individual{us_individual},
             us_inverse_individual{us_inverse_individual},
             hs{hs},
             u0_inverse_u_last{u0_inverse_u_last}
           {};
           
       DMatrix<dim, 1> propagate(DMatrix<Dynamic, n_ctrl> ctrl_amp,
                                 DMatrix<dim, 1> state,
                                 double dt) {
           size_t steps = ctrl_amp.rows();
           DMatrix<dim, 1> exp_d0 = (d0*dt).exp().matrix();
   
           ctrl_amp *= dt;
           state = u0_inverse * state;
           state = state.cwiseProduct(exp_d0);
           for (size_t j = 0; j < length; j++) {
               state = us[j] * state;
               state = state.cwiseProduct((ds[j]*ctrl_amp(0, j)).exp().matrix());
           }
           for (size_t i = 1; i < steps; i++) {
               state = u0_inverse_u_last * state;
               state = state.cwiseProduct(exp_d0);
               for (size_t j = 0; j < length; j++) {
                   state = us[j] * state;
                   state = state.cwiseProduct(
                       (ds[j]*ctrl_amp(i, j)).exp().matrix());
               }
           }
           state = us[length] * state;
           return state;
       };
   
       template<int l = Dynamic>
       DMatrix<dim, l> propagate_collection(DMatrix<Dynamic, n_ctrl> ctrl_amp,
                                            DMatrix<dim, l> states,
                                            double dt) {
           size_t steps = ctrl_amp.rows();
           Eigen::Array<complex<double>, dim, 1> exp_d0 = (d0*dt).exp();
   
           ctrl_amp *= dt;
           states = u0_inverse*states;
           states = states.array().colwise() * exp_d0;
           for (size_t j = 0; j < length; j++) {
               states = us[j] * states;
               states = states.array().colwise() * (ds[j]*ctrl_amp(0, j)).exp();
           }
           for (size_t i = 1; i < steps; i++) {
               states = u0_inverse_u_last * states;
               states = states.array().colwise() * exp_d0;
               for (size_t j = 0; j < length; j++) {
                   states = us[j] * states;
                   states = states.array().colwise()
                           *(ds[j]*ctrl_amp(i, j)).exp();
               }
           }
           states = us[length] * states;
           return states;
       };
   
       DMatrix<dim, Dynamic> propagate_all(DMatrix<Dynamic, n_ctrl> ctrl_amp,
                                         DMatrix<dim, 1> state,
                                         double dt) {
           size_t steps = ctrl_amp.rows();
           DMatrix<dim, 1> exp_d0 = (d0*dt).exp().matrix();
           DMatrix<dim, Dynamic> phi(state.rows(), steps+1);
           ctrl_amp *= dt;
           phi.col(0) = state;
           phi.col(1) = u0_inverse * state;
           phi.col(1) = phi.col(1).cwiseProduct(exp_d0);
           for (size_t j = 0; j < length; j++) {
               phi.col(1) = us[j] * phi.col(1);
               phi.col(1) = phi.col(1).cwiseProduct(
                   (ds[j]*ctrl_amp(0, j)).exp().matrix());
           }
           phi.col(1) = us[length] * phi.col(1);
           for (size_t i = 1, k = 2; i < steps; i++, k++) {
               phi.col(k) = u0_inverse * phi.col(i);
               phi.col(k) = phi.col(k).cwiseProduct(exp_d0);
               for (size_t j = 0; j < length; j++) {
                   phi.col(k) = us[j] * phi.col(k);
                   phi.col(k) = phi.col(k).cwiseProduct(
                       (ds[j]*ctrl_amp(i, j)).exp().matrix());
               }
               phi.col(k) = us[length] * phi.col(k);
           }
           return phi;
       };
   
       complex<double> evolved_expectation_value(DMatrix<Dynamic, n_ctrl> ctrl_amp,
                                                 DMatrix<dim, 1> state,
                                                 double dt,
                                                 DMatrix<dim, dim> observable) {
           state = propagate(ctrl_amp, state, dt);
           return state.dot(observable * state); // dot conjugates and transposes
       };
   
       DMatrix<Dynamic, 1> evolved_expectation_value_all(
                                                 DMatrix<Dynamic, n_ctrl> ctrl_amp,
                                                 DMatrix<dim, 1> state,
                                                 double dt,
                                                 DMatrix<dim, dim> observable) {
           DMatrix<dim, Dynamic> phi = propagate_all(ctrl_amp, state, dt);
           return (phi.adjoint() * observable * phi).diagonal();
       };
   
       std::tuple<complex<double>, Eigen::Matrix<double, Dynamic, n_ctrl>>
       switching_function(DMatrix<Dynamic, n_ctrl> ctrl_amp,
                          DMatrix<dim, 1> state,
                          double dt,
                          DMatrix<dim, dim> cost) {
           // forward propagation
           size_t steps = ctrl_amp.rows();
           DMatrix<dim, 1> exp_d0 = (d0*dt).exp().matrix();
           ctrl_amp *= dt;
           state = u0_inverse * state;
           state = state.cwiseProduct(exp_d0);
           for (size_t j = 0; j < length; j++) {
               state = us[j] * state;
               state = state.cwiseProduct((ds[j]*ctrl_amp(0, j)).exp().matrix());
           }
           for (size_t i = 1; i < steps; i++) {
               state = u0_inverse_u_last * state;
               state = state.cwiseProduct(exp_d0);
               for (size_t j = 0; j < length; j++) {
                   state = us[j] * state;
                   state = state.cwiseProduct(
                                            (ds[j]*ctrl_amp(i, j)).exp().matrix());
               }
           }
           state = us[length] * state;
           // Lambda
           DMatrix<dim, 1> lambda = cost * state;
           // Energy
           complex<double> E = state.dot(lambda);
           // Back propagation
           DMatrix<dim, 1> exp_d0_c = exp_d0.conjugate();
           ctrl_amp *= -1;
           DMatrix<Dynamic, n_ctrl> out = DMatrix<Dynamic, n_ctrl>::Zero(steps,
                                                                         length);
           for (size_t k = steps-1; true;) {
               for (size_t j = length-1; true;) {
                   out(k, j) = lambda.dot(hs[j] * state);
                   state = us_inverse_individual[j] * state;
                   state = state.cwiseProduct(
                                            (ds[j]*ctrl_amp(k, j)).exp().matrix());
                   state = us_individual[j] * state;
                   lambda = us_inverse_individual[j] * lambda;
                   lambda = lambda.cwiseProduct(
                       (ds[j]*ctrl_amp(k, j)).exp().matrix());
                   lambda = us_individual[j] * lambda;
                   if (j == 0) {
                       break;
                   }
                   j--;
               }
               if (k == 0) {
                   break;
               }
               state = u0_inverse * state;
               state = state.cwiseProduct(exp_d0_c);
               state = u0 * state;
               lambda = u0_inverse * lambda;
               lambda = lambda.cwiseProduct(exp_d0_c);
               lambda = u0 * lambda;
               k--;
           }
           return
               std::tuple<complex<double>, Eigen::Matrix<double, Dynamic, n_ctrl>>(
                   E,
                   2*out.imag()
           );
       };
   };
   }