Benchmarks

benchmarks-sourcecode/generate_html.py

@@ -0,0 +1,21 @@
+import os
+import subprocess
+
+sourcedir = "."
+builddir = "."
+
+names = os.listdir(sourcedir)
+names = filter(lambda x: not x == os.path.basename(__file__) and not x.endswith(".html"), names)
+# names = filter(lambda x: "raman" in x, names)
+
+for name in names:
+    print("Converting", name)
+    barename, extension = os.path.splitext(name)
+    sourcepath = os.path.join(sourcedir, name)
+    targetpath = os.path.join(builddir, barename + "_" + extension[1:] + ".html")
+    subprocess.run(["pygmentize",
+                        "-o", targetpath,
+                        sourcepath])
+
+f = open("../css/sourcecode.css", "w")
+subprocess.run(["pygmentize", "-f", "html", "-a", ".highlight", "-S", "default"], stdout=f)

benchmarks-sourcecode/timeevolution_master.jl

@@ -0,0 +1,21 @@
+κ = 1.
+η = 4κ
+Δ = 0
+tmax = 100
+tsteps = 201
+T = Vector(linspace(0, tmax, tsteps))
+
+b = FockBasis(Ncutoff-1)
+a = destroy(b)
+ad = dagger(a)
+n = number(b)
+H = Δ*ad*a + η*(a + ad)
+J = [destroy(b)]
+Γ = [2κ]
+
+Ψ₀ = fockstate(b, 0)
+ρ₀ = Ψ₀ ⊗ dagger(Ψ₀)
+exp_n = Float64[]
+fout(t, ρ) = push!(exp_n, real(expect(n, ρ)))
+timeevolution.master(T, ρ₀, H, J; Gamma=Γ, fout=fout, reltol=1e-6, abstol=1e-8)
+return exp_n

benchmarks-sourcecode/timeevolution_master.m

