|
| 1 | +""" |
| 2 | +@author Kevin Li, Michael Schenk |
| 3 | +@date 31. January 2014 |
| 4 | +@brief Collection of cython functions to calculate statistics |
| 5 | + of bunch and slice_set data. |
| 6 | +@copyright CERN |
| 7 | +""" |
| 8 | +import numpy as np |
| 9 | +cimport numpy as np |
| 10 | +cimport libc.math as cmath |
| 11 | + |
| 12 | +cimport cython.boundscheck |
| 13 | +cimport cython.cdivision |
| 14 | + |
| 15 | + |
| 16 | +@cython.boundscheck(False) |
| 17 | +cpdef double mean(double[::1] u): |
| 18 | + """ Cython function to calculate the mean value of dataset u. """ |
| 19 | + cdef double mean_u = 0 |
| 20 | + cdef unsigned int n = u.shape[0] |
| 21 | + |
| 22 | + cdef unsigned int i |
| 23 | + for i in xrange(n): |
| 24 | + mean_u += u[i] |
| 25 | + if n: |
| 26 | + mean_u /= n |
| 27 | + |
| 28 | + return mean_u |
| 29 | + |
| 30 | +@cython.boundscheck(False) |
| 31 | +cpdef double std(double[::1] u): |
| 32 | + """ Cython function to calculate the standard deviation of |
| 33 | + dataset u. """ |
| 34 | + cdef double mean_u = mean(u) |
| 35 | + cdef double std_u = 0 |
| 36 | + cdef double du = 0 |
| 37 | + |
| 38 | + cdef unsigned int n = u.shape[0] |
| 39 | + cdef unsigned int i |
| 40 | + for i in xrange(n): |
| 41 | + du = u[i] - mean_u |
| 42 | + std_u += du * du |
| 43 | + if n: |
| 44 | + std_u /= n |
| 45 | + |
| 46 | + return cmath.sqrt(std_u) |
| 47 | + |
| 48 | +@cython.boundscheck(False) |
| 49 | +cpdef double emittance(double[::1] u, double[::1] up): |
| 50 | + """ Cython function to calculate the emittance of datasets |
| 51 | + u and up, i.e. a coordinate-momentum pair. To calculate the |
| 52 | + emittance, one needs the mean values of quantities u and |
| 53 | + up. """ |
| 54 | + cdef double mean_u = mean(u) |
| 55 | + cdef double mean_up = mean(up) |
| 56 | + |
| 57 | + cdef double u2 = 0 |
| 58 | + cdef double up2 = 0 |
| 59 | + cdef double uup = 0 |
| 60 | + cdef double du = 0 |
| 61 | + cdef double dup = 0 |
| 62 | + |
| 63 | + cdef unsigned int n = u.shape[0] |
| 64 | + cdef unsigned int i |
| 65 | + for i in xrange(n): |
| 66 | + du = u[i] - mean_u |
| 67 | + dup = up[i] - mean_up |
| 68 | + |
| 69 | + u2 += du * du |
| 70 | + up2 += dup * dup |
| 71 | + uup += du * dup |
| 72 | + if n: |
| 73 | + u2 /= n |
| 74 | + up2 /= n |
| 75 | + uup /= n |
| 76 | + |
| 77 | + return cmath.sqrt(u2*up2 - uup*uup) |
| 78 | + |
| 79 | + |
| 80 | +''' |
| 81 | +Cython statistics functions for an instance of a SliceSet class. |
| 82 | +''' |
| 83 | + |
| 84 | +@cython.boundscheck(False) |
| 85 | +@cython.cdivision(True) |
| 86 | +cpdef count_macroparticles_per_slice(int[::1] slice_index_of_particle, |
| 87 | + int[::1] particles_within_cuts, |
| 88 | + int[::1] n_macroparticles): |
| 89 | + """ Cython function to count the number of macroparticles in |
| 90 | + each slice. """ |
| 91 | + cdef unsigned int n_particles_within_cuts = particles_within_cuts.shape[0] |
| 92 | + cdef unsigned int s_idx, i |
| 93 | + |
| 94 | + for i in xrange(n_particles_within_cuts): |
| 95 | + s_idx = slice_index_of_particle[particles_within_cuts[i]] |
| 96 | + n_macroparticles[s_idx] += 1 |
| 97 | + |
| 98 | + |
| 99 | +@cython.boundscheck(False) |
| 100 | +@cython.cdivision(True) |
| 101 | +cpdef sort_particle_indices_by_slice(int[::1] slice_index_of_particle, |
| 102 | + int[::1] particles_within_cuts, |
| 103 | + int[::1] slice_positions, |
| 104 | + int[::1] particle_indices_by_slice): |
| 105 | + """ Iterate once through all the particles within the slicing |
| 106 | + region and assign their position in the bunch.z array to the |
| 107 | + respective slice they are in. |
| 108 | + This is to provide a method to the user that allows to see |
| 109 | + which particles are in a specific slice (see |
| 110 | + particle_indices_of_slice in SliceSet class). """ |
| 111 | + cdef unsigned int n_part_in_cuts = particles_within_cuts.shape[0] |
| 112 | + cdef unsigned int n_slices = slice_positions.shape[0] - 1 |
| 113 | + |
| 114 | + cdef unsigned int[::1] pos_ctr = np.zeros(n_slices, dtype=np.uint32) |
| 115 | + cdef unsigned int i, p_idx, s_idx |
| 116 | + cdef unsigned int pos |
| 117 | + |
| 118 | + for i in xrange(n_part_in_cuts): |
| 119 | + p_idx = particles_within_cuts[i] |
| 120 | + s_idx = slice_index_of_particle[p_idx] |
| 121 | + |
| 122 | + pos = slice_positions[s_idx] + pos_ctr[s_idx] |
| 123 | + particle_indices_by_slice[pos] = p_idx |
| 124 | + pos_ctr[s_idx] += 1 |
| 125 | + |
| 126 | + |
| 127 | +@cython.boundscheck(False) |
| 128 | +@cython.cdivision(True) |
| 129 | +cpdef mean_per_slice(int[::1] slice_index_of_particle, |
| 130 | + int[::1] particles_within_cuts, |
| 131 | + int[::1] n_macroparticles, |
| 132 | + double[::1] u, double[::1] mean_u): |
| 133 | + """ Iterate once through all the particles within the |
| 134 | + slicing region and calculate simultaneously the mean |
| 135 | + value of quantity u for each slice separately. """ |
| 136 | + cdef unsigned int n_part_in_cuts = particles_within_cuts.shape[0] |
| 137 | + cdef unsigned int n_slices = mean_u.shape[0] |
| 138 | + cdef unsigned int p_idx, s_idx, i |
| 139 | + |
| 140 | + for i in xrange(n_part_in_cuts): |
| 141 | + p_idx = particles_within_cuts[i] |
| 142 | + s_idx = slice_index_of_particle[p_idx] |
| 143 | + mean_u[s_idx] += u[p_idx] |
| 144 | + |
| 145 | + for i in xrange(n_slices): |
| 146 | + if n_macroparticles[i]: |
| 147 | + mean_u[i] /= n_macroparticles[i] |
| 148 | + |
| 149 | + |
| 150 | +@cython.boundscheck(False) |
| 151 | +@cython.cdivision(True) |
| 152 | +cpdef std_per_slice(int[::1] slice_index_of_particle, |
| 153 | + int[::1] particles_within_cuts, |
| 154 | + int[::1] n_macroparticles, |
| 155 | + double[::1] u, double[::1] std_u): |
| 156 | + """ Iterate once through all the particles within the |
| 157 | + slicing region and calculate simultaneously the |
| 158 | + standard deviation of quantity u for each slice |
| 159 | + separately. """ |
| 160 | + cdef unsigned int n_part_in_cuts = particles_within_cuts.shape[0] |
| 161 | + cdef unsigned int n_slices = std_u.shape[0] |
| 162 | + cdef unsigned int p_idx, s_idx, i |
| 163 | + cdef double du |
| 164 | + |
| 165 | + cdef double[::1] mean_u = np.zeros(n_slices, dtype=np.double) |
| 166 | + mean_per_slice(slice_index_of_particle, particles_within_cuts, |
| 167 | + n_macroparticles, u, mean_u) |
| 168 | + |
| 169 | + for i in xrange(n_part_in_cuts): |
| 170 | + p_idx = particles_within_cuts[i] |
| 171 | + s_idx = slice_index_of_particle[p_idx] |
| 172 | + |
| 173 | + du = u[p_idx] - mean_u[s_idx] |
| 174 | + std_u[s_idx] += du * du |
| 175 | + |
| 176 | + for i in xrange(n_slices): |
| 177 | + if n_macroparticles[i]: |
| 178 | + std_u[i] /= n_macroparticles[i] |
| 179 | + |
| 180 | + std_u[i] = cmath.sqrt(std_u[i]) |
| 181 | + |
| 182 | + |
| 183 | +@cython.boundscheck(False) |
| 184 | +@cython.cdivision(True) |
| 185 | +cpdef emittance_per_slice(int[::1] slice_index_of_particle, |
| 186 | + int[::1] particles_within_cuts, |
| 187 | + int[::1] n_macroparticles, |
| 188 | + double[::1] u, double[::1] up, double[::1] epsn_u): |
| 189 | + """ Iterate once through all the particles within the |
| 190 | + slicing region and calculate simultaneously the emittance |
| 191 | + of quantities u and up, i.e. a coordinate-momentum pair, |
| 192 | + for each slice separately. To calculate the emittance per |
| 193 | + slice, one needs the mean values of quantities u and up |
| 194 | + for each slice. """ |
| 195 | + cdef unsigned int n_part_in_cuts = particles_within_cuts.shape[0] |
| 196 | + cdef unsigned int n_slices = epsn_u.shape[0] |
| 197 | + cdef unsigned int p_idx, s_idx, i |
| 198 | + |
| 199 | + # Determine mean values of u and up for each slice. |
| 200 | + cdef double[::1] mean_u = np.zeros(n_slices, dtype=np.double) |
| 201 | + cdef double[::1] mean_up = np.zeros(n_slices, dtype=np.double) |
| 202 | + mean_per_slice(slice_index_of_particle, particles_within_cuts, |
| 203 | + n_macroparticles, u, mean_u) |
| 204 | + mean_per_slice(slice_index_of_particle, particles_within_cuts, |
| 205 | + n_macroparticles, up, mean_up) |
| 206 | + |
| 207 | + cdef double du, dup |
| 208 | + cdef double[::1] u2 = np.zeros(n_slices, dtype=np.double) |
| 209 | + cdef double[::1] up2 = np.zeros(n_slices, dtype=np.double) |
| 210 | + cdef double[::1] uup = np.zeros(n_slices, dtype=np.double) |
| 211 | + |
| 212 | + for i in xrange(n_part_in_cuts): |
| 213 | + p_idx = particles_within_cuts[i] |
| 214 | + s_idx = slice_index_of_particle[p_idx] |
| 215 | + |
| 216 | + du = u[p_idx] - mean_u[s_idx] |
| 217 | + dup = up[p_idx] - mean_up[s_idx] |
| 218 | + |
| 219 | + u2[s_idx] += du * du |
| 220 | + up2[s_idx] += dup * dup |
| 221 | + uup[s_idx] += du * dup |
| 222 | + |
| 223 | + for i in xrange(n_slices): |
| 224 | + if n_macroparticles[i]: |
| 225 | + u2[i] /= n_macroparticles[i] |
| 226 | + up2[i] /= n_macroparticles[i] |
| 227 | + uup[i] /= n_macroparticles[i] |
| 228 | + |
| 229 | + epsn_u[i] = cmath.sqrt(u2[i]*up2[i] - uup[i]*uup[i]) |
0 commit comments