release-script: Merge branch 'release/v1.13.0'

LICENSE.txt

@@ -0,0 +1,33 @@
+BSD 4-Clause License
+
+Copyright (c) 2018, CERN
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+
+* Redistributions of source code must retain the above copyright notice, this
+  list of conditions and the following disclaimer.
+
+* Redistributions in binary form must reproduce the above copyright notice,
+  this list of conditions and the following disclaimer in the documentation
+  and/or other materials provided with the distribution.
+
+* Neither the name of the copyright holder nor the names of its
+  contributors may be used to endorse or promote products derived from
+  this software without specific prior written permission.
+
+* Neither the name of the copyright holder nor the names of its contributors
+  may be used to endorse or promote products derived from this software
+  without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
+FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

PyHEADTAIL/_version.py

@@ -1 +1 @@
-__version__ = '1.12.4'
+__version__ = '1.13.0'

PyHEADTAIL/cobra_functions/stats.pyx

@@ -418,8 +418,11 @@ cpdef emittance_per_slice(int[::1] slice_index_of_particle,
 @cython.boundscheck(False)
 @cython.cdivision(True)
 @cython.wraparound(False)
-cpdef calc_cell_stats(bunch, double beta_z, double radial_cut,
-                      int n_rings, int n_azim_slices):
+cpdef calc_cell_stats(
+        double[::1] x, double[::1] xp, double[::1] y,
+        double[::1] yp, double[::1] z, double[::1] dp,
+        double beta_z, double radial_cut,
+        int n_rings, int n_azim_slices):
 
     # Prepare arrays to store cell statistics.
     cdef int[:,::1] n_particles_cell = np.zeros((n_azim_slices, n_rings),
@@ -438,12 +441,12 @@ cpdef calc_cell_stats(bunch, double beta_z, double radial_cut,
                                                dtype=np.double)
 
     # Declare datatypes of bunch coords.
-    cdef double[::1] x = bunch.x
-    cdef double[::1] xp = bunch.xp
-    cdef double[::1] y = bunch.y
-    cdef double[::1] yp = bunch.yp
-    cdef double[::1] z = bunch.z
-    cdef double[::1] dp = bunch.dp
+    # cdef double[::1] x = bunch.x
+    # cdef double[::1] xp = bunch.xp
+    # cdef double[::1] y = bunch.y
+    # cdef double[::1] yp = bunch.yp
+    # cdef double[::1] z = bunch.z
+    # cdef double[::1] dp = bunch.dp
     cdef unsigned int n_particles = x.shape[0]
 
     cdef double ring_width = radial_cut / <double>n_rings

PyHEADTAIL/feedback/transverse_damper.py

@@ -14,42 +14,81 @@
 
 class TransverseDamper(Element):
 
-    def __init__(self, dampingrate_x, dampingrate_y, *args, **kwargs):
+    def __init__(self, dampingrate_x, dampingrate_y, phase=90,
+                 local_beta_function=None, *args, **kwargs):
+        '''Ideal transverse damper with an in-place "measurement"
+        (transverse "pick-up") of the transverse dipole moment.
+        Note: a single bunch in the beam is assumed, i.e. this works on
+        the entire beam's moments.
+
+        Arguments:
+            - dampingrate_x, dampingrate_y: horizontal and vertical
+                damping rates in turns (e.g. 50 turns for a typical 2018
+                LHC ADT set-up)
+            - phase: phase of the damper kick in degrees with respect to
+                the transverse position "pick-up". The default value of
+                90 degrees corresponds to a typical resistive damper.
+            - local_beta_function: the optics beta function at the
+                transverse position "pick-up" (e.g. in the local place
+                of this Element). This is required if the damper is not
+                a purely resistive damper (or exciter), i.e. if the
+                phase is not 90 (or 270) degrees. The beta function is
+                assumed to be the same for both transverse planes,
+                otherwise use two instances of the TransverseDamper.
+        '''
 
         if dampingrate_x and not dampingrate_y:
             self.gain_x = 2/dampingrate_x
             self.track = self.track_horizontal
-            self.prints('Damper in V active')
+            self.prints('Damper in horizontal plane active')
         elif not dampingrate_x and dampingrate_y:
             self.gain_y = 2/dampingrate_y
             self.track = self.track_vertical
-            self.prints('Damper in Y active')
+            self.prints('Damper in vertical plane active')
         elif not dampingrate_x and not dampingrate_y:
             self.prints('Dampers not active')
         else:
             self.gain_x = 2/dampingrate_x
             self.gain_y = 2/dampingrate_y
             self.track = self.track_all
             self.prints('Dampers active')
+        if phase != 90 and phase != 270 and not local_beta_function:
+            raise TypeError(
+                'TransverseDamper: numeric local_beta_function value at '
+                'position of damper missing! (Required because of non-zero '
+                'reactive damper component.)')
+        self.phase_in_2pi = phase / 360. * 2*np.pi
+        self.local_beta_function = local_beta_function
 
     # will be overwritten at initialisation
     def track(self, beam):
         pass
 
     def track_horizontal(self, beam):
-        beam.xp -= self.gain_x * beam.mean_xp()
+        beam.xp -= self.gain_x * np.sin(self.phase_in_2pi) * beam.mean_xp()
+        if self.local_beta_function:
+            beam.xp -= (self.gain_x * np.cos(self.phase_in_2pi) *
+                        beam.mean_x() / self.local_beta_function)
 
     def track_vertical(self, beam):
-        beam.yp -= self.gain_y * beam.mean_yp()
+        beam.yp -= self.gain_y * np.sin(self.phase_in_2pi) * beam.mean_yp()
+        if self.local_beta_function:
+            beam.yp -= (self.gain_y * np.cos(self.phase_in_2pi) *
+                        beam.mean_y() / self.local_beta_function)
 
     def track_all(self, beam):
-        beam.xp -= self.gain_x * beam.mean_xp()
-        beam.yp -= self.gain_y * beam.mean_yp()
+        beam.xp -= self.gain_x * np.sin(self.phase_in_2pi) * beam.mean_xp()
+        beam.yp -= self.gain_y * np.sin(self.phase_in_2pi) * beam.mean_yp()
+        if self.local_beta_function:
+            beam.xp -= (self.gain_x * np.cos(self.phase_in_2pi) *
+                        beam.mean_x() / self.local_beta_function)
+            beam.yp -= (self.gain_y * np.cos(self.phase_in_2pi) *
+                        beam.mean_y() / self.local_beta_function)
 
     @classmethod
-    def horizontal(cls, dampingrate_x):
-        return cls(dampingrate_x, 0)
+    def horizontal(cls, dampingrate_x, *args, **kwargs):
+        return cls(dampingrate_x, 0, *args, **kwargs)
 
     @classmethod
-    def vertical(cls, dampingrate_y):
-        return cls(0, dampingrate_y)
+    def vertical(cls, dampingrate_y, *args, **kwargs):
+        return cls(0, dampingrate_y, *args, **kwargs)

PyHEADTAIL/gpu/gpu_wrap.py

@@ -47,6 +47,11 @@
     #        'No GPU capabilities available')
     has_pycuda = False
 
+def _empty_like(gpuarray):
+    return pycuda.gpuarray.empty(
+        shape=gpuarray.shape, dtype=gpuarray.dtype,
+        allocator=gpu_utils.memory_pool.allocate)
+
 
 if has_pycuda:
     # define all compilation depending functions (e.g. ElementwiseKernel)
@@ -57,7 +62,7 @@
     )
     def sub_scalar(gpuarr, scalar, out=None, stream=None):
         if out is None:
-            out = pycuda.gpuarray.empty_like(gpuarr)
+            out = _empty_like(gpuarr)
         _sub_1dgpuarr(out, gpuarr, scalar, stream=stream)
         return out
 
@@ -71,14 +76,14 @@ def _mul_scalar(gpuarr, scalar, out=None, stream=None):
         to specify a stream
         '''
         if out is None:
-            out = pycuda.gpuarray.empty_like(gpuarr)
+            out = _empty_like(gpuarr)
         _mul_with_factor(out, gpuarr, scalar, stream=stream)
 
     def _multiply(a, b, out=None, stream=None):
         '''Elementwise multiply of two gpuarray specifying a stream
         Required because gpuarray.__mul__ has no stream argument'''
         if out is None:
-            out = pycuda.gpuarray.empty_like(a)
+            out = _empty_like(a)
         func = pycuda.elementwise.get_binary_op_kernel(a.dtype, b.dtype,
             out.dtype, "*")
         func.prepared_async_call(a._grid, a._block, stream, a.gpudata,
@@ -136,7 +141,7 @@ def _compute_sigma(a, b, c, d, out=None, stream=None):
         '''Computes elementwise a - b*c/d as required in compute sigma for
         the emittance '''
         if out is None:
-            out = pycuda.gpuarray.empty_like(a)
+            out = _empty_like(a)
         _comp_sigma(out, a, b, c, d, stream=stream)
         return out
 
@@ -160,7 +165,7 @@ def _emittance_dispersion(
             n, cov_u2, cov_u_up, cov_up2, cov_u_dp, cov_up_dp,
             cov_dp2, out=None, stream=None):
         if out is None:
-            out = pycuda.gpuarray.empty_like(cov_u2)
+            out = _empty_like(cov_u2)
         _emitt_disp(out, cov_u2, cov_u_up, cov_up2, cov_u_dp, cov_up_dp,
                     cov_dp2, np.float64(n), stream=stream)
         return out
@@ -176,7 +181,7 @@ def _emittance_dispersion(
     def _emittance_no_dispersion(
             n, cov_u2, cov_u_up, cov_up2, out=None, stream=None):
         if out is None:
-            out = pycuda.gpuarray.empty_like(cov_u2)
+            out = _empty_like(cov_u2)
         _emitt_nodisp(out, cov_u2, cov_u_up, cov_up2, np.float64(n),
                       stream=stream)
         return out
@@ -194,9 +199,9 @@ def wofz(in_real, in_imag, out_real=None, out_imag=None, stream=None):
         part of z.
         '''
         if out_real is None:
-            out_real = pycuda.gpuarray.empty_like(in_real)
+            out_real = _empty_like(in_real)
         if out_imag is None:
-            out_imag = pycuda.gpuarray.empty_like(in_imag)
+            out_imag = _empty_like(in_imag)
         _wofz(in_real, in_imag, out_real, out_imag, stream=stream)
         return out_real, out_imag
 
@@ -207,7 +212,7 @@ def wofz(in_real, in_imag, out_real=None, out_imag=None, stream=None):
     )
     def sign(array, out=None, stream=None):
         if out is None:
