Untitled

//#include "/home/lipegal/Documents/OIST/code/Modular/kp_lib.h"
#include "/flash/BuschU/andres/Transfer_protocol/kp_lib.h"
#include <vector>
#include <functional>
#include <iostream>
#include <fstream>
#include <cmath>
#include <boost/numeric/odeint.hpp>
using namespace boost::numeric::odeint;
using namespace std;

int in_trunc_dim = 1600;
double d_tol = 1e-7;

using HamKernel = std::function<double(int, int)>;

// ------------------------------------------------------------------
// Helpers
// ------------------------------------------------------------------
// Diagonalize hamiltonian with given parameters and store data
Ham_data H_diag_spec (const std::function<double(int, int)>& H, int trunc_d_init, const double tol_overlap, const int num_comp) {
    int trunc_d = trunc_d_init, last_trunc_d = 0;
    std::vector<double> converged_eigvals, current_eigvals;
    std::vector<double> last_eigvecs;
    std::vector<double> A; // Flattened hamiltonian matrix
    bool converged = false;

    // Check for convergence in diagonalization
    while (!converged) {
        A.resize (trunc_d * trunc_d);
        current_eigvals.resize(trunc_d);
        // Flatten Hamiltonian matrix
        for(int n = 0; n < trunc_d; n++) {
            for(int m = 0; m < trunc_d; m++) {
                A [n * trunc_d + m] = H (n + 1,m + 1);
                // There is a distintion between physical (n,m) and array index
            }
        }

        // Define working variables
        int n = trunc_d, lda = trunc_d, info, lwork = -1;
        double work_test;
        char jobz = 'V', uplo = 'U';

        // Diagonalization (is not being done twice, first just assigns workspace sizes)
        dsyev_(&jobz, &uplo, &n, A.data(), &lda, current_eigvals.data(), &work_test, &lwork, &info);
        lwork = static_cast<int>(work_test);
        std::vector<double> work_opt(lwork);
        dsyev_(&jobz, &uplo, &n, A.data(), &lda, current_eigvals.data(), work_opt.data(), &lwork, &info);

        //if (info != 0) {
          //  throw std::runtime_error("Diagonalization failed : " + std::to_string(info));
        //}

        // Convergence check (similar to old version)
        if (last_eigvecs.empty()) {
            last_eigvecs = A;
            last_trunc_d = trunc_d;
            trunc_d += 20;
            continue; // Skip convergence check on first iteration
        }

        converged = true;
        int max_check = std::min(num_comp, last_trunc_d);
        for (int i = 0; i < max_check; i++) {
            double overlap = 0.0;

            // Eigenvectors are column-wise in LAPACK
            for (int k = 0; k < last_trunc_d; k++) {
                overlap += last_eigvecs[k + i * (last_trunc_d)]
                    * A[k + i * trunc_d];
            }
            overlap = std::abs(overlap);

            //std::cout << std::setprecision(16) << std::scientific;
            //std::cout << "1 - ovlp = " << 1 - overlap << endl;

            if (1.0 - overlap > tol_overlap) {
                converged = false;
                break;
            }
        }

        if (!converged) {
            last_eigvecs = A;
            last_trunc_d = trunc_d;
            trunc_d += 500;
        }
    }

    // Convert to output format
    RealVector eigenvals(trunc_d);
    for(int i = 0; i < trunc_d; i++) {
        eigenvals(i) = current_eigvals[i];
    }

    ComplexVector eigenvecs(trunc_d * trunc_d);
    for(int i = 0; i < trunc_d * trunc_d; i++) {
        eigenvecs(i) = A[i];
    }

    return {eigenvals, eigenvecs, trunc_d};
}

