Changeset 862

Timestamp:

04/18/16 14:14:14 (4 weeks ago)

Author:

mmckerns

Message:

add mystic.math.Distribution; use it for InitialPoints?, ensemble Distribution

Location:

mystic

Files:

• _math/__init__.py (modified) (2 diffs)
• _math/grid.py (modified) (3 diffs)
• _math/samples.py (modified) (6 diffs)
• mystic/abstract_ensemble_solver.py (modified) (3 diffs)
• mystic/abstract_solver.py (modified) (2 diffs)
• mystic/ensemble.py (modified) (7 diffs)

mystic/_math/__init__.py

@@ -29 +29 @@
     samplepts    -- generate a set of randomly sampled points 
     almostEqual  -- test if equal within some absolute or relative tolerance
+    Distribution -- generate a sampling distribution instance
 
 
@@ -43 +44 @@
 import distance as paramtrans
 
+
+# distribution object
+class Distribution(object):
+    """
+Sampling distribution for mystic optimizers
+    """
+    def __init__(self, generator=None, *args, **kwds):
+        """
+generate a sampling distribution with interface dist(size=None)
+
+input::
+    - generator: a 'distribution' method from scipy.stats or numpy.random
+    - args: positional arguments for the distribtution object
+    - kwds: keyword arguments for the distribution object
+
+note::
+    this method only accepts numpy.random methods with the keyword 'size'
+        """
+        from mystic.tools import random_state
+        rng = random_state(module='numpy.random')
+        if generator is None: generator = rng.random
+        if getattr(generator, 'rvs', False): 
+            d = generator(*args, **kwds)
+            self.rvs = lambda size=None: d.rvs(size=size, random_state=rng)
+        else:
+            d = getattr(rng, generator.__name__)
+            self.rvs = lambda size=None: d(size=size, *args, **kwds)
+        return
+    def __call__(self, size=None):
+        """generate a sample of given size (tuple) from the distribution"""
+        return self.rvs(size)
+
 # end of file

mystic/_math/grid.py

@@ -9 +9 @@
 """
 
-def gridpts(q):
+def gridpts(q, dist=None):
     """
 takes a list of lists of arbitrary length q = [[1,2],[3,4]]
 and produces a list of gridpoints g = [[1,3],[1,4],[2,3],[2,4]]
+
+Note:
+    if a mystic.math.Distribution is provided, use it to inject randomness
     """
     w = [[] for i in range(len(q[-1]))]
@@ -20 +23 @@
           w[l].append(q[j][k])
       if j: w += [i[:] for i in w[:]*(len(q[j-1])-1)]
-    return [list(reversed(w[i])) for i in range(len(w))]
+    pts = [list(reversed(w[i])) for i in range(len(w))]
+    # inject some randomness
+    if dist is None: return pts
+    if not len(pts): return pts
+    pts += dist((len(pts),len(pts[0])))
+    return pts.tolist()
 
 
-def samplepts(lb,ub,npts):
+def samplepts(lb,ub,npts,dist=None):
     """
 takes lower and upper bounds (e.g. lb = [0,3], ub = [2,4])
@@ -32 +40 @@
     ub  --  a list of the upper bounds
     npts  --  number of sample points
+    dist -- a mystic.math.Distribution instance
     """
     from mystic.math.samples import random_samples
-    q = random_samples(lb,ub,npts)
-    q = [list(i) for i in q]
-    q = zip(*q)
-    return [list(i) for i in q]
+    q = random_samples(lb,ub,npts,dist)
+    return q.T.tolist()
+   #q = [list(i) for i in q]
+   #q = zip(*q)
+   #return [list(i) for i in q]
 
 

mystic/_math/samples.py

@@ -14 +14 @@
 
 # SAMPLING #
-def random_samples(lb,ub,npts=10000):
+def _random_samples(lb,ub,npts=10000):
   """
 generate npts random samples between given lb & ub
