<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Development notebook
Step1:  Photo-count detection
Theory
Stochastic Schrödinger equation for photocurrent detection, in Milburn's formulation, is
$$
d|\psi(t)\rangle = -iH|\psi(t)\rangle dt + \left(\frac{a}{\sqrt{\langle a^\dagger a \rangle_\psi(t)}} - 1\right) |\psi(t)\rangle dN(t) - \frac{1}{2}\gamma\left(\langle a^\dagger a \rangle_\psi(t) - a^\dagger a \right) |\psi(t)\rangle dt
$$
and $dN(t)$ is a Poisson distributed increment with $E[dN(t)] = \gamma \langle a^\dagger a\rangle (t)$.
Formulation in QuTiP
In QuTiP we write the stochastic Schrödinger equation on the form
Step2:  Solve using stochastic master equation
$\displaystyle D_{1}[a, |\psi\rangle] = 
\frac{1}{2} ( \langle \psi|a^\dagger a |\psi\rangle  - a^\dagger a )|\psi\rangle
$
$\displaystyle D_{2}[a, |\psi\rangle] = 
\frac{a|\psi\rangle}{\sqrt{\langle \psi| a^\dagger a|\psi\rangle}} - |\psi\rangle
$
Using QuTiP built-in photo-current detection functions for $D_1$ and $D_2$
Step3:  Homodyne detection
Step4:  Theory
Stochastic master equation for homodyne can be written as (see Gardiner and Zoller, Quantum Noise)
$$
d|\psi(t)\rangle = -iH|\psi(t)\rangle dt + 
\frac{1}{2}\gamma \left(
\langle a + a^\dagger \rangle_\psi a  - a^\dagger a   - \frac{1}{4}\langle a + a^\dagger \rangle_\psi^2  \right)|\psi(t)\rangle dt
+ 
\sqrt{\gamma}\left(
a - \frac{1}{2} \langle a + a^\dagger \rangle_\psi
\right)|\psi(t)\rangle dW
$$
where $dW(t)$ is a normal distributed increment with $E[dW(t)] = \sqrt{dt}$.
In QuTiP format we have
Step5:  $$
D_{2}[c,|\psi\rangle] = 
\left(c  - \frac{1}{2} \langle c + c^\dagger \rangle_\psi \right)|\psi(t)\rangle
$$
Step6:  Solve problem again, this time with a specified noise (from previous run)
Step7:  Using QuTiP built-in homodyne detection functions for $D_1$ and $D_2$
Step8:  Solve problem again, this time with a specified noise (from previous run)
Step9:  Heterodyne detection
Stochastic Schrödinger equation for heterodyne detection can be written as
$$
d|\psi(t)\rangle = -iH|\psi(t)\rangle dt 
-\frac{1}{2}\gamma\left(a^\dagger a - \langle a^\dagger \rangle a +  \frac{1}{2}\langle a \rangle\langle a^\dagger \rangle)\right)|\psi(t)\rangle dt
+ 
\sqrt{\gamma/2}\left(
a - \frac{1}{2} \langle a + a^\dagger \rangle_\psi
\right)|\psi(t)\rangle dW_1
-i
\sqrt{\gamma/2}\left(
a - \frac{1}{2} \langle a - a^\dagger \rangle_\psi
\right)|\psi(t)\rangle dW_2
$$
where $dW_i(t)$ is a normal distributed increment with $E[dW_i(t)] = \sqrt{dt}$.
In QuTiP format we have
Step10:  $$D_{2}^{(1)}[c, |\psi(t)\rangle] = \sqrt{1/2} (c - \langle c + c^\dagger \rangle / 2) \psi$$
$$D_{2}^{(1)}[c, |\psi(t)\rangle] = -i\sqrt{1/2} (c - \langle c - c^\dagger \rangle /2) \psi$$
Step11:  Using QuTiP built-in heterodyne detection functions for $D_1$ and $D_2$
Step12:  Solve problem again, this time with a specified noise (from previous run)
Step13:  Software version
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
import matplotlib.pyplot as plt
from qutip import *
import numpy as np
Explanation: Development notebook: Tests for QuTiP's stochastic Schrödinger equation solver
Copyright (C) 2011 and later, Paul D. Nation & Robert J. Johansson
In this notebook we test the qutip stochastic Schrödinger equation solver (ssesolve) with a few examples. The same examples are solved using the stochastic master equation solver in the notebook development-smesolve-tests.
<style>
.rendered_html {
    font-family: Liberation Serif;
}
.rendered_html h1 {
    font-family: Liberation Sans;
    margin: 0 0;
}
.rendered_html h2 {
    font-family: Liberation Sans;
    margin: 0 0;
}
</style>
End of explanation
N = 10
w0 = 0.5 * 2 * np.pi
times = np.linspace(0, 15, 150)
dt = times[1] - times[0]
gamma = 0.25
A = 2.5
ntraj = 50
nsubsteps = 100
a = destroy(N)
x = a + a.dag()
H = w0 * a.dag() * a
psi0 = fock(N, 5)
sc_ops = [np.sqrt(gamma) * a]
e_ops = [a.dag() * a, a + a.dag(), (-1j)*(a - a.dag())]
result_ref = mesolve(H, psi0, times, sc_ops, e_ops)
plot_expectation_values(result_ref);
Explanation: Photo-count detection
Theory
Stochastic Schrödinger equation for photocurrent detection, in Milburn's formulation, is
$$
d|\psi(t)\rangle = -iH|\psi(t)\rangle dt + \left(\frac{a}{\sqrt{\langle a^\dagger a \rangle_\psi(t)}} - 1\right) |\psi(t)\rangle dN(t) - \frac{1}{2}\gamma\left(\langle a^\dagger a \rangle_\psi(t) - a^\dagger a \right) |\psi(t)\rangle dt
$$
and $dN(t)$ is a Poisson distributed increment with $E[dN(t)] = \gamma \langle a^\dagger a\rangle (t)$.
Formulation in QuTiP
In QuTiP we write the stochastic Schrödinger equation on the form:
$$
d|\psi(t)\rangle = -iH|\psi(t)\rangle dt + D_{1}[c] |\psi(t)\rangle dt  + D_{2}[c] |\psi(t)\rangle dW
$$
where $c = \sqrt{\gamma} a$ is the collapse operator including the rate of the process as a coefficient in the operator. We can identify
$$
D_{1}[c] |\psi(t)\rangle = - \frac{1}{2}\left(\langle c^\dagger c \rangle_\psi(t) - c^\dagger c \right)
$$
$$
D_{2}[c] |\psi(t)\rangle = \frac{c}{\sqrt{\langle c^\dagger c \rangle_\psi(t)}} - 1
$$
$$dW = dN(t)$$
Reference solution: deterministic master equation
End of explanation
result = photocurrent_sesolve(H, psi0, times, sc_ops, e_ops,
                  ntraj=ntraj*5, nsubsteps=nsubsteps, store_measurement=True, normalize=0)
