Source Code for Module parsimony.datasets.simulate.utils

  1  # -*- coding: utf-8 -*- 
  2  """ 
  3  Created on Thu Sep 26 10:50:17 2013 
  4   
  5  Copyright (c) 2013-2014, CEA/DSV/I2BM/Neurospin. All rights reserved. 
  6   
  7  @author:  Tommy Löfstedt 
  8  @email:   tommy.loefstedt@cea.fr 
  9  @license: BSD 3-clause. 
 10  """ 
 11  import numpy as np 
 12   
 13  __all__ = ['TOLERANCE', 'RandomUniform', 'ConstantValue', 
 14             'norm2', 'bisection_method'] 
 15   
 16  TOLERANCE = 5e-8 
 17   
 18   
 19 -class RandomUniform(object): 
 20      """ 
 21      Example 
 22      ------- 
 23      >>> rnd = RandomUniform(-1, 1) 
 24      >>> rnd(3) #doctest: +ELLIPSIS 
 25      array([...]) 
 26      >>> rnd(2, 2) #doctest: +ELLIPSIS 
 27      array([[..., ...], 
 28             [..., ...]]) 
 29      """ 
 30 -    def __init__(self, a=0, b=1): 
 31   
 32          self.a = float(a) 
 33          self.b = float(b) 
 34   
 35 -    def rand(self, *d): 
 36   
 37          R = np.random.rand(*d) 
 38          R = R * (self.b - self.a) + self.a 
 39   
 40          return R 
 41   
 42 -    def __call__(self, *d): 
 43   
 44          return self.rand(*d) 
 45   
 46   
 47 -class ConstantValue(object): 
 48      """Random-like number generator that returns a constant value. 
 49   
 50      Example 
 51      ------- 
 52      >>> rnd = ConstantValue(5.) 
 53      >>> rnd(3) 
 54      array([ 5.,  5.,  5.]) 
 55      >>> rnd(2, 2) 
 56      array([[ 5.,  5.], 
 57             [ 5.,  5.]]) 
 58      """ 
 59 -    def __init__(self, val): 
 60   
 61          self.val = val 
 62   
 63 -    def __call__(self, *shape): 
 64   
 65          return np.repeat(self.val, np.prod(shape)).reshape(shape) 
 66   
 67   
 68  #def U(a, b): 
 69  # 
 70  #    t = max(a, b) 
 71  #    a = float(min(a, b)) 
 72  #    b = float(t) 
 73  #    return (np.random.rand() * (b - a)) + a 
 74   
 75   
 76 -def norm2(x): 
 77   
 78      return np.sqrt(np.sum(x ** 2.0)) 
 79   
 80   
 81 -def bisection_method(f, low=0.0, high=1.0, maxiter=30, eps=TOLERANCE): 
 82      """ Finds the value of x such that |f(x)|_2 < eps, for x > 0 and where 
 83      |.|_2 is the 2-norm. 
 84   
 85      Parameters 
 86      ---------- 
 87      f : The function for which a root is found. The function must be increasing 
 88              for increasing x, and decresing for decreasing x. 
 89   
 90      low : A value for which f(low) < 0. If f(low) >= 0, a lower value, low', 
 91              will be found such that f(low') < 0 and used instead of low. 
 92   
 93      high : A value for which f(high) > 0. If f(high) <= 0, a higher value, 
 94              high', will be found such that f(high') < 0 and used instead of 
 95              high. 
 96   
 97      maxiter : The maximum number of iterations. 
 98   
 99      eps : A positive value such that |f(x)|_2 < eps. Only guaranteed if 
100              |f(x)|_2 < eps in less than maxiter iterations. 
101   
102      Returns 
103      ------- 
104      x : The value such that |f(x)|_2 < eps. 
105      """ 
106      # Find start values. If the low and high 
107      # values are feasible this will just break 
108      for i in xrange(maxiter * 2): 
109          l = f(low) 
110          h = f(high) 
111   
112          if l < 0 and h > 0: 
113              break 
114   
115          if l >= 0: 
116              low /= 2.0 
117          if h <= 0: 
118              high *= 2.0 
119   
120      # Use the bisection method to find where |f(x)|_2 < eps. 
121      for i in xrange(maxiter): 
122          mid = (low + high) / 2.0 
123          val = f(mid) 
124          if val < 0: 
125              low = mid 
126          if val > 0: 
127              high = mid 
128   
129          if np.sqrt(np.sum((high - low) ** 2.0)) <= eps: 
130              break 
131   
132      return (low + high) / 2.0 
133   
134  if __name__ == "__main__": 
135      import doctest 
136      doctest.testmod() 
137