nifty_utilities.py 4.92 KB
 Ultima committed Aug 25, 2015 1 2 3 4 ``````# -*- coding: utf-8 -*- import numpy as np `````` Ultima committed Nov 12, 2015 5 6 7 `````` def hermitianize_gaussian(x): # make the point inversions `````` Ultima committed Aug 25, 2015 8 9 `````` flipped_x = _hermitianize_inverter(x) flipped_x = flipped_x.conjugate() `````` Ultima committed Nov 12, 2015 10 `````` # check if x was already hermitian `````` Ultima committed Aug 25, 2015 11 12 `````` if (x == flipped_x).all(): return x `````` Ultima committed Nov 12, 2015 13 14 `````` # average x and flipped_x. # Correct the variance by multiplying sqrt(0.5) `````` Ultima committed Aug 25, 2015 15 `````` x = (x + flipped_x) * np.sqrt(0.5) `````` Ultima committed Nov 12, 2015 16 17 18 `````` # The fixed points of the point inversion must not be avaraged. # Hence one must multiply them again with sqrt(0.5) # -> Get the middle index of the array `````` Ultima committed Aug 25, 2015 19 20 `````` mid_index = np.array(x.shape, dtype=np.int)//2 dimensions = mid_index.size `````` Ultima committed Nov 12, 2015 21 22 `````` # Use ndindex to iterate over all combinations of zeros and the # mid_index in order to correct all fixed points. `````` Ultima committed Aug 25, 2015 23 24 25 26 27 28 29 `````` for i in np.ndindex((2,)*dimensions): temp_index = tuple(i*mid_index) x[temp_index] *= np.sqrt(0.5) try: x.hermitian = True except(AttributeError): pass `````` Ultima committed Sep 10, 2015 30 `````` `````` Ultima committed Aug 25, 2015 31 32 33 `````` return x `````` Ultima committed Nov 12, 2015 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 ``````def hermitianize(x): # make the point inversions flipped_x = _hermitianize_inverter(x) flipped_x = flipped_x.conjugate() # check if x was already hermitian if (x == flipped_x).all(): return x # average x and flipped_x. # Correct the variance by multiplying sqrt(0.5) x = (x + flipped_x) / 2. try: x.hermitian = True except(AttributeError): pass return x `````` Ultima committed Aug 25, 2015 51 52 `````` def _hermitianize_inverter(x): `````` Ultima committed Nov 12, 2015 53 `````` # calculate the number of dimensions the input array has `````` Ultima committed Aug 25, 2015 54 `````` dimensions = len(x.shape) `````` Ultima committed Nov 12, 2015 55 56 57 `````` # prepare the slicing object which will be used for mirroring slice_primitive = [slice(None), ]*dimensions # copy the input data `````` Ultima committed Aug 25, 2015 58 `````` y = x.copy() `````` Ultima committed Nov 12, 2015 59 `````` # flip in every direction `````` Ultima committed Aug 25, 2015 60 61 `````` for i in xrange(dimensions): slice_picker = slice_primitive[:] `````` Ultima committed Nov 12, 2015 62 `````` slice_picker[i] = slice(1, None, None) `````` ultimanet committed Mar 03, 2016 63 64 `````` slice_picker = tuple(slice_picker) `````` Ultima committed Aug 25, 2015 65 66 `````` slice_inverter = slice_primitive[:] slice_inverter[i] = slice(None, 0, -1) `````` ultimanet committed Mar 03, 2016 67 `````` slice_inverter = tuple(slice_inverter) `````` Ultima committed Aug 25, 2015 68 69 `````` try: `````` ultimanet committed Jan 25, 2016 70 71 `````` y.set_data(to_key=slice_picker, data=y, from_key=slice_inverter) `````` Ultima committed Aug 25, 2015 72 73 74 `````` except(AttributeError): y[slice_picker] = y[slice_inverter] return y `````` Ultima committed Sep 10, 2015 75 76 `````` `````` Ultima committed Dec 11, 2015 77 ``````def direct_vdot(x, y): `````` Ultima committed Nov 12, 2015 78 `````` # the input could be fields. Try to extract the data `````` Ultima committed Aug 25, 2015 79 80 81 82 83 84 85 86 `````` try: x = x.get_val() except(AttributeError): pass try: y = y.get_val() except(AttributeError): pass `````` Ultima committed Nov 12, 2015 87 `````` # try to make a direct vdot `````` Ultima committed Aug 25, 2015 88 89 90 91 `````` try: return x.vdot(y) except(AttributeError): pass `````` Ultima committed Sep 10, 2015 92 `````` `````` Ultima committed Aug 25, 2015 93 94 95 `````` try: return y.vdot(x) except(AttributeError): `````` Ultima committed Sep 10, 2015 96 97 `````` pass `````` Ultima committed Nov 12, 2015 98 `````` # fallback to numpy `````` Ultima committed Sep 10, 2015 99 100 `````` return np.vdot(x, y) `````` Ultima committed Aug 25, 2015 101 102 `````` def convert_nested_list_to_object_array(x): `````` Ultima committed Nov 12, 2015 103 104 105 `````` # if x is a nested_list full of ndarrays all having the same size, # np.shape returns the shape of the ndarrays, too, i.e. too many # dimensions `````` Ultima committed Aug 25, 2015 106 `````` possible_shape = np.shape(x) `````` Ultima committed Nov 12, 2015 107 `````` # Check if possible_shape goes too deep. `````` Ultima committed Sep 10, 2015 108 `````` dimension_counter = 0 `````` Ultima committed Aug 25, 2015 109 110 111 112 113 114 115 116 `````` current_extract = x for i in xrange(len(possible_shape)): if isinstance(current_extract, list) == False and\ isinstance(current_extract, tuple) == False: break current_extract = current_extract[0] dimension_counter += 1 real_shape = possible_shape[:dimension_counter] `````` Ultima committed Nov 12, 2015 117 `````` # if the numpy array was not encapsulated at all, return x directly `````` Ultima committed Aug 25, 2015 118 119 `````` if real_shape == (): return x `````` Ultima committed Nov 12, 2015 120 121 `````` # Prepare the carrier-object carrier = np.empty(real_shape, dtype=np.object) `````` Ultima committed Aug 25, 2015 122 123 `````` for i in xrange(np.prod(real_shape)): ii = np.unravel_index(i, real_shape) `````` Ultima committed Sep 10, 2015 124 `````` try: `````` Ultima committed Aug 25, 2015 125 126 127 `````` carrier[ii] = x[ii] except(TypeError): extracted = x `````` Ultima committed Sep 10, 2015 128 `````` for j in xrange(len(ii)): `````` Ultima committed Aug 25, 2015 129 `````` extracted = extracted[ii[j]] `````` Ultima committed Sep 10, 2015 130 `````` carrier[ii] = extracted `````` Ultima committed Aug 25, 2015 131 132 133 134 135 136 137 138 139 140 141 142 143 144 `````` return carrier def field_map(ishape, function, *args): if ishape == (): return function(*args) else: if args == (): result = np.empty(ishape, dtype=np.object) for i in xrange(np.prod(ishape)): ii = np.unravel_index(i, ishape) result[ii] = function() return result else: `````` Ultima committed Nov 12, 2015 145 146 147 148 `````` # define a helper function in order to clip the get-indices # to be suitable for the foreign arrays in args. # This allows you to do operations, like adding to fields # with ishape (3,4,3) and (3,4,1) `````` Ultima committed Aug 25, 2015 149 150 151 152 153 154 155 156 157 158 `````` def get_clipped(w, ind): w_shape = np.array(np.shape(w)) get_tuple = tuple(np.clip(ind, 0, w_shape-1)) return w[get_tuple] result = np.empty_like(args[0]) for i in xrange(np.prod(result.shape)): ii = np.unravel_index(i, result.shape) result[ii] = function(*map( lambda z: get_clipped(z, ii), args) ) `````` Ultima committed Nov 12, 2015 159 160 `````` # result[ii] = function(*map(lambda z: z[ii], args)) return result``````