Sqrt fix (still not working)

calculations/ECCC.jl

@@ -0,0 +1,72 @@
+using Roots, LinearAlgebra, Plots
+
+include("../EC.jl")
+include("../common.jl")
+include("../p_space.jl")
+
+μ = 0.5
+V_system(c) = (p, q) -> c*(-5*g0(sqrt(3), p, q) + 2*g0(sqrt(10), p, q)) # ResonanceEC: Eq. (20)
+
+# determining c0 with EC
+temp_c = range(1.1, 0.9, 3)
+p, w = get_mesh([0, 8 - 3im], [512])
+H0 = get_T_matrix(p, μ)
+dim = size(H0, 1)
+V = get_V_matrix(V_system(1), p, w)
+EC = affine_EC(H0, V, w)
+train!(EC, temp_c; ref_eval=-0.2, CAEC=false, verbose=false)
+quick_extrapolate(c) = minimum(abs2, get_extrapolated_evals(EC.H0_EC, EC.H1_EC, EC.N_EC, c, 0))
+c0 = find_zero(quick_extrapolate, 0.85)
+
+# Calculation of training vectors
+training_c = range(1.2, 0.9, 9) # original: range(1.35, 0.9, 5)
+ref_E = -0.3
+order::Int = ceil((length(training_c) - 1) / 2) # order of the Pade approximant
+
+training_E = ComplexF64[]
+training_vecs = Vector{ComplexF64}[]
+
+for c in training_c
+    println("Training for c = $c")
+    H = H0 + c .* V
+    evals, evecs = eigs(H, sigma=ref_E, maxiter=5000, ritzvec=true, check=1)
+    val = nearest(evals, ref_E)
+    push!(training_E, val)
+    vec = evecs[:, nearestIndex(evals, val)]
+    vec ./= sqrt(only(transpose(vec) * vec)) # normalize
+    push!(training_vecs, vec)
+end
+
+training_vecs = hcat(training_vecs...)
+
+# Calculation of target E
+target_c = range(0.78, 0.45, 7) # original: range(0.75, 0.40, 8)
+target_E = [quick_pole_E(V_system(c)) for c in target_c]
+
+# Solve coefficients as a linear system
+M_left_element(c, i) = alt_sqrt(c - c0)^i
+M_left = M_left_element.(training_c, (0:order)')
+a = zeros(ComplexF64, dim, order+1)
+b = zeros(ComplexF64, dim, order+1)
+for i in 1:dim
+    println("Fitting coefficients $i/$dim")
+    M_right = -training_vecs[i, :] .* M_left[:, 2:end] # remove the first column
+    M = hcat(M_left, M_right) # M = [M_left | M_right]
+    sol = M \ training_vecs[i, :]
+    a[i, :] .= sol[1:order+1]
+    b[i, :] .= [1; sol[order+2:end]]
+end
+
+# Pade approximant
+polynomial(a, c) = sum(i -> a[:, i+1] .* alt_sqrt(c - c0)^i, 0:order)
+pade_approx(c) = polynomial(a, c) ./ polynomial(b, c)
+
+# Extrapolate
+extrapolating_c = [training_c; target_c]
+extrapolated_vec = pade_approx.(extrapolating_c)
+extrapolated_E = [only(transpose(vec) * (H0 + c .* V) * vec) for (vec, c) in zip(extrapolated_vec, extrapolating_c)]
+
+# Plotting
+scatter(real.(training_E), imag.(training_E), label="training")
+scatter!(real.(target_E), imag.(target_E), label="target")
+scatter!(real.(extrapolated_E), imag.(extrapolated_E), label="extrapolated", m=:star5)