plot_expectation_values([result, result_ref]);
for m in result.measurement:
    plt.step(times, dt * m.real)
Explanation: Solve using stochastic master equation
$\displaystyle D_{1}[a, |\psi\rangle] = 
\frac{1}{2} ( \langle \psi|a^\dagger a |\psi\rangle  - a^\dagger a )|\psi\rangle
$
$\displaystyle D_{2}[a, |\psi\rangle] = 
\frac{a|\psi\rangle}{\sqrt{\langle \psi| a^\dagger a|\psi\rangle}} - |\psi\rangle
$
Using QuTiP built-in photo-current detection functions for $D_1$ and $D_2$
End of explanation
ntraj = 50
nsubsteps = 200
H = w0 * a.dag() * a + A * (a + a.dag())
result_ref = mesolve(H, psi0, times, sc_ops, e_ops)
Explanation: Homodyne detection
End of explanation
op = sc_ops[0]
opp = (op + op.dag()).data
opn = (op.dag()*op).data
op = op.data
Hd = H.data * -1j
def d1_psi_func(t, psi):
    e_x = cy.cy_expect_psi(opp, psi, 0)
    return  cy.spmv(Hd, psi) + 0.5 * (e_x * cy.spmv(op, psi) - cy.spmv(opn, psi) - 0.25 * e_x ** 2 * psi)
