Merge pull request #122 from SotaYoshida/dev

Dev
SotaYoshida · Jul 2, 2024 · 6c39820 · 6c39820
2 parents e66e293 + 0ccda1a
commit 6c39820
Show file tree

Hide file tree

Showing 18 changed files with 794 additions and 187 deletions.
diff --git a/Manifest.toml b/Manifest.toml
@@ -2,7 +2,7 @@
 
 julia_version = "1.10.0"
 manifest_format = "2.0"
-project_hash = "0f8ce17359409f799c28940c96570dda0c8ecfc0"
+project_hash = "2659f8769473273827cbb29dd4cb2724ffd6c600"
 
 [[deps.Adapt]]
 deps = ["LinearAlgebra", "Requires"]
@@ -90,16 +90,6 @@ git-tree-sha1 = "4b859a208b2397a7a623a03449e4636bdb17bcf2"
 uuid = "83423d85-b0ee-5818-9007-b63ccbeb887a"
 version = "1.16.1+1"
 
-[[deps.ChainRulesCore]]
-deps = ["Compat", "LinearAlgebra"]
-git-tree-sha1 = "2118cb2765f8197b08e5958cdd17c165427425ee"
-uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
-version = "1.19.0"
-weakdeps = ["SparseArrays"]
-
-    [deps.ChainRulesCore.extensions]
-    ChainRulesCoreSparseArraysExt = "SparseArrays"
-
 [[deps.CodecZlib]]
 deps = ["TranscodingStreams", "Zlib_jll"]
 git-tree-sha1 = "59939d8a997469ee05c4b4944560a820f9ba0d73"
@@ -332,12 +322,6 @@ git-tree-sha1 = "ff38ba61beff76b8f4acad8ab0c97ef73bb670cb"
 uuid = "0656b61e-2033-5cc2-a64a-77c0f6c09b89"
 version = "3.3.9+0"
 
-[[deps.GPUArraysCore]]
-deps = ["Adapt"]
-git-tree-sha1 = "2d6ca471a6c7b536127afccfa7564b5b39227fe0"
-uuid = "46192b85-c4d5-4398-a991-12ede77f4527"
-version = "0.1.5"
-
 [[deps.GR]]
 deps = ["Artifacts", "Base64", "DelimitedFiles", "Downloads", "GR_jll", "HTTP", "JSON", "Libdl", "LinearAlgebra", "Pkg", "Preferences", "Printf", "Random", "Serialization", "Sockets", "TOML", "Tar", "Test", "UUIDs", "p7zip_jll"]
 git-tree-sha1 = "27442171f28c952804dede8ff72828a96f2bfc1f"
@@ -471,12 +455,6 @@ git-tree-sha1 = "49fb3cb53362ddadb4415e9b73926d6b40709e70"
 uuid = "b14d175d-62b4-44ba-8fb7-3064adc8c3ec"
 version = "0.2.4"
 
