Cuda (#19)

* Add hamiltonian to __init__ and target.py

* Add runge-kutta propagator

* structure

* Make the scan run

* Profiling traces

* refactor efp

* trive.ini fix

* hamiltonian folder structure

also fixed unhelpful settings of target.py (set while debugging)

* default name

* TRSF hamiltonian

* Initial cuda implementation, Work in Progress

* More cuda device functions

* Kernel and memory transfer -- Compiles, but has memory access issue

* silence prints, Hamiltoinian POINTER, datatypes

* datatypes, remove memoryview, hardcoded mu (temp)

* datatypes, fix dot, statically allocated hamiltonian matrices, hard coded rho_0

* Remove extraneous prints, use sys.argv for number of points, use perf_counter for timing

* remove multiprocessing

* Remove import from propagate

* Scaling results

* Report

* README

Author: ksunden

README.rst

@@ -1,3 +1,66 @@
 WrightSim
 =========
 A simulation package for multidimensional spectroscopy.
+
+
+Installation
+------------
+
+.. code-block:: bash
+    
+    $ git clone https://github.com/wright-group/WrightSim
+    $ cd WrightSim
+    $ python setup.py develop
+
+Note: This will install all required dependencies.
+PyCUDA is not a required dependency in general, but is required if GPU simulations are desired.
+
+.. code-block:: bash
+
+    $ pip install pycuda
+
+PyCUDA requires an Nvida graphics card and drivers, and the CUDA libraries installed.
+
+Usage
+-----
+
+An example script is provided at ``./scripts/target.py``
+
+This script can be modified to suit an individual simulation.
+
+The basic steps are:
+
+#. Select an experiment
+#. Set up the axes of the scan
+#. Set the time interfal and buffers
+#. Create a Hamiltonian object
+#. Run the scan
+#. (optional) review the results
+
+Level of parallelism is selected by the ``mp`` parameter of teh ``exp.run`` method.
+
+- ``True`` or ``"cpu"`` enables CPU multiprocessing
+- ``"gpu"`` enables CUDA 
+- ``False`` or ``""`` runs in single threaded mode
+
+
+The script is set up to read the dimensions from arguments.
+it can be run like so:
+
+.. code-block:: bash
+    
+    $ ./target.py 32 16
+
+
+This will run a 3D simulation of 32x32x16 Freq-Freq-Delay.
+
+It will print the time it took to compute.
+
+An example of how to visualize results is provided at the end of the script, but is not active, as it only works for 2D scans.
+
+To convert this script to a 2D scan and visualize, uncomment the lines which set ``exp.d2.points`` and ``exp.d2.active`` and remove the triple quotes around the plotting routine.
+There it will ignore the second argument, producing a 32x32 2D frequency spectrum. (The second argument is still parsed, so it is required, just ignored.)
+
+
+
+

WrightSim/__init__.py

@@ -1,6 +1,7 @@
 #  flake8: noqa
 
 from . import experiment
+from . import hamiltonian
 from . import integration
 from . import measure
 from . import response

WrightSim/experiment/__init__.py

@@ -2,12 +2,11 @@
 
 
 import os
-import collections
 
 import WrightTools as wt
 
-from . import experiment
-from . import pulse
+from . import _pulse
+from ._experiment import Experiment
 from ._axis import Axis
 
 
@@ -27,7 +26,10 @@ def builtin(name):
 
 def from_ini(p):
     ini = wt.kit.INI(p)
-    # get axes
+    # pulse
+    pulse_class_name = ini.read('main', 'pulse class name')
+    pulse_class = _pulse.__dict__[pulse_class_name]
+    # axes
     axis_names = ini.sections
     axis_names.remove('main')
     axes = []
@@ -37,12 +39,11 @@ def from_ini(p):
         pulses = ini.read(name, 'pulses')
         parameter = ini.read(name, 'parameter')
         axis = Axis(points=points, units=None, name=name, active=active, pulses=pulses,
-                    parameter=parameter)
+                    parameter=parameter, cols=pulse_class.cols)
         axes.append(axis)
     # construct experiment object
     name = ini.read('main', 'name')
     pm = ini.read('main', 'pm')
