speed up (::PiecewiseLegendreFTArray)(n)

Project.toml

@@ -18,9 +18,9 @@ SpecialFunctions = "2"
 julia = "1"
 
 [extras]
+AssociatedLegendrePolynomials = "2119f1ac-fb78-50f5-8cc0-dda848ebdb19"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
-AssociatedLegendrePolynomials = "2119f1ac-fb78-50f5-8cc0-dda848ebdb19"
 
 [targets]
 test = ["Random", "Test", "AssociatedLegendrePolynomials"]

playground_pluto.jl

@@ -1,24 +1,9 @@
 ### A Pluto.jl notebook ###
-# v0.19.0
+# v0.19.3
 
 using Markdown
 using InteractiveUtils
 
-# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
-macro bind(def, element)
-    quote
-        local iv = try
-            Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"),
-                                           "AbstractPlutoDingetjes")].Bonds.initial_value
-        catch
-            b -> missing
-        end
-        local el = $(esc(element))
-        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
-        el
-    end
-end
-
 # ╔═╡ 842c7aec-b17a-11ec-127b-3f91b4687627
 begin
     using Revise
@@ -33,147 +18,124 @@ using LinearAlgebra, Random
 # ╔═╡ 36915eda-1b70-4fea-95f4-f252df56c060
 using Plots, PlutoUI
 
-# ╔═╡ 8d00dcf9-e24e-41d7-af8a-a5e7e4ab3939
-using PyCall
-
-# ╔═╡ b8c2c8c7-9a42-4394-a390-aecf9a1c6f34
-#basis = DimensionlessBasis(boson, 10)
-
-# ╔═╡ 19d59ed7-aaea-4e47-9708-606acb03cdbb
+# ╔═╡ 227eb6c4-d35f-4d6c-b79e-dc494088fc97
+using BenchmarkTools
+
+# ╔═╡ 953d4311-7ce9-4359-a7a2-b5f5f3192c65
+md"""
+```python
+import sparse_ir as ir, numpy as np
+basis = ir.FiniteTempBasis('F', 10, 4, 1e-6)
+U = 0.5
+def rhow(w):
+  return np.sqrt(1 - w.clip(-1,1)**2) / np.pi
+rho0l = basis.v.overlap(rhow)
+G0l = -basis.s * rho0l
+Gl = G0l
+stau = ir.TauSampling(basis)
+siw = ir.MatsubaraSampling(basis)
+for _ in range(30):
+  Gtau = stau.evaluate(Gl)
+  Sigmatau = U**2 * Gtau**3
+  Sigmal = stau.fit(Sigmatau)
+  Sigmaiw = siw.evaluate(Sigmal)
+  G0iw = siw.evaluate(G0l)
+  Giw = 1/(1/G0iw - Sigmaiw)
+  Gl = siw.fit(Giw)
+```
+"""
+
+# ╔═╡ 0af86ef4-0740-4101-b209-2eb3d23002e3
+basis = FiniteTempBasis(fermion, 10000.0, 4.0, 1e-6)
+
+# ╔═╡ f636e647-c321-4625-9923-3ff0b8a37a06
+U = 0.5
+
+# ╔═╡ d6e782e5-5ec9-4cdf-b308-647c4ae37ce0
+rhow(w) = √(1 - clamp(w, -1, 1)^2) / π
+
+# ╔═╡ 353eedbc-664b-463d-b795-6f9887a81dc0
+rho0l = overlap.(basis.v, rhow)
+
+# ╔═╡ 1a2b33e9-6d3a-4c7b-b05a-253a13509dfe
+G0l = - basis.s .* rho0l
+
+# ╔═╡ a12f1f1e-1f21-4cbd-ad8a-be0f2ed579bc
+size(basis)
+
+# ╔═╡ f4f0a07d-a579-4ed7-ada5-7e882cdaed55
 begin
