Changeset 569


Ignore:
Timestamp:
09/24/12 15:59:25 (4 years ago)
Author:
mmckerns
Message:

merge of dirac_measure plus compressed, paramtrans, and legacydata to _math

Files:
14 deleted
14 edited
4 moved

Legend:

Unmodified
Added
Removed
  • branches/UQ/math/legacy/MINMAX_StAlData.py

    r556 r569  
    1717# the dataset 
    1818####################################################################### 
    19 from legacydata import load_dataset 
     19from mystic.math.legacydata import load_dataset 
    2020data = load_dataset('StAlDataset.txt') 
    2121 
     
    6767def maximize(params,npts,bounds): 
    6868 
    69   from dirac_measure import scenario 
     69  from mystic.math.dirac_measure import scenario 
    7070  from numpy import inf 
    7171  target,error = params 
     
    9696    ##################### begin function-specific ##################### 
    9797    # impose norm on the weights of the discrete measures 
    98     from dirac_measure import norm_wts_constraintsFactory as factory 
     98    from mystic.math.dirac_measure import norm_wts_constraintsFactory as factory 
    9999    constrain = factory(npts) 
    100100 
     
    158158  nz = 1  #NOTE: SET THE NUMBER OF 'v' POINTS HERE! 
    159159  npts = (nx,ny,nz) 
    160   from paramtrans import _npts 
     160  from mystic.math.paramtrans import _npts 
    161161  _n = _npts(npts) 
    162162 
     
    219219 
    220220  from numpy import array 
    221   from dirac_measure import scenario 
     221  from mystic.math.dirac_measure import scenario 
    222222  c = scenario() 
    223223  c.load(solved,npts) 
  • branches/UQ/math/legacy/MM_OUQ_StAlData.py

    r556 r569  
    2323# the dataset 
    2424####################################################################### 
    25 from legacydata import load_dataset 
     25from mystic.math.legacydata import load_dataset 
    2626data = load_dataset('StAlDataset.txt') 
    2727 
     
    7575def maximize(params,npts,bounds): 
    7676 
    77   from dirac_measure import scenario 
     77  from mystic.math.dirac_measure import scenario 
    7878  from numpy import inf 
    7979  target,error = params 
     
    104104    ##################### begin function-specific ##################### 
    105105    # impose mean on the values of the product measure 
    106     from dirac_measure import mean_y_norm_wts_constraintsFactory as factory 
     106    from mystic.math.dirac_measure import mean_y_norm_wts_constraintsFactory as factory 
    107107    constrain = factory((target[0],error[0]), npts) 
    108108 
     
    187187  nz = 2  #NOTE: SET THE NUMBER OF 'v' POINTS HERE! 
    188188  npts = (nx,ny,nz) 
    189   from paramtrans import _npts 
     189  from mystic.math.paramtrans import _npts 
    190190  _n = _npts(npts) 
    191191 
     
    250250 
    251251  from numpy import array 
    252   from dirac_measure import scenario 
     252  from mystic.math.dirac_measure import scenario 
    253253  c = scenario() 
    254254  c.load(solved,npts) 
  • branches/UQ/math/legacy/TEST_OUQ_1dData.py

    r556 r569  
    2121# the dataset 
    2222####################################################################### 
    23 from legacydata import load_dataset 
     23from mystic.math.legacydata import load_dataset 
    2424data = load_dataset('ExampleDataset.txt', range(0,1)) 
    2525 
     
    7171def maximize(params,npts,bounds): 
    7272 
    73   from dirac_measure import scenario 
     73  from mystic.math.dirac_measure import scenario 
    7474  from numpy import inf 
    7575  target,error = params 
     
    100100    ##################### begin function-specific ##################### 
    101101    # impose mean on the values of the product measure 
    102     from dirac_measure import mean_y_norm_wts_constraintsFactory as factory 
     102    from mystic.math.dirac_measure import mean_y_norm_wts_constraintsFactory as factory 
    103103    constrain = factory((target[0],error[0]), npts) 
    104104 
     
    183183  nz = 1  #NOTE: SET THE NUMBER OF 'v' POINTS HERE! 
    184184  npts = (nx,ny,nz) 
    185   from paramtrans import _npts 
     185  from mystic.math.paramtrans import _npts 
    186186  _n = _npts(npts) 
    187187 
     
    246246 
    247247  from numpy import array 
    248   from dirac_measure import scenario 
     248  from mystic.math.dirac_measure import scenario 
    249249  c = scenario() 
    250250  c.load(solved,npts) 
  • branches/UQ/math/legacy/TEST_OUQ_1dSurr_CxCy.py

    r557 r569  
    7070def maximize(params,npts,bounds): 
    7171 
    72   from dirac_measure import scenario 
     72  from mystic.math.dirac_measure import scenario 
    7373  from numpy import inf 
    7474  target,error = params 
     
    9999    ##################### begin function-specific ##################### 
    100100    # impose mean on the values of the product measure 
    101     from dirac_measure import mean_y_norm_wts_constraintsFactory as factory 
     101    from mystic.math.dirac_measure import mean_y_norm_wts_constraintsFactory as factory 
    102102    constrain = factory((target[0],error[0]), npts) 
    103103 
     
    185185  nz = 1  #NOTE: SET THE NUMBER OF 'v' POINTS HERE! 
    186186  npts = (nx,ny,nz) 
    187   from paramtrans import _npts 
     187  from mystic.math.paramtrans import _npts 
    188188  _n = _npts(npts) 
    189189 
     
    247247 
    248248  from numpy import array 
    249   from dirac_measure import scenario 
     249  from mystic.math.dirac_measure import scenario 
    250250  c = scenario() 
    251251  c.load(solved,npts) 
     
    261261  print "sum_wts: %s == 1.0" % [sum(w) for w in c.wts] 
    262262 
    263   from paramtrans import graphical_distance 
     263  from mystic.math.paramtrans import graphical_distance 
    264264  Ry = graphical_distance(model, c, ytol=Cy, xtol=Cx, cutoff=0.0, imax=0) 
    265265  print "vertical_distance: %s <= %s" % (Ry, Cy + max(Cx)) 
  • branches/UQ/math/legacy/TEST_OUQ_1dSurr_Cy.py

    r557 r569  
    7070def maximize(params,npts,bounds): 
    7171 
    72   from dirac_measure import scenario 
     72  from mystic.math.dirac_measure import scenario 
    7373  from numpy import inf 
    7474  target,error = params 
     
    9999    ##################### begin function-specific ##################### 
    100100    # impose mean on the values of the product measure 
    101     from dirac_measure import mean_y_norm_wts_constraintsFactory as factory 
     101    from mystic.math.dirac_measure import mean_y_norm_wts_constraintsFactory as factory 
    102102    constrain = factory((target[0],error[0]), npts) 
    103103 
     
    185185  nz = 1  #NOTE: SET THE NUMBER OF 'v' POINTS HERE! 
    186186  npts = (nx,ny,nz) 
    187   from paramtrans import _npts 
     187  from mystic.math.paramtrans import _npts 
    188188  _n = _npts(npts) 
    189189 
     
    247247 
    248248  from numpy import array 
    249   from dirac_measure import scenario 
     249  from mystic.math.dirac_measure import scenario 
    250250  c = scenario() 
    251251  c.load(solved,npts) 
     
    261261  print "sum_wts: %s == 1.0" % [sum(w) for w in c.wts] 
    262262 
    263   from paramtrans import graphical_distance 
     263  from mystic.math.paramtrans import graphical_distance 
    264264  Ry = graphical_distance(model, c, ytol=Cy, xtol=Cx, cutoff=0.0, imax=0) 
    265265  print "vertical_distance: %s <= %s" % (Ry, Cy + max(Cx)) 
  • branches/UQ/math/legacy/TEST_OUQ_StAlData.py

    r558 r569  
    2121# the dataset 
    2222####################################################################### 
    23 from legacydata import load_dataset 
     23from mystic.math.legacydata import load_dataset 
    2424data = load_dataset('StAlDataset.txt') 
    2525 
     
    7171def maximize(params,npts,bounds): 
    7272 
    73   from dirac_measure import scenario 
     73  from mystic.math.dirac_measure import scenario 
    7474  from numpy import inf 
    7575  target,error = params 
     
    100100    ##################### begin function-specific ##################### 
    101101    # impose mean on the values of the product measure 
    102     from dirac_measure import mean_y_norm_wts_constraintsFactory as factory 
     102    from mystic.math.dirac_measure import mean_y_norm_wts_constraintsFactory as factory 
    103103    constrain = factory((target[0],error[0]), npts) 
    104104 
     
    184184  nz = 2  #NOTE: SET THE NUMBER OF 'v' POINTS HERE! 
    185185  npts = (nx,ny,nz) 
    186   from paramtrans import _npts 
     186  from mystic.math.paramtrans import _npts 
    187187  _n = _npts(npts) 
    188188 
     
    247247 
    248248  from numpy import array 
    249   from dirac_measure import scenario 
     249  from mystic.math.dirac_measure import scenario 
    250250  c = scenario() 
    251251  c.load(solved,npts) 
  • branches/UQ/math/legacy/TEST_OUQ_StStSurr_Cy.py

    r557 r569  
    7070def maximize(params,npts,bounds): 
    7171 
    72   from dirac_measure import scenario 
     72  from mystic.math.dirac_measure import scenario 
    7373  from numpy import inf 
    7474  target,error = params 
     
    9999    ##################### begin function-specific ##################### 
    100100    # impose mean on the values of the product measure 
    101     from dirac_measure import mean_y_norm_wts_constraintsFactory as factory 
     101    from mystic.math.dirac_measure import mean_y_norm_wts_constraintsFactory as factory 
    102102    constrain = factory((target[0],error[0]), npts) 
    103103 
     
    185185  nz = 1  #NOTE: SET THE NUMBER OF 'v' POINTS HERE! 
    186186  npts = (nx,ny,nz) 
    187   from paramtrans import _npts 
     187  from mystic.math.paramtrans import _npts 
    188188  _n = _npts(npts) 
    189189 
     
    247247 
    248248  from numpy import array 
    249   from dirac_measure import scenario 
     249  from mystic.math.dirac_measure import scenario 
    250250  c = scenario() 
    251251  c.load(solved,npts) 
     
    261261  print "sum_wts: %s == 1.0" % [sum(w) for w in c.wts] 
    262262 
    263   from paramtrans import graphical_distance 
     263  from mystic.math.paramtrans import graphical_distance 
    264264  Ry = graphical_distance(model, c, ytol=Cy, xtol=Cx, cutoff=0.0, imax=0) 
    265265  print "vertical_distance: %s <= %s" % (Ry, Cy + max(Cx)) 
  • branches/UQ/math/legacy/test_ExampleDataset.py

    r557 r569  
    77Mathematical Modeling and Numerical Analysis (submitted 2012). 
    88""" 
    9 from legacydata import load_dataset 
     9from mystic.math.legacydata import load_dataset 
    1010datafile = 'ExampleDataset.txt' 
    1111 
     
    3131 
    3232  # Check shortness 
    33   from paramtrans import lipschitz_distance, graphical_distance 
     33  from mystic.math.paramtrans import lipschitz_distance, graphical_distance 
    3434  print("\nshort: %s" % ex1d_data.short()) 
    3535  L = ex1d_data.lipschitz; 
  • branches/UQ/math/legacy/test_ModeledDataset.py

    r557 r569  
    1111 
    1212  # Build a initial dataset 
    13   from legacydata import datapoint, dataset 
     13  from mystic.math.legacydata import datapoint, dataset 
    1414  x = [[0,0,0],[1,1,1],[2,2,2],[3,3,3]] 
    1515  y = [0,2,4,9] 
  • branches/UQ/math/legacy/test_StAlDataset.py

    r557 r569  
    77Mathematical Modeling and Numerical Analysis (submitted 2012). 
    88""" 
    9 from legacydata import load_dataset 
     9from mystic.math.legacydata import load_dataset 
    1010datafile = 'StAlDataset.txt' 
    1111st_al_data = load_dataset(datafile) 
     
    2121 
    2222  # Check shortness 
    23   from paramtrans import lipschitz_distance, graphical_distance 
     23  from mystic.math.paramtrans import lipschitz_distance, graphical_distance 
    2424  print("\nshort: %s" % st_al_data.short()) 
    2525  L = st_al_data.lipschitz; 
  • branches/UQ/math/sausage/TEST_OUQ_1dSurr_diam.py

    r553 r569  
    2222# the model function and the dataset 
    2323####################################################################### 
    24 #from StAlSurrogate import st_al_surr as model 
    25 #from legacydata import load_dataset 
    26 #data = load_dataset('ExampleDataset.txt', range(0,1)) 
    27 #L = data.lipschitz 
    28  
    2924def model(x): 
    3025  return x[0]; 
     26 
    3127 
    3228####################################################################### 
     
    7672def maximize(params,npts,bounds): 
    7773 
    78   from dirac_measure import scenario 
     74  from mystic.math.dirac_measure import scenario 
    7975  from mystic.math import almostEqual 
    8076  from numpy import inf 
     
    189185  nz = 1  #NOTE: SET THE NUMBER OF 'v' POINTS HERE! 
    190186  npts = (nx,ny,nz) 
    191   from paramtrans import _npts 
     187  from mystic.math.paramtrans import _npts 
    192188  _n = _npts(npts) 
    193189 
     
    252248 
    253249  from numpy import array 
    254   from dirac_measure import scenario 
     250  from mystic.math.dirac_measure import scenario 
    255251  c = scenario() 
    256252  c.load(solved,npts) 
  • branches/UQ/math/sausage/TEST_OUQ_StStSurr.py

    r555 r569  
    2323####################################################################### 
    2424from StStSurrogate import marc_surr as model 
    25 #from legacydata import load_dataset 
    26 #data = load_dataset('ExampleDataset.txt', range(0,1)) 
    27 #L = data.lipschitz 
    2825 
    2926 
     
    7471def maximize(params,npts,bounds): 
    7572 
    76   from dirac_measure import scenario 
     73  from mystic.math.dirac_measure import scenario 
    7774  from mystic.math import almostEqual 
    7875  from numpy import inf 
     
    200197  nz = 1  #NOTE: SET THE NUMBER OF 'v' POINTS HERE! 
    201198  npts = (nx,ny,nz) 
    202   from paramtrans import _npts 
     199  from mystic.math.paramtrans import _npts 
    203200  _n = _npts(npts) 
    204201 
     
    263260 
    264261  from numpy import array 
    265   from dirac_measure import scenario 
     262  from mystic.math.dirac_measure import scenario 
    266263  c = scenario() 
    267264  c.load(solved,npts) 
  • mystic/_math/dirac_measure.py

    r508 r569  
    22""" 
    33Classes for Dirac measure data objects. 
    4 Includes point, dirac_measure, and product_measure classes. 
     4Includes point, dirac_measure, product_measure, and scenario classes. 
    55""" 
    66# Adapted from seesaw2d.py in branches/UQ/math/examples2/  
     
    1212 
    1313class point(object): 
    14   """ 1-d object with weight and position """ 
     14  """ 1-d object with weight and position 
     15 
     16 queries: 
     17  p.weight   --  returns weight 
     18  p.position  --  returns position 
     19  p.rms  --  returns the square root of sum of squared position 
     20 
     21 settings: 
     22  p.weight = w1  --  set the weight 
     23  p.position = x1  --  set the position 
     24""" 
    1525 
    1626  def __init__(self, position, weight): 
     
    2131  def __repr__(self): 
    2232    return "(%s @%s)" % (self.weight, self.position) 
     33 
     34  def __rms(self): # square root of sum of squared positions 
     35    from math import sqrt 
     36    return sqrt(sum([i**2 for i in self.position])) 
     37 
     38  # interface 
     39  rms = property(__rms) 
    2340 
    2441  pass 
     
    7289  def __mean(self): 
    7390    from mystic.math.measures import mean 
    74     return mean(self.coords, self.weights) 
     91    return mean(self.coords, self.weights)  
    7592 
    7693  def __range(self): 
     
    174191 methods: 
    175192  c.pof(f)  --  calculate the probability of failure 
     193  c.sampled_pof(f, npts) -- calculate the pof using sampled points 
    176194  c.get_expect(f)  --  calculate the expectation 
    177195  c.set_expect((center,delta), f)  --  impose expectation by adjusting positions 
    178196  c.flatten()  --  convert measure to a flat list of parameters 
    179197  c.load(params, pts)  --  'fill' the measure from a flat list of parameters 
     198  c.update(params) -- 'update' the measure from a flat list of parameters 
    180199 
    181200 notes: 
     
    185204  - weight wxi should be same for each (yj,zk) at xi; similarly for wyi & wzi 
    186205""" 
     206  def __val(self): 
     207    raise NotImplementedError, "'value' is undefined in a dirac_measure" 
    187208 
    188209  def __pts(self): 
     
    276297  #     u += self.weights[i] 
    277298  # return u  #XXX: does this need to be normalized? 
    278      
     299 
    279300  def sampled_pof(self, f, npts=10000): 
    280301    """calculate probability of failure over a given function, f, 
     
    307328    from numpy import transpose 
    308329    return transpose(pts)  #XXX: assumes 'coords' is a list of floats 
     330 
     331  def update(self, params): 
     332    """update the product measure from a list of parameters 
     333 
     334The dimensions of the product measure will not change""" 
     335    pts = self.pts 
     336    _len = 2 * sum(pts) 
     337 
     338    if len(params)  >  _len:  # if Y-values are appended to params 
     339      params, values  =  params[:_len], params[_len:] 
     340 
     341    pm = unflatten(params, pts) 
     342    zo = pm.count([]) 
     343    self[:] = pm[:len(self) - zo] + self[len(pm) - zo:] 
     344    return 
    309345 
    310346  def load(self, params, pts): 
     
    327363    pattern followed for each new dimension desired for the product measure. 
    328364""" 
    329     c = unflatten(params, pts) 
    330     for i in range(len(c)): 
    331       self.append(c[i])  
     365    _len = 2 * sum(pts) 
     366    if len(params)  >  _len:  # if Y-values are appended to params 
     367      params, values  =  params[:_len], params[_len:] 
     368 
     369    self.extend( unflatten(params, pts) ) 
    332370    return 
    333371 
     
    349387    pattern followed for each dimension of the product measure. 
    350388""" 
    351     return flatten(self) 
     389    params = flatten(self) 
     390    return params 
     391 
     392  #XXX: name stinks... better as "non_redundant"? ...is really a helper 
     393  def differs_by_one(self, ith, all=True, index=True): 
     394    """get the product measure coordinates where the associated binary 
     395string differs by exactly one index 
     396 
     397  Inputs: 
     398    ith   = the target index 
     399    all   = if False, return only the results for indicies < i 
     400    index = if True, return the index of the results (not results themselves) 
     401""" 
     402    from mystic.math.compressed import index2binary, differs_by_one 
     403    b = index2binary(range(self.npts), self.npts) 
     404    return differs_by_one(ith, b, all, index)  
     405 
     406  def select(self, *index, **kwds): 
     407    """generator for product measure coords due to selected position indicies 
     408 (NOTE: only works for product measures of dimension 2^K) 
     409 
     410  >>> r 
     411  [[9, 8], [1, 3], [4, 2]] 
     412  >>> r.select(*range(r.npts)) 
     413  [(9, 1, 4), (8, 1, 4), (9, 3, 4), (8, 3, 4), (9, 1, 2), (8, 1, 2), (9, 3, 2), (8, 3, 2)] 
     414  >>> 
     415  >>> _pack(r) 
     416  [(9, 1, 4), (8, 1, 4), (9, 3, 4), (8, 3, 4), (9, 1, 2), (8, 1, 2), (9, 3, 2), (8, 3, 2)] 
     417""" 
     418    from mystic.math.compressed import index2binary, binary2coords 
     419    v = index2binary(list(index), self.npts) 
     420    return binary2coords(v, self.pos, **kwds) 
     421    #XXX: '_pack' requires resorting ([::-1]) so that indexing is wrong. 
     422    #     Better if modify mystic's pack to match sorting of binary strings ? 
    352423 
    353424 #__center = None 
     
    368439 
    369440 
     441class scenario(product_measure):  #FIXME: meant to only accept sets... 
     442  """ a N-d product measure (collection of dirac measures) with values 
     443  s = scenario(product_measure, [value1, value2, ..., valueN])   
     444    where each point in the product measure is paried with a value 
     445    (essentially, a dataset in product_measure representation) 
     446 
     447 queries: 
     448  s.npts  --  returns total number of points 
     449  s.weights   --  returns list of weights 
     450  s.coords  --  returns list of position tuples 
     451  s.values  --  returns list of values 
     452  s.mass  --  returns list of weight norms 
     453  s.pts  --  returns number of points for each discrete measure 
     454  s.wts  --  returns list of weights for each discrete measure 
     455  s.pos  --  returns list of positions for each discrete measure 
     456 
     457 settings: 
     458  s.coords = [(x1,y1,z1),...]  --  set the positions (tuples in product measure) 
     459  s.values = [v1,v2,v3,...]  --  set the values (correspond to position tuples) 
     460 
     461 methods: 
     462  s.pof(f)  --  calculate the probability of failure 
     463  s.pof_value(f)  --  calculate the probability of failure using the values 
     464  s.sampled_pof(f, npts) -- calculate the pof using sampled points 
     465  s.get_expect(f)  --  calculate the expectation 
     466  s.set_expect((center,delta), f)  --  impose expectation by adjusting positions 
     467  s.get_mean_value()  --  calculate the mean values for a scenario 
     468  s.set_mean_value(m)  --  impose mean value by adjusting values 
     469  s.set_feasible(data)  --  impose shortness by adjusting positions and values 
     470  s.short_wrt_data(data) -- check for shortness with respect to data 
     471  s.short_wrt_self(L) -- check for shortness with respect to self 
     472  s.set_valid(model) -- impose validity by adjusting positions and values 
     473  s.valid_wrt_model(model) -- check for validity with respect to the model 
     474  s.flatten()  --  convert measure to a flat list of parameters 
     475  s.load(params, pts)  --  'fill' the measure from a flat list of parameters 
     476  s.update(params) -- 'update' the measure from a flat list of parameters 
     477 
     478 notes: 
     479  - constraints impose expect (center - delta) <= E <= (center + delta) 
     480  - constraints impose sum(weights) == 1.0 for each set 
     481  - assumes that s.npts = len(s.coords) == len(s.weights) 
     482  - weight wxi should be same for each (yj,zk) at xi; similarly for wyi & wzi 
     483""" 
     484  def __init__(self, pm=None, values=None): 
     485    super(product_measure,self).__init__() 
     486    if pm:  
     487      pm = product_measure(pm) 
     488      self.load(pm.flatten(), pm.pts) 
     489    if not values: values = [] 
     490    self.__Y = values # storage for values of s.coords 
     491    return 
     492 
     493  def __values(self): 
     494    return self.__Y 
     495 
     496  def __set_values(self, values): 
     497    self.__Y = values[:] 
     498    return 
     499 
     500  def get_mean_value(self):  # get mean of y's 
     501    """calculate the mean of the associated values for a scenario""" 
     502    from mystic.math.measures import mean 
     503    return mean(self.values, self.weights) 
     504 
     505  def set_mean_value(self, m):  # set mean of y's 
     506    """set the mean for the associated values of a scenario""" 
     507    from mystic.math.measures import impose_mean 
     508    self.values = impose_mean(m, self.values, self.weights) 
     509    return 
     510 
     511  def valid_wrt_model(self, model, blamelist=False, pairs=True, \ 
     512                                   all=False, raw=False, **kwds): 
     513    """check for scenario validity with respect to the model 
     514 
     515Inputs: 
     516    model -- the model function, y' = F(x') 
     517    blamelist -- if True, report which points are infeasible 
     518    pairs -- if True, report indicies of infeasible points 
     519    all -- if True, report results for each point (opposed to all points) 
     520    raw -- if True, report numerical results (opposed to boolean results) 
     521 
     522Additional Inputs: 
     523    ytol -- maximum acceptable difference |y - F(x')|; a single value 
     524    xtol -- maximum acceptable difference |x - x'|; an iterable or single value 
     525    cutoff -- zero out distances less than cutoff; typically: ytol, 0.0, or None 
     526 
     527Notes: 
     528    xtol defines the n-dimensional base of a pilar of height ytol, centered at 
     529    each point. The region inside the pilar defines the space where a "valid" 
     530    model must intersect. If xtol is not specified, then the base of the pilar 
     531    will be a dirac at x' = x. This function performs an optimization for each 
     532    x to find an appropriate x'. While cutoff and ytol are very tightly related, 
     533    they play a distinct role; ytol is used to set the optimization termination 
     534    for an acceptable |y - F(x')|, while cutoff is applied post-optimization. 
     535""" 
     536    from mystic.math.legacydata import dataset  
     537    data = dataset()  
     538    data.load(self.coords, self.values) 
     539   #data.lipschitz = L 
     540    for i in range(len(data)): 
     541      data[i].id = i 
     542    return data.valid(model, blamelist=blamelist, pairs=pairs, \ 
     543                                       all=all, raw=raw, **kwds) 
     544 
     545  def short_wrt_self(self, L, blamelist=False, pairs=True, \ 
     546                              all=False, raw=False, **kwds): 
     547    """check for shortness with respect to the scenario itself 
     548 
     549Inputs: 
     550    L -- the lipschitz constant 
     551    blamelist -- if True, report which points are infeasible 
     552    pairs -- if True, report indicies of infeasible points 
     553    all -- if True, report results for each point (opposed to all points) 
     554    raw -- if True, report numerical results (opposed to boolean results) 
     555 
     556Additional Inputs: 
     557    tol -- maximum acceptable deviation from shortness 
     558    cutoff -- zero out distances less than cutoff; typically: tol, 0.0, or None 
     559 
     560Notes: 
     561    Each point x,y can be thought to have an associated double-cone with slope 
     562    equal to the lipschitz constant. Shortness with respect to another point is 
     563    defined by the first point not being inside the cone of the second. We can 
     564    allow for some error in shortness, a short tolerance 'tol', for which the 
     565    point x,y is some acceptable y-distance inside the cone. While very tightly 
     566    related, cutoff and tol play distinct roles; tol is subtracted from 
     567    calculation of the lipschitz_distance, while cutoff zeros out the value 
     568    of any element less than the cutoff. 
     569""" 
     570    from mystic.math.legacydata import dataset  
     571    data = dataset()  
     572    data.load(self.coords, self.values) 
     573    data.lipschitz = L 
     574    for i in range(len(data)): 
     575      data[i].id = i 
     576    return data.short(blamelist=blamelist, pairs=pairs, \ 
     577                                           all=all, raw=raw, **kwds) 
     578 
     579  def short_wrt_data(self, data, L=None, blamelist=False, pairs=True, \ 
     580                                         all=False, raw=False, **kwds): 
     581    """check for shortness with respect to the given data 
     582 
     583Inputs: 
     584    data -- a collection of data points 
     585    L -- the lipschitz constant, if different from that provided with data 
     586    blamelist -- if True, report which points are infeasible 
     587    pairs -- if True, report indicies of infeasible points 
     588    all -- if True, report results for each point (opposed to all points) 
     589    raw -- if True, report numerical results (opposed to boolean results) 
     590 
     591Additional Inputs: 
     592    tol -- maximum acceptable deviation from shortness 
     593    cutoff -- zero out distances less than cutoff; typically cutoff = tol or 0.0 
     594 
     595Notes: 
     596    Each point x,y can be thought to have an associated double-cone with slope 
     597    equal to the lipschitz constant. Shortness with respect to another point is 
     598    defined by the first point not being inside the cone of the second. We can 
     599    allow for some error in shortness, a short tolerance 'tol', for which the 
     600    point x,y is some acceptable y-distance inside the cone. While very tightly 
     601    related, cutoff and tol play distinct roles; tol is subtracted from 
     602    calculation of the lipschitz_distance, while cutoff zeros out the value 
     603    of any element less than the cutoff. 
     604""" 
     605    from mystic.math.legacydata import dataset  
     606    _self = dataset()  
     607    _self.load(self.coords, self.values) 
     608    _self.lipschitz = data.lipschitz 
     609    for i in range(len(_self)): 
     610      _self[i].id = i 
     611    return _self.short(data, L=L, blamelist=blamelist, pairs=pairs, \ 
     612                                               all=all, raw=raw, **kwds) 
     613 
     614  def set_feasible(self, data, cutoff=0.0, bounds=None, constraints=None, \ 
     615                                                  with_self=True, **kwds): 
     616    """impose shortness on a scenario with respect to given data points 
     617 
     618Inputs: 
     619    data -- a collection of data points 
     620    cutoff -- acceptable deviation from shortness 
     621 
     622Additional Inputs: 
     623    with_self -- if True, shortness will also be imposed with respect to self 
     624    tol -- acceptable optimizer termination before sum(infeasibility) = 0. 
     625    bounds -- a tuple of sample bounds:   bounds = (lower_bounds, upper_bounds) 
     626    constraints -- a function that takes a flat list parameters 
     627        x' = constraints(x) 
     628""" 
     629    # imposes: is_short(x, x'), is_short(x, z ) 
     630    # use additional 'constraints' kwds to impose: y >= m, norm(wi) = 1.0 
     631    pm = impose_feasible(cutoff, data, guess=self.pts, bounds=bounds, \ 
     632                         constraints=constraints, with_self=with_self, **kwds) 
     633    self.update( pm.flatten(all=True) ) 
     634    return 
     635 
     636  def set_valid(self, model, cutoff=0.0, bounds=None, constraints=None, **kwds): 
     637    """impose validity on a scenario with respect to given data points 
     638 
     639Inputs: 
     640    model -- the model function, y' = F(x'), that approximates reality, y = G(x) 
     641    cutoff -- acceptable model invalidity |y - F(x')| 
     642 
     643Additional Inputs: 
     644    xtol -- acceptable pointwise graphical distance of model from reality 
     645    tol -- acceptable optimizer termination before sum(infeasibility) = 0. 
     646    bounds -- a tuple of sample bounds:   bounds = (lower_bounds, upper_bounds) 
     647    constraints -- a function that takes a flat list parameters 
     648        x' = constraints(x) 
     649 
     650Notes: 
     651    xtol defines the n-dimensional base of a pilar of height cutoff, centered at 
     652    each point. The region inside the pilar defines the space where a "valid" 
     653    model must intersect. If xtol is not specified, then the base of the pilar 
     654    will be a dirac at x' = x. This function performs an optimization to find 
     655    a set of points where the model is valid. Here, tol is used to set the 
     656    optimization termination for the sum(graphical_distances), while cutoff is 
     657    used in defining the graphical_distance between x,y and x',F(x'). 
     658""" 
     659    # imposes is_feasible(R, Cv), where R = graphical_distance(model, pts) 
     660    # use additional 'constraints' kwds to impose: y >= m, norm(wi) = 1.0 
     661    pm = impose_valid(cutoff, model, guess=self, \ 
     662                      bounds=bounds, constraints=constraints, **kwds) 
     663    self.update( pm.flatten(all=True) ) 
     664    return 
     665  
     666  def pof_value(self, f): 
     667    """calculate probability of failure over a given function, f, 
     668where f takes a list of (scenario) values and returns a single value 
     669 
     670Inputs: 
     671    f -- a function that returns True for 'success' and False for 'failure' 
     672""" 
     673    u = 0 
     674    set = zip(self.values, self.weights) 
     675    for x in set: 
     676      if f(x[0]) <= 0.0: 
     677        u += x[1] 
     678    return u 
     679 
     680  def update(self, params): #XXX: overwritten.  create standalone instead ? 
     681    """update the scenario from a list of parameters 
     682 
     683The dimensions of the scenario will not change""" 
     684    pts = self.pts 
     685    _len = 2 * sum(pts) 
     686 
     687    if len(params)  >  _len:  # if Y-values are appended to params 
     688      params, values  =  params[:_len], params[_len:] 
     689      self.values = values[:len(self.values)] + self.values[len(values):]  
     690 
     691    pm = unflatten(params, pts) 
     692    zo = pm.count([]) 
     693    self[:] = pm[:len(self) - zo] + self[len(pm) - zo:] 
     694    return 
     695 
     696  def load(self, params, pts): #XXX: overwritten.  create standalone instead ? 
     697    """load a list of parameters corresponding to N x 1D discrete measures 
     698 
     699Inputs: 
     700    params -- a list of parameters (see 'notes') 
     701    pts -- number of points in each of the underlying discrete measures 
     702 
     703Notes: 
     704    To append len(pts) new discrete measures to scenario c, where 
     705    pts = (M, N, ...) 
     706    params = [wt_x1, ..., wt_xM, \ 
     707                 x1, ..., xM,    \ 
     708              wt_y1, ..., wt_yN, \ 
     709                 y1, ..., yN,    \ 
     710                     ...] 
     711    Thus, the provided list is M weights and the corresponding M positions, 
     712    followed by N weights and the corresponding N positions, with this 
     713    pattern followed for each new dimension desired for the scenario. 
     714""" 
     715    _len = 2 * sum(pts) 
     716    if len(params)  >  _len:  # if Y-values are appended to params 
     717      params, self.values  =  params[:_len], params[_len:] 
     718 
     719    self.extend( unflatten(params, pts) ) 
     720    return 
     721 
     722  def flatten(self, all=True): #XXX: overwritten.  create standalone instead ? 
     723    """flatten the scenario into a list of parameters 
     724 
     725Returns: 
     726    params -- a list of parameters (see 'notes') 
     727 
     728Notes: 
     729    For a scenario c where c.pts = (M, N, ...), then 
     730    params = [wt_x1, ..., wt_xM, \ 
     731                 x1, ..., xM,    \ 
     732              wt_y1, ..., wt_yN, \ 
     733                 y1, ..., yN,    \ 
     734                     ...] 
     735    Thus, the returned list is M weights and the corresponding M positions, 
     736    followed by N weights and the corresponding N positions, with this 
     737    pattern followed for each dimension of the scenario. 
     738""" 
     739    params = flatten(self) 
     740    if all: params.extend(self.values) # if Y-values, return those as well 
     741    return params 
     742 
     743  # interface 
     744  values = property(__values, __set_values ) 
     745  pass 
     746 
     747 
    370748#--------------------------------------------- 
    371749# creators and destructors from parameter list 
    372750 
    373751def _mimic(samples, weights): 
    374   """Generate a product_measure object from a list N product measure 
     752  """Generate a product_measure object from a list of N product measure 
    375753positions and a list of N weights. The resulting product measure will 
    376754mimic the original product measure's statistics, but be larger in size. 
     
    395773 
    396774 
    397 def _list_of_measures(samples, weights): 
     775def _uniform_weights(samples): 
     776  """generate a nested list of N x 1D weights from a nested list of N x 1D 
     777discrete measure positions, where the weights have norm 1.0 and are uniform. 
     778 
     779>>> c.pos 
     780[[1, 2, 3], [4, 5], [6]] 
     781>>> _uniform_weights(c.pos) 
     782[[0.333333333333333, 0.333333333333333, 0.333333333333333], [0.5, 0.5], [1.0]] 
     783""" 
     784  from mystic.math.measures import normalize 
     785  return [normalize([1.]*len(xi)) for xi in samples] 
     786 
     787 
     788def _list_of_measures(samples, weights=None): 
    398789  """generate a list of N x 1D discrete measures from a nested list of N x 1D 
    399 discrete measure positions and a nested list of N x 1D weights.""" 
     790discrete measure positions and a nested list of N x 1D weights. 
     791 
     792Note this function does not return a product measure, it returns a list.""" 
    400793  total = [] 
     794  if not weights: weights = _uniform_weights(samples) 
    401795  for i in range(len(samples)): 
    402796    next = dirac_measure() 
     
    407801 
    408802 
    409 def compose(samples, weights): 
     803def compose(samples, weights=None): 
    410804  """Generate a product_measure object from a nested list of N x 1D 
    411 discrete measure positions and a nested list of N x 1D weights.""" 
     805discrete measure positions and a nested list of N x 1D weights. If weights 
     806are not provided, a uniform distribution with norm = 1.0 will be used.""" 
     807  if not weights: weights = _uniform_weights(samples) 
    412808  total = _list_of_measures(samples, weights) 
    413809  c = product_measure(total) 
     
    423819 
    424820 
     821#def expand(data, npts): 
     822#  """Generate a scenario object from a dataset. The scenario will have 
     823#uniformly distributed weights and have dimensions given by pts.""" 
     824#  coords,values = data.fetch() 
     825#  from mystic.math.measures import _unpack 
     826#  pm = compose( _unpack(coords, npts) ) 
     827#  return scenario(pm, values[:pm.npts]) 
     828 
     829 
    425830def unflatten(params, npts): 
    426831  """Map a list of random variables to N x 1D discrete measures 
     
    431836 
    432837 
     838from itertools import chain #XXX: faster, but sloppy to have as importable 
    433839def flatten(c): 
    434840  """Flattens a product_measure object into a list.""" 
    435   rv = [] 
    436   for i in range(len(c)): 
    437     rv.append(c[i].weights) 
    438     rv.append(c[i].coords) 
     841  rv = [(i.weights,i.coords) for i in c] 
    439842  # now flatten list of lists into just a list 
    440   from itertools import chain 
    441   rv = list(chain(*rv)) 
    442   return rv 
     843  return list(chain(*chain(*rv))) # faster than mystic.tools.flatten 
     844 
     845 
     846##### bounds-conserving-mean: borrowed from seismic/seismic.py ##### 
     847def bounded_mean(mean_x, samples, xmin, xmax, wts=None): 
     848  from mystic.math.measures import impose_mean, impose_spread 
     849  from mystic.math.measures import spread, mean 
     850  from numpy import asarray 
     851  a = impose_mean(mean_x, samples, wts) 
     852  if min(a) < xmin:   # maintain the bound 
     853    #print "violate lo(a)" 
     854    s = spread(a) - 2*(xmin - min(a)) #XXX: needs compensation (as below) ? 
     855    a = impose_mean(mean_x, impose_spread(s, samples, wts), wts) 
     856  if max(a) > xmax:   # maintain the bound 
     857    #print "violate hi(a)" 
     858    s = spread(a) + 2*(xmax - max(a)) #XXX: needs compensation (as below) ? 
     859    a = impose_mean(mean_x, impose_spread(s, samples, wts), wts) 
     860  return asarray(a) 
     861##################################################################### 
     862 
     863 
     864#-------------------------------------------------- 
     865# constraints solvers and factories for feasibility 
     866 
     867# used in self-consistent constraints function c(x) for 
     868#   is_short(x, x') and is_short(x, z) 
     869def norm_wts_constraintsFactory(pts): 
     870  """factory for a constraints function that: 
     871  - normalizes weights 
     872""" 
     873 #from dirac_measure import scenario 
     874  def constrain(rv): 
     875    "constrain:  sum(wi)_{k} = 1 for each k in K" 
     876    pm = scenario() 
     877    pm.load(rv, pts)      # here rv is param: w,x,y 
     878    #impose: sum(wi)_{k} = 1 for each k in K 
     879    norm = 1.0 
     880    for i in range(len(pm)): 
     881      w = pm[i].weights 
     882      w[-1] = norm - sum(w[:-1]) 
     883      pm[i].weights = w 
     884    rv = pm.flatten(all=True) 
     885    return rv 
     886  return constrain 
     887 
     888# used in self-consistent constraints function c(x) for 
     889#   is_short(x, x'), is_short(x, z), and y >= m 
     890def mean_y_norm_wts_constraintsFactory(target, pts): 
     891  """factory for a constraints function that: 
     892  - imposes a mean on scenario values 
     893  - normalizes weights 
     894""" 
     895 #from dirac_measure import scenario 
     896  from mystic.math.measures import mean, impose_mean 
     897 #target[0] is target mean 
     898 #target[1] is acceptable deviation 
     899  def constrain(rv): 
     900    "constrain:  y >= m  and  sum(wi)_{k} = 1 for each k in K" 
     901    pm = scenario() 
     902    pm.load(rv, pts)      # here rv is param: w,x,y 
     903    #impose: sum(wi)_{k} = 1 for each k in K 
     904    norm = 1.0 
     905    for i in range(len(pm)): 
     906      w = pm[i].weights 
     907      w[-1] = norm - sum(w[:-1]) 
     908      pm[i].weights = w 
     909    #impose: y >= m  
     910    values, weights = pm.values, pm.weights 
     911    y = float(mean(values, weights)) 
     912    if not (y >= float(target[0])): 
     913      pm.values = impose_mean(target[0]+target[1], values, weights) 
     914    rv = pm.flatten(all=True)  
     915    return rv 
     916  return constrain 
     917 
     918def impose_feasible(cutoff, data, guess=None, **kwds): 
     919  """impose shortness on a given list of parameters w,x,y. 
     920 
     921Optimization on w,x,y over the given bounds seeks sum(infeasibility) = 0. 
     922  (this function is not ???-preserving) 
     923 
     924Inputs: 
     925    cutoff -- maximum acceptable deviation from shortness 
     926    data -- a dataset of observed points (these points are 'static') 
     927    guess -- the scenario providing an initial guess at feasibility, 
     928        or a tuple of dimensions of the target scenario 
     929 
     930Additional Inputs: 
     931    tol -- acceptable optimizer termination before sum(infeasibility) = 0. 
     932    bounds -- a tuple of sample bounds:   bounds = (lower_bounds, upper_bounds) 
     933    constraints -- a function that takes a flat list parameters 
     934        x' = constraints(x) 
     935 
     936Outputs: 
     937    pm -- a scenario with desired shortness 
     938""" 
     939  from numpy import sum, asarray 
     940  from mystic.math.legacydata import dataset 
     941  from mystic.math.paramtrans import lipschitz_distance, infeasibility, _npts 
     942  if guess is None: 
     943    message = "Requires a guess scenario, or a tuple of scenario dimensions." 
     944    raise TypeError, message 
     945  # get initial guess 
     946  if hasattr(guess, 'pts'): # guess is a scenario 
     947    pts = guess.pts    # number of x 
     948    guess = guess.flatten(all=True) 
     949  else: 
     950    pts = guess        # guess is given as a tuple of 'pts' 
     951    guess = None 
     952  npts = _npts(pts)    # number of Y 
     953  long_form = len(pts) - list(pts).count(2) # can use '2^K compressed format' 
     954 
     955  # prepare bounds for solver 
     956  bounds = kwds.pop('bounds', None) 
     957  # if bounds are not set, use the default optimizer bounds 
     958  if bounds is None: 
     959    lower_bounds = []; upper_bounds = [] 
     960    for n in pts:  # bounds for n*x in each dimension  (x2 due to weights) 
     961      lower_bounds += [None]*n * 2 
     962      upper_bounds += [None]*n * 2 
     963    # also need bounds for npts*y values 
     964    lower_bounds += [None]*npts 
     965    upper_bounds += [None]*npts 
     966    bounds = lower_bounds, upper_bounds 
     967  bounds = asarray(bounds).T 
     968 
     969  # plug in the 'constraints' function:  param' = constraints(param) 
     970  # constraints should impose_mean(y,w), and possibly sum(weights) 
     971  constraints = kwds.pop('constraints', None) # default is no constraints 
     972  if not constraints:  # if None (default), there are no constraints 
     973    constraints = lambda x: x 
     974 
     975  _self = kwds.pop('with_self', True) # default includes self in shortness 
     976  if _self is not False: _self = True 
     977  # tolerance for optimization on sum(y) 
     978  tol = kwds.pop('tol', 0.0) # default 
     979  npop = kwds.pop('npop', 20) #XXX: tune npop? 
     980  maxiter = kwds.pop('maxiter', 1000) #XXX: tune maxiter? 
     981 
     982  # if no guess was made, then use bounds constraints 
     983  if guess is None: 
     984    if npop: 
     985      guess = bounds 
     986    else:  # fmin_powell needs a list params (not bounds) 
     987      guess = [(a + b)/2. for (a,b) in bounds] 
     988 
     989  # construct cost function to reduce sum(lipschitz_distance) 
     990  def cost(rv): 
     991    """compute cost from a 1-d array of model parameters, 
     992    where:  cost = | sum(lipschitz_distance) | """ 
     993    _data = dataset() 
     994    _pm = scenario() 
     995    _pm.load(rv, pts)      # here rv is param: w,x,y 
     996    if not long_form: 
     997      coords = _pm.select(*range(npts)) 
     998    else: coords = _pm.coords 
     999    _data.load( data.coords, data.values )                        # LOAD static 
     1000    if _self: 
     1001      _data.load( coords, _pm.values )                            # LOAD dynamic 
     1002    _data.lipschitz = data.lipschitz                              # LOAD L 
     1003    Rv = lipschitz_distance(_data.lipschitz, _pm, _data, tol=cutoff, **kwds) 
     1004    v = infeasibility(Rv, cutoff) 
     1005    return abs(sum(v)) 
     1006 
     1007  # construct and configure optimizer 
     1008  debug = False  #!!! 
     1009  maxfun = 1e+6 
     1010  crossover = 0.9; percent_change = 0.9 
     1011  ftol = abs(tol); gtol = None 
     1012 
     1013  if debug: 
     1014    print "lower bounds: %s" % bounds.T[0] 
     1015    print "upper bounds: %s" % bounds.T[1] 
     1016  # print "initial value: %s" % guess 
     1017  # use optimization to get feasible points 
     1018  from mystic.solvers import diffev2, fmin_powell 
     1019  from mystic.monitors import Monitor, VerboseMonitor 
     1020  from mystic.strategy import Best1Bin, Best1Exp 
     1021  evalmon = Monitor();  stepmon = Monitor(); strategy = Best1Exp 
     1022  if debug: stepmon = VerboseMonitor(10)  #!!! 
     1023  if npop: # use VTR 
     1024    results = diffev2(cost, guess, npop, ftol=ftol, gtol=gtol, bounds=bounds,\ 
     1025                      maxiter=maxiter, maxfun=maxfun, constraints=constraints,\ 
     1026                      cross=crossover, scale=percent_change, strategy=strategy,\ 
     1027                      evalmon=evalmon, itermon=stepmon,\ 
     1028                      full_output=1, disp=0, handler=False) 
     1029  else: # use VTR 
     1030    results = fmin_powell(cost, guess, ftol=ftol, gtol=gtol, bounds=bounds,\ 
     1031                      maxiter=maxiter, maxfun=maxfun, constraints=constraints,\ 
     1032                      evalmon=evalmon, itermon=stepmon,\ 
     1033                      full_output=1, disp=0, handler=False) 
     1034  # repack the results 
     1035  pm = scenario() 
     1036  pm.load(results[0], pts)            # params: w,x,y 
     1037 #if debug: print "final cost: %s" % results[1] 
     1038  if debug and results[2] >= maxiter: # iterations 
     1039    print "Warning: constraints solver terminated at maximum iterations" 
     1040 #func_evals = results[3]           # evaluation 
     1041  return pm 
     1042 
     1043 
     1044def impose_valid(cutoff, model, guess=None, **kwds): 
     1045  """impose model validity on a given list of parameters w,x,y 
     1046 
     1047Optimization on w,x,y over the given bounds seeks sum(infeasibility) = 0. 
     1048  (this function is not ???-preserving) 
     1049 
     1050Inputs: 
     1051    cutoff -- maximum acceptable model invalidity |y - F(x')|; a single value 
     1052    model -- the model function, y' = F(x'), that approximates reality, y = G(x) 
     1053    guess -- the scenario providing an initial guess at validity, 
     1054        or a tuple of dimensions of the target scenario 
     1055 
     1056Additional Inputs: 
     1057    xtol -- acceptable pointwise graphical distance of model from reality 
     1058    tol -- acceptable optimizer termination before sum(infeasibility) = 0. 
     1059    bounds -- a tuple of sample bounds:   bounds = (lower_bounds, upper_bounds) 
     1060    constraints -- a function that takes a flat list parameters 
     1061        x' = constraints(x) 
     1062 
     1063Outputs: 
     1064    pm -- a scenario with desired model validity 
     1065 
     1066Notes: 
     1067    xtol defines the n-dimensional base of a pilar of height cutoff, centered at 
     1068    each point. The region inside the pilar defines the space where a "valid" 
     1069    model must intersect. If xtol is not specified, then the base of the pilar 
     1070    will be a dirac at x' = x. This function performs an optimization to find 
     1071    a set of points where the model is valid. Here, tol is used to set the 
     1072    optimization termination for the sum(graphical_distances), while cutoff is 
     1073    used in defining the graphical_distance between x,y and x',F(x'). 
     1074""" 
     1075  from numpy import sum as _sum, asarray 
     1076  from mystic.math.paramtrans import graphical_distance, infeasibility, _npts 
     1077  if guess is None: 
     1078    message = "Requires a guess scenario, or a tuple of scenario dimensions." 
     1079    raise TypeError, message 
     1080  # get initial guess 
     1081  if hasattr(guess, 'pts'): # guess is a scenario 
     1082    pts = guess.pts    # number of x 
     1083    guess = guess.flatten(all=True) 
     1084  else: 
     1085    pts = guess        # guess is given as a tuple of 'pts' 
     1086    guess = None 
     1087  npts = _npts(pts)    # number of Y 
     1088 
     1089  # prepare bounds for solver 
     1090  bounds = kwds.pop('bounds', None) 
     1091  # if bounds are not set, use the default optimizer bounds 
     1092  if bounds is None: 
     1093    lower_bounds = []; upper_bounds = [] 
     1094    for n in pts:  # bounds for n*x in each dimension  (x2 due to weights) 
     1095      lower_bounds += [None]*n * 2 
     1096      upper_bounds += [None]*n * 2 
     1097    # also need bounds for npts*y values 
     1098    lower_bounds += [None]*npts 
     1099    upper_bounds += [None]*npts 
     1100    bounds = lower_bounds, upper_bounds 
     1101  bounds = asarray(bounds).T 
     1102 
     1103  # plug in the 'constraints' function:  param' = constraints(param) 
     1104  constraints = kwds.pop('constraints', None) # default is no constraints 
     1105  if not constraints:  # if None (default), there are no constraints 
     1106    constraints = lambda x: x 
     1107 
     1108  # 'wiggle room' tolerances 
     1109  ipop = kwds.pop('ipop', 10) #XXX: tune ipop (inner optimization)? 
     1110  imax = kwds.pop('imax', 10) #XXX: tune imax (inner optimization)? 
     1111  # tolerance for optimization on sum(y) 
     1112  tol = kwds.pop('tol', 0.0) # default 
     1113  npop = kwds.pop('npop', 20) #XXX: tune npop (outer optimization)? 
     1114  maxiter = kwds.pop('maxiter', 1000) #XXX: tune maxiter (outer optimization)? 
     1115 
     1116  # if no guess was made, then use bounds constraints 
     1117  if guess is None: 
     1118    if npop: 
     1119      guess = bounds 
     1120    else:  # fmin_powell needs a list params (not bounds) 
     1121      guess = [(a + b)/2. for (a,b) in bounds] 
     1122 
     1123  # construct cost function to reduce sum(infeasibility) 
     1124  def cost(rv): 
     1125    """compute cost from a 1-d array of model parameters, 
     1126    where: cost = | sum( infeasibility ) | """ 
     1127    # converting rv to scenario 
     1128    points = scenario() 
     1129    points.load(rv, pts) 
     1130    # calculate infeasibility 
     1131    Rv = graphical_distance(model, points, ytol=cutoff, ipop=ipop, \ 
     1132                                                        imax=imax, **kwds) 
     1133    v = infeasibility(Rv, cutoff) 
     1134    # converting v to E 
     1135    return _sum(v) #XXX: abs ? 
     1136 
     1137  # construct and configure optimizer 
     1138  debug = False  #!!! 
     1139  maxfun = 1e+6 
     1140  crossover = 0.9; percent_change = 0.8 
     1141  ftol = abs(tol); gtol = None #XXX: optimally, should be VTRCOG... 
     1142 
     1143  if debug: 
     1144    print "lower bounds: %s" % bounds.T[0] 
     1145    print "upper bounds: %s" % bounds.T[1] 
     1146  # print "initial value: %s" % guess 
     1147  # use optimization to get model-valid points 
     1148  from mystic.solvers import diffev2, fmin_powell 
     1149  from mystic.monitors import Monitor, VerboseMonitor 
     1150  from mystic.strategy import Best1Bin, Best1Exp 
     1151  evalmon = Monitor();  stepmon = Monitor(); strategy = Best1Exp 
     1152  if debug: stepmon = VerboseMonitor(2)  #!!! 
     1153  if npop: # use VTR 
     1154    results = diffev2(cost, guess, npop, ftol=ftol, gtol=gtol, bounds=bounds,\ 
     1155                      maxiter=maxiter, maxfun=maxfun, constraints=constraints,\ 
     1156                      cross=crossover, scale=percent_change, strategy=strategy,\ 
     1157                      evalmon=evalmon, itermon=stepmon,\ 
     1158                      full_output=1, disp=0, handler=False) 
     1159  else: # use VTR 
     1160    results = fmin_powell(cost, guess, ftol=ftol, gtol=gtol, bounds=bounds,\ 
     1161                      maxiter=maxiter, maxfun=maxfun, constraints=constraints,\ 
     1162                      evalmon=evalmon, itermon=stepmon,\ 
     1163                      full_output=1, disp=0, handler=False) 
     1164  # repack the results 
     1165  pm = scenario() 
     1166  pm.load(results[0], pts)            # params: w,x,y 
     1167 #if debug: print "final cost: %s" % results[1] 
     1168  if debug and results[2] >= maxiter: # iterations 
     1169    print "Warning: constraints solver terminated at maximum iterations" 
     1170 #func_evals = results[3]           # evaluation 
     1171  return pm 
     1172 
     1173 
     1174if __name__ == '__main__': 
     1175  from mystic.math.paramtrans import * 
     1176  model = lambda x:sum(x) 
     1177  a = [0,1,9,8, 1,0,4,6, 1,0,1,2, 0,1,2,3,4,5,6,7] 
     1178  feasability = 0.0; deviation = 0.01 
     1179  validity = 5.0; wiggle = 1.0 
     1180  y_mean = 5.0; y_buffer = 0.0 
     1181  L = [.75,.5,.25] 
     1182  bc = [(0,7,2),(3,0,2),(2,0,3),(1,0,3),(2,4,2)] 
     1183  bv = [5,3,1,4,8] 
     1184  pts = (2,2,2) 
     1185  from mystic.math.legacydata import dataset 
     1186  data = dataset() 
     1187  data.load(bc, bv) 
     1188  data.lipschitz = L 
     1189  pm = scenario() 
     1190  pm.load(a, pts) 
     1191  pc = pm.coords 
     1192  pv = pm.values 
     1193  #--- 
     1194  _data = dataset() 
     1195  _data.load(bc, bv) 
     1196  _data.load(pc, pv) 
     1197  _data.lipschitz = data.lipschitz 
     1198  from numpy import sum 
     1199  ans = sum(lipschitz_distance(L, pm, _data)) 
     1200  print "original: %s @ %s\n" % (ans, a) 
     1201 #print "pm: %s" % pm 
     1202 #print "data: %s" % data 
     1203  #--- 
     1204  lb = [0,.5,-100,-100,  0,.5,-100,-100,  0,.5,-100,-100,   0,0,0,0,0,0,0,0] 
     1205  ub = [.5,1, 100, 100,  .5,1, 100, 100,  .5,1, 100, 100,   9,9,9,9,9,9,9,9] 
     1206  bounds = (lb,ub) 
     1207 
     1208  _constrain = mean_y_norm_wts_constraintsFactory((y_mean,y_buffer), pts) 
     1209  results = impose_feasible(feasability, data, guess=pts, tol=deviation, \ 
     1210                            bounds=bounds, constraints=_constrain) 
     1211  from mystic.math.measures import mean 
     1212  print "solved: %s" % results.flatten(all=True) 
     1213  print "mean(y): %s >= %s" % (mean(results.values, results.weights), y_mean) 
     1214  print "sum(wi): %s == 1.0" % [sum(w) for w in results.wts] 
     1215 
     1216  print "\n---------------------------------------------------\n" 
     1217 
     1218  bc = bc[:-2] 
     1219  ids = ['1','2','3'] 
     1220  t = dataset() 
     1221  t.load(bc, map(model, bc), ids) 
     1222  t.update(t.coords, map(model, t.coords)) 
     1223# r = dataset() 
     1224# r.load(t.coords, t.values) 
     1225# L = [0.1, 0.0, 0.0] 
     1226  print "%s" % t 
     1227  print "L: %s" % L 
     1228  print "shortness:" 
     1229  print lipschitz_distance(L, t, t, tol=0.0) 
     1230   
     1231  print "\n---------------------------------------------------\n" 
     1232 
     1233  print "Y: %s" % str(results.values) 
     1234  print "sum(wi): %s == 1.0" % [sum(w) for w in results.wts] 
    4431235 
    4441236 
  • mystic/_math/legacydata.py

    r558 r569  
    328328    if L is None: L = self.lipschitz 
    329329    if data is None: data = self 
    330     from paramtrans import lipschitz_distance, is_feasible 
     330    from mystic.math.paramtrans import lipschitz_distance, is_feasible 
    331331    # calculate the shortness 
    332332    Rv = lipschitz_distance(L, self, data, **kwds) 
     
    372372    elif cutoff is False: cutoff = None 
    373373 
    374     from paramtrans import graphical_distance, is_feasible 
     374    from mystic.math.paramtrans import graphical_distance, is_feasible 
    375375    # calculate the model validity 
    376376    Rv = graphical_distance(model, self, **kwds) 
  • mystic/_math/paramtrans.py

    r558 r569  
    305305 #L = list of lipschitz constants, for use when lipschitz metric is desired 
    306306 #constraints = constraints function for finding minimum distance 
    307   from legacydata import dataset 
     307  from mystic.math.legacydata import dataset 
    308308  from numpy import asarray 
    309309  from mystic.solvers import diffev2, fmin_powell 
     
    422422  print "\nbuilding a scenario from the params..." 
    423423  # [store Y as 'values' OR register(F) for Y=F(X) OR points store y as 'val' ?] 
    424   from dirac_measure import scenario 
     424  from mystic.math.dirac_measure import scenario 
    425425  pm = scenario() 
    426426  pm.load(param1, pts) 
     
    437437  # build a dataset (using X,Y) 
    438438  # [store W as 'weights' ?] 
    439   from legacydata import dataset 
     439  from mystic.math.legacydata import dataset 
    440440  d = dataset() 
    441441  d.load(X, Y)  
  • mystic/scripts/support_hypercube_scenario.py

    r566 r569  
    266266    # would be nice to use meta = ['wx','wx2','x','x2','wy',...] 
    267267 
    268   from dirac_measure import scenario 
    269   from legacydata import dataset 
     268  from mystic.math.dirac_measure import scenario 
     269  from mystic.math.legacydata import dataset 
    270270  try: # select whether to plot the cones 
    271271    cones = parsed_opts.cones 
     
    280280  try:  # get the name of the dataset file 
    281281    file = parsed_args[1] 
    282     from legacydata import load_dataset 
     282    from mystic.math.legacydata import load_dataset 
    283283    data = load_dataset(file) 
    284284  except: 
  • mystic/setup.py

    r530 r569  
    316316             'scripts/support_convergence.py', 
    317317             'scripts/support_hypercube.py', 
    318              'scripts/support_hypercube_measures.py']) 
     318             'scripts/support_hypercube_measures.py', 
     319             'scripts/support_hypercube_scenario.py']) 
    319320""" 
    320321 
Note: See TracChangeset for help on using the changeset viewer.