// Derivative of Hamiltonian
double Ham_dd(int n, int m, double d_2,
              const std::vector<Domain>& domains, const RealVector& offsets, double L_total) {
    double sum = 0.;
    double k = PI / L_total;
    const Domain &d = domains[1]; // second domain
    for (int i = 1; i <= d.M; i++) { // BARRIER POSITION GOES FROM 1 TO M
        double ym = y_(1, i, d_2, domains, offsets, L_total) - xc + ( L_total / 2.) ;
        sum += k * (n + m) * sin( k * (n + m) * ym) - k * (n - m) * sin( k * (n - m) * ym);
    }
    return sum * ( (d.h * d.L ) / (2 * L_total * d.M) )  ;
}

double Ham_dd_num (int n, int m, double dd,
            const vector<Domain>& domains, const RealVector& offsets, double L) {
    
    vector<Domain> domains_p = domains;
    vector<Domain> domains_m = domains;

    // Only vary the height parameter
    domains_p[1].D += dd;
    domains_m[1].D -= dd;

    double Lp = syst_L(domains_p);
    double Lm = syst_L(domains_m);

    RealVector offsets_p = offsets_(domains_p);
    RealVector offsets_m = offsets_(domains_m);

    double Hp = Ham_uni(n, m, domains_p, offsets_p, Lp);
    double Hm = Ham_uni(n, m, domains_m, offsets_m, Lm);

    return (Hp - Hm) / (2.0 * dd);
}
// ------------------------------------------------------------------
// FAQUAD single-Delta computation
// ------------------------------------------------------------------
void d_FAQUAD(const RealVector& d_vals, int d_idx, double d_sep, double h,
              int idx_1, int idx_2, ofstream& out) {
    double d = d_vals(d_idx);

    double E1 = 0., E2 = 0.;
    Complex M_an = 0., M_num = 0.;

    int doms_num = 2;
    std::vector<Domain> domains;
    double L_[] = {20., 20.};
    int    M_[] = {20, 20};
    double B_[] = {10., 0.};
    double W_[] = {1., 0.};
    double D_[] = {0., d};
    double h_[] = {h, h};

    dom_load(domains, L_, M_, B_, W_, D_, h_, doms_num);
    double L = syst_L(domains);
    RealVector offsets = offsets_(domains);

    // Construct Hamiltonian
    auto H_in = [&](int n, int m) {
        return Ham_uni(n, m, domains, offsets, L);
    };

    auto H_dh_fun_an = [&](int n, int m) {
        return Ham_dd(n, m, d, domains, offsets, L);
    };

    auto H_dh_fun_num = [&](int n, int m) {
        return Ham_dd_num(n, m, 1e-5, domains, offsets, L);
    };

    Ham_data h_data = H_diag_spec(H_in, 4000, 1e-06, 6000);

    E1 = h_data.eigvals(idx_1);
    E2 = h_data.eigvals(idx_2);

    //static ComplexVector v1_prev, v2_prev;
    //static bool first = true;

    ComplexVector v1 = eigen_extract(h_data, idx_1);
    ComplexVector v2 = eigen_extract(h_data, idx_2);

    // Compute derivative matrix element M
    M_an = 0.0;
    M_num = 0.0;
    int N = v1.size();
    for (int n = 0; n < N; ++n) {
        for (int m = 0; m < N; ++m) {
            M_an += std::conj(v1(n)) * H_dh_fun_an(n, m) * v2(m);
            M_num += std::conj(v1(n)) * H_dh_fun_num(n, m) * v2(m);
        }
    }

    out << d << " " << E1 << " " << E2 << " "
        << real(M_an) << " " << imag(M_an) << " "
        << real(M_num) << " " << imag(M_num) << endl;
}