-    ε = 1e-15
-    stat = boson
-
-    Λ = 1e4
-    wmax = 1.0
-    pole = 0.1 * wmax
-    β = Λ / wmax
-end
-
-# ╔═╡ 1cb616f4-7957-4f9d-8942-c881b87b94bc
-basis = FiniteTempBasis(stat, β, wmax, ε)
-
-# ╔═╡ 8886718d-414b-4d87-a85d-0c5638e440c7
-@bind n Slider(eachindex(basis.u))
-
-# ╔═╡ 48dc95f8-fa62-4857-9d4a-b336832677d0
-# ╠═╡ disabled = true
-#=╠═╡
-let
-	x = -1:0.0001:1
-	p = plot()
-	poly = basis.u[n]
-	plot!(p, x, poly, label="poly")
-	for i in 1:6#basis.u.polyorder
-		poly = deriv(poly)
-		plot!(p, x, x -> poly(x) / (2basis.kernel.Λ)^i, label="derivative $i")
+	Gl = G0l
+	stau = TauSampling(basis);
+	siw = MatsubaraSampling(basis);
+	for _ in 1:30
+		Gtau = evaluate(stau, Gl)
+		Sigmatau = U^2 * Gtau.^3
+		Sigmal = fit(stau, Sigmatau)
+		Sigmaiw = evaluate(siw, Sigmal)
+		G0iw = evaluate(siw, G0l)
+		Giw = 1 ./ (1 ./ G0iw - Sigmaiw)
+		Glold = copy(Gl)
+		Gl = fit(siw, Giw)
+		println(norm(Gl - Glold))
 	end
-	p
-end
-  ╠═╡ =#
-
-# ╔═╡ a9ce84bb-7f44-40f0-9558-6ba930d9fb6d
-# ╠═╡ disabled = true
-#=╠═╡
-@bind d Slider(0:basis.u.polyorder, show_value=true)
-  ╠═╡ =#
-
-# ╔═╡ 60508545-94c1-4af0-b5ed-a128c8ef8aef
-# ╠═╡ disabled = true
-#=╠═╡
-let
-	x = -1:0.0001:1
-	poly = basis.u[n]
-	plot(x, deriv(poly, d))
-end
-  ╠═╡ =#
-
-# ╔═╡ 2cd01891-3cdf-48f4-a60c-e2052d2c5a59
-basis_py = pyimport("sparse_ir").FiniteTempBasis("B", β, wmax, ε)
-
-# ╔═╡ 00bb8304-c307-4a2b-b5dc-9fb4991116a0
-for i in 1:104
-    #basis.uhat[i].model.moments .= basis_py.uhat[i-1]._model.moments
+	real.(Gl)
 end
 
-# ╔═╡ 6fa0439c-198f-4802-9e81-e25ac6fc580b
-basis.uhat[2].model.moments
-
-# ╔═╡ 23a0ff60-e431-40c9-b066-bbd787238f8f
-basis_py.uhat.__getitem__(0:1)._model.moments
+# ╔═╡ d1eed5b7-4fdb-4a4e-a08f-4beb7f709dc9
+moments(polyFTs::SparseIR.PiecewiseLegendreFTArray) = collect(eachrow(reduce(hcat, p.model.moments for p in polyFTs)))
 
-# ╔═╡ 53b9baba-4651-4d96-adf5-6122bdb41e97
-vec(basis_py.uhat[0]._model.moments)
-
-# ╔═╡ d7e3a562-8f93-4791-9ab8-d57d6bb6bc6a
-begin
-    stat_shift = (stat == fermion) ? 1 : 0
-    weight = (stat == fermion) ? 1 : 1 / tanh(0.5 * Λ * pole / wmax)
-    gl = -basis.s .* basis.v(pole) * weight
-    func_G(n) = 1 / (im * (2n + stat_shift) * π / β - pole)
+# ╔═╡ 88ea79e3-03bd-4edc-9f04-7682a5e9d36d
+function moments2(polyFTs::SparseIR.PiecewiseLegendreFTArray)
+	n = length(first(polyFTs).model.moments)
+	[[p.model.moments[i] for p in polyFTs] for i in 1:n]
 end
 