-[[deps.KrylovKit]]
-deps = ["ChainRulesCore", "GPUArraysCore", "LinearAlgebra", "Printf"]
-git-tree-sha1 = "1a5e1d9941c783b0119897d29f2eb665d876ecf3"
-uuid = "0b1a1467-8014-51b9-945f-bf0ae24f4b77"
-version = "0.6.0"
-
 [[deps.LAME_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
 git-tree-sha1 = "f6250b16881adf048549549fba48b1161acdac8c"
@@ -995,11 +973,13 @@ deps = ["IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_j
 git-tree-sha1 = "e2cfc4012a19088254b3950b85c3c1d8882d864d"
 uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
 version = "2.3.1"
-weakdeps = ["ChainRulesCore"]
 
     [deps.SpecialFunctions.extensions]
     SpecialFunctionsChainRulesCoreExt = "ChainRulesCore"
 
+    [deps.SpecialFunctions.weakdeps]
+    ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
+
 [[deps.SplittablesBase]]
 deps = ["Setfield", "Test"]
 git-tree-sha1 = "e08a62abc517eb79667d0a29dc08a3b589516bb5"

diff --git a/Project.toml b/Project.toml
@@ -1,7 +1,7 @@
 name = "NuclearToolkit"
 uuid = "89bb3bae-bcec-43ae-87b7-9dd181dc6334"
 authors = ["SotaYoshida <[email protected]> and contributors"]
-version = "0.4.2"
+version = "0.4.3"
 
 [deps]
 Arpack = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97"
@@ -12,7 +12,6 @@ DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 FLoops = "cc61a311-1640-44b5-9fba-1b764f453329"
 Glob = "c27321d9-0574-5035-807b-f59d2c89b15c"
 HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
-KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77"
 LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
 LatinHypercubeSampling = "a5e1c1ea-c99a-51d3-a14d-a9a37257b02d"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -38,7 +37,6 @@ DocStringExtensions = "0.9"
 FLoops = "0.2"
 Glob = "1"
 HDF5 = "0.16, 0.17"
-KrylovKit = "0.6, 0.7, 0.8"
 LaTeXStrings = "1"
 LatinHypercubeSampling = "1.8"
 MKL = "0.6, 0.7"

diff --git a/example/log_sample_script.txt b/example/log_sample_script.txt
@@ -36,12 +36,12 @@ option in optional_parameters.jl will be used.
 size of dWS (jmax 9 lmax 40 e2max 8 Nnmax 20):
      dtri    4.46 MB    dcgm0    1.11 MB  d6j_int    1.11 MB   d6j_lj    0.28 MB
   d9j_lsj    1.76 MB  dictHOB    0.44 MB
-target: O16 Ref. => Z=8 N=8 E:    -130.997774E1b:     337.500000E2b:    -468.497774 E3b        0.000000
-E:  -132.733407  = E1b 363.64899  + E2b -496.38240   (   -55.149 -367.063  -74.170), + E3b   0.00000
+target: O16 Ref. => Z=8 N=8 E:  -132.733407  = E1b 363.64899  + E2b -496.38240   (   -55.149 -367.063  -74.170), + E3b   0.00000
 EMP2    -23.70204 1b   -0.00000 pp   -1.43732 pn  -20.80963 nn   -1.45510
 EMP3      0.37877 pp    -0.110 hh    -1.895 ph     2.384
 E_HF   -132.73341  E_MBPT(3) =    -156.0567  Eexp:     -127.619
 
+optional_parameters.jl is used for IMSRG
 parameters in optional_parameters.jl will be used.
 def-by-run d6j_lj done! 4191
  step:        s             E0   ||Omega_1||   ||Omega_2||     ||Eta_1||     ||Eta_2||    Ncomm.  nwritten

diff --git a/src/IMSRG.jl/emulator_imsrg.jl b/src/IMSRG.jl/emulator_imsrg.jl
@@ -250,7 +250,6 @@ function svd_Op_twobody(s,Op::Operator,Chan2b;verbose=true,max_rank=20)
         SVD = LinearAlgebra.svd(mat)
         if verbose 
             print("ch ",@sprintf("%5i",ch), " JPT=($J,$P,$Tz)\t fullrank ", @sprintf("%5i",fullrank), "/ dim= ",@sprintf("%5i",dim),"  ")
-            #print_vec("singular values", SVD.S)
             tol = 1.e-8
             hit = 0
             dimfull += div(dim*(dim+1),2)
@@ -292,16 +291,41 @@ function constructor_snapshot_matrix(fns)
     return parse(Float64,s_end),X,Y
 end
 
-function get_DMD_operator(X,Y,r)
-    Z = svds(X; nsv=r)[1]
-    U_r, S_r, Vt_r = Z.U, Z.S, Z.Vt
+function call_SVD(X, r, method)
+    if method == "Arpack"
+        Z = svds(X; nsv=r)[1]
+        return Z.U, Z.S, Z.Vt
+    elseif method == "KrylovKit"
+        m, n = size(X)
+        vals, lvecs, rvevs, info  = svdsolve(X, m, r)
+        U_r = zeros(Float64, m, r)
+        S_r = zeros(Float64, r)
+        Vt_r = zeros(Float64, r, n)
+        S_r .= vals[1:r]
+        for i = 1:r
+            U_r[:,i] .= lvecs[i]
+            Vt_r[i,:] .= rvevs[i]
+        end
+        return U_r, S_r, Vt_r
+    else
+        if method != "full"
+            println("$method is not supported! svd in LinearAlgebra is called instead.")
+        end
+        Z = svd(X)
+        if r <= rank(X) 
+            return Z.U[:,1:r], Z.S[1:r], Z.Vt[1:r,:]
+        else
+            return Z.U, Z.S, Z.Vt            
+        end
+    end
+end
 
+function get_DMD_operator(U_r, S_r, Vt_r, Y)
     S_r_inv = diagm( 1.0 ./ S_r)
     Z = BLAS.gemm('T', 'N', 1.0, Vt_r, S_r_inv)
     YZ = BLAS.gemm('N','N',1.0, Y, Z)
     Atilde = BLAS.gemm('T', 'N', 1.0, U_r, YZ)
-
-    return U_r, Atilde
+    return Atilde
 end
 
 function check_DMD_norm(X, Y, r, U_r, Atilde; verbose=false)
@@ -313,8 +337,7 @@ function check_DMD_norm(X, Y, r, U_r, Atilde; verbose=false)
     for k = 1:size(Y)[2]
         x_k .= x_new
         BLAS.gemv!('N', 1.0, Atilde, x_k, 0.0, x_new)
-        Y_latent[:,k] .= x_new
-
+        Y_latent[:,k] .= x_new        
     end
     Yapprox = BLAS.gemm('N', 'N', 1.0, U_r, Y_latent)
     if verbose 
@@ -325,6 +348,8 @@ function check_DMD_norm(X, Y, r, U_r, Atilde; verbose=false)
         end
     end
     println("norm(Y-Yapprox,Inf) = ", @sprintf("%10.4e",norm(Y - Yapprox,Inf)), " Fro. ", @sprintf("%10.4e",norm(Y - Yapprox,2)))
+    tnorm = max(norm(Y - Yapprox,Inf), norm(Y - Yapprox,2))
+    return tnorm
 end
 
 function check_stationarity(z, z_k, z_pred,  Atilde, itnum=1000)
@@ -359,11 +384,10 @@ function extrapolate_DMD(x_start, U_r, Atilde, s_pred, fn_exact, s_end, ds, nuc,
             z_k .= 0.0
             z_new .= z1_r
             s = s_end
-            while true
+            while s < s_target
                 z_k .= z_new
                 BLAS.gemv!('N', 1.0, Atilde, z_k, 0.0, z_new)
                 s += ds
-                if s >= s_target; break; end
             end
             x_pred = BLAS.gemv('N', 1.0, U_r, z_k)
             E_imsrg = 0.0
@@ -403,6 +427,57 @@ function read_dmdvec_hdf5(fn)
     return vec
 end
 
+function util_SVD(X, Y, r_max, tol_svd, method, allow_fullSVD, to; dont_care_stationarity=false)
+    @assert method == "full" || method == "Arpack" || method == "KrylovKit" "method must be full, KrylovKit, or Arpack"
+    fullrank = rank(X)
+    r_used = min(r_max,fullrank)
+    println("fullrank $fullrank r_used $r_used")
+    if allow_fullSVD || method == "full"
+        @timeit to "full SVD" U,Sigma,Vt = call_SVD(X, r_used, "full")
+        r_used = min(r_max, fullrank, sum(Sigma .> tol_svd))
+        print_vec("checking full Sigma[1:$r_used]...", Sigma[1:r_used];ine=true)
+    end
+    @timeit to "SVD($method)" U_r, S_r, Vt_r = call_SVD(X,r_used,method)
+
+    redright = true
+    u = @view U_r[:,1:r_used]
+    s = @view S_r[1:r_used]
+    vt = @view Vt_r[1:r_used,:]
+    Atilde = get_DMD_operator(u,s,vt,Y)
+    while redright && !dont_care_stationarity 
+        Atilde = get_DMD_operator(u,s,vt,Y)
+        evals_Atilde = eigvals(Atilde)
+        redright = !check_stationarity_Atilde(evals_Atilde)
+        print_vec("checking Sigma...", S_r[1:r_used];ine=true)
+        if redright == false; break; end
+        println("truncation at rank <= $r_used does not give converged results")
+        r_used -= 1
+        @assert r_used > 0 "Aborted because the snapshots don't give converged resutls. Please modify smin/smax of snapshots."
+        u = @view U_r[:,1:r_used]
+        s = @view S_r[1:r_used]
+        vt = @view Vt_r[1:r_used,:]
+    end
+    print_vec("singular values r <= $r_used ", S_r[1:r_used];ine=true)
+    return fullrank, r_used, u, Atilde
+end
+
+function check_packages_SVD(X)
+    r = min(10,rank(X))
+    m, n = size(X)
+    SVD = svd(X)
+    s = SVD.S    
+    s_K = svdsolve(X, m, r)[1]
+    s_A = svds(X; nsv=r)[1].S
+
+    @assert issorted(s,rev=true) "s must be sorted in descending order"
+
+    println("Comparing the singular values by different packages....")
+    print_vec("full SVD (LinearAlgebra)", s[1:r];ine=true)
+    print_vec("trun.SVD (KrylovKit)    ", s_K[1:r];ine=true)
+    print_vec("trun.SVD (Arpack)       ", s_A[1:r];ine=true)
+    return nothing
+end
+
 """
 main API for DMD
 
@@ -418,11 +493,22 @@ main API for DMD
 - `s_pred::Vector{Float64}`: values of `s` for the extrapolation
 - `fn_exact::Vector{String}`: filenames of the exact data for the extrapolation, which must have the same length as `s_pred`
 - `allow_fullSVD::Bool`: if `true`, the full SVD is performed and the rank is determined by `trank` and the tolerance `tol_svd`
-- `tol_svd::Float64`: tolerance for the singular values for the truncated SVD
+- `tol_svd::Float64`: tolerance for the singular values for the truncated SVD, too small value gives noisy behavior
 - `inttype::String`: interaction type, which is used for the filename of the output data
+- `methodSVD::String`: method for the truncated SVD.
+    - `"full"`: full SVD (LinearAlgebra)
+    - `"Arpack"`: truncated SVD svds Arpack.jl
+    - `"KrylovKit"`: truncated SVD using `svdsolve[1]` in KrylovKit.jl, this may be not appropriate for the current problem, so do not use it here.
+- `oupdir::String`: output directory for the emulated dmdvec data
+- `is_show::Bool`: if `true`, the TimerOutput result is shown
+- `debugmode::Bool`: if `true`, the packages for the SVD are checked
+- `dont_care_stationarity::Bool`: if `true`, the stationarity of the Atilde is not checked
 """
-function dmd_main(emax, nuc, fns, trank, smin, ds;s_pred=Float64[],fn_exact=String[],
-                  allow_fullSVD=true,tol_svd=1e-7,inttype="",oupdir="flowOmega/")
+function dmd_main(emax, nuc, fns, r_max, smin, ds;s_pred=Float64[],fn_exact=String[],
+                  allow_fullSVD=true,tol_svd=1e-6,inttype="",
+                  methodSVD="Arpack", oupdir="flowOmega/",is_show=true, debugmode=false,
+                  dont_care_stationarity=true)
+    to = TimerOutput()
     if !isdir("flowOmega")
         println("dir. flowOmega is created!")
         mkdir("flowOmega")
@@ -433,21 +519,55 @@ function dmd_main(emax, nuc, fns, trank, smin, ds;s_pred=Float64[],fn_exact=Stri
     s_end, X,Y = constructor_snapshot_matrix(fns)
     println("Snapshot from smin $smin s_end $s_end ds $ds")
 
-    # truncated SVD using Arpack.jl and construct tilde(A)
-    fullrank = rank(X)
-    r = trank
-    if allow_fullSVD
-        SVD = svd(X)
-        sigma_full = SVD.S
-        r = min(r, fullrank, sum(sigma_full .> tol_svd))
-        println("fullrank $(fullrank) rank $r")
-        print_vec("singular values", sigma_full[1:r];ine=true)
+    if debugmode
+        check_packages_SVD(X)
     end
-    U_r, Atilde = get_DMD_operator(X,Y,r)    
-    check_DMD_norm(X, Y, r, U_r, Atilde)
+
+    fullrank, r, U_r, Atilde = util_SVD(X, Y, r_max, tol_svd, methodSVD, allow_fullSVD, to; dont_care_stationarity=dont_care_stationarity) 
+    evals_Atilde = eigvals(Atilde)
+    check_stationarity_Atilde(evals_Atilde)
+    plot_Atilde_eigvals(evals_Atilde, emax, nuc, smin, s_end, ds, inttype, fullrank, r)
 
+    if is_show
+        show(to); println()
+    end
     # extrapolation
     extrapolate_DMD(Y[:,end],U_r, Atilde, s_pred, fn_exact, s_end, ds, nuc, inttype, emax, oupdir)
 
     return nothing
 end
+
+function check_stationarity_Atilde(evals; tol=1.e-2)
+    tf = all(abs.(evals) .<= 1.0+tol)
+    println("eigenvalues of Atilde ", evals)
+    println("stationary?  may be... ", tf)
+    return tf
+end
+
+function plot_Atilde_eigvals(evals, emax, nuc, smin, s_end, ds, inttype, fullrank, r)
+    if !isdir("pic") 
+        mkdir("pic")
+    end
+    fn = "pic/Atilde_eigvals_e$(emax)_$(nuc)_smin$(smin)_send$(s_end)_ds$(ds)_rank$(r)_of_$(fullrank).pdf"
+    # make figure square
+    p = plot(size=(300,300))
+
+    # draw unit circle
+    θ = LinRange(0,2π,100)
+    plot!(cos.(θ),sin.(θ),lw=2,alpha=0.8,label=:none)
+
+    # plot eigenvalues .<= 1.0
+    inon = evals[findall(abs.(evals) .<= 1.0)]
+    plot!(real(inon),imag(inon),st=:scatter,label=L"|z| \leq 1", alpha=0.8)
+
+    out = evals[findall(abs.(evals) .> 1.0)]
+    if length(out) > 0
+        plot!(real(out),imag(out),st=:scatter,label=L"|z| > 1",alpha=0.8)
+    end
+
+    xlabel!(L"\mathrm{Re} (z)")
+    ylabel!(L"\mathrm{Im} (z)")
+    savefig(p,fn)
+
+    return nothing
+end
diff --git a/src/IMSRG.jl/imsrg_util.jl b/src/IMSRG.jl/imsrg_util.jl
@@ -862,7 +862,7 @@ Function to write temporary binary files of Operator matrix elements, when split
 """
 function write_omega_bin(binfo::basedat,Chan2b::Vector{chan2b},n_written::Int,Omega::Operator,s::Float64,E0::Float64;Oplabel="Omega")
     if !isdir("flowOmega")  && n_written==0
-        run(`mkdir flowOmega`)
+        mkdir("flowOmega")
     end    
     pid = getpid()
     nw = n_written + 1
@@ -1076,6 +1076,7 @@ function init_IMSRGobject(HFobj,filename;smax=500.0,dsmax=0.5,maxnormOmega=0.25,
         println("Since $filename is not found, the default parameters will be used.")
         return IMSRGobj
     else
+        println("$filename is used for IMSRG")
         read_imsrg_parameter!(filename,IMSRGobj)
         return IMSRGobj
     end

diff --git a/src/NuclearToolkit.jl b/src/NuclearToolkit.jl
@@ -9,7 +9,6 @@ using DocStringExtensions
 using FLoops
 using Glob
 using HDF5
-using KrylovKit
 using LatinHypercubeSampling
 using LaTeXStrings
 using LinearAlgebra

diff --git a/src/ShellModel/KSHELL.jl b/src/ShellModel/KSHELL.jl
@@ -6,7 +6,7 @@ struct ksl_state
 end
 
 struct kshell_nuc
-    nuc::nuclei
+    nuc::nucleus
     sntf::String
     Egs::Float64
     Egs_target::Float64