@@ -0,0 +1,25 @@
+kap = 1.;
+eta = 4*kap;
+delta = 0;
+tmax = 100;
+tsteps = 201;
+tlist = linspace(0, tmax, tsteps);
+a = destroy(Ncutoff);
+ad = a';
+H = delta*ad*a + eta*(a + ad);
+C = sqrt(2*kap)*a;
+LH = -1i * (spre(H) - spost(H));
+L1 = spre(C)*spost(C')-0.5*spre(C'*C)-0.5*spost(C'*C);
+L = LH+L1;
+
+
+psi0 = basis(Ncutoff, 1);
+rho0 = psi0 * psi0';
+
+% Set up options, if required
+options.reltol = 1e-6;
+options.abstol = 1e-6;
+% Write out the data file
+ode2file('file1.dat', L, rho0, tlist, options);
+% Call the equation solver
+odesolve('file1.dat','file2.dat');

benchmarks-sourcecode/timeevolution_master.py

@@ -0,0 +1,15 @@
+kap = 1.
+eta = 4*kap
+delta = 0
+tmax = 100
+tsteps = 201
+tlist = np.linspace(0, tmax, tsteps)
+
+a = qt.destroy(Ncutoff)
+ad = a.dag()
+H = delta*ad*a + eta*(a + ad)
+c_ops = [np.sqrt(2*kap)*a]
+
+psi0 = qt.basis(Ncutoff, 0)
+n = qt.mesolve(H, psi0, tlist, c_ops, [ad*a], options=options).expect[0]
+return n

benchmarks-sourcecode/timeevolution_master_jl.html

@@ -0,0 +1,22 @@
+<div class="highlight"><pre><span></span><span class="n">κ</span> <span class="o">=</span> <span class="mf">1.</span>
+<span class="n">η</span> <span class="o">=</span> <span class="mi">4</span><span class="n">κ</span>
+<span class="n">Δ</span> <span class="o">=</span> <span class="mi">0</span>
+<span class="n">tmax</span> <span class="o">=</span> <span class="mi">100</span>
+<span class="n">tsteps</span> <span class="o">=</span> <span class="mi">201</span>
+<span class="n">T</span> <span class="o">=</span> <span class="kt">Vector</span><span class="p">(</span><span class="n">linspace</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span> <span class="n">tmax</span><span class="p">,</span> <span class="n">tsteps</span><span class="p">))</span>
+
+<span class="n">b</span> <span class="o">=</span> <span class="n">FockBasis</span><span class="p">(</span><span class="n">Ncutoff</span><span class="o">-</span><span class="mi">1</span><span class="p">)</span>
+<span class="n">a</span> <span class="o">=</span> <span class="n">destroy</span><span class="p">(</span><span class="n">b</span><span class="p">)</span>
+<span class="n">ad</span> <span class="o">=</span> <span class="n">dagger</span><span class="p">(</span><span class="n">a</span><span class="p">)</span>
+<span class="n">n</span> <span class="o">=</span> <span class="n">number</span><span class="p">(</span><span class="n">b</span><span class="p">)</span>
+<span class="n">H</span> <span class="o">=</span> <span class="n">Δ</span><span class="o">*</span><span class="n">ad</span><span class="o">*</span><span class="n">a</span> <span class="o">+</span> <span class="n">η</span><span class="o">*</span><span class="p">(</span><span class="n">a</span> <span class="o">+</span> <span class="n">ad</span><span class="p">)</span>
+<span class="n">J</span> <span class="o">=</span> <span class="p">[</span><span class="n">destroy</span><span class="p">(</span><span class="n">b</span><span class="p">)]</span>
+<span class="n">Γ</span> <span class="o">=</span> <span class="p">[</span><span class="mi">2</span><span class="n">κ</span><span class="p">]</span>
+
+<span class="n">Ψ₀</span> <span class="o">=</span> <span class="n">fockstate</span><span class="p">(</span><span class="n">b</span><span class="p">,</span> <span class="mi">0</span><span class="p">)</span>
+<span class="n">ρ₀</span> <span class="o">=</span> <span class="n">Ψ₀</span> <span class="n">⊗</span> <span class="n">dagger</span><span class="p">(</span><span class="n">Ψ₀</span><span class="p">)</span>
+<span class="n">exp_n</span> <span class="o">=</span> <span class="kt">Float64</span><span class="p">[]</span>
+<span class="n">fout</span><span class="p">(</span><span class="n">t</span><span class="p">,</span> <span class="n">ρ</span><span class="p">)</span> <span class="o">=</span> <span class="n">push!</span><span class="p">(</span><span class="n">exp_n</span><span class="p">,</span> <span class="n">real</span><span class="p">(</span><span class="n">expect</span><span class="p">(</span><span class="n">n</span><span class="p">,</span> <span class="n">ρ</span><span class="p">)))</span>
+<span class="n">timeevolution</span><span class="o">.</span><span class="n">master</span><span class="p">(</span><span class="n">T</span><span class="p">,</span> <span class="n">ρ₀</span><span class="p">,</span> <span class="n">H</span><span class="p">,</span> <span class="n">J</span><span class="p">;</span> <span class="n">Gamma</span><span class="o">=</span><span class="n">Γ</span><span class="p">,</span> <span class="n">fout</span><span class="o">=</span><span class="n">fout</span><span class="p">,</span> <span class="n">reltol</span><span class="o">=</span><span class="mf">1e-6</span><span class="p">,</span> <span class="n">abstol</span><span class="o">=</span><span class="mf">1e-8</span><span class="p">)</span>
+<span class="k">return</span> <span class="n">exp_n</span>
+</pre></div>

benchmarks-sourcecode/timeevolution_master_m.html

