sign error (#39)

* sign error

cf. Kohler 2017 Eqns. S9-14

* units

* Update travis for using services xvfb, python version

* Propegate changes to cuda implementation

* Have a target.py that actually runs, even if it isn't as clean as I'd like

Co-authored-by: ddkohler <[email protected]>
Co-authored-by: Kyle Sunden <[email protected]>

Author: ddkohler

.travis.yml

@@ -1,7 +1,7 @@
 language: python
 python:
-  - "3.5"
-  - "3.6"
+  - "3.7"
+  - "3.8"
 addons:
     apt:
         packages:
@@ -12,8 +12,6 @@ addons:
 install:
   - pip install -e .
   - pip install -U pytest
-before_script:
-  - "export DISPLAY=:99.0"
-  - "sh -e /etc/init.d/xvfb start"
-  - sleep 3 # give xvfb some time to start
+services:
+  - xvfb
 script: pytest

WrightSim/experiment/_pulse.py

@@ -113,7 +113,7 @@ def pulse(cls, eparams, pm=None):
             cc = np.ones((eparams.shape[-1]))
         else:
             cc = np.sign(pm)
-        x = np.exp(-1j*(cc[:,None]*(freq[:,None]*(t[None,:] - mu[:,None])+p[:,None])))
+        x = np.exp(1j*(cc[:,None]*(freq[:,None]*(t[None,:] - mu[:,None])+p[:,None])))
         x*= y[:,None] * np.exp(-(t[None,:] - mu[:,None])**2 / (2*sigma[:,None]**2) )
         return t, x
 

WrightSim/hamiltonian/default.py

@@ -1,8 +1,9 @@
-
 import numpy as np
 
 from ..mixed import propagate
 
+wn_to_omega = 2 * np.pi * 3 * 10 ** -5
+
 
 class Hamiltonian:
     cuda_struct = """
@@ -23,14 +24,32 @@ class Hamiltonian:
         int* recorded_indices;
     };
     """
-    cuda_mem_size = 4*4 + np.intp(0).nbytes*6
-
-
-    def __init__(self, rho=None, tau=None, mu=None,
-                        omega=None, w_central=7000., coupling=0,
-                        propagator=None, phase_cycle=False,
-                        labels=['00','01 -2','10 2\'','10 1','20 1+2\'','11 1-2','11 2\'-2', '10 1-2+2\'', '21 1-2+2\''],
-                        time_orderings=list(range(1,7)), recorded_indices = [7, 8]):
+    cuda_mem_size = 4 * 4 + np.intp(0).nbytes * 6
+
+    def __init__(
+        self,
+        rho=None,
+        tau=None,
+        mu=None,
+        omega=None,
+        w_central=7000.0,
+        coupling=0,
+        propagator=None,
+        phase_cycle=False,
+        labels=[
+            "00",
+            "01 -2",
+            "10 2'",
+            "10 1",
+            "20 1+2'",
+            "11 1-2",
+            "11 2'-2",
+            "10 1-2+2'",
+            "21 1-2+2'",
+        ],
+        time_orderings=list(range(1, 7)),
+        recorded_indices=[7, 8],
+    ):
         """Create a Hamiltonian object.
 
         Parameters
@@ -80,32 +99,33 @@ def __init__(self, rho=None, tau=None, mu=None,
         """
         if rho is None:
             self.rho = np.zeros(len(labels), dtype=np.complex128)
-            self.rho[0] = 1.
+            self.rho[0] = 1.0
         else:
             self.rho = np.array(rho, dtype=np.complex128)
 
         if tau is None:
-            tau = 50. #fs
+            tau = 50.0  # fs
         if np.isscalar(tau):
             self.tau = np.array([np.inf, tau, tau, tau, tau, np.inf, np.inf, tau, tau])
         else:
             self.tau = tau
 
-        #TODO: Think about dictionaries or some other way of labeling Mu values
+        # TODO: Think about dictionaries or some other way of labeling Mu values
         if mu is None:
-            self.mu = np.array([1., 1.], dtype=np.complex128)
+            self.mu = np.array([1.0, 1.0], dtype=np.complex128)
         else:
             self.mu = np.array(mu, dtype=np.complex128)
 
         if omega is None:
             w_ag = w_central
             w_2aa = w_ag - coupling
-            w_2ag = 2*w_ag - coupling
-            w_gg = 0.
+            w_2ag = 2 * w_ag - coupling
+            w_gg = 0.0
             w_aa = w_gg
             self.omega = np.array([w_gg, -w_ag, w_ag, w_ag, w_2ag, w_aa, w_aa, w_ag, w_2aa])
+            self.omega *= wn_to_omega
         else:
-            self.omega = w_0
+            self.omega = omega
 
         if propagator is None:
             self.propagator = propagate.runge_kutta
@@ -116,23 +136,28 @@ def __init__(self, rho=None, tau=None, mu=None,
         self.recorded_indices = recorded_indices
 
         self.time_orderings = time_orderings
-        self.Gamma = 1./self.tau
+        self.Gamma = 1.0 / self.tau
+
+    @property
+    def omega_wn(self):
+        return self.omega / wn_to_omega
 
     def to_device(self, pointer):
         """Transfer the Hamiltonian to a C struct in CUDA device memory.
 
         Currently expects a pointer to an already allocated chunk of memory.
         """
         import pycuda.driver as cuda
-        #TODO: Reorganize to allocate here and return the pointer, this is more friendly
+
+        # TODO: Reorganize to allocate here and return the pointer, this is more friendly
         # Transfer the arrays which make up the hamiltonian
         rho = cuda.to_device(self.rho)
         mu = cuda.to_device(self.mu)
         omega = cuda.to_device(self.omega)
         Gamma = cuda.to_device(self.Gamma)
 
         # Convert time orderings into a C boolean array of 1 and 0, offset by one
-        tos = [1 if i in self.time_orderings else 0 for i in range(1,7)]
+        tos = [1 if i in self.time_orderings else 0 for i in range(1, 7)]
 
         # Transfer time orderings and recorded indices
         time_orderings = cuda.to_device(np.array(tos, dtype=np.int8))
@@ -141,7 +166,7 @@ def to_device(self, pointer):
         # Transfer metadata about the lengths of feilds
         cuda.memcpy_htod(pointer, np.array([len(self.rho)], dtype=np.int32))
         cuda.memcpy_htod(int(pointer) + 4, np.array([len(self.mu)], dtype=np.int32))
-        #TODO: generalize nTimeOrderings
+        # TODO: generalize nTimeOrderings
         cuda.memcpy_htod(int(pointer) + 8, np.array([6], dtype=np.int32))
         cuda.memcpy_htod(int(pointer) + 12, np.array([len(self.recorded_indices)], dtype=np.int32))
 
@@ -153,8 +178,6 @@ def to_device(self, pointer):
         cuda.memcpy_htod(int(pointer) + 48, np.intp(int(time_orderings)))
         cuda.memcpy_htod(int(pointer) + 56, np.intp(int(recorded_indices)))
 
-        
-
     def matrix(self, efields, time):
         """Generate the time dependant Hamiltonian Coupling Matrix.
 
@@ -172,9 +195,9 @@ def matrix(self, efields, time):
             Shape T x N x N array with the full Hamiltonian at each time step.
             N is the number of states in the Density vector.
         """
-        #TODO: Just put the body of this method in here, rather than calling _gen_matrix
-        E1,E2,E3 = efields[0:3]
-        return  self._gen_matrix(E1, E2, E3, time, self.omega)
+        # TODO: Just put the body of this method in here, rather than calling _gen_matrix
+        E1, E2, E3 = efields[0:3]
+        return self._gen_matrix(E1, E2, E3, time, self.omega)
 
     def _gen_matrix(self, E1, E2, E3, time, energies):
         """
@@ -187,13 +210,13 @@ def _gen_matrix(self, E1, E2, E3, time, energies):
         outside of the matrix
         """
         # Define transition energies
-        wag  = energies[1]
+        wag = energies[1]
         w2aa = energies[-1]
-        
+
         # Define dipole moments
         mu_ag = self.mu[0]
         mu_2aa = self.mu[-1]
-    
+
         # Define helpful variables
         A_1 = 0.5j * mu_ag * E1 * np.exp(-1j * wag * time)
         A_2 = 0.5j * mu_ag * E2 * np.exp(1j * wag * time)
@@ -207,46 +230,46 @@ def _gen_matrix(self, E1, E2, E3, time, energies):
 
         # Add appropriate array elements, according to the time orderings
         if 3 in self.time_orderings or 5 in self.time_orderings:
-            out[:,1,0] = -A_2
+            out[:, 1, 0] = -A_2
         if 4 in self.time_orderings or 6 in self.time_orderings:
-            out[:,2,0] = A_2prime
+            out[:, 2, 0] = A_2prime
         if 1 in self.time_orderings or 2 in self.time_orderings:
-            out[:,3,0] = A_1
+            out[:, 3, 0] = A_1
         if 3 in self.time_orderings:
-            out[:,5,1] = A_1
+            out[:, 5, 1] = A_1
         if 5 in self.time_orderings:
-            out[:,6,1] = A_2prime
+            out[:, 6, 1] = A_2prime
         if 4 in self.time_orderings:
-            out[:,4,2] = B_1
+            out[:, 4, 2] = B_1
         if 6 in self.time_orderings:
-            out[:,6,2] = -A_2
+            out[:, 6, 2] = -A_2
         if 2 in self.time_orderings:
-            out[:,4,3] = B_2prime
+            out[:, 4, 3] = B_2prime
         if 1 in self.time_orderings:
-            out[:,5,3] = -A_2
+            out[:, 5, 3] = -A_2
         if 2 in self.time_orderings or 4 in self.time_orderings:
-            out[:,7,4] = B_2
-            out[:,8,4] = -A_2
+            out[:, 7, 4] = B_2
+            out[:, 8, 4] = -A_2
         if 1 in self.time_orderings or 3 in self.time_orderings:
-            out[:,7,5] = -2 * A_2prime
-            out[:,8,5] = B_2prime
+            out[:, 7, 5] = -2 * A_2prime
+            out[:, 8, 5] = B_2prime
         if 5 in self.time_orderings or 6 in self.time_orderings:
-            out[:,7,6] = -2 * A_1
-            out[:,8,6] = B_1
+            out[:, 7, 6] = -2 * A_1
+            out[:, 8, 6] = B_1
 
         # Add Gamma along the diagonal
         for i in range(len(self.Gamma)):
-            out[:,i,i] = -1 * self.Gamma[i]
+            out[:, i, i] = -1 * self.Gamma[i]
 
-        #NOTE: NISE multiplied outputs by the approriate mu in here
+        # NOTE: NISE multiplied outputs by the approriate mu in here
         #      This mu factors out, remember to use it where needed later
         #      Removed for clarity and aligning with Equation S15 of Kohler2017
 
         return out
 
     cuda_matrix_source = """
     /**
-     *  Hamiltonian_matrix: Computes the Hamiltonian matrix for an indevidual time step.
+     *  Hamiltonian_matrix: Computes the Hamiltonian matrix for an individual time step.
      *  NOTE: This differs from the Python implementation, which computes the full time 
      *          dependant hamiltonian, this only computes for a single time step
      *          (to conserve memory).
@@ -330,5 +353,4 @@ def _gen_matrix(self, E1, E2, E3, time, energies):
         // Put Gamma along the diagonal
         for(int i=0; i<ham.nStates; i++) out[i*ham.nStates + i] = -1. * ham.Gamma[i];
     }
-"""    
-
+"""

WrightSim/mixed/propagate.py

@@ -1,5 +1,6 @@
 import numpy as np
 
+
 def runge_kutta(t, efields, n_recorded, hamiltonian):
     """ Evolves the hamiltonian in time using the runge_kutta method.
 
@@ -23,7 +24,7 @@ def runge_kutta(t, efields, n_recorded, hamiltonian):
         2-D array of recorded density vector elements for each time step in n_recorded.
     """
     # can only call on n_recorded and t after efield_object.E is called
-    dt = np.abs(t[1]-t[0])
+    dt = np.abs(t[1] - t[0])
     # extract attributes of the system
     rho_emitted = np.empty((len(hamiltonian.recorded_indices), n_recorded), dtype=np.complex128)
 
@@ -32,27 +33,28 @@ def runge_kutta(t, efields, n_recorded, hamiltonian):
     # index to keep track of elements of rho_emitted
     emitted_index = 0
     rho_i = hamiltonian.rho.copy()
-    for k in range(len(t)-1):
+    for k in range(len(t) - 1):
         # calculate delta rho based on previous rho values
         temp_delta_rho = np.dot(H[k], rho_i)
-        temp_rho_i = rho_i + temp_delta_rho*dt
+        temp_rho_i = rho_i + temp_delta_rho * dt
         # second order term
-        delta_rho = np.dot(H[k+1], temp_rho_i)
-        rho_i += dt/2 * (temp_delta_rho + delta_rho)
-        # if we are close enough to final coherence emission, start 
+        delta_rho = np.dot(H[k + 1], temp_rho_i)
+        rho_i += dt / 2 * (temp_delta_rho + delta_rho)
+        # if we are close enough to final coherence emission, start
         # storing these values
         if k >= len(t) - n_recorded:
-            for rec,native in enumerate(hamiltonian.recorded_indices):
+            for rec, native in enumerate(hamiltonian.recorded_indices):
                 rho_emitted[rec, emitted_index] = rho_i[native]
             emitted_index += 1
     # Last timestep
     temp_delta_rho = np.dot(H[-1], rho_i)
-    rho_i += temp_delta_rho*dt
-    for rec,native in enumerate(hamiltonian.recorded_indices):
+    rho_i += temp_delta_rho * dt
+    for rec, native in enumerate(hamiltonian.recorded_indices):
         rho_emitted[rec, emitted_index] = rho_i[native]
 
     return rho_emitted
 
+
 muladd_cuda_source = """
 /*
  * muladd: Linear combination of two vectors with two scalar multiples
@@ -139,7 +141,7 @@ def runge_kutta(t, efields, n_recorded, hamiltonian):
     for(int i=0; i < n; i++)
     {
         // Complex phase and magnitude
-        out[i] = pycuda::exp(-1. * I * ((double)(phase_matching[i]) *
+        out[i] = pycuda::exp(1. * I * ((double)(phase_matching[i]) *
                                         (params[3 + i*5] * (t - params[2 + i*5]) + params[4 + i*5])));
         // Gaussian envelope
         out[i] *= params[0 + i*5] * exp(-1 * (t-params[2 + i*5]) * (t-params[2 + i*5])

scripts/target.py

@@ -27,8 +27,8 @@
 dt = 50.  # pulse duration (fs)
 slitwidth = 120.  # mono resolution (wn)
 
-nw = sys.argv[1]#16  # number of frequency points (w1 and w2)
-nt = sys.argv[2]#1376  # number of delay points (d2)
+nw = 256  # number of frequency points (w1 and w2)
+nt = 256  # number of delay points (d2)
 
 
 # --- workspace -----------------------------------------------------------------------------------
@@ -41,26 +41,25 @@
 #exp.w2.points = 0.
 exp.d2.points = np.linspace(-2 * dt, 8 * dt, nt)
 exp.w1.active = exp.w2.active = exp.d2.active = True
-#exp.d2.points = 4 * dt
-#exp.d2.active = False
+exp.d2.points = 4 * dt
+exp.d2.active = False
 exp.timestep = 2.
 exp.early_buffer = 100.0
 exp.late_buffer  = 400.0
 
 # create hamiltonian
-ham = ws.hamiltonian.Hamiltonian(w_central=0.)
+ham = ws.hamiltonian.Hamiltonian(w_central=300.)
 #ham.time_orderings = [5]
 ham.recorded_elements = [7,8]
 
 # do scan
 begin = time.perf_counter()
-scan = exp.run(ham, mp='')
+scan = exp.run(ham, mp='cpu')
 print(time.perf_counter()-begin)
 gpuSig = scan.sig.copy()
 #with wt.kit.Timer():
 #    scan2 = exp.run(ham, mp=None)
 #cpuSig = scan2.sig.copy()
-'''
 plt.close('all')
 # measure and plot
 fig, gs = wt.artists.create_figure(cols=[1, 'cbar'])
@@ -78,4 +77,3 @@
 wt.artists.plot_colorbar(label='ampiltude')
 # finish
 plt.show()
-'''