Explanation: Theory
Stochastic master equation for homodyne can be written as (see Gardiner and Zoller, Quantum Noise)
$$
d|\psi(t)\rangle = -iH|\psi(t)\rangle dt + 
\frac{1}{2}\gamma \left(
\langle a + a^\dagger \rangle_\psi a  - a^\dagger a   - \frac{1}{4}\langle a + a^\dagger \rangle_\psi^2  \right)|\psi(t)\rangle dt
+ 
\sqrt{\gamma}\left(
a - \frac{1}{2} \langle a + a^\dagger \rangle_\psi
\right)|\psi(t)\rangle dW
$$
where $dW(t)$ is a normal distributed increment with $E[dW(t)] = \sqrt{dt}$.
In QuTiP format we have:
$$
d|\psi(t)\rangle = -iH|\psi(t)\rangle dt + D_{1}[c, |\psi(t)\rangle] dt  + D_{2}[c, |\psi(t)\rangle] dW
$$
where $c = \sqrt{\gamma} a$, so we can identify
$$
D_{1}[c, |\psi\rangle] = 
\frac{1}{2}\left(
\langle c + c^\dagger \rangle_\psi c   - c^\dagger c  - \frac{1}{4}\langle c + c^\dagger \rangle_\psi^2 
\right) |\psi(t)\rangle 
$$
End of explanation
def d2_psi_func(t, psi):
    out = np.zeros((1,len(psi)), dtype=complex)
    e_x = cy.cy_expect_psi(opp, psi, 0)
    out[0,:] = cy.spmv(op,psi)
    out -= 0.5 * e_x * psi
    return out
result = general_stochastic(psi0, times, d1=d1_psi_func, d2=d2_psi_func, 
                            e_ops=e_ops, ntraj=ntraj, 
                            m_ops=[a + a.dag()], dW_factors=[1/np.sqrt(gamma)],
                            nsubsteps=nsubsteps, store_measurement=True)
plot_expectation_values([result, result_ref]);
for m in result.measurement:
    plt.plot(times, m[:, 0].real, 'b', alpha=0.1)
plt.plot(times,  np.array(result.measurement).mean(axis=0)[:,0].real, 'r', lw=2)
    
plt.plot(times, result_ref.expect[1], 'k', lw=2)
plt.ylim(-15, 15);
Explanation: $$
D_{2}[c,|\psi\rangle] = 
\left(c  - \frac{1}{2} \langle c + c^\dagger \rangle_\psi \right)|\psi(t)\rangle
$$
End of explanation
result = general_stochastic(psi0, times, d1=d1_psi_func, d2=d2_psi_func, 
                            e_ops=e_ops, ntraj=ntraj, noise=result.noise,
                            m_ops=[a + a.dag()], dW_factors=[1/np.sqrt(gamma)],
                            nsubsteps=nsubsteps, store_measurement=True)
plot_expectation_values([result, result_ref]);
for m in result.measurement:
    plt.plot(times, m[:, 0].real, 'b', alpha=0.1)
plt.plot(times,  np.array(result.measurement).mean(axis=0)[:,0].real, 'r', lw=2)
    