-            out = pycuda.gpuarray.empty_like(array)
+            out = _empty_like(array)
         _sign(array, out, stream=stream)
         return out
 
@@ -260,13 +265,9 @@ def thrust_mean_and_std_per_slice(sliceset, u, stream=None):
         p_sids = sliceset.slice_index_of_particle
         # slice_index_of_particle may have slice indices outside of slicing area,
         # the following arrays therefore can comprise non valid slice entries
-        slice_ids_noncontained = pycuda.gpuarray.empty(
-            p_sids.shape, dtype=p_sids.dtype,
-            allocator=gpu_utils.memory_pool.allocate)
-        slice_means_noncontained = pycuda.gpuarray.empty(
-            u.shape, dtype=u.dtype, allocator=gpu_utils.memory_pool.allocate)
-        slice_stds_noncontained = pycuda.gpuarray.empty(
-            u.shape, dtype=u.dtype, allocator=gpu_utils.memory_pool.allocate)
+        slice_ids_noncontained = _empty_like(p_sids)
+        slice_means_noncontained = _empty_like(u)
+        slice_stds_noncontained = _empty_like(u)
 
         (_, _, _, new_end) = thrust.thrust_stats_per_slice(
             p_sids, u, slice_ids_noncontained, slice_means_noncontained,
@@ -368,8 +369,8 @@ def covariance(a,b, stream=None):
         b: pycuda.GPUArray
     '''
     n = len(a)
-    x = pycuda.gpuarray.empty_like(a)
-    y = pycuda.gpuarray.empty_like(b)
+    x = _empty_like(a)
+    y = _empty_like(b)
     mean_a = skcuda.misc.mean(a)
     #x -= mean_a
     _sub_1dgpuarr(x, a, mean_a, stream=stream)
@@ -397,7 +398,7 @@ def std(a, stream=None):
     #return skcuda.misc.std(a, ddof=1)
     n = len(a)
     #mean_a = skcuda.misc.mean(a)
-    x = pycuda.gpuarray.empty_like(a)
+    x = _empty_like(a)
     mean_a = mean(a, stream=stream)
     _sub_1dgpuarr(x, a, mean_a, stream=stream)
     _inplace_pow(x, 2, stream=stream)
@@ -523,7 +524,7 @@ def emittance(u, up, dp, stream=None):
     n = len(u)
     mean_u = mean(u, stream=stream)
     mean_up = mean(up, stream=stream)
-    out = pycuda.gpuarray.empty_like(mean_u)
+    out = _empty_like(mean_u)
     tmp_u = sub_scalar(u, mean_u, stream=stream)
     tmp_up = sub_scalar(up, mean_up, stream=stream)
     tmp_space = _multiply(tmp_u, tmp_u, stream=stream)
@@ -568,7 +569,7 @@ def emittance_multistream(u, up, dp, stream=None):
     tmp_u = sub_scalar(u, mean_u, stream=streams[0])
     tmp_space = _multiply(tmp_u, tmp_u, stream=streams[0])
     cov_u2 = pycuda.gpuarray.sum(tmp_space, stream=streams[0])
-    out = pycuda.gpuarray.empty_like(mean_u)
+    out = _empty_like(mean_u)
     tmp_up = sub_scalar(up, mean_up, stream=streams[1])
     streams[0].synchronize()
     streams[1].synchronize()
@@ -601,11 +602,11 @@ def cumsum(array, dest=None):
     '''
     if array.dtype == np.int32:
         if dest is None:
-            dest = pycuda.gpuarray.empty_like(array)
+            dest = _empty_like(array)
         thrust_interface.thrust_cumsum_int(array, dest)
     elif array.dtype == np.float64:
         if dest is None:
-            dest = pycuda.gpuarray.empty_like(array)
+            dest = _empty_like(array)
         thrust_interface.thrust_cumsum_double(array, dest)
     else:
         dest = array.copy()
@@ -660,7 +661,7 @@ def apply_permutation(array, permutation):
         permutation permutation array: must be np.int32 (or int32), is asserted
     '''
     assert(permutation.dtype.itemsize == 4 and permutation.dtype.kind is 'i')
-    tmp = pycuda.gpuarray.empty_like(array)
+    tmp = _empty_like(array)
     dtype = array.dtype
     if dtype.itemsize == 8 and dtype.kind is 'f':
         thrust.apply_sort_perm_double(array, tmp, permutation)
@@ -832,7 +833,7 @@ def sorted_emittance_per_slice(sliceset, u, up, dp=None, stream=None):
     cov_u2 = sorted_cov_per_slice(sliceset, u, u, stream=streams[0])
     cov_up2= sorted_cov_per_slice(sliceset, up, up, stream=streams[1])
     cov_u_up = sorted_cov_per_slice(sliceset, u, up, stream=streams[2])
-    out = pycuda.gpuarray.empty_like(cov_u2)
+    out = _empty_like(cov_u2)
     # use this factor in emitt_disp: the code has a 1/(n*n+n) factor which is not
     # required here since the scaling is done in the cov_per_slice
     # --> 1/(n*n + n) must be 1. ==> n = sqrt(5)/2 -0.5

PyHEADTAIL/impedances/wake_kicks.py

@@ -105,7 +105,7 @@ def _convolution_numpy(self, target_times, source_times,
         try:
             source_times = source_times.get()
         except AttributeError:
-            pass #is already on GPU
+            pass #is already on CPU
         dt_to_target_slice = np.concatenate(
             (target_times - source_times[-1],
             (target_times - source_times[0])[1:]))