Changeset 862


Ignore:
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:
6 edited

Legend:

Unmodified
Added
Removed
  • mystic/_math/__init__.py

    r855 r862  
    2929    samplepts    -- generate a set of randomly sampled points  
    3030    almostEqual  -- test if equal within some absolute or relative tolerance 
     31    Distribution -- generate a sampling distribution instance 
    3132 
    3233 
     
    4344import distance as paramtrans 
    4445 
     46 
     47# distribution object 
     48class Distribution(object): 
     49    """ 
     50Sampling distribution for mystic optimizers 
     51    """ 
     52    def __init__(self, generator=None, *args, **kwds): 
     53        """ 
     54generate a sampling distribution with interface dist(size=None) 
     55 
     56input:: 
     57    - generator: a 'distribution' method from scipy.stats or numpy.random 
     58    - args: positional arguments for the distribtution object 
     59    - kwds: keyword arguments for the distribution object 
     60 
     61note:: 
     62    this method only accepts numpy.random methods with the keyword 'size' 
     63        """ 
     64        from mystic.tools import random_state 
     65        rng = random_state(module='numpy.random') 
     66        if generator is None: generator = rng.random 
     67        if getattr(generator, 'rvs', False):  
     68            d = generator(*args, **kwds) 
     69            self.rvs = lambda size=None: d.rvs(size=size, random_state=rng) 
     70        else: 
     71            d = getattr(rng, generator.__name__) 
     72            self.rvs = lambda size=None: d(size=size, *args, **kwds) 
     73        return 
     74    def __call__(self, size=None): 
     75        """generate a sample of given size (tuple) from the distribution""" 
     76        return self.rvs(size) 
     77 
    4578# end of file 
  • mystic/_math/grid.py

    r855 r862  
    99""" 
    1010 
    11 def gridpts(q): 
     11def gridpts(q, dist=None): 
    1212    """ 
    1313takes a list of lists of arbitrary length q = [[1,2],[3,4]] 
    1414and produces a list of gridpoints g = [[1,3],[1,4],[2,3],[2,4]] 
     15 
     16Note: 
     17    if a mystic.math.Distribution is provided, use it to inject randomness 
    1518    """ 
    1619    w = [[] for i in range(len(q[-1]))] 
     
    2023          w[l].append(q[j][k]) 
    2124      if j: w += [i[:] for i in w[:]*(len(q[j-1])-1)] 
    22     return [list(reversed(w[i])) for i in range(len(w))] 
     25    pts = [list(reversed(w[i])) for i in range(len(w))] 
     26    # inject some randomness 
     27    if dist is None: return pts 
     28    if not len(pts): return pts 
     29    pts += dist((len(pts),len(pts[0]))) 
     30    return pts.tolist() 
    2331 
    2432 
    25 def samplepts(lb,ub,npts): 
     33def samplepts(lb,ub,npts,dist=None): 
    2634    """ 
    2735takes lower and upper bounds (e.g. lb = [0,3], ub = [2,4]) 
     
    3240    ub  --  a list of the upper bounds 
    3341    npts  --  number of sample points 
     42    dist -- a mystic.math.Distribution instance 
    3443    """ 
    3544    from mystic.math.samples import random_samples 
    36     q = random_samples(lb,ub,npts) 
    37     q = [list(i) for i in q] 
    38     q = zip(*q) 
    39     return [list(i) for i in q] 
     45    q = random_samples(lb,ub,npts,dist) 
     46    return q.T.tolist() 
     47   #q = [list(i) for i in q] 
     48   #q = zip(*q) 
     49   #return [list(i) for i in q] 
    4050 
    4151 
  • mystic/_math/samples.py

    r855 r862  
    1414 
    1515# SAMPLING # 
    16 def random_samples(lb,ub,npts=10000): 
     16def _random_samples(lb,ub,npts=10000): 
    1717  """ 
    1818generate npts random samples between given lb & ub 
     
    3232 
    3333 
     34def random_samples(lb,ub, npts=10000, dist=None, clip=False): 
     35  """ 
     36generate npts samples from the given distribution between given lb & ub 
     37 
     38Inputs: 
     39    dist  --  a mystic.tools.Distribution instance 
     40    lower bounds  --  a list of the lower bounds 
     41    upper bounds  --  a list of the upper bounds 
     42    npts  --  number of sample points [default = 10000] 
     43    clip  --  if True, clip at bounds, else resample [default = False] 
     44""" 
     45  if dist is None: 
     46    return _random_samples(lb,ub, npts) 
     47  import numpy as np 
     48  dim = len(lb) 
     49  pts = dist((npts,dim)) # transpose of desired shape 
     50  pts = np.clip(pts, lb, ub).T 
     51  if clip: return pts  #XXX: returns a numpy.array 
     52  bad = ((pts.T == lb) + (pts.T == ub)).T 
     53  new = bad.sum() 
     54  _n, n = 1, 1000 #FIXME: fixed number of max tries 
     55  while new: 
     56    if _n == n: #XXX: slows the while loop... 
     57      raise RuntimeError('bounds could not be applied in %s iterations' % n) 
     58    pts[bad] = dist(new) 
     59    pts = np.clip(pts.T, lb, ub).T 
     60    bad = ((pts.T == lb) + (pts.T == ub)).T 
     61    new = bad.sum() 
     62    _n += 1 
     63  return pts  #XXX: returns a numpy.array 
     64 
     65 
    3466def sample(f,lb,ub,npts=10000): 
    3567  """ 
     
    4375""" 
    4476  from numpy import transpose 
    45   pts = random_samples(lb, ub, npts) 
     77  pts = _random_samples(lb, ub, npts) 
    4678 
    4779  failure = 0; success = 0 
     
    68100  from numpy import inf, transpose 
    69101  from mystic.tools import wrap_bounds 
    70   pts = random_samples(lb, ub, npts) 
     102  pts = _random_samples(lb, ub, npts) 
    71103  f = wrap_bounds(f,lb,ub) 
    72104  ave = 0; count = 0 
     
    111143    npts -- the number of points to sample [Default is npts=10000] 
    112144""" 
    113   pts = random_samples(lb, ub, npts) 
     145  pts = _random_samples(lb, ub, npts) 
    114146  return _pof_given_samples(f, pts) 
    115147 
     
    197229    upper = [100.0, 30.0, 2.8] 
    198230 
    199     pts = random_samples(lower,upper,num_sample_points) 
     231    pts = _random_samples(lower,upper,num_sample_points) 
    200232    print "randomly sampled points\nbetween %s and %s" % (lower, upper) 
    201233    print pts 
  • mystic/mystic/abstract_ensemble_solver.py

    r860 r862  
    118118        # default settings for nested optimization 
    119119        #XXX: move nbins and npts to _InitialPoints? 
     120        self._dist = None #kwds['dist'] if 'dist' in kwds else None 
    120121        nbins = kwds['nbins'] if 'nbins' in kwds else [1]*dim 
    121122        if isinstance(nbins, int): 
     
    206207        raise NotImplementedError, "must be overwritten..." 
    207208 
    208     def SetDistributionInitialPoints(self, dist): 
     209    def SetSampledInitialPoints(self, dist=None): 
    209210        """Generate Random Initial Points from Distribution (dist) 
    210211 
    211212input:: 
    212     - dist: a scipy.stats distribution instance 
     213    - dist: a mystic.math.Distribution instance 
    213214 
    214215*** this method must be overwritten ***""" 
     
    272273        return 
    273274 
     275    def SetDistribution(self, dist=None): 
     276        """Set the distribution used for determining solver starting points 
     277 
     278Inputs: 
     279    - dist: a mystic.math.Distribution instance 
     280""" 
     281        from mystic.math import Distribution 
     282        if dist and Distribution not in dist.__class__.mro(): 
     283            dist = Distribution(dist) #XXX: or throw error? 
     284        self._dist = dist 
     285        return 
     286 
    274287    def _InitialPoints(self): 
    275288        """Generate a grid of starting points for the ensemble of optimizers 
  • mystic/mystic/abstract_solver.py

    r861 r862  
    411411        """ 
    412412        from mystic.tools import random_state 
    413         prng = random_state(module='numpy.random') 
     413        rng = random_state(module='numpy.random') 
    414414        assert(len(mean) == self.nDim) 
    415415        if var is None: 
     
    423423                var = var * numpy.eye(self.nDim) 
    424424        for i in range(len(self.population)): 
    425             self.population[i] = prng.multivariate_normal(mean, var).tolist() 
    426         return 
    427  
    428     def SetDistributionInitialPoints(self, dist, *args): 
     425            self.population[i] = rng.multivariate_normal(mean, var).tolist() 
     426        return 
     427 
     428    def SetSampledInitialPoints(self, dist=None): 
    429429        """Generate Random Initial Points from Distribution (dist) 
    430430 
    431431input:: 
    432     - dist: a scipy.stats distribution instance (or a numpy.random method) 
    433  
    434 note:: 
    435     this method only accepts numpy.random methods with the keyword 'size'""" 
    436         from mystic.tools import random_state 
    437         prng = random_state(module='numpy.random') 
    438         if not args and type(dist).__name__ is 'rv_frozen': #from scipy.stats 
    439             for i in range(self.nPop): 
    440                 self.population[i] = dist.rvs(self.nDim, random_state=prng).tolist() 
    441         else: #from numpy.random with *args and size in kwds 
    442             for i in range(self.nPop):  
    443                 self.population[i] = dist(*args, size=self.nDim).tolist() 
     432    - dist: a mystic.math.Distribution instance 
     433""" 
     434        from mystic.math import Distribution 
     435        if dist is None: 
     436            dist = Distribution() 
     437        elif type(Distribution) not in dist.__class__.mro(): 
     438            dist = Distribution(dist) #XXX: or throw error? 
     439        for i in range(self.nPop): 
     440            self.population[i] = dist(self.nDim) 
    444441        return 
    445442 
  • mystic/mystic/ensemble.py

    r860 r862  
    3838class LatticeSolver(AbstractEnsembleSolver): 
    3939    """ 
    40 parallel mapped optimization starting from the center of N grid points 
     40parallel mapped optimization starting from the centers of N grid points 
    4141    """ 
    4242    def __init__(self, dim, nbins=8): 
     
    7272        # build a grid of starting points 
    7373        from mystic.math import gridpts 
    74         return gridpts(bins) 
     74        return gridpts(bins, self._dist) 
    7575 
    7676 
    7777class BuckshotSolver(AbstractEnsembleSolver): 
    7878    """ 
    79 parallel mapped optimization starting from the N random points 
     79parallel mapped optimization starting from N uniform randomly sampled points 
    8080    """ 
    8181    def __init__(self, dim, npts=8): 
     
    104104        # build a grid of starting points 
    105105        from mystic.math import samplepts 
    106         return samplepts(lower,upper,npts) 
     106        return samplepts(lower,upper,npts, self._dist) 
    107107 
    108108 
     
    150150        This function should return y', with y' == 0 when the encoded 
    151151        constraints are satisfied, and y' > 0 otherwise. 
     152    dist -- an optional mystic.math.Distribution instance.  If provided, 
     153        this distribution generates randomness in ensemble starting position. 
    152154 
    153155Returns: (xopt, {fopt, iter, funcalls, warnflag, allfuncalls}, {allvecs}) 
     
    186188    solver.SetEvaluationMonitor(evalmon) 
    187189    solver.SetGenerationMonitor(stepmon) 
     190    if 'dist' in kwds: 
     191        solver.SetDistribution(kwds['dist']) 
    188192    if 'penalty' in kwds: 
    189193        solver.SetPenalty(kwds['penalty']) 
     
    271275        This function should return y', with y' == 0 when the encoded 
    272276        constraints are satisfied, and y' > 0 otherwise. 
     277    dist -- an optional mystic.math.Distribution instance.  If provided, 
     278        this distribution generates randomness in ensemble starting position. 
    273279 
    274280Returns: (xopt, {fopt, iter, funcalls, warnflag, allfuncalls}, {allvecs}) 
     
    307313    solver.SetEvaluationMonitor(evalmon) 
    308314    solver.SetGenerationMonitor(stepmon) 
     315    if 'dist' in kwds: 
     316        solver.SetDistribution(kwds['dist']) 
    309317    if 'penalty' in kwds: 
    310318        solver.SetPenalty(kwds['penalty']) 
Note: See TracChangeset for help on using the changeset viewer.