-# ╔═╡ af735714-2fff-41de-aa2e-51c563313bb5
-begin
-    # Compute G(iwn) using unl
-    matsu_test = Int[-1, 0, 1, 1e2, 1e4, 1e6, 1e8, 1e10, 1e12]
-    prj_w = transpose(basis.uhat(2matsu_test .+ stat_shift))
-    Giwn_t = prj_w * gl
+# ╔═╡ fcc7578b-29b6-44dc-8ea0-cb2a0b96c6b3
+function moments3(polyFTs::SparseIR.PiecewiseLegendreFTArray)
+	n = length(first(polyFTs).model.moments)
+	#[[p.model.moments[i] for p in polyFTs] for i in 1:n]
+	mat = reduce(hcat, p.model.moments for p in polyFTs)
+	eachrow(mat)
 end
 
-# ╔═╡ e7c4cc78-57d8-4a9a-8282-051779375abe
-# Compute G(iwn) from analytic expression
-Giwn_ref = func_G.(matsu_test)
-
-# ╔═╡ 682d942c-6442-44d8-99af-60324cdd85ae
-basis.uhat[1:2](Int(4e5))
-
-# ╔═╡ 7f17ea66-1c6c-49f0-87c6-dc226f15af59
-magnitude = maximum(abs, Giwn_ref)
-
-# ╔═╡ 7140a21f-ee18-4e33-86da-d622f69e9956
-diff = abs.(Giwn_t - Giwn_ref)
+# ╔═╡ 9666b21b-4ef5-4588-9d41-1d43ccc1eadf
+size(basis.uhat)
 
-# ╔═╡ 0f1b8564-b891-4910-a3df-4decb0ae7d96
-tol = max(10 * last(basis.s) / first(basis.s), 1e-10)
+# ╔═╡ 1773939a-4ea4-4856-a769-e66e9018b3e6
+@btime moments(basis.uhat)
 
-# ╔═╡ e2898e3b-0aac-4909-a4bb-fbc6161fcc41
-diff ./ magnitude
+# ╔═╡ b934f8e1-1ba1-4281-8e68-b927e10c2973
+@btime moments2(basis.uhat)
 
-# ╔═╡ 27a6ffed-2e5f-4ee2-8368-f49e9226a13e
-maximum(diff ./ magnitude)
+# ╔═╡ 951a8aa2-0e42-4822-afa0-b5746842f3df
+@btime moments3(basis.uhat)
 
-# ╔═╡ 6ac780bd-d5dc-4230-96f3-8a6b71cf15cb
-diff ./ Giwn_ref
+# ╔═╡ e95e6459-83d2-4769-936c-c368e00dcc65
+@time basis.uhat(10000003)
 
-# ╔═╡ ffb9a88b-7b65-47cb-9ee6-9bb31e433e2c
-maximum(abs, diff ./ Giwn_ref)
+# ╔═╡ 2be2f599-6152-4a8f-bb7c-ef4293e2de9e
+@time [û(10000003) for û in basis.uhat]
 
 # ╔═╡ Cell order:
 # ╠═842c7aec-b17a-11ec-127b-3f91b4687627
 # ╠═c92e8377-29e9-4b29-ae11-e0423459722d
 # ╠═36915eda-1b70-4fea-95f4-f252df56c060
