/*! Type II quasisolitons in Lieb-Liniger: amplitudes_IDeltax

  Purpose:
  Produce the .amplitudes file
  containing the (complex) amplitudes of each basis state
  for amplitudes optimizing depletion density integral I(Δx)
  for (multi-)type II hole wavepackets in Lieb-Liniger.

  This executable requires Abacus version 2.

  See README for compilation instructions.

  Copyright © Jean-Sébastien Caux, Anahita Sarvi and Cesare Vianello.

  This program is free software: you can redistribute it and/or modify it under the terms of the GNU Affero General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

  This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Affero General Public License for more details.

  You should have received a copy of the GNU Affero General Public License along with this program. If not, see <https://www.gnu.org/licenses/>.
*/


import std;

import abacus;


int main(int argc, char* argv[])
{
  using namespace std::complex_literals;

  if (argc != 8) {
    std::cout << "Executable amplitudes_IDeltax\n"
	      << " part of the Type II quasisolitons in Lieb-Liniger software suite\n"
	      << " copyright © Jean-Sébastien Caux, Anahita Sarvi and Cesare Vianello.\n";
    std::cout << "\nPurpose:\n"
	      << " Produce the .amplitudes file\n"
	      << " containing the (complex) amplitudes of each basis state\n"
	      << " for amplitudes optimizing depletion density integral I(Δx)\n"
	      << " for (multi-)type II hole wavepackets in Lieb-Liniger.\n";

    std::cout << "\nPrerequirements:\n----------------\n"
	      << " - .states and .rho files produced by executable rho\n";

    std::cout << "\nUsage:\n------\n";
    std::cout << "holeIx2 <c> <L> <N> <nr holes> <width> <offset>\n\n";

    int warg { 16 }, wtype { 10 }, wcons { 26 };
    std::cout << std::left << std::setw(warg) << "Argument" << std::setw(wtype) << "Type"
	      << std::setw(wcons) << "Constraints" << "Description\n";
    std::cout << std::left << std::setw(warg) << "--------" << std::setw(wtype) << "----"
	      << std::setw(wcons) << "-----------" << "-----------\n";
    std::cout << std::left << std::setw(warg) << "c" << std::setw(wtype) << "Real"
	      << std::setw(wcons) << "> 0"
	      << "Value of the interaction parameter\n";
    std::cout << std::left << std::setw(warg) << "L" << std::setw(wtype) << "Real"
	      << std::setw(wcons) << "> 0"
	      << "System size\n";
    std::cout << std::left << std::setw(warg) << "N" << std::setw(wtype) << "int"
	      << std::setw(wcons) << "> 0"
	      << "Number of particles\n";
    std::cout << std::left << std::setw(warg) << "nr holes" << std::setw(wtype) << "int"
	      << std::setw(wcons) << "1 <= nr holes <= N"
	      << "Number of holes (Type II modes)\n";
    std::cout << std::left << std::setw(warg) << "width" << std::setw(wtype) << "int"
	      << std::setw(wcons) << "nr holes < width <= N"
	      << "Width of the hole window\n";
    std::cout << std::left << std::setw(warg) << "offset" << std::setw(wtype) << "int"
	      << std::setw(wcons) << "0 <= offset <= N-width"
	      << "Offset of the hole window w/r to the right Fermi edge\n";
    std::cout << std::left << std::setw(warg) << "Δx" << std::setw(wtype) << "Real"
	      << std::setw(wcons) << "0 < Δx <= L"
	      << "Spatial width of depletion density to optimize\n";

    return 0;
  }

  std::cout << std::setprecision(std::numeric_limits<Real>::digits10 + 1);

  Real c { std::stold(argv[1]) };
  Real L { std::stold(argv[2]) };
  int N { std::stoi(argv[3]) };
  int nholes { std::stoi(argv[4]) };
  int width { std::stoi(argv[5]) };
  int offset { std::stoi(argv[6]) };
  Real Δx {std::stold(argv[7]) };


  // Define the input files

  std::stringstream filename_base;
  filename_base << "c_" << c << "_N_" << N << "_L_" << L << "_nholes_" << nholes;
  filename_base << "_width_" << width;
  filename_base << "_offset_" << offset;

  std::stringstream states_filename;
  states_filename << filename_base.str() << ".states";
  std::ifstream states_file;
  states_file.open(states_filename.str());
  states_file >> std::setprecision(std::numeric_limits<Real>::digits10 + 1);

  std::stringstream rho_ME_filename;
  rho_ME_filename << filename_base.str() << ".rho";
  std::ifstream rho_ME_file;
  rho_ME_file.open(rho_ME_filename.str());
  rho_ME_file >> std::setprecision(std::numeric_limits<Real>::digits10 + 1);


  // Input the states info

  std::vector<std::string> label;
  std::vector<int> iK;
  std::vector<Real> E;
  std::string tmp_label;
  int tmp_iK;
  Real tmp_E;
  states_file >> tmp_label;
  do {
    states_file >> tmp_iK >> tmp_E;
    label.push_back(tmp_label);
    iK.push_back(tmp_iK);
    E.push_back(tmp_E);
  } while (states_file >> tmp_label);
  states_file.close();

  int nr_states { int(label.size()) };


  // Input the density operator matrix elements

  std::map<std::string, std::map<std::string, std::complex<Real>>> rho_ME;

  for (int ibra { 0 }; ibra < nr_states; ++ibra) {
    for (int iket { 0 }; iket <= ibra; ++iket) {
      rho_ME_file >> rho_ME[label[ibra]][label[iket]];
    }
  }
  rho_ME_file.close();



  // Define the M matrix for the given Δx

  Matrix<Real> M (nr_states);
  Real K;
  M.setZero();

  for (int ibra { 0 }; ibra < nr_states; ++ibra) {
    for (int iket { 0 }; iket < ibra; ++iket) {
      K = pi_r * (iK[ibra] - iK[iket]) * Δx/L; // K_{ab}*Δx/2
      // To balance the matrix numerically, it is multiplied by sqrt(N)
      M(ibra, iket) = std::sqrt(N) *
	std::real(rho_ME[label[ibra]][label[iket]]) * (K != Real(0) ? std::sin(K)/K : 1);
      M(iket, ibra) = M(ibra, iket); // matrix is symmetric
    }
  }


  // Householder reduction and QL factorization

  std::vector<Real> d(nr_states);
  std::vector<Real> e(nr_states);
  tred2(M, d, e);
  tqli(d, e, M);
  // std::cout << "Eigenvalues according to QL: " << d << std::endl;

  // Find minimal eigenvalue; the optimal wavefunction is then the corresponding eigenvector
  int min_red_ev_index { 0 };
  Real min_red_ev { d[0] };
  for (int ibra { 0 }; ibra < nr_states; ++ibra) {
    if (d[ibra] < min_red_ev) {
      min_red_ev_index = ibra;
      min_red_ev = d[ibra];
    }
  }


  // Define the output files

  std::stringstream amplitudes_filename;
  amplitudes_filename << filename_base.str() << "_IDeltax_" << Δx << ".amplitudes";
  std::ofstream amplitudes_file;
  amplitudes_file.open(amplitudes_filename.str(), std::ios::out | std::ios::trunc);
  amplitudes_file << std::setprecision(std::numeric_limits<Real>::digits10 + 1);


  // Output the state amplitudes

  for (int ibra { 0 }; ibra < nr_states; ++ibra) {
    amplitudes_file << M(ibra, min_red_ev_index) << "\t";
  }
  amplitudes_file.close();


  return 0;
}