@@ -32 +32 @@
 
 
+def random_samples(lb,ub, npts=10000, dist=None, clip=False):
+  """
+generate npts samples from the given distribution between given lb & ub
+
+Inputs:
+    dist  --  a mystic.tools.Distribution instance
+    lower bounds  --  a list of the lower bounds
+    upper bounds  --  a list of the upper bounds
+    npts  --  number of sample points [default = 10000]
+    clip  --  if True, clip at bounds, else resample [default = False]
+"""
+  if dist is None:
+    return _random_samples(lb,ub, npts)
+  import numpy as np
+  dim = len(lb)
+  pts = dist((npts,dim)) # transpose of desired shape
+  pts = np.clip(pts, lb, ub).T
+  if clip: return pts  #XXX: returns a numpy.array
+  bad = ((pts.T == lb) + (pts.T == ub)).T
+  new = bad.sum()
+  _n, n = 1, 1000 #FIXME: fixed number of max tries
+  while new:
+    if _n == n: #XXX: slows the while loop...
+      raise RuntimeError('bounds could not be applied in %s iterations' % n)
+    pts[bad] = dist(new)
+    pts = np.clip(pts.T, lb, ub).T
+    bad = ((pts.T == lb) + (pts.T == ub)).T
+    new = bad.sum()
+    _n += 1
+  return pts  #XXX: returns a numpy.array
+
+
 def sample(f,lb,ub,npts=10000):
   """
@@ -43 +75 @@
 """
   from numpy import transpose
-  pts = random_samples(lb, ub, npts)
+  pts = _random_samples(lb, ub, npts)
 
   failure = 0; success = 0
@@ -68 +100 @@
   from numpy import inf, transpose
   from mystic.tools import wrap_bounds
-  pts = random_samples(lb, ub, npts)
+  pts = _random_samples(lb, ub, npts)
   f = wrap_bounds(f,lb,ub)
   ave = 0; count = 0
@@ -111 +143 @@
     npts -- the number of points to sample [Default is npts=10000]
 """
-  pts = random_samples(lb, ub, npts)
+  pts = _random_samples(lb, ub, npts)
   return _pof_given_samples(f, pts)
 
@@ -197 +229 @@
     upper = [100.0, 30.0, 2.8]
 
-    pts = random_samples(lower,upper,num_sample_points)
+    pts = _random_samples(lower,upper,num_sample_points)
     print "randomly sampled points\nbetween %s and %s" % (lower, upper)
     print pts

mystic/mystic/abstract_ensemble_solver.py

@@ -118 +118 @@
         # default settings for nested optimization
         #XXX: move nbins and npts to _InitialPoints?
+        self._dist = None #kwds['dist'] if 'dist' in kwds else None
         nbins = kwds['nbins'] if 'nbins' in kwds else [1]*dim
         if isinstance(nbins, int):
@@ -206 +207 @@
         raise NotImplementedError, "must be overwritten..."
 
-    def SetDistributionInitialPoints(self, dist):
+    def SetSampledInitialPoints(self, dist=None):
         """Generate Random Initial Points from Distribution (dist)
 
 input::
-    - dist: a scipy.stats distribution instance
+    - dist: a mystic.math.Distribution instance
 
 *** this method must be overwritten ***"""
@@ -272 +273 @@
         return
 
+    def SetDistribution(self, dist=None):
+        """Set the distribution used for determining solver starting points
+
+Inputs:
+    - dist: a mystic.math.Distribution instance
+"""
+        from mystic.math import Distribution
+        if dist and Distribution not in dist.__class__.mro():
+            dist = Distribution(dist) #XXX: or throw error?
+        self._dist = dist
+        return
+
     def _InitialPoints(self):
         """Generate a grid of starting points for the ensemble of optimizers

mystic/mystic/abstract_solver.py