-# ╠═b8c2c8c7-9a42-4394-a390-aecf9a1c6f34
-# ╠═8886718d-414b-4d87-a85d-0c5638e440c7
-# ╠═48dc95f8-fa62-4857-9d4a-b336832677d0
-# ╠═a9ce84bb-7f44-40f0-9558-6ba930d9fb6d
-# ╠═60508545-94c1-4af0-b5ed-a128c8ef8aef
-# ╠═19d59ed7-aaea-4e47-9708-606acb03cdbb
-# ╠═1cb616f4-7957-4f9d-8942-c881b87b94bc
-# ╠═2cd01891-3cdf-48f4-a60c-e2052d2c5a59
-# ╠═00bb8304-c307-4a2b-b5dc-9fb4991116a0
-# ╠═6fa0439c-198f-4802-9e81-e25ac6fc580b
-# ╠═23a0ff60-e431-40c9-b066-bbd787238f8f
-# ╠═53b9baba-4651-4d96-adf5-6122bdb41e97
-# ╠═d7e3a562-8f93-4791-9ab8-d57d6bb6bc6a
-# ╠═af735714-2fff-41de-aa2e-51c563313bb5
-# ╠═e7c4cc78-57d8-4a9a-8282-051779375abe
-# ╠═682d942c-6442-44d8-99af-60324cdd85ae
-# ╠═7f17ea66-1c6c-49f0-87c6-dc226f15af59
-# ╠═7140a21f-ee18-4e33-86da-d622f69e9956
-# ╠═0f1b8564-b891-4910-a3df-4decb0ae7d96
-# ╠═e2898e3b-0aac-4909-a4bb-fbc6161fcc41
-# ╠═27a6ffed-2e5f-4ee2-8368-f49e9226a13e
-# ╠═6ac780bd-d5dc-4230-96f3-8a6b71cf15cb
-# ╠═ffb9a88b-7b65-47cb-9ee6-9bb31e433e2c
-# ╠═8d00dcf9-e24e-41d7-af8a-a5e7e4ab3939
+# ╟─953d4311-7ce9-4359-a7a2-b5f5f3192c65
+# ╠═0af86ef4-0740-4101-b209-2eb3d23002e3
+# ╠═f636e647-c321-4625-9923-3ff0b8a37a06
+# ╠═d6e782e5-5ec9-4cdf-b308-647c4ae37ce0
+# ╠═353eedbc-664b-463d-b795-6f9887a81dc0
+# ╠═1a2b33e9-6d3a-4c7b-b05a-253a13509dfe
+# ╠═a12f1f1e-1f21-4cbd-ad8a-be0f2ed579bc
+# ╠═f4f0a07d-a579-4ed7-ada5-7e882cdaed55
+# ╠═227eb6c4-d35f-4d6c-b79e-dc494088fc97
+# ╠═d1eed5b7-4fdb-4a4e-a08f-4beb7f709dc9
+# ╠═88ea79e3-03bd-4edc-9f04-7682a5e9d36d
+# ╠═fcc7578b-29b6-44dc-8ea0-cb2a0b96c6b3
+# ╠═9666b21b-4ef5-4588-9d41-1d43ccc1eadf
+# ╠═1773939a-4ea4-4856-a769-e66e9018b3e6
+# ╠═b934f8e1-1ba1-4281-8e68-b927e10c2973
+# ╠═951a8aa2-0e42-4822-afa0-b5746842f3df
+# ╠═e95e6459-83d2-4769-936c-c368e00dcc65
+# ╠═2be2f599-6152-4a8f-bb7c-ef4293e2de9e

src/SparseIR.jl

@@ -2,7 +2,7 @@
 module SparseIR
 
 using IntervalRootFinding: roots as roots_irf, Interval, isunique, interval, mid, Newton
-using LinearAlgebra: svd, SVD, QRIteration
+using LinearAlgebra: dot, svd, SVD, QRIteration
 using LowRankApprox: psvd
 using QuadGK: gauss, kronrod, quadgk
 using SpecialFunctions: sphericalbesselj as sphericalbesselj_sf

src/_specfuncs.jl

@@ -54,7 +54,7 @@ legvander(x, deg)
 
 Pseudo-Vandermonde matrix of degree `deg`.
 """
-function legvander(x::Array{T,N}, deg::Integer) where {T,N}
+function legvander(x::AbstractArray{T,N}, deg::Integer) where {T,N}
     deg ≥ 0 || throw(DomainError(deg, "legvander needs a non-negative degree"))
 
     # leading dimensions

src/kernel.jl

@@ -52,7 +52,7 @@ end
 
 function LogisticKernel(Λ)
     Λ ≥ 0 || throw(DomainError(Λ, "Kernel cutoff Λ must be non-negative"))
-    LogisticKernel(float(Λ))
+    return LogisticKernel(float(Λ))
 end
 
 @doc raw"""