using LinearAlgebra, Printf, PyPlot
Base.show(io::IO, f::Float64) = @printf(io, "%1.2f", f)

function analyze_operator(A)
    S, U = eigen(A); u1, u2 = U[:, 1], U[:, 2]
    println("A="); show(stdout, "text/plain", A);
    println(" "); println(" ")
    println("norm(A-A†)= ", norm(A-A'), ", eigvalues of A: ", S, ", and (u₁, u₂) = ", dot(u1, u2));
    println(" ");
end

A = [-3/2 -1/2; -1/2 -3/2]; analyze_operator(A)
A = [0 1; 1 0]; analyze_operator(A)
A = rand(2, 2) + im*rand(2,2); A = (A+A')/2; # this is to convert the random A to a self-adjoint one
analyze_operator(A)

A = rand(2, 2) + im*rand(2,2); 

S, U = eigen(A); u1 = U[:, 1]; u2 =U[:, 2]
Sadj, V = eigen(copy(A')); v1 = V[:, 1]; v2 =V[:, 2]

println("A="); show(stdout, "text/plain", A);
println(" "); println(" ")
println("eigvalues of A: λ₁=", S[1], " , λ₂=", S[2]); println(" "); 

println("A†="); show(stdout, "text/plain", A');
println(" "); println(" ")
println("eigvalues of A†: μ₁", Sadj[1], " , μ₂=", Sadj[2]); println(" ");

println("inner products between eigenvectors of A and of A†");
println("(u₁, v₁) = ", dot(u1, v1), ", (u₁, v₂) = ", dot(u1, v2))
println("(u₂, v₁) = ", dot(u2, v1), ", (u₂, v₂) = ", dot(u2, v2))

# A function that integrates (1) and computes E(t)

function integrateforward(x0; tfinal=1)
    t = range(0, stop=tfinal, length=1000)
    energy = zeros(length(t))
    upperbound = zeros(length(t))
    dt = t[2]-t[1]
    xt = x0
    P = Matrix(I, 2, 2);
    expAdt = exp(A*dt)
    for j=1:length(t)
        energy[j] = dot(xt, xt)
        upperbound[j] = maximum(eigvals(copy(P')*P))
        P = expAdt*P
        xt = expAdt*xt
    end
    return t, energy, upperbound
end;

A = [-1 5; 0 -2]; 
S, U = eigen(A); u1 = U[:, 1]; u2 =U[:, 2]
Sadj, V = eigen(copy(A')); v1 = V[:, 2]; v2 =V[:, 1]; #to make sure v1 <--> u1
Sinst, Uinst = eigen(A + copy(A')); u01 = Uinst[:, 1]; u02 =Uinst[:, 2];
topt = 0.3; Sopt, Uopt = eigen(exp(A'*topt)*exp(A*topt)); w1 = Uopt[:, 1]; w2 =Uopt[:, 2];
tfinal = 2
figure(figsize=(10, 3.5))
t, energy, upperbound = integrateforward(u1; tfinal=tfinal)
semilogy(t, upperbound, "--k", label="upper bound")
semilogy(t, energy, label="\$ A\$")
t, energy = integrateforward(v2; tfinal=tfinal)
semilogy(t, energy, label="\$A^\\dagger\$")
t, energy = integrateforward(u02; tfinal=tfinal)
semilogy(t, energy, label="\$A +A^\\dagger\$")
t, energy = integrateforward(w2; tfinal=tfinal)
semilogy(t, energy, label="\$e^{A^\\dagger t_*}e^{A t_*}\$ for \$t_*=0.3\$")
xlabel(L"t"); ylabel(L"E(t)/E(0)"); legend();

tfinal = 1
figure(figsize=(10, 3.5))
t, energy, upperbound = integrateforward(u1; tfinal=tfinal)
semilogy(t, upperbound, "--k", label="upper bound")
semilogy(t, energy, label="\$ A\$")
t, energy = integrateforward(v2; tfinal=tfinal)
semilogy(t, energy, label="\$A^\\dagger\$")
t, energy = integrateforward(u02; tfinal=tfinal)
semilogy(t, energy, label="\$A +A^\\dagger\$")
t, energy = integrateforward(w2; tfinal=tfinal)
semilogy(t, energy, label="\$e^{A^\\dagger t_*}e^{A t_*}\$ for \$t_*=0.3\$")
title("zooming in near \$t=0\$")
ylim(0.9, 2); xlabel(L"t"); ylabel(L"E(t)/E(0)"); legend();

A = [-1 1e3; 0 -2]; 
Sadj, V = eigen(copy(A')); v1 = V[:, 2]; v2 =V[:, 1]; #to make sure v1 <--> u1
tfinal = 5
figure(figsize=(10, 3.5))
t, energy, upperbound = integrateforward(v2; tfinal=tfinal)
semilogy(t, upperbound, "--k", label="upper bound")
semilogy(t, energy, label="\$A^\\dagger\$")
t, energy = integrateforward(u02; tfinal=tfinal)
semilogy(t, energy, label="\$A +A^\\dagger\$")
semilogy(t, 1e6exp.(-2*t), ":k", label="\$\\propto e^{-2\\lambda_1 t}\$")
xlabel(L"t"); ylabel(L"E(t)/E(0)"); legend();

A = rand(2, 2) + im*rand(2,2); 

S, U = eigen(A); u1 = U[:, 1]; u2 =U[:, 2]
Sadj, V = eigen(copy(A')); v1 = V[:, 1]; v2 =V[:, 2]

D = Diagonal([dot(u1, v1), dot(u2, v2)])

println("   V = ", V)
println("U⁻¹D = ", inv(U)'*D)


A = [-1 5; 0 -2]; 
S, U = eigen(A); u1 = U[:, 1]; u2 =U[:, 2]
Sadj, V = eigen(copy(A')); v1 = V[:, 2]; v2 =V[:, 1]; #to make sure v1 <--> u1
tfinal = 5
figure(figsize=(10, 3.5))
t, energy = integrateforward(u1; tfinal=tfinal)
semilogy(t, energy, label="\$ A\$")
t, energy = integrateforward(v1; tfinal=tfinal)
semilogy(t, energy, label="\$A^\\dagger\$")
semilogy(t, exp.(-2*t)/abs2(dot(u1, v1)), "-.k", label="eq. (23)")
xlabel(L"t"); ylabel(L"E(t)/E(0)"); legend();