@@ -411 +411 @@
         """
         from mystic.tools import random_state
-        prng = random_state(module='numpy.random')
+        rng = random_state(module='numpy.random')
         assert(len(mean) == self.nDim)
         if var is None:
@@ -423 +423 @@
                 var = var * numpy.eye(self.nDim)
         for i in range(len(self.population)):
-            self.population[i] = prng.multivariate_normal(mean, var).tolist()
-        return
-
-    def SetDistributionInitialPoints(self, dist, *args):
+            self.population[i] = rng.multivariate_normal(mean, var).tolist()
+        return
+
+    def SetSampledInitialPoints(self, dist=None):
         """Generate Random Initial Points from Distribution (dist)
 
 input::
-    - dist: a scipy.stats distribution instance (or a numpy.random method)
-
-note::
-    this method only accepts numpy.random methods with the keyword 'size'"""
-        from mystic.tools import random_state
-        prng = random_state(module='numpy.random')
-        if not args and type(dist).__name__ is 'rv_frozen': #from scipy.stats
-            for i in range(self.nPop):
-                self.population[i] = dist.rvs(self.nDim, random_state=prng).tolist()
-        else: #from numpy.random with *args and size in kwds
-            for i in range(self.nPop): 
-                self.population[i] = dist(*args, size=self.nDim).tolist()
+    - dist: a mystic.math.Distribution instance
+"""
+        from mystic.math import Distribution
+        if dist is None:
+            dist = Distribution()
+        elif type(Distribution) not in dist.__class__.mro():
+            dist = Distribution(dist) #XXX: or throw error?
+        for i in range(self.nPop):
+            self.population[i] = dist(self.nDim)
         return
 

mystic/mystic/ensemble.py

@@ -38 +38 @@
 class LatticeSolver(AbstractEnsembleSolver):
     """
-parallel mapped optimization starting from the center of N grid points
+parallel mapped optimization starting from the centers of N grid points
     """
     def __init__(self, dim, nbins=8):
@@ -72 +72 @@
         # build a grid of starting points
         from mystic.math import gridpts
-        return gridpts(bins)
+        return gridpts(bins, self._dist)
 
 
 class BuckshotSolver(AbstractEnsembleSolver):
     """
-parallel mapped optimization starting from the N random points
+parallel mapped optimization starting from N uniform randomly sampled points
     """
     def __init__(self, dim, npts=8):
@@ -104 +104 @@
         # build a grid of starting points
         from mystic.math import samplepts
-        return samplepts(lower,upper,npts)
+        return samplepts(lower,upper,npts, self._dist)
 
 
@@ -150 +150 @@
         This function should return y', with y' == 0 when the encoded
         constraints are satisfied, and y' > 0 otherwise.
+    dist -- an optional mystic.math.Distribution instance.  If provided,
+        this distribution generates randomness in ensemble starting position.
 
 Returns: (xopt, {fopt, iter, funcalls, warnflag, allfuncalls}, {allvecs})
@@ -186 +188 @@
     solver.SetEvaluationMonitor(evalmon)
     solver.SetGenerationMonitor(stepmon)
+    if 'dist' in kwds:
+        solver.SetDistribution(kwds['dist'])
     if 'penalty' in kwds:
         solver.SetPenalty(kwds['penalty'])
@@ -271 +275 @@
         This function should return y', with y' == 0 when the encoded
         constraints are satisfied, and y' > 0 otherwise.
+    dist -- an optional mystic.math.Distribution instance.  If provided,
+        this distribution generates randomness in ensemble starting position.
 
 Returns: (xopt, {fopt, iter, funcalls, warnflag, allfuncalls}, {allvecs})
@@ -307 +313 @@
     solver.SetEvaluationMonitor(evalmon)
     solver.SetGenerationMonitor(stepmon)
+    if 'dist' in kwds:
+        solver.SetDistribution(kwds['dist'])
     if 'penalty' in kwds:
         solver.SetPenalty(kwds['penalty'])