@@ -0,0 +1,26 @@
+<div class="highlight"><pre><span></span><span class="n">kap</span> <span class="p">=</span> <span class="mf">1.</span><span class="p">;</span>
+<span class="n">eta</span> <span class="p">=</span> <span class="mi">4</span><span class="o">*</span><span class="n">kap</span><span class="p">;</span>
+<span class="n">delta</span> <span class="p">=</span> <span class="mi">0</span><span class="p">;</span>
+<span class="n">tmax</span> <span class="p">=</span> <span class="mi">100</span><span class="p">;</span>
+<span class="n">tsteps</span> <span class="p">=</span> <span class="mi">201</span><span class="p">;</span>
+<span class="n">tlist</span> <span class="p">=</span> <span class="nb">linspace</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span> <span class="n">tmax</span><span class="p">,</span> <span class="n">tsteps</span><span class="p">);</span>
+<span class="n">a</span> <span class="p">=</span> <span class="n">destroy</span><span class="p">(</span><span class="n">Ncutoff</span><span class="p">);</span>
+<span class="n">ad</span> <span class="p">=</span> <span class="n">a</span><span class="o">&#39;</span><span class="p">;</span>
+<span class="n">H</span> <span class="p">=</span> <span class="n">delta</span><span class="o">*</span><span class="n">ad</span><span class="o">*</span><span class="n">a</span> <span class="o">+</span> <span class="n">eta</span><span class="o">*</span><span class="p">(</span><span class="n">a</span> <span class="o">+</span> <span class="n">ad</span><span class="p">);</span>
+<span class="n">C</span> <span class="p">=</span> <span class="nb">sqrt</span><span class="p">(</span><span class="mi">2</span><span class="o">*</span><span class="n">kap</span><span class="p">)</span><span class="o">*</span><span class="n">a</span><span class="p">;</span>
+<span class="n">LH</span> <span class="p">=</span> <span class="o">-</span><span class="mi">1</span><span class="n">i</span> <span class="o">*</span> <span class="p">(</span><span class="n">spre</span><span class="p">(</span><span class="n">H</span><span class="p">)</span> <span class="o">-</span> <span class="n">spost</span><span class="p">(</span><span class="n">H</span><span class="p">));</span>
+<span class="n">L1</span> <span class="p">=</span> <span class="n">spre</span><span class="p">(</span><span class="n">C</span><span class="p">)</span><span class="o">*</span><span class="n">spost</span><span class="p">(</span><span class="n">C</span><span class="o">&#39;</span><span class="p">)</span><span class="o">-</span><span class="mf">0.5</span><span class="o">*</span><span class="n">spre</span><span class="p">(</span><span class="n">C</span><span class="o">&#39;*</span><span class="n">C</span><span class="p">)</span><span class="o">-</span><span class="mf">0.5</span><span class="o">*</span><span class="n">spost</span><span class="p">(</span><span class="n">C</span><span class="o">&#39;*</span><span class="n">C</span><span class="p">);</span>
+<span class="n">L</span> <span class="p">=</span> <span class="n">LH</span><span class="o">+</span><span class="n">L1</span><span class="p">;</span>
+
+
+<span class="n">psi0</span> <span class="p">=</span> <span class="n">basis</span><span class="p">(</span><span class="n">Ncutoff</span><span class="p">,</span> <span class="mi">1</span><span class="p">);</span>
+<span class="n">rho0</span> <span class="p">=</span> <span class="n">psi0</span> <span class="o">*</span> <span class="n">psi0</span><span class="o">&#39;</span><span class="p">;</span>
+
+<span class="c">% Set up options, if required</span>
+<span class="n">options</span><span class="p">.</span><span class="n">reltol</span> <span class="p">=</span> <span class="mf">1e-6</span><span class="p">;</span>
+<span class="n">options</span><span class="p">.</span><span class="n">abstol</span> <span class="p">=</span> <span class="mf">1e-6</span><span class="p">;</span>
+<span class="c">% Write out the data file</span>
+<span class="n">ode2file</span><span class="p">(</span><span class="s">&#39;file1.dat&#39;</span><span class="p">,</span> <span class="n">L</span><span class="p">,</span> <span class="n">rho0</span><span class="p">,</span> <span class="n">tlist</span><span class="p">,</span> <span class="n">options</span><span class="p">);</span>
+<span class="c">% Call the equation solver</span>
+<span class="n">odesolve</span><span class="p">(</span><span class="s">&#39;file1.dat&#39;</span><span class="p">,</span><span class="s">&#39;file2.dat&#39;</span><span class="p">);</span>
+</pre></div>

benchmarks-sourcecode/timeevolution_master_py.html

@@ -0,0 +1,16 @@
+<div class="highlight"><pre><span></span><span class="n">kap</span> <span class="o">=</span> <span class="mf">1.</span>
+<span class="n">eta</span> <span class="o">=</span> <span class="mi">4</span><span class="o">*</span><span class="n">kap</span>
+<span class="n">delta</span> <span class="o">=</span> <span class="mi">0</span>
+<span class="n">tmax</span> <span class="o">=</span> <span class="mi">100</span>
+<span class="n">tsteps</span> <span class="o">=</span> <span class="mi">201</span>
+<span class="n">tlist</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">linspace</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span> <span class="n">tmax</span><span class="p">,</span> <span class="n">tsteps</span><span class="p">)</span>
+
+<span class="n">a</span> <span class="o">=</span> <span class="n">qt</span><span class="o">.</span><span class="n">destroy</span><span class="p">(</span><span class="n">Ncutoff</span><span class="p">)</span>
+<span class="n">ad</span> <span class="o">=</span> <span class="n">a</span><span class="o">.</span><span class="n">dag</span><span class="p">()</span>
+<span class="n">H</span> <span class="o">=</span> <span class="n">delta</span><span class="o">*</span><span class="n">ad</span><span class="o">*</span><span class="n">a</span> <span class="o">+</span> <span class="n">eta</span><span class="o">*</span><span class="p">(</span><span class="n">a</span> <span class="o">+</span> <span class="n">ad</span><span class="p">)</span>
+<span class="n">c_ops</span> <span class="o">=</span> <span class="p">[</span><span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="mi">2</span><span class="o">*</span><span class="n">kap</span><span class="p">)</span><span class="o">*</span><span class="n">a</span><span class="p">]</span>
+
+<span class="n">psi0</span> <span class="o">=</span> <span class="n">qt</span><span class="o">.</span><span class="n">basis</span><span class="p">(</span><span class="n">Ncutoff</span><span class="p">,</span> <span class="mi">0</span><span class="p">)</span>
+<span class="n">n</span> <span class="o">=</span> <span class="n">qt</span><span class="o">.</span><span class="n">mesolve</span><span class="p">(</span><span class="n">H</span><span class="p">,</span> <span class="n">psi0</span><span class="p">,</span> <span class="n">tlist</span><span class="p">,</span> <span class="n">c_ops</span><span class="p">,</span> <span class="p">[</span><span class="n">ad</span><span class="o">*</span><span class="n">a</span><span class="p">],</span> <span class="n">options</span><span class="o">=</span><span class="n">options</span><span class="p">)</span><span class="o">.</span><span class="n">expect</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span>
+<span class="k">return</span> <span class="n">n</span>
+</pre></div>

benchmarks-sourcecode/timeevolution_particle.jl

@@ -0,0 +1,21 @@
+xmin = -10
+xmax = 10
+
+x0 = 2
+p0 = 1
+sigma0 = 1
+
+
+# Setup
+bx = PositionBasis(xmin, xmax, Npoints)
+x = position(bx)
+p = momentum(bx)
+H = p^2 + full(2*x^2)
+psi0 = gaussianstate(bx, x0, p0, sigma0)
+
+
+# Benchmark
+T = [0:1.:10;]
+exp_x = Float64[]
+fout(t, psi) = push!(exp_x, real(expect(x, psi)))
+timeevolution.schroedinger(T, psi0, H; fout=fout, reltol=1e-6, abstol=1e-8)

benchmarks-sourcecode/timeevolution_particle.py

@@ -0,0 +1,41 @@
+xmin = -10
+xmax = 10
+
+x0 = 2
+p0 = 1
+sigma0 = 1
+
+
+# Setup
+dx = (xmax - xmin)/Npoints
+pmin = -np.pi/dx
+pmax = np.pi/dx
+dp = (pmax - pmin)/Npoints
+
+samplepoints_x = np.linspace(xmin, xmax, Npoints, endpoint=False)
+samplepoints_p = np.linspace(pmin, pmax, Npoints, endpoint=False)
+
+x = qt.Qobj(np.diag(samplepoints_x))
+row0 = [sum([p*np.exp(-1j*p*dxji) for p in samplepoints_p])*dp*dx/(2*np.pi) for dxji in samplepoints_x - xmin]
+row0 = np.array(row0)
+col0 = row0.conj()
+
+a = np.zeros([Npoints, Npoints], dtype=complex)
+for i in range(Npoints):
+    a[i, i:] = row0[:Npoints-i]
+    a[i:, i] = col0[:Npoints-i]
+p = qt.Qobj(a)
+
+H = p**2 + 2*x**2
+
+def gaussianstate(x0, p0, sigma):
+    alpha = 1./(np.pi**(1/4)*np.sqrt(sigma))*np.sqrt(dx)
+    data = alpha*np.exp(1j*p0*(samplepoints_x-x0) - (samplepoints_x-x0)**2/(2*sigma**2))
+    return qt.Qobj(data)
+
+psi0 = gaussianstate(2., 1., 1.)
+
+
+# Benchmark
+tlist = np.linspace(0, 10, 11)
+exp_x = qt.mesolve(H, psi0, tlist, [], [x], options=options).expect[0]

benchmarks-sourcecode/timeevolution_particle_jl.html

@@ -0,0 +1,22 @@
+<div class="highlight"><pre><span></span><span class="n">xmin</span> <span class="o">=</span> <span class="o">-</span><span class="mi">10</span>
+<span class="n">xmax</span> <span class="o">=</span> <span class="mi">10</span>
+
+<span class="n">x0</span> <span class="o">=</span> <span class="mi">2</span>
+<span class="n">p0</span> <span class="o">=</span> <span class="mi">1</span>
+<span class="n">sigma0</span> <span class="o">=</span> <span class="mi">1</span>
+
+
+<span class="c"># Setup</span>
+<span class="n">bx</span> <span class="o">=</span> <span class="n">PositionBasis</span><span class="p">(</span><span class="n">xmin</span><span class="p">,</span> <span class="n">xmax</span><span class="p">,</span> <span class="n">Npoints</span><span class="p">)</span>
+<span class="n">x</span> <span class="o">=</span> <span class="n">position</span><span class="p">(</span><span class="n">bx</span><span class="p">)</span>
+<span class="n">p</span> <span class="o">=</span> <span class="n">momentum</span><span class="p">(</span><span class="n">bx</span><span class="p">)</span>
+<span class="n">H</span> <span class="o">=</span> <span class="n">p</span><span class="o">^</span><span class="mi">2</span> <span class="o">+</span> <span class="n">full</span><span class="p">(</span><span class="mi">2</span><span class="o">*</span><span class="n">x</span><span class="o">^</span><span class="mi">2</span><span class="p">)</span>
+<span class="n">psi0</span> <span class="o">=</span> <span class="n">gaussianstate</span><span class="p">(</span><span class="n">bx</span><span class="p">,</span> <span class="n">x0</span><span class="p">,</span> <span class="n">p0</span><span class="p">,</span> <span class="n">sigma0</span><span class="p">)</span>
+
+
+<span class="c"># Benchmark</span>
+<span class="n">T</span> <span class="o">=</span> <span class="p">[</span><span class="mi">0</span><span class="o">:</span><span class="mf">1.</span><span class="o">:</span><span class="mi">10</span><span class="p">;]</span>
+<span class="n">exp_x</span> <span class="o">=</span> <span class="kt">Float64</span><span class="p">[]</span>
+<span class="n">fout</span><span class="p">(</span><span class="n">t</span><span class="p">,</span> <span class="n">psi</span><span class="p">)</span> <span class="o">=</span> <span class="n">push!</span><span class="p">(</span><span class="n">exp_x</span><span class="p">,</span> <span class="n">real</span><span class="p">(</span><span class="n">expect</span><span class="p">(</span><span class="n">x</span><span class="p">,</span> <span class="n">psi</span><span class="p">)))</span>
+<span class="n">timeevolution</span><span class="o">.</span><span class="n">schroedinger</span><span class="p">(</span><span class="n">T</span><span class="p">,</span> <span class="n">psi0</span><span class="p">,</span> <span class="n">H</span><span class="p">;</span> <span class="n">fout</span><span class="o">=</span><span class="n">fout</span><span class="p">,</span> <span class="n">reltol</span><span class="o">=</span><span class="mf">1e-6</span><span class="p">,</span> <span class="n">abstol</span><span class="o">=</span><span class="mf">1e-8</span><span class="p">)</span>
+</pre></div>