-    pulse_class_name = ini.read('main', 'pulse class name')
-    e = experiment.Experiment(axes=axes, name=name, pm=pm, pulse_class_name=pulse_class_name)
+    e = Experiment(axes=axes, name=name, pm=pm, pulse_class=pulse_class)
     # finish
     return e

WrightSim/experiment/_axis.py

@@ -4,6 +4,9 @@
 # --- import --------------------------------------------------------------------------------------
 
 
+import numpy as np
+
+
 # --- define --------------------------------------------------------------------------------------
 
 
@@ -15,10 +18,28 @@
 
 class Axis(object):
 
-    def __init__(self, name, units=None, points=None, active=False, parameter=None, pulses=None):
+    def __init__(self, name, units=None, points=None, active=False, parameter=None, pulses=None,
+                 cols=None):
         self.name = name
         self.units = units
         self.points = points
         self.active = active
         self.parameter = parameter
         self.pulses = pulses
+        self.cols = cols
+        self._coords = None
+
+    def __repr__(self):
+        return "<WrightSim.Axis '{0}' at {1}>".format(self.name, id(self))
+
+    @property
+    def coords(self):
+        if self._coords is not None:
+            return self._coords
+        parameter_index = self.cols.index(self.parameter)
+        coords = np.zeros((len(self.pulses), 2), dtype=np.int)
+        for i, pulse in enumerate(self.pulses):
+            # specify pulse and then efield param of pulse
+            coords[i] = [pulse, parameter_index]
+        self._coords = coords
+        return coords

WrightSim/experiment/experiment.py WrightSim/experiment/_experiment.pyWrightSim/experiment/experiment.py renamed to WrightSim/experiment/_experiment.py

@@ -7,7 +7,8 @@
 
 import WrightTools as wt
 
-from . import pulse
+from . import _pulse
+from ._scan import Scan
 
 
 # --- define --------------------------------------------------------------------------------------
@@ -28,7 +29,7 @@
 class Experiment:
     """Experiment."""
 
-    def __init__(self, axes, name, pm, pulse_class_name):
+    def __init__(self, axes, name, pm, pulse_class):
         # basic attributes
         self.axes = axes
         for a in self.axes:
@@ -40,11 +41,11 @@ def __init__(self, axes, name, pm, pulse_class_name):
         self.early_buffer = early_buffer
         self.late_buffer = late_buffer
         # pulse
-        self.pulse_class_name = pulse_class_name
-        self.pulses = [pulse.__dict__[self.pulse_class_name]() for _ in self.pm]
+        self.pulse_class = pulse_class
+        self.pulses = [self.pulse_class() for _ in self.pm]
 
     def __repr__(self):
-        return 'WrightSim.Experiment object \'{0}\' at {1}'.format(self.name, str(id(self)))
+        return '<WrightSim.Experiment object \'{0}\' at {1}>'.format(self.name, str(id(self)))
 
     @property
     def active_axes(self):
@@ -54,6 +55,25 @@ def active_axes(self):
     def axis_names(self):
         return [a.name for a in self.axes]
 
