dft_python/Kernel.py at master · mathisrichter/dft_python · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
import math
import random
from scipy import ndimage
import numpy
import copy

def convolve(input, kernel):
    convolution_result = copy.copy(input)
    for dimension_index in range(kernel.get_dimensionality()):
        ndimage.convolve1d(convolution_result, \
                           kernel.get_separated_kernel_part(dimension_index), \
                           axis = dimension_index, \
                           output = convolution_result, \
                           mode = 'wrap')

    return convolution_result


class Kernel:
    "n-dimensional kernel"

    def __init__(self, amplitude):
        self._dimensionality = None
        self._amplitude = amplitude

    def get_amplitude(self):
        return self._amplitude

    def get_separated_kernel_part(self, dimension_index):
        self._check_dimension_index(dimension_index)
        return self._separated_kernel_parts[dimension_index]

    def get_separated_kernel_parts(self):
        return self._separated_kernel_parts

    def get_dimension_size(self, dimension_index):
        return self._dimension_sizes[dimension_index]

    def get_dimension_sizes(self):
        return self._dimension_sizes

    def get_dimensionality(self):
        return self._dimensionality

    def set_amplitude(self, amplitude):
        self._amplitude = amplitude
        self._calculate_kernel()

    def _check_dimension_index(self, dimension_index):
        if not (dimension_index >= 0 and dimension_index < self._dimensionality):
            print("Error. Kernel only has", self._dimensionality, "dimensions.")

    def _calculate_kernel(self):
        pass


class GaussKernel(Kernel):
    "n-dimensional Gauss kernel"

    def __init__(self, amplitude, widths, shifts=None):
        Kernel.__init__(self, amplitude)
        self._widths = widths
        self._shifts = None
        self._dimensionality = len(self._widths)
        self._dimension_sizes = None
        self._limit = 0.1
        self._separated_kernel_parts = None

        if (shifts is None):
            self._shifts = [0.0] * self._dimensionality
        else:
            if (len(self._widths) != len(self._shifts)):
                print("Error. Number of shift and width values does not match.")

        self._calculate_separated_kernel_parts()

    def get_width(self, dimension_index):
        self._check_dimension_index(dimension_index)
        return self._widths[dimension_index]

    def get_widths(self):
        return self._widths

    def get_shift(self, dimension_index):
        self._check_dimension_index(dimension_index)
        return self._shifts[dimension_index]

    def get_shifts(self):
        return self._shifts

    def set_width(self, width, dimension_index):
        self._check_dimension_index(dimension_index)
        self._widths[dimension_index] = width
        self._calculate_separated_kernel_parts()

    def set_shift(self, shift, dimension_index):
        self._check_dimension_index(dimension_index)
        self._shifts = shifts
        self._calculate_separated_kernel_parts()

    def _calculate_dimension_size(self, dimension_index):
        if not (dimension_index >= 0 and dimension_index < self._dimensionality):
            print("Error. Kernel only supports ", self._dimensionality, " dimensions.")

        dimension_width = self._widths[dimension_index]

        if (dimension_width < 10000 and dimension_width > 0):
            dimension_size = int(round(math.sqrt(2.0 * math.pow(dimension_width, 2.0) \
                                 * math.log(math.fabs(self._amplitude) / self._limit))) + 1)
        else:
            print("Error. Selected mode with is not in the proper bounds (0 < width < 10000).")

        if (dimension_size % 2 == 0):
            dimension_size += 1

        return dimension_size

    def _calculate_kernel(self):
        self._calculate_separated_kernel_parts()

    def _calculate_separated_kernel_parts(self):
        if (self._separated_kernel_parts is not None):
            del(self._separated_kernel_parts[:])
        else:
            self._separated_kernel_parts = []

        for dimension_index in range(self._dimensionality):
            dimension_size = self._calculate_dimension_size(dimension_index)

            center = (dimension_size / 2.0) + self._shifts[dimension_index]

            kernel_part = numpy.zeros(shape=dimension_size)
            ramp = numpy.linspace(0, dimension_size, dimension_size)
            for i in range(dimension_size):
                kernel_part[i] = math.exp(-math.pow(ramp[i] - center, 2.0) / \
                                         (2.0 * math.pow(self._widths[dimension_index], 2.0)))

            # normalize kernel part
            kernel_part *= 1.0 / kernel_part.sum()

            # multiply the first kernel part with the amplitude.
            # when convolving with all separated kernel parts, this will lead
            # to the correct amplitude value for the "whole kernel"
            if (dimension_index == 0):
                kernel_part *= self._amplitude

            self._separated_kernel_parts.append(kernel_part)

class BoxKernel(Kernel):
    "n-dimensional box kernel"

    def __init__(self, amplitude = 5.0):
        Kernel.__init__(self, amplitude)
        self._dimensionality = 1
        self._calculate_kernel()

    def get_separated_kernel_part(self, dimension_index):
        return self._kernel

    def _calculate_kernel(self):
        self._kernel = numpy.ones(shape=(1)) * self._amplitude