1 """
2 Framework for performing a Proper Orthogonal Decomposition (POD).
3
4 Useful references:
5 - http://en.wikipedia.org/wiki/Homogeneous_coordinates
6 - http://en.wikipedia.org/wiki/Transformation_matrix#Affine_transformations
7
8 Usage example:
9
10 >>> import pod
11 >>> refs = [ [-4.0, -1.0],
12 >>> [-2.0, -1.0],
13 >>> [3.0, 4.0] ]
14 >>> target = [-2.0, 1.5]
15 >>> decomposition = pod.decompose(target, refs, epsilon=1E-6, max_iter=90)
16 >>> print decomposition.get_decomposition()
17 [-1.9999991745134178, 1.4999993808850638]
18 >>> print decomposition.get_reference_weights()
19 [0.96153806466991254, 0.0, 0.61538436138874408]
20
21 The example above shows the reconstruction of the target using 3 reference
22 signals, from which only reference 1 and reference 3 are useful (reference 2
23 is assigned a weight of 0).
24
25 @author: Christophe Alexandre <ch.alexandre at bluewin dot ch>
26
27 @license:
28 Copyright(C) 2010 Christophe Alexandre
29
30 This program is free software: you can redistribute it and/or modify
31 it under the terms of the GNU Lesser General Public License as published by
32 the Free Software Foundation, either version 3 of the License, or
33 (at your option) any later version.
34
35 This program is distributed in the hope that it will be useful,
36 but WITHOUT ANY WARRANTY; without even the implied warranty of
37 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
38 GNU General Public License for more details.
39
40 You should have received a copy of the GNU Lesser General Public License
41 along with this program. If not, see <http://www.gnu.org/licenses/lgpl.txt>.
42 """
43 __all__ = ['decompose', 'combined_distance', 'DecompositionBasic']
44
45 import os
46 import logging
47
48 import vecspace
49 from util import NullHandler
50
51 _h = NullHandler()
52 _logger = logging.getLogger('pod')
53 _logger.addHandler(_h)
54
56 """
57 Distance function used for ordering the projections.
58
59 A weight of 0.0 defines the distance to the projected point, while a weight
60 of 1.0 defines the distance relative to the point generating the line.
61
62 At each iteration step the current point is projected onto the closest of
63 all lines.
64
65 @param generator_weight: how much weight is assigned to the generator point
66 @type generator_weight: float usually in range [0.0, 1.0]
67 """
68
69 def func(start_point, projected_point,
70 generator_point, w=generator_weight):
71 d1 = start_point.sub(generator_point).norm()
72 d2 = start_point.sub(projected_point).norm()
73 return w * d1 + (1.0 - w) * d2
74
75 return func
76
77 -def decompose(source, references,
78 epsilon=1E-10, max_iter=20,
79 max_factors=None, max_weight=None,
80 distance=combined_distance(0.0)
81 ):
82 """
83 Decomposing the source using the proposed reference points.
84
85 @param source: input point
86 @type source: list
87 @param references: list of reference points
88 @type references: list
89 @param epsilon: limit of the error sequence for stopping iteration
90 @type epsilon: float
91 @param max_iter: safeguard for stopping iteration
92 @type max_iter: int
93 @param max_factors: limit for the number of reference points, None allowing to use all of them
94 @type max_factors: int
95 @param distance: function used for finding the closest line to project on
96 @type distance: a function of start point, projected point, generator point
97 @return: decomposition details
98 @rtype: IterativeDecomposition
99 """
100 r = BaseDecomposition(references, epsilon, max_iter, max_factors, max_weight)
101 r.resolve(source)
102 return r
103
105 """
106 Decomposition interface definition.
107 """
108
109 - def __init__(self, references, epsilon=1E-10, max_iter=20, max_factors=None):
110 assert len(references) > 0, 'at least one reference is required'
111 dim = len(references[0])
112 for r in references:
113 assert len(r) == dim
114 self._vector_space = vecspace.VectorSpace(dim)
115 self._reference_points = []
116 for r in references:
117 ref = self._vector_space.define_point(*r)
118 self._reference_points.append(ref)
119 self._weights = dict()
120 for p in self._reference_points:
121 self._weights[p] = 0.0
122 self._epsilon = epsilon
123 self._start = None
124 self._max_iter = max_iter
125 if max_factors is None:
126 self._max_factors = len(self._reference_points)
127 else:
128 self._max_factors = max_factors
129 self._error_norm = None
130
132 """
133 Computes current decomposition result on the fly.
134
135 @return: the current relicate
136 @rtype: linalg.Point
137 """
138 decomposition = self._vector_space.origin
139 for d in self._weights.keys():
140 w_d = self._weights[d]
141 decomposition = decomposition.add(d.scale(w_d))
142 return decomposition
143
145 """
146 Iterates decomposition until convergence or max iteration is reached.
147
148 @param point: coordinates of the point to be decompositiond
149 @type point: list
150 @return: coordinates of decompositiond point
151 @rtype: list
152 """
153 pass
154
156 """
157 Returns the weight assigned to the reference provided in the constructor
158 at the indicated position.
159
160 @param position: position of the reference in the list provided in the constructor
161 @type position: int
162 """
163 ref = self._reference_points[position]
164 return self._weights[ref]
165
167 """
168 Returns the weights assigned to the references in order to construct the
169 proposed input.
170 """
171 return [self._weights[self._reference_points[i]]
172 for i in range(len(self._reference_points))]
173
175 """
176 Returns the result of the decomposition process.
177
178 @return: decomposition result
179 @rtype: list
180 """
181 return self._compute_decomposition().get_data()
182
184 """
185 Returns a measure of the decomposition error.
186
187 @return: length of the difference between the result and the initial point
188 @rtype: float
189 """
190 return self._compute_decomposition().sub(self._start).norm()
191
193 out = 'reference points:' + os.linesep
194 for p in self._reference_points:
195 out += str(p) + os.linesep
196 out += 'weightings:' + os.linesep
197 out += str(self._weights)
198 return out
199
201 """
202 Decomposition on a set of reference points.
203 """
204
205 - def __init__(self, references,
206 epsilon=1E-10,
207 max_iter=20,
208 max_factors=None,
209 max_weight=None,
210 distance=combined_distance(0.0)):
211 """
212 @param distance: function of start point, projected point and generator point
213 """
214 IterativeDecomposition.__init__(self, references, epsilon, max_iter, max_factors)
215 self._max_weight = max_weight
216 self._distance = distance
217
219 """ Projects onto the closest of the straight lines defined by
220 the reference points.
221 """
222
223 origin = self._vector_space.origin
224 if point.sub(origin).norm() <= self._epsilon:
225
226 return point
227
228 _logger.debug('projecting ' + str(point))
229 projections = dict()
230 distances = dict()
231 for ref in reference_points:
232 line = self._vector_space.define_line(ref, origin)
233 _logger.debug('computing projection onto ' + str(line))
234 ref_proj = line.project(point)
235 projections[ref] = ref_proj.projected
236 distances[ref] = self._distance(point, ref_proj.projected, ref)
237 _logger.debug('distance ' + str(distances[ref]))
238
239 ref_points = []
240 for p in reference_points:
241 if self._max_weight is None:
242 ref_points.append(p)
243 else:
244 additional_weight = projections[p].units(p)
245 suggested_weight = self._weights[p] + additional_weight
246 if abs(suggested_weight) < self._max_weight:
247 ref_points.append(p)
248
249 if len(ref_points) == 0:
250
251 return point
252
253
254 def by_dist(ref1, ref2, d=distances):
255 return cmp(d[ref1], d[ref2])
256 ref_points.sort(by_dist)
257
258 closest = ref_points[0]
259 additional_weight = projections[closest].units(closest)
260 self._weights[closest] += additional_weight
261 _logger.debug('closest driver: %s, weight=%f' % (str(closest), self._weights[closest]))
262
263 return point.sub(projections[closest])
264
266 """
267 Iterates decomposition until convergence or max iteration is reached.
268
269 @param point: coordinates of the point to be decompositiond
270 @type point: list
271 @return: coordinates of decompositiond point
272 @rtype: list
273 """
274 IterativeDecomposition.resolve(self, point)
275 _logger.debug(' ------------- STARTING PROCESS -------------')
276 target = self._vector_space.define_point(*point)
277 self._start = target
278 reference_points = self._reference_points
279 projector = self._project_point(target, reference_points)
280 diff = None
281 _logger.debug('distance to projection: %f' % projector.norm())
282 _logger.debug('drivers: ' + str(self._weights))
283 i = 0
284 while (diff is None) or (diff > self._epsilon and i < self._max_iter):
285 i = i + 1
286 previous = self._compute_decomposition()
287 _logger.debug(' ------------- ITERATION %d -------------' % i)
288 projector = self._project_point(projector, reference_points)
289 enabled_drivers = [p for p in self._weights.keys() if abs(self._weights[p]) > 0.0]
290 drivers_count = len(enabled_drivers)
291 _logger.debug('number of drivers: %d' % drivers_count)
292 if drivers_count >= self._max_factors:
293
294 _logger.debug('count limit reached for drivers: %d' % self._max_factors)
295 reference_points = enabled_drivers
296
297 decomposition = self._compute_decomposition()
298 diff = decomposition.sub(previous).norm()
299 _logger.debug('improvement: %f' % diff)
300 _logger.debug('distance to projection: %f' % projector.norm())
301 _logger.debug('drivers: ' + str(self._weights))
302 _logger.debug('decomposition: ' + str(decomposition))
303
304 _logger.debug('start:' + str(target))
305 _logger.debug('diff:' + str(decomposition.sub(target)))
306 return decomposition.get_data()
307
309
310 - def __init__(self, ref1, ref2, max_iter):
311 BaseDecomposition.__init__(self, [ref1, ref2], max_iter=max_iter)
312