Package pod
[frames] | no frames]

Source Code for Package pod

  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           
55 -def combined_distance(generator_weight):
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
104 -class IterativeDecomposition(object):
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
131 - def _compute_decomposition(self):
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
144 - def resolve(self, point):
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
155 - def get_reference_weight(self, position):
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
166 - def get_reference_weights(self):
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
174 - def get_decomposition(self):
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
183 - def get_error_norm(self):
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
192 - def __repr__(self):
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
200 -class BaseDecomposition(IterativeDecomposition):
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
218 - def _project_point(self, point, reference_points):
219 """ Projects onto the closest of the straight lines defined by 220 the reference points. 221 """ 222 # computes projection of source point to each subspace defined by ref points 223 origin = self._vector_space.origin 224 if point.sub(origin).norm() <= self._epsilon: 225 # already matched: do nothing 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 # no eligible point left 251 return point 252 253 # finds main driver (shortest distance to ref line) 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
265 - def resolve(self, point):
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 # limits number of drivers 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
308 -class DecompositionBasic(BaseDecomposition):
309
310 - def __init__(self, ref1, ref2, max_iter):
311 BaseDecomposition.__init__(self, [ref1, ref2], max_iter=max_iter)
312