plt.plot(times, result_ref.expect[1], 'k', lw=2)
plt.ylim(-15, 15);
Explanation: Solve problem again, this time with a specified noise (from previous run)
End of explanation
result = ssesolve(H, psi0, times, sc_ops, e_ops, ntraj=ntraj, nsubsteps=nsubsteps,
                  method='homodyne', store_measurement=True, dW_factors=[1])
plot_expectation_values([result, result_ref]);
for m in result.measurement:
    plt.plot(times, m[:, 0].real, 'b', alpha=0.1)
plt.plot(times,  np.array(result.measurement).mean(axis=0)[:,0].real/np.sqrt(gamma), 'r', lw=2)
    
plt.plot(times, result_ref.expect[1], 'k', lw=2)
plt.ylim(-15, 15);
plt.plot(times,  np.array(result.measurement).mean(axis=0)[:,0].real/np.sqrt(gamma), 'r', lw=2)
plt.plot(times, result_ref.expect[1], 'k', lw=2)
plt.plot(times, result.expect[1], 'b', lw=2)
Explanation: Using QuTiP built-in homodyne detection functions for $D_1$ and $D_2$
End of explanation
result = ssesolve(H, psi0, times, sc_ops, e_ops, ntraj=ntraj, nsubsteps=nsubsteps,
                  method='homodyne', store_measurement=True, noise=result.noise)
plot_expectation_values([result, result_ref]);
for m in result.measurement:
    plt.plot(times, m[:, 0].real, 'b', alpha=0.1)
plt.plot(times,  np.array(result.measurement).mean(axis=0)[:,0].real/np.sqrt(gamma), 'r', lw=2)
    
plt.plot(times, result_ref.expect[1], 'k', lw=2)
plt.ylim(-15, 15);
Explanation: Solve problem again, this time with a specified noise (from previous run)
End of explanation
op = sc_ops[0]
opd = (op.dag()).data
opp = (op + op.dag()).data
opm = (op + op.dag()).data
opn = (op.dag()*op).data
op = op.data
Hd = H.data * -1j
def d1_psi_func(t, psi):
    e_xd = cy.cy_expect_psi(opd, psi, 0)
    e_x = cy.cy_expect_psi(op, psi, 0)
    return  cy.spmv(Hd, psi) - 0.5 * (cy.spmv(opn, psi) - e_xd * cy.spmv(op, psi) + 0.5 * e_x * e_xd * psi)
Explanation: Heterodyne detection
Stochastic Schrödinger equation for heterodyne detection can be written as
$$
d|\psi(t)\rangle = -iH|\psi(t)\rangle dt 
-\frac{1}{2}\gamma\left(a^\dagger a - \langle a^\dagger \rangle a +  \frac{1}{2}\langle a \rangle\langle a^\dagger \rangle)\right)|\psi(t)\rangle dt
+ 
\sqrt{\gamma/2}\left(
a - \frac{1}{2} \langle a + a^\dagger \rangle_\psi
\right)|\psi(t)\rangle dW_1
-i
\sqrt{\gamma/2}\left(
a - \frac{1}{2} \langle a - a^\dagger \rangle_\psi
\right)|\psi(t)\rangle dW_2
$$
where $dW_i(t)$ is a normal distributed increment with $E[dW_i(t)] = \sqrt{dt}$.
In QuTiP format we have:
$$
d|\psi(t)\rangle = -iH|\psi(t)\rangle dt + D_{1}[c, |\psi(t)\rangle] dt  + \sum_n D_{2}^{(n)}[c, |\psi(t)\rangle] dW_n
$$
where $c = \sqrt{\gamma} a$, so we can identify
$$
D_1[c, |\psi(t)\rangle] =
-\frac{1}{2}\left(c^\dagger c - \langle c^\dagger \rangle c +  \frac{1}{2}\langle c \rangle\langle c^\dagger \rangle \right) |\psi(t)\rangle
$$
End of explanation
sqrt2 = np.sqrt(0.5)
def d2_psi_func(t, psi):
    out = np.zeros((2,len(psi)), dtype=complex)
    e_p = cy.cy_expect_psi(opp, psi, 0)
    e_m = cy.cy_expect_psi(opm, psi, 0)
    out[0,:] = (cy.spmv(op,psi) - e_p * 0.5 * psi)*sqrt2
    out[1,:] = (cy.spmv(op,psi) - e_m * 0.5 * psi)*sqrt2*-1j
    return out