void d_FAQUAD_conv(const RealVector& d_vals, int d_idx, double d_sep, double h,
              int idx_1, int idx_2, ofstream& out) {
    double d = d_vals(d_idx);

    int num_comp_local = 2000;
    int step = 500;
    double eps_abs = 1e-10;
    double eps_rel = 1e-6;

    Complex M_prev = 0.0;
    Complex M_curr = 0.0;
    double E1 = 0., E2 = 0.;
    Complex M_an = 0., M_num = 0.;

    int doms_num = 2;
    std::vector<Domain> domains;
    double L_[] = {20., 20.};
    int    M_[] = {20, 20};
    double B_[] = {10., 0.};
    double W_[] = {1., 0.};
    double D_[] = {0., d};
    double h_[] = {h, h};

    dom_load(domains, L_, M_, B_, W_, D_, h_, doms_num);
    double L = syst_L(domains);
    RealVector offsets = offsets_(domains);

    // Construct Hamiltonian
    auto H_in = [&](int n, int m) {
        return Ham_uni(n, m, domains, offsets, L);
    };

    auto H_dh_fun_an = [&](int n, int m) {
        return Ham_dd(n, m, d, domains, offsets, L);
    };

    auto H_dh_fun_num = [&](int n, int m) {
        return Ham_dd_num(n, m, 1e-5, domains, offsets, L);
    };

    bool first = true;
    while (true) {
        Ham_data h_data = H_diag_spec(H_in, in_trunc_dim, d_tol, num_comp_local);

        E1 = h_data.eigvals(idx_1);
        E2 = h_data.eigvals(idx_2);

        //static ComplexVector v1_prev, v2_prev;
        //static bool first = true;

        ComplexVector v1 = eigen_extract(h_data, idx_1);
        ComplexVector v2 = eigen_extract(h_data, idx_2);

        // Compute derivative matrix element M
        M_an = 0.0;
        M_num = 0.0;
        int N = v1.size();
        for (int n = 0; n < N; ++n) {
            for (int m = 0; m < N; ++m) {
                M_an += std::conj(v1(n)) * H_dh_fun_an(n, m) * v2(m);
                M_num += std::conj(v1(n)) * H_dh_fun_num(n, m) * v2(m);
            }
        }

        cout << std::setprecision(16) << "M = " << real(M_an) << " - " << real(M_num) << endl;

        if (!first) {
            double diff = abs(M_an - M_prev);
            double scale = max(abs(M_an), abs(M_prev));

            if (diff < eps_abs || diff / scale < eps_rel) {
                M_curr = M_an;
                break;
            }
        }

        M_prev = M_an;
        first = false;
        num_comp_local += step;

        if (num_comp_local > 7000) {
            cout << "WARNING: did not converge up to max num_comp\n";
            M_curr = M_an;
            break;
        }
    }

    out << d << " " << E1 << " " << E2 << " "
        << real(M_curr) << " " << imag(M_curr) << " "
        << real(M_num) << " " << imag(M_num) << endl;
}

// ------------------------------------------------------------------
// Main
// ------------------------------------------------------------------
int main(int argc, char* argv[]) {
    if (argc < 5) {
        cerr << "Usage: ./run [d_in][d_f][d_sep][d_idx][JOB_ID]" << endl;
        return 1;
    }

    double d_sep = stod(argv[3]);
    int d_idx = stoi(argv[4]);
    int job_id = stoi(argv[5]);
    // -------------------- System setup ------------------------------
    double d_i = stod(argv[1]), d_f = stod(argv[2]), h = 1.6;
    int N_steps = static_cast<int>((d_f - d_i) / d_sep) + 1;
    RealVector d_vals(N_steps);
    for (int i = 0; i < N_steps; ++i)
        d_vals(i) = d_i + i * d_sep;

    // Output file for this specific index
    string fname = "data/gen_data_" + to_string(job_id) + ".dat";
    // can change job_id -> d_idx if single d per job and ios::trunc -> ios::app
    ofstream gen_File(fname, ios::app);

    //for(int d_idx_i = 0; d_idx_i <= N_steps; d_idx_i++) {
    d_FAQUAD(d_vals, d_idx, d_sep, h, 38, 39, gen_File);
    cout << " completed d = " << d_idx << endl;
    //}
    //cout << "FLAG: completed d_idx = " << d_idx << endl;
    return 0;
}