In [5]:
import Potapov_Code.Roots as Roots
import Potapov_Code.Potapov as Potapov
import Potapov_Code.Time_Delay_Network as Time_Delay_Network
import Potapov_Code.Time_Sims as Time_Sims
import Potapov_Code.functions as functions
import Potapov_Code.tests as tests
import numpy as np
import numpy.linalg as la
from scipy.integrate import ode
import matplotlib.pyplot as plt
%matplotlib inline
In [6]:
def contour_plot(Mat,func = abs):
'''
Make a simple plot to view a matrix
Args:
Mat (compelx-valued matrix): a matrix to view
func (optional[function]): a function to apply to each component
Generates a plot of the matrix.
'''
fig = plt.figure()
ax = fig.add_subplot(111)
cax = ax.matshow(func(Mat), interpolation='nearest')
fig.colorbar(cax)
plt.show()
In [400]:
Ex = Time_Delay_Network.Example1(r = 0.7, max_linewidth=35.)
Ex.run_Potapov()
E = Ex.E
roots = Ex.roots
M1 = Ex.M1
delays = Ex.delays
modes = functions.spatial_modes(roots,M1,E)
In [401]:
vecs = Ex.vecs
In [402]:
roots.sort(key = lambda z: z.imag)
In [403]:
roots
Out[403]:
In [463]:
root_z = lambda z: 1j*z ## a fake root we will vary
In [471]:
roots_to_use = lambda z: [root_z(z).imag,roots[0].imag,roots[-1].imag]
In [472]:
modes_to_use = modes
In [473]:
delay_indices = 0
In [474]:
plus_or_minus_arr = [1,1,1]
In [488]:
x = np.linspace(-10000, 10000,4000)
In [489]:
reload(Potapov_Code.functions)
Out[489]:
In [490]:
f = lambda z: functions.make_nonlinear_interaction(roots_to_use(z),modes_to_use,Ex.delays,0,
0,0.01,plus_or_minus_arr)
In [491]:
plt.plot(x, [abs(f(z))**2 for z in x])
Out[491]:
In [492]:
Ex = Time_Delay_Network.Example3(r1 = 0.7, r3 = 0.7, max_linewidth=35.)
Ex.run_Potapov()
E = Ex.E
roots = Ex.roots
M1 = Ex.M1
delays = Ex.delays
modes = functions.spatial_modes(roots,M1,E)
In [493]:
vecs = Ex.vecs
In [494]:
roots.sort(key = lambda z: z.imag)
In [495]:
roots
Out[495]:
In [496]:
roots[-1]+roots[-2]
Out[496]:
In [497]:
root_z = lambda z: 1j*z ## a fake root we will vary
In [563]:
roots_to_use = lambda z: [root_z(z).imag,roots[0].imag,roots[1].imag]
In [564]:
modes_to_use = [modes[0],modes[0], modes[0]]
In [565]:
delay_indices = 0
In [566]:
plus_or_minus_arr = [1,1,1]
In [582]:
x = np.linspace(-800,1000,4000)
In [583]:
reload(Potapov_Code.functions)
Out[583]:
In [584]:
f = lambda z: functions.make_nonlinear_interaction(roots_to_use(z),modes_to_use,Ex.delays,0,
0,0.01,plus_or_minus_arr)
In [585]:
### by default, indices_of_refraction=[1,1,1]
We will vary the fake root we introduced to obtain the phase-mismatch diagram. That is, the phase mismatch $\delta k$ is going to be some linear function of $z$.
In [586]:
plt.plot(x, [abs(f(z))**2 for z in x])
Out[586]:
What happens when we change the indices of refraction for the different modes? The phase-mismatch will shift depending on where the new $\delta k = 0$ occurs. The width of the peak may also change if the indices of refraction are large.
In [605]:
indices_of_refraction=[3.,5.,10.]
In [606]:
f = lambda z: functions.make_nonlinear_interaction(roots_to_use(z),
modes_to_use,Ex.delays,0,0,0.1,plus_or_minus_arr,indices_of_refraction=indices_of_refraction)
In [607]:
plt.plot(x, [abs(f(z))**2 for z in x])
Out[607]:
In this section we will use example 3 to generate a Hamiltonian with nonlinaer coefficients resulting when inserting a nonlinearity in a circuit. We will assume that the nonlinearity is inserting at the delay line of index 0 corresponding to $\tau_1$.
In [611]:
import sympy as sp
import itertools
from qnet.algebra.circuit_algebra import *
In [612]:
Ex = Time_Delay_Network.Example3(r1 = 0.9, r3 = 0.9, max_linewidth=35.,max_freq=25.)
Ex.run_Potapov()
E = Ex.E
roots = Ex.roots
M1 = Ex.M1
delays = Ex.delays
modes = functions.spatial_modes(roots,M1,E)
In [613]:
roots
Out[613]:
In [614]:
## nonlinearity information
delay_index = 0
start_nonlin = 0.
duration_nonlin = .1
In [806]:
NONLIN_WEIGHT = 10.
In [615]:
m = len(roots)
In [616]:
indices = range(m)
In [617]:
chi_order = 3 ## i.e. chi-3 nonlinearity
In [618]:
plus_minus_combinations = list(itertools.combinations(range(chi_order + 1), 2)) ## pick which fields are annihilated
In [619]:
list_of_pm_arr = []
for tup in plus_minus_combinations:
ls = [1]*(chi_order+1)
for i in tup:
ls[i]=-1
list_of_pm_arr.append(ls)
In [620]:
a = [sp.symbols('a_'+str(i)) for i in range(m)]
a_H = [sp.symbols('a^H_'+str(i)) for i in range(m)]
In [ ]:
In [621]:
A,B,C,D = Potapov.get_Potapov_ABCD(Ex.roots,Ex.vecs,Ex.T,z=0.)
In [622]:
#Omega = (A-A.H)/(2j) #### closed dynamics only. i.e. not damping
In [623]:
Omega = -1j*A ## full dynamics
In [624]:
H_lin_sp = 0
## with sympy only
for i in range(m):
for j in range(m):
H_lin_sp += a_H[i]*a[j]*Omega[i,j]
In [625]:
def make_nonlin_term_sp(combination,pm_arr):
'''
Make symbolic term
With sympy only
'''
r = 1
for index,sign in zip(combination,pm_arr):
if sign == 1:
r*= a_H[index]
else:
r *= a[index]
return r
Let's impose a large 'index of refraction'. In the future we will replaces this by better conditions for phase-mismatch, including realistic values. For now, this will narrow the gain versus $\Delta k$ function so that few interaction terms remain.
In [626]:
def weight(combination,pm_arr):
roots_to_use = np.array([roots[i].imag for i in combination])
modes_to_use = [modes[i] for i in combination]
return functions.make_nonlinear_interaction(roots_to_use, modes_to_use, delays, delay_indices,
start_nonlin,duration_nonlin,pm_arr,
indices_of_refraction = [1000.]*len(combination),
eps=1e-12,)
In [627]:
## TODO: add a priori check to restrict exponential growth
weights = {}
In [628]:
count = 0
for pm_arr in list_of_pm_arr:
field_combinations = itertools.combinations_with_replacement(range(m), chi_order+1)
for combination in field_combinations:
count += 1
weights[tuple(combination),tuple(pm_arr)] = weight(combination,pm_arr)
print count
In [629]:
plt.hist([abs(x) for x in [weights[key] for key in weights] ],bins=100);
As we see above, most of the interactions are negligible. Let's drop them out.
In [630]:
significant_weight_keys = [key for key in weights if abs(weights[key]) > 1e-4]
In [631]:
significant_weights = dict((key,weights[key]) for key in significant_weight_keys)
In [632]:
significant_weights = {k:v for k,v in weights.iteritems() if abs(v) > 1e-4} ## more elegant
In [633]:
len(significant_weights)
Out[633]:
In [782]:
H_nonlin_sp = 0 ## with sympy only
for combination,pm_arr in significant_weights:
H_nonlin_sp += make_nonlin_term_sp(combination,pm_arr)*significant_weights[combination,pm_arr]
In [807]:
H_sp = H_lin_sp + H_nonlin_sp*NONLIN_WEIGHT
In [808]:
def make_sp_conj(A):
'''
Returns the symbolic conjugate of A.
Args:
A (symbolic expression in symbols a[i] and a_H[i])
Returns:
The complex conjugate of A
'''
A_H = sp.conjugate(H_sp)
for i in range(len(a)):
A_H = A_H.subs(sp.conjugate(a[i]),a_H[i])
A_H = A_H.subs(sp.conjugate(a_H[i]),a[i])
return A_H
In [809]:
def make_eq_motion(H_sp):
'''
Input is a tuple or list, output is a matrix vector
'''
A_H = make_sp_conj(H_sp)
diff_ls = [1j*sp.diff(H_sp,var) for var in a_H] + [-1j*sp.diff(A_H,var) for var in a]
fs = [sp.lambdify( tuple(a+a_H),expression) for expression in diff_ls ]
return lambda arr: (np.asmatrix([ f(* arr ) for f in fs])).T
In [810]:
A_H = sp.conjugate(H_sp)
for i in range(len(a)):
A_H = A_H.subs(sp.conjugate(a[i]),a_H[i])
A_H = A_H.subs(sp.conjugate(a_H[i]),a[i])
In [811]:
eq_mot = make_eq_motion(H_sp)
In [812]:
def double_up(M1,M2=None):
if M2 == None:
M2 = np.zeros_like(M1)
top = np.hstack([M1,M2])
bottom = np.hstack([np.conj(M2),np.conj(M1)])
return np.vstack([top,bottom])
In [813]:
A_d,C_d,D_d = map(double_up,(A,C,D))
In [814]:
B_d = -double_up(C.H)
In [815]:
def make_f(eq_mot,B,a_in):
'''
Nonlinear equations of motion
'''
return lambda t,a: np.asarray(eq_mot(a)+B*a_in(t)).T[0]
In [816]:
def make_f_lin(A,B,a_in):
'''
Linear equations of motion
'''
return lambda t,a: np.asarray(A*np.asmatrix(a).T+B*a_in(t)).T[0]
In [817]:
a_in = lambda t: np.asmatrix([1.]*4).T
In [818]:
f = make_f(eq_mot,B_d,a_in)
In [819]:
f_lin = make_f_lin(A_d,B_d,a_in)
In [820]:
eq_res = eq_mot([1.]*10)
print eq_res
In [821]:
mat_res = A_d*np.asmatrix([1.]*10).T
print mat_res
In [822]:
## compute L2 error between
np.sqrt(sum(np.asarray(abs(eq_res - mat_res))**2))
Out[822]:
In [845]:
r = ode(f).set_integrator('zvode', method='bdf')
r_lin = ode(f_lin).set_integrator('zvode', method='bdf')
In [846]:
y0 = np.asmatrix([0.]*10).T
t0=0.
In [847]:
r.set_initial_value(y0, t0)
r_lin.set_initial_value(y0, t0)
Out[847]:
In [848]:
t1 = 100
dt = 0.01
In [849]:
Y = []
while r.successful() and r.t < t1:
r.integrate(r.t+dt)
u = a_in(r.t)
Y.append(C_d*r.y+D_d*u)
In [850]:
Y_lin = []
while r_lin.successful() and r_lin.t < t1:
r_lin.integrate(r_lin.t+dt)
u = a_in(r_lin.t)
Y_lin.append(C_d*r_lin.y+D_d*u)
In [851]:
for i in range(4):
plt.plot([(y).real[i][0,0] for y in Y ])
for i in range(4):
plt.plot([(y).real[i][0,0] for y in Y_lin ])
Using the symbolic packages is kind of slow. For classical simulations maybe we can avoid that. We just need to extract the equations of motion, which should end up being sparse in the interaction terms.
TODO: implement without sympy, e.g. with Julia
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [206]:
roots_to_use = np.array([roots[i] for i in combination])
modes_to_use = [modes[i] for i in combination]
In [207]:
def call_make_non_lin():
return functions.make_nonlinear_interaction(roots_to_use, modes_to_use, delays, delay_indices,
start_nonlin,duration_nonlin,pm_arr,
indices_of_refraction,
eps=1e-12,func=lambda z : z.imag)
In [208]:
call_make_non_lin()
Out[208]:
In [209]:
indices_of_refraction = 1000.
call_make_non_lin()
Out[209]:
In [210]:
start_nonlin = -1 ## this shouldn't happen
call_make_non_lin()
In [211]:
start_nonlin = [1]*len(roots_to_use)
call_make_non_lin()
In [212]:
start_nonlin = 0.00001
duration_nonlin = .099
call_make_non_lin()
Out[212]:
In [213]:
start_nonlin = [0.00001]*len(roots_to_use)
duration_nonlin = .099
call_make_non_lin()
Out[213]:
In [131]:
## consolidated weights do not take into account which modes are createad or annihilated.
consolidated_weightes = {}
for key in significant_weights:
if not key[0] in consolidated_weightes:
consolidated_weightes[key[0]] = significant_weights[key]
else:
consolidated_weightes[key[0]] += significant_weights[key]
Out[131]:
In [ ]:
## QNET annihilation and creation operators
a_ = [Destroy(local_space('fock', namespace = str(i))) for i in range(m)]
In [ ]:
## Make linear Hamiltonian with QNET
H_lin = sum([a_[i].dag()*a_[i]*Omega[i,i] for i in range(m)]) ## with QNET
def make_nonlin_term(combination,pm_arr):
'''
Make symbolic term
With QNET
'''
r = 1
for index,sign in zip(combination,pm_arr):
if sign == 1:
r*= a_[index].dag()
else:
r *= a_[index]
return r
In [ ]:
## Make nonlinear Hamiltonian in QNET
H_nonlin = 0 ## with QNET
for combination,pm_arr in significant_weights:
H_nonlin += make_nonlin_term(combination,pm_arr)*significant_weights[combination,pm_arr]
H_qnet = H_lin+H_nonlin