+    def run(self, hamiltonian, mp=True):
+        """Run the experiment.
+
+        Parameters
+        ----------
+        hamiltonian : WrightSim Hamiltonian
+            Hamiltonian.
+        mp : boolean (optional)
+            Toggle CPU multiprocessing. Default is True.
+
+        Returns
+        -------
+        WrightSim Scan
+            Scan that was run."""
+        out = Scan(self, hamiltonian)
+        out.run(mp=mp)
+        # finish
+        return out
+
     def set_axis(self, axis_name, points):
         '''
         Activate and define points for one of the experimental axes.

WrightSim/experiment/_pulse.py

@@ -0,0 +1,122 @@
+# --- import --------------------------------------------------------------------------------------
+
+
+import numpy as np
+
+
+# --- define --------------------------------------------------------------------------------------
+
+
+wn_to_omega = 2 * np.pi * 3e-5  # omega is radians / fs
+
+
+# factor used to convert FWHM to stdev for function definition
+# gaussian defined such that FWHM,intensity scale = sigma * 2 * sqrt(ln(2))
+# (i.e. FWHM, intensity scale ~ 1.67 * sigma, amplitude scale)
+FWHM_to_sigma = 1. / (2 * np.sqrt(np.log(2)))
+
+
+# TODO: these should probably not be hard-coded
+timestep = 4.0
+early_buffer = 100.0
+late_buffer  = 400.0
+
+# --- functions -----------------------------------------------------------------------------------
+
+
+def _get_t(obj, d):
+    """
+        returns the t array that is appropriate for the given pulse
+        parameters; needs the appropriate buffers/timestep as well
+    """
+    if obj.fixed_bounds:
+        t_min = obj.fixed_bounds_min
+        t_max = obj.fixed_bounds_max
+    else:
+        t_min = d.min() - obj.early_buffer
+        t_max = d.max() + obj.late_buffer
+    # span up to and including t_max now
+    t = np.arange(t_min, t_max + obj.timestep, obj.timestep)
+    # alternate form:
+    return t
+
+
+# --- class ---------------------------------------------------------------------------------------
+
+
+class GaussRWA:
+    """
+    returns an array of gaussian values with arguments of
+        input array vec,
+        mean mu,
+        standard deviation sigma
+        amplitude amp
+        frequency freq
+        phase offset p
+
+    we are in rw approx, so if interaction is supposed to be negative, give
+    frequency a negative sign
+    """
+    # dictionary to associate array position to pulse attribute
+    cols = ['A',  # amplitude, a.u.
+            's',  # pulse FWHM (fs),
+            'd',  # pulse center delay (fs),
+            'w',  # frequency (wn),
+            'p']  # phase shift (radians)
+    # initial vars come from misc module, just as with scan module
+    timestep = timestep
+    early_buffer = early_buffer
+    late_buffer = late_buffer
+    # fixed bounds set to true if you want fixed time indices for the pulses
+    fixed_bounds = False
+    fixed_bounds_min = None
+    fixed_bounds_max = None
+
+    @classmethod
+    def pulse(cls, eparams, pm=None):
+        """Get efields for each pulse.
+
+        Parameters
+        ----------
+        eparams : 2D numpy ndarray
+            Array in (pulse, parameter). Refer to cols attribute for list
+            of parameters.
+        pm : iterable of integers
+            Phase matching. If None, fully positive phase matching is
+            assumed.
+
+        Returns
+        -------
+        (1D array, 1D array)
+            Tuple of 1D numpy ndarray time (fs) and complex efield (a.u.).
+        """
+        # import if the size is right
+        area  = eparams[:, 0].copy().astype(float)
+        sigma = eparams[:, 1].copy().astype(float)
+        mu    = eparams[:, 2].copy().astype(float)
+        freq  = eparams[:, 3].copy().astype(float)
+        p     = eparams[:, 4].copy().astype(float)
+        # proper unit conversions
+        sigma *= FWHM_to_sigma
+        freq *= wn_to_omega
+        # redefine delays s.t. always relative to the first pulse listed
+        offset = mu[0]
+        # subtract off the value
+        mu -= offset
+        # normalize amplitdue to stdev
+        y = area / (sigma * np.sqrt(2 * np.pi))
+        #print y, sigma, mu, freq, p
+        # calculate t
+        t = cls.get_t(mu)
+        # incorporate complex conjugates if necessary
+        if pm is 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*= y[:,None] * np.exp(-(t[None,:] - mu[:,None])**2 / (2*sigma[:,None]**2) )
+        return t, x
+
+    @classmethod
+    def get_t(cls, d):
+        return _get_t(cls, d)