nifty_utilities.py 4.92 KB
Newer Older
Ultima's avatar
Ultima committed
1
2
3
4
# -*- coding: utf-8 -*-

import numpy as np

Ultima's avatar
Ultima committed
5
6
7

def hermitianize_gaussian(x):
    # make the point inversions
Ultima's avatar
Ultima committed
8
9
    flipped_x = _hermitianize_inverter(x)
    flipped_x = flipped_x.conjugate()
Ultima's avatar
Ultima committed
10
    # check if x was already hermitian
Ultima's avatar
Ultima committed
11
12
    if (x == flipped_x).all():
        return x
Ultima's avatar
Ultima committed
13
14
    # average x and flipped_x.
    # Correct the variance by multiplying sqrt(0.5)
Ultima's avatar
Ultima committed
15
    x = (x + flipped_x) * np.sqrt(0.5)
Ultima's avatar
Ultima committed
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's avatar
Ultima committed
19
20
    mid_index = np.array(x.shape, dtype=np.int)//2
    dimensions = mid_index.size
Ultima's avatar
Ultima committed
21
22
    # Use ndindex to iterate over all combinations of zeros and the
    # mid_index in order to correct all fixed points.
Ultima's avatar
Ultima committed
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
30

Ultima's avatar
Ultima committed
31
32
33
    return x


Ultima's avatar
Ultima committed
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's avatar
Ultima committed
51
52

def _hermitianize_inverter(x):
Ultima's avatar
Ultima committed
53
    # calculate the number of dimensions the input array has
Ultima's avatar
Ultima committed
54
    dimensions = len(x.shape)
Ultima's avatar
Ultima committed
55
56
57
    # prepare the slicing object which will be used for mirroring
    slice_primitive = [slice(None), ]*dimensions
    # copy the input data
Ultima's avatar
Ultima committed
58
    y = x.copy()
Ultima's avatar
Ultima committed
59
    # flip in every direction
Ultima's avatar
Ultima committed
60
61
    for i in xrange(dimensions):
        slice_picker = slice_primitive[:]
Ultima's avatar
Ultima committed
62
        slice_picker[i] = slice(1, None, None)
ultimanet's avatar
ultimanet committed
63
64
        slice_picker = tuple(slice_picker)

Ultima's avatar
Ultima committed
65
66
        slice_inverter = slice_primitive[:]
        slice_inverter[i] = slice(None, 0, -1)
ultimanet's avatar
ultimanet committed
67
        slice_inverter = tuple(slice_inverter)
Ultima's avatar
Ultima committed
68
69

        try:
70
71
            y.set_data(to_key=slice_picker, data=y,
                       from_key=slice_inverter)
Ultima's avatar
Ultima committed
72
73
74
        except(AttributeError):
            y[slice_picker] = y[slice_inverter]
    return y
75
76


77
def direct_vdot(x, y):
Ultima's avatar
Ultima committed
78
    # the input could be fields. Try to extract the data
Ultima's avatar
Ultima committed
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's avatar
Ultima committed
87
    # try to make a direct vdot
Ultima's avatar
Ultima committed
88
89
90
91
    try:
        return x.vdot(y)
    except(AttributeError):
        pass
92

Ultima's avatar
Ultima committed
93
94
95
    try:
        return y.vdot(x)
    except(AttributeError):
96
97
        pass

Ultima's avatar
Ultima committed
98
    # fallback to numpy
99
100
    return np.vdot(x, y)

Ultima's avatar
Ultima committed
101
102

def convert_nested_list_to_object_array(x):
Ultima's avatar
Ultima committed
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's avatar
Ultima committed
106
    possible_shape = np.shape(x)
Ultima's avatar
Ultima committed
107
    # Check if possible_shape goes too deep.
108
    dimension_counter = 0
Ultima's avatar
Ultima committed
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's avatar
Ultima committed
117
    # if the numpy array was not encapsulated at all, return x directly
Ultima's avatar
Ultima committed
118
119
    if real_shape == ():
        return x
Ultima's avatar
Ultima committed
120
121
    # Prepare the carrier-object
    carrier = np.empty(real_shape, dtype=np.object)
Ultima's avatar
Ultima committed
122
123
    for i in xrange(np.prod(real_shape)):
        ii = np.unravel_index(i, real_shape)
124
        try:
Ultima's avatar
Ultima committed
125
126
127
            carrier[ii] = x[ii]
        except(TypeError):
            extracted = x
128
            for j in xrange(len(ii)):
Ultima's avatar
Ultima committed
129
                extracted = extracted[ii[j]]
130
            carrier[ii] = extracted
Ultima's avatar
Ultima committed
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's avatar
Ultima committed
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's avatar
Ultima committed
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's avatar
Ultima committed
159
160
                # result[ii] = function(*map(lambda z: z[ii], args))
            return result