result = general_stochastic(psi0, times, d1=d1_psi_func, d2=d2_psi_func, 
                            e_ops=e_ops, ntraj=ntraj, len_d2=2,
                            m_ops=[(a + a.dag()), (-1j)*(a - a.dag())], dW_factors=[2/np.sqrt(gamma), 2/np.sqrt(gamma)],
                            nsubsteps=nsubsteps, store_measurement=True)
plot_expectation_values([result, result_ref]);
#fig, ax = subplots()
for m in result.measurement:
    plt.plot(times, m[:, 0].real, 'r', alpha=0.025)
    plt.plot(times, m[:, 1].real, 'b', alpha=0.025)
    
plt.plot(times, result_ref.expect[1], 'k', lw=2);
plt.plot(times, result_ref.expect[2], 'k', lw=2);
plt.ylim(-10, 10)
plt.plot(times, np.array(result.measurement).mean(axis=0)[:,0].real, 'r', lw=2);
plt.plot(times, np.array(result.measurement).mean(axis=0)[:,1].real, 'b', lw=2);
Explanation: $$D_{2}^{(1)}[c, |\psi(t)\rangle] = \sqrt{1/2} (c - \langle c + c^\dagger \rangle / 2) \psi$$
$$D_{2}^{(1)}[c, |\psi(t)\rangle] = -i\sqrt{1/2} (c - \langle c - c^\dagger \rangle /2) \psi$$
End of explanation
result = ssesolve(H, psi0, times, sc_ops, e_ops, ntraj=ntraj, nsubsteps=nsubsteps,
                  method='heterodyne', store_measurement=True)
plot_expectation_values([result, result_ref]);
for m in result.measurement:
    plt.plot(times, m[:, 0, 0].real, 'r', alpha=0.025)
    plt.plot(times, m[:, 0, 1].real, 'b', alpha=0.025)
    
plt.plot(times, result_ref.expect[1], 'k', lw=2)
plt.plot(times, result_ref.expect[2], 'k', lw=2)
plt.plot(times, np.array(result.measurement).mean(axis=0)[:,0,0].real/np.sqrt(gamma), 'r', lw=2)
plt.plot(times, np.array(result.measurement).mean(axis=0)[:,0,1].real/np.sqrt(gamma), 'b', lw=2)
Explanation: Using QuTiP built-in heterodyne detection functions for $D_1$ and $D_2$
End of explanation
result = ssesolve(H, psi0, times, sc_ops, e_ops, ntraj=ntraj, nsubsteps=nsubsteps,
                  method='heterodyne', store_measurement=True, noise=result.noise)
plot_expectation_values([result, result_ref]);
for m in result.measurement:
    plt.plot(times, m[:, 0, 0].real, 'r', alpha=0.025)
    plt.plot(times, m[:, 0, 1].real, 'b', alpha=0.025)
    
plt.plot(times, result_ref.expect[1], 'k', lw=2);
plt.plot(times, result_ref.expect[2], 'k', lw=2);
plt.plot(times, np.array(result.measurement).mean(axis=0)[:,0,0].real/np.sqrt(gamma), 'r', lw=2);
plt.plot(times, np.array(result.measurement).mean(axis=0)[:,0,1].real/np.sqrt(gamma), 'b', lw=2);
Explanation: Solve problem again, this time with a specified noise (from previous run)
End of explanation
from qutip.ipynbtools import version_table
version_table()
Explanation: Software version
End of explanation
<END_TASK>