Changeset 556 for branches


Ignore:
Timestamp:
09/10/12 07:35:06 (4 years ago)
Author:
mmckerns
Message:

fix: graphical_distance optimization is short-circuited when Cx = 0;
tuned impose_valid for Cx = 0 and Cx != 0
exposed npop and maxiter to user in dataset.valid, scenario.set_valid, etc
fix: account for the non-standard return of direc from fmin_powell
add code for StStSurrogate? tested for validity with Cy (no Cx)

Location:
branches/UQ/math/legacy
Files:
1 added
1 deleted
6 edited
2 copied
4 moved

Legend:

Unmodified
Added
Removed

  • branches/UQ/math/legacy/MINMAX_StAlData.py

    r546 r556  
    44MINMAX = -1  ## NOTE: sup = maximize = -1; inf = minimize = 1 
    55####################################################################### 
    6 # scaling and mpi info; also optimizer configuration parameters 
    7 # hard-wired: use DE solver, don't use mpi, F-F' calculation 
     6# optimizer configuration parameters 
    87####################################################################### 
    98npop = 20 #40  #!!! 
    … …  
    1615 
    1716####################################################################### 
    18 # the model function and the dataset 
    19 ####################################################################### 
    20 #from StAlSurrogate import st_al_surr as model 
     17# the dataset 
     18####################################################################### 
    2119from legacydata import load_dataset 
    2220data = load_dataset('StAlDataset.txt') 
    … …  
    150148####################################################################### 
    151149if __name__ == '__main__': 
    152  #function_name = model.__name__ 
    153150  dataset_name = data.__name__ 
    154151 
    … …  
    164161  _n = _npts(npts) 
    165162 
    166   #FIXME: check units versus those in model and in dataset 
     163  #XXX: check units versus those in dataset 
    167164  w_lower = [0.0]; Y_lower = [0.0] 
    168165  w_upper = [1.0]; Y_upper = [100.0] 
    … …  
    189186  print "..............\n" 
    190187 
    191  #print " model: f(x) = %s(x)" % function_name 
    192188  print " data: z = {%s}" % dataset_name 
    193189  print " target: %s" % str(target) 

  • branches/UQ/math/legacy/MM_OUQ_StAlData.py

    r546 r556  
    1010MINMAX = -1  ## NOTE: sup = maximize = -1; inf = minimize = 1 
    1111####################################################################### 
    12 # scaling and mpi info; also optimizer configuration parameters 
    13 # hard-wired: use DE solver, don't use mpi, F-F' calculation 
     12# optimizer configuration parameters 
    1413####################################################################### 
    1514npop = 32  #!!! 
    … …  
    2221 
    2322####################################################################### 
    24 # the model function and the dataset 
    25 ####################################################################### 
    26 #from StAlSurrogate import st_al_surr as model 
     23# the dataset 
     24####################################################################### 
    2725from legacydata import load_dataset 
    2826data = load_dataset('StAlDataset.txt') 
    … …  
    111109    # check mean value, and if necessary use constrain to set mean value 
    112110    y = float(mean(c.values, c.weights)) 
    113     if not (y >= float(target[0])):   #XXX: or target[0]-error[0] ? 
    114      #c.values = impose_mean(target[0]+target[1], c.values, c.weights) 
     111    if not (y >= float(target[0] - error[0])): 
    115112      c.update( constrain( c.flatten(all=True) ) ) 
    116113 
    … …  
    127124      if not [sum(w) for w in c.wts] == [1.0] * len(c.wts): 
    128125        print "norm(w) == 1.0: False" 
    129       if not c.get_mean_value() >= target[0]: 
    130         print "mean(y) >= %s: False"#%s" % (target[0], c.get_mean_value() \ 
    131         #                                                   >= target[0]) 
     126      if not c.get_mean_value() >= float(target[0] - error[0]): 
     127        print "mean(y) >= %s: False" % str(target[0] - error[0]) 
    132128    ###################### end function-specific ###################### 
    133129    # extract weights and positions and values 
    … …  
    146142    ##################### begin function-specific ##################### 
    147143    y = float(mean(c.values, c.weights)) 
    148     if not (y >= float(target[0])): 
     144    if not (y >= float(target[0] - error[0])): 
    149145      if debug: print "skipping mean: %s" % y 
    150146      return inf  #XXX: FORCE TO SATISFY E CONSTRAINTS 
    … …  
    176172####################################################################### 
    177173if __name__ == '__main__': 
    178  #function_name = model.__name__ 
    179174  dataset_name = data.__name__ 
    180  
    181  #model_tol = 2.0  #expect difference is in [-1.677, 1.507], where 2 > |-1.677| 
    182175 
    183176  from mystic.math.measures import mean 
    184177# y_mean = mean(data.values) 
    185178  y_mean = 11.0    #NOTE: SET THE 'mean' HERE! 
    186   y_error = 0.0001 #NOTE: SET THE 'error' HERE! 
     179  y_mean_error = 0.0001 #NOTE: SET THE 'error' HERE! 
    187180  theta = 7.0      #NOTE: SET THE 'failure criteria' HERE! 
    188181  short_tol = 1e-9 #NOTE: SET THE 'short tolerance' HERE! 
    189   target = (y_mean,theta,None); error = (y_error,None,short_tol) 
     182  target = (y_mean,theta,None); error = (y_mean_error,None,short_tol) 
    190183  pars = (target,error) 
    191184 
    … …  
    197190  _n = _npts(npts) 
    198191 
    199   #FIXME: check units versus those in model and in dataset 
     192  #XXX: check units versus those in dataset 
    200193  w_lower = [0.0]; Y_lower = [0.0] 
    201194  w_upper = [1.0]; Y_upper = [100.0] 
    … …  
    205198 
    206199 #XXX XXX: EDITED TO USE w_split 
    207  #lower_bounds = (nx * w_lower) + (nx * h_lower) \ 
    208  #             + (ny * w_lower) + (ny * a_lower) \ 
    209  #             + (nz * w_lower) + (nz * v_lower) \ 
    210200  lower_bounds = (w_lower) + (w_split) + (nx * h_lower) \ 
    211201               + (w_lower) + (w_split) + (ny * a_lower) \ 
    212202               + (w_lower) + (w_split) + (nz * v_lower) \ 
    213203               + (_n * Y_lower)  
    214  #upper_bounds = (nx * w_upper) + (nx * h_upper) \ 
    215  #             + (ny * w_upper) + (ny * a_upper) \ 
    216  #             + (nz * w_upper) + (nz * v_upper) \ 
    217204  upper_bounds = (w_split) + (w_upper) + (nx * h_upper) \ 
    218205               + (w_split) + (w_upper) + (ny * a_upper) \ 
    … …  
    230217  print "..............\n" 
    231218 
    232  #print " model: f(x) = %s(x)" % function_name 
    233219  print " data: z = {%s}" % dataset_name 
    234220  print " target: %s" % str(target) 
    … …  
    277263  print "fails short wrt data:\n%s" % c.short_wrt_data(data, tol=short_tol,\ 
    278264                                                             blamelist=True) 
    279   print "mean(y): %s >= %s" % (str( c.get_mean_value() ), y_mean) 
     265  print "mean(y): %s >= %s" % (str( c.get_mean_value() ), y_mean - y_mean_error) 
    280266  print "sum_wts: %s == 1.0" % [sum(w) for w in c.wts] 
    281267 

  • branches/UQ/math/legacy/TEST_OUQ_1dData.py

    r546 r556  
    88MINMAX = -1  ## NOTE: sup = maximize = -1; inf = minimize = 1 
    99####################################################################### 
    10 # scaling and mpi info; also optimizer configuration parameters 
    11 # hard-wired: use DE solver, don't use mpi, F-F' calculation 
     10# optimizer configuration parameters 
    1211####################################################################### 
    1312npop = 32  #!!! 
    … …  
    2019 
    2120####################################################################### 
    22 # the model function and the dataset 
    23 ####################################################################### 
    24 #from StAlSurrogate import st_al_surr as model 
     21# the dataset 
     22####################################################################### 
    2523from legacydata import load_dataset 
    2624data = load_dataset('ExampleDataset.txt', range(0,1)) 
    … …  
    107105    # check mean value, and if necessary use constrain to set mean value 
    108106    y = float(mean(c.values, c.weights)) 
    109     if not (y >= float(target[0])):   #XXX: or target[0]-error[0] ? 
    110      #c.values = impose_mean(target[0]+target[1], c.values, c.weights) 
     107    if not (y >= float(target[0] - error[0])): 
    111108      c.update( constrain( c.flatten(all=True) ) ) 
    112109 
    … …  
    123120      if not [sum(w) for w in c.wts] == [1.0] * len(c.wts): 
    124121        print "norm(w) == 1.0: False" 
    125       if not c.get_mean_value() >= target[0]: 
    126         print "mean(y) >= %s: False"#%s" % (target[0], c.get_mean_value() \ 
    127         #                                                   >= target[0]) 
     122      if not c.get_mean_value() >= float(target[0] - error[0]): 
     123        print "mean(y) >= %s: False" % str(target[0] - error[0]) 
    128124    ###################### end function-specific ###################### 
    129125    # extract weights and positions and values 
    … …  
    142138    ##################### begin function-specific ##################### 
    143139    y = float(mean(c.values, c.weights)) 
    144     if not (y >= float(target[0])): 
     140    if not (y >= float(target[0] - error[0])): 
    145141      if debug: print "skipping mean: %s" % y 
    146142      return inf  #XXX: FORCE TO SATISFY E CONSTRAINTS 
    … …  
    172168####################################################################### 
    173169if __name__ == '__main__': 
    174  #function_name = model.__name__ 
    175170  dataset_name = data.__name__ 
    176  
    177  #model_tol = 2.0  #expect difference is in [-1.677, 1.507], where 2 > |-1.677| 
    178171 
    179172  from mystic.math.measures import mean 
    180173# y_mean = mean(data.values) 
    181174  y_mean = 0.5     #NOTE: SET THE 'mean' HERE! 
    182   y_error = 0.0001 #NOTE: SET THE 'error' HERE! 
     175  y_mean_error = 0.0001 #NOTE: SET THE 'error' HERE! 
    183176  theta = 0.0      #NOTE: SET THE 'failure criteria' HERE! 
    184177  short_tol = 1e-9 #NOTE: SET THE 'short tolerance' HERE! 
    185   target = (y_mean,theta,None); error = (y_error,None,short_tol) 
     178  target = (y_mean,theta,None); error = (y_mean_error,None,short_tol) 
    186179  pars = (target,error) 
    187180 
    … …  
    193186  _n = _npts(npts) 
    194187 
    195   #FIXME: check units versus those in model and in dataset 
     188  #XXX: check units versus those in dataset 
    196189  w_lower = [0.0]; Y_lower = [0.0] 
    197190  w_upper = [1.0]; Y_upper = [1.0] 
    … …  
    201194 
    202195 #XXX XXX: EDITED TO USE w_split *AND* npts=(2,1,1) *AND* FIX Y1 == Y_lower 
    203  #lower_bounds = (nx * w_lower) + (nx * h_lower) \ 
    204  #             + (ny * w_lower) + (ny * a_lower) \ 
    205  #             + (nz * w_lower) + (nz * v_lower) \ 
    206196  lower_bounds = (w_lower) + (w_split) + (nx * h_lower) \ 
    207197               + (ny * w_upper) + (ny * a_lower) \ 
    208198               + (nz * w_upper) + (nz * v_lower) \ 
    209199               + (_n * Y_lower)  
    210  #upper_bounds = (nx * w_upper) + (nx * h_upper) \ 
    211200  upper_bounds = (w_split) + (w_upper) + (nx * h_upper) \ 
    212201               + (ny * w_upper) + (ny * a_upper) \ 
    213202               + (nz * w_upper) + (nz * v_upper) \ 
    214203               + (Y_lower) + (Y_upper)  
    215  #             + (_n * Y_upper)  
    216204  bounds = (lower_bounds,upper_bounds) 
    217205 
    … …  
    225213  print "..............\n" 
    226214 
    227  #print " model: f(x) = %s(x)" % function_name 
    228215  print " data: z = {%s}" % dataset_name 
    229216  print " target: %s" % str(target) 
    … …  
    272259  print "fails short wrt data:\n%s" % c.short_wrt_data(data, tol=short_tol,\ 
    273260                                                             blamelist=True) 
    274   print "mean(y): %s >= %s" % (str( c.get_mean_value() ), y_mean) 
     261  print "mean(y): %s >= %s" % (str( c.get_mean_value() ), y_mean - y_mean_error) 
    275262  print "sum_wts: %s == 1.0" % [sum(w) for w in c.wts] 
    276263 

  • branches/UQ/math/legacy/TEST_OUQ_1dSurr_CxCy.py

    r552 r556  
    22""" 
    33OUQ with model sausage: valid_wrt_model, mean(y), norm(wts) 
    4   with lipschitz and bounds from ExampleDataset 
     4  with 'wiggle room' from full sausage around the model 
    55""" 
    66 
    7 debug = True#False 
     7debug = False 
    88MINMAX = -1  ## NOTE: sup = maximize = -1; inf = minimize = 1 
    99####################################################################### 
    10 # scaling and mpi info; also optimizer configuration parameters 
    11 # hard-wired: use DE solver, don't use mpi, F-F' calculation 
    12 ####################################################################### 
    13 npop = 20 #20 #40  #!!! 
     10# optimizer configuration parameters 
     11####################################################################### 
     12npop = 10 #15 #20 #32 #40  #!!! 
    1413maxiter = 1000 
    1514maxfun = 1e+6 
    … …  
    2019 
    2120####################################################################### 
    22 # the model function and the dataset 
    23 ####################################################################### 
    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  
    29 def model(x): 
    30   return x[0]; 
     21# the model function 
     22####################################################################### 
     23def model(x): return x[0] 
     24 
    3125 
    3226####################################################################### 
    … …  
    110104    # check mean value, and if necessary use constrain to set mean value 
    111105    y = float(mean(c.values, c.weights)) 
    112     if not (y >= float(target[0])):   #XXX: or target[0]-error[0] ? 
     106    if not (y >= float(target[0] - error[0])): 
    113107      c.update( constrain( c.flatten(all=True) ) ) 
    114108 
    115109    # then test if valid... then impose model validity on product measure 
    116     if not c.valid_wrt_model(model, tol=target[2], xtol=target[3]): 
     110    if not c.valid_wrt_model(model, tol=target[2], xtol=target[3], \ 
     111                                                   maxiter=error[4]): 
    117112      c.set_valid(model, cutoff=target[2], bounds=bounds, tol=error[2], \ 
    118                                  constraints=constrain, xtol=target[3]) 
     113                         constraints=constrain, xtol=target[3], imax=target[4]) 
    119114    ###################### more function-specific ##################### 
    120115    if debug: 
    121       if not c.valid_wrt_model(model, tol=target[2], xtol=target[3]): 
     116      if not c.valid_wrt_model(model, tol=target[2], xtol=target[3], \ 
     117                                                     maxiter=error[4]): 
    122118        print "valid_wrt_model: False" 
    123119      if not [sum(w) for w in c.wts] == [1.0] * len(c.wts): 
    124120        print "norm(w) == 1.0: False" 
    125       if not c.get_mean_value() >= target[0]: 
    126         print "mean(y) >= %s: False" 
     121      if not c.get_mean_value() >= float(target[0] - error[0]): 
     122        print "mean(y) >= %s: False" % str(target[0] - error[0]) 
    127123    ###################### end function-specific ###################### 
    128124    # extract weights and positions and values 
    … …  
    141137    ##################### begin function-specific ##################### 
    142138    y = float(mean(c.values, c.weights)) 
    143     if not (y >= float(target[0])): 
     139    if not (y >= float(target[0] - error[0])): 
    144140      if debug: print "skipping mean: %s" % y 
    145141      return inf  #XXX: FORCE TO SATISFY E CONSTRAINTS 
    146142 
    147     if not c.valid_wrt_model(model, tol=target[2], xtol=target[3]): 
     143    if not c.valid_wrt_model(model, tol=target[2], xtol=target[3], \ 
     144                                                   maxiter=error[4]): 
    148145      if debug: print "skipping model-invalidity" 
    149146      return inf  #XXX: FORCE TO SATISFY E CONSTRAINTS 
    … …  
    168165if __name__ == '__main__': 
    169166  function_name = model.__name__ 
    170  #dataset_name = data.__name__ 
    171167 
    172168  from mystic.math.measures import mean 
    … …  
    174170  y_mean_error = 0.0 #NOTE: SET THE 'error' HERE! 
    175171  theta = 0.0      #NOTE: SET THE 'failure criteria' HERE! 
    176   Cy = 1.0          #NOTE: SET THE 'cutoff' HERE! 
    177   Cx = (0.0,0.0,0.0) #NOTE: SET THE 'wiggle' HERE! 
    178   valid_tol = 1.e-6  #NOTE: SET THE 'model tolerance' HERE! 
    179   target = (y_mean,theta,Cy,Cx); error = (y_mean_error,None,valid_tol,None) 
     172  Cy = 0.1           #NOTE: SET THE 'cutoff' HERE! 
     173  Cx = (2.0,0.0,0.0) #NOTE: SET THE 'wiggle' HERE! 
     174  valid_tol = 0.0  #NOTE: SET THE 'model tolerance' HERE! 
     175  iter = 200       #NOTE: SET THE 'Cy-validation maxiter' HERE! 
     176  imax = 10        #NOTE: SET THE 'Cx-validation maxiter' HERE! 
     177  target = (y_mean,theta,Cy,Cx,imax) 
     178  error = (y_mean_error,None,valid_tol,None,iter) 
    180179  pars = (target,error) 
    181180 
    … …  
    187186  _n = _npts(npts) 
    188187 
    189   #FIXME: check units versus those in model and in dataset 
     188  #XXX: check units versus those in model 
    190189  w_lower = [0.0]; Y_lower = [theta] 
    191   w_upper = [1.0]; Y_upper = [5.0] 
     190  w_upper = [1.0]; Y_upper = [10.0] 
    192191  h_lower = [0.0]; a_lower = [0.0];  v_lower = [0.0] 
    193192  h_upper = [1.0]; a_upper = [0.0];  v_upper = [0.0] 
    194   w_split = [0.5] 
    195  
    196  #XXX XXX: EDITED TO USE w_split *AND* npts=(2,1,1) *AND* FIX Y1 == Y_lower 
    197  #lower_bounds = (nx * w_lower) + (nx * h_lower) \ 
    198  #             + (ny * w_lower) + (ny * a_lower) \ 
    199  #             + (nz * w_lower) + (nz * v_lower) \ 
    200   lower_bounds = (w_lower) + (w_split) + (nx * h_lower) \ 
    201                + (ny * w_upper) + (ny * a_lower) \ 
    202                + (nz * w_upper) + (nz * v_lower) \ 
    203                + (_n * Y_lower)  
    204  #upper_bounds = (nx * w_upper) + (nx * h_upper) \ 
    205   upper_bounds = (w_split) + (w_upper) + (nx * h_upper) \ 
     193 
     194 #XXX XXX: EDITED TO USE npts=(2,1,1) *AND* FIX Y1 == Y_lower 
     195  lower_bounds = (nx * w_lower) + (nx * h_lower) \ 
     196               + (ny * w_lower) + (ny * a_lower) \ 
     197               + (nz * w_lower) + (nz * v_lower) \ 
     198               + (_n * Y_lower) 
     199  upper_bounds = (nx * w_upper) + (nx * h_upper) \ 
    206200               + (ny * w_upper) + (ny * a_upper) \ 
    207201               + (nz * w_upper) + (nz * v_upper) \ 
    208                + (Y_lower) + (Y_upper)  
    209  #             + (_n * Y_upper)  
     202               + (Y_lower) + (Y_upper) 
    210203  bounds = (lower_bounds,upper_bounds) 
    211204 
    … …  
    220213 
    221214  print " model: f(x) = %s(x)" % function_name 
    222   #print " lipschitz constant: L = %s" % L 
    223215  print " target: %s" % str(target) 
    224216  print " error: %s" % str(error) 
    … …  
    262254 #print "solved: %s" % str( c.flatten(all=True) ) 
    263255 
    264   print "fails valid wrt model:\n%s" % c.valid_wrt_model(model, xtol=Cx, \ 
    265                                                          tol=Cy, blamelist=True) 
    266   print "mean(y): %s >= %s" % (str( c.get_mean_value() ), y_mean) 
     256  print "fails valid wrt model:\n%s" % \ 
     257        c.valid_wrt_model(model, xtol=Cx, tol=Cy, blamelist=True, maxiter=iter) 
     258  print "mean(y): %s >= %s" % (str( c.get_mean_value() ), y_mean - y_mean_error) 
    267259  print "sum_wts: %s == 1.0" % [sum(w) for w in c.wts] 
    268260 
     261  from paramtrans import graphical_distance 
     262  Ry = graphical_distance(model, c, tol=Cy, xtol=Cx, cutoff=0.0, maxiter=0) 
     263  print "vertical_distance: %s <= %s" % (Ry, Cy + max(Cx)) 
     264  Rv = graphical_distance(model, c, tol=Cy, xtol=Cx, cutoff=0.0, maxiter=iter) 
     265  print "graphical_distance: %s <= %s" % (Rv, Cy) 
     266 
    269267# EOF 

  • branches/UQ/math/legacy/TEST_OUQ_1dSurr_Cy.py

    r552 r556  
    22""" 
    33OUQ with model sausage: valid_wrt_model, mean(y), norm(wts) 
    4   with lipschitz and bounds from ExampleDataset 
     4  with 'wiggle room' around the model in Y 
    55""" 
    66 
    7 debug = True#False 
     7debug = False 
    88MINMAX = -1  ## NOTE: sup = maximize = -1; inf = minimize = 1 
    99####################################################################### 
    10 # scaling and mpi info; also optimizer configuration parameters 
    11 # hard-wired: use DE solver, don't use mpi, F-F' calculation 
    12 ####################################################################### 
    13 npop = 20 #20 #40  #!!! 
     10# optimizer configuration parameters 
     11####################################################################### 
     12npop = 40  #!!! 
    1413maxiter = 1000 
    1514maxfun = 1e+6 
    … …  
    2019 
    2120####################################################################### 
    22 # the model function and the dataset 
    23 ####################################################################### 
    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  
    29 def model(x): 
    30   return x[0]; 
     21# the model function 
     22####################################################################### 
     23def model(x): return x[0] 
     24 
    3125 
    3226####################################################################### 
    … …  
    110104    # check mean value, and if necessary use constrain to set mean value 
    111105    y = float(mean(c.values, c.weights)) 
    112     if not (y >= float(target[0])):   #XXX: or target[0]-error[0] ? 
     106    if not (y >= float(target[0] - error[0])): 
    113107      c.update( constrain( c.flatten(all=True) ) ) 
    114108 
    115109    # then test if valid... then impose model validity on product measure 
    116     if not c.valid_wrt_model(model, tol=target[2], xtol=target[3]): 
     110    if not c.valid_wrt_model(model, tol=target[2], xtol=target[3], \ 
     111                                                   maxiter=error[4]): 
    117112      c.set_valid(model, cutoff=target[2], bounds=bounds, tol=error[2], \ 
    118                                  constraints=constrain, xtol=target[3]) 
     113                         constraints=constrain, xtol=target[3], imax=target[4]) 
    119114    ###################### more function-specific ##################### 
    120115    if debug: 
    121       if not c.valid_wrt_model(model, tol=target[2], xtol=target[3]): 
     116      if not c.valid_wrt_model(model, tol=target[2], xtol=target[3], \ 
     117                                                     maxiter=error[4]): 
    122118        print "valid_wrt_model: False" 
    123119      if not [sum(w) for w in c.wts] == [1.0] * len(c.wts): 
    124120        print "norm(w) == 1.0: False" 
    125       if not c.get_mean_value() >= target[0]: 
    126         print "mean(y) >= %s: False" 
     121      if not c.get_mean_value() >= float(target[0] - error[0]): 
     122        print "mean(y) >= %s: False" % str(target[0] - error[0]) 
    127123    ###################### end function-specific ###################### 
    128124    # extract weights and positions and values 
    … …  
    141137    ##################### begin function-specific ##################### 
    142138    y = float(mean(c.values, c.weights)) 
    143     if not (y >= float(target[0])): 
     139    if not (y >= float(target[0] - error[0])): 
    144140      if debug: print "skipping mean: %s" % y 
    145141      return inf  #XXX: FORCE TO SATISFY E CONSTRAINTS 
    146142 
    147     if not c.valid_wrt_model(model, tol=target[2], xtol=target[3]): 
     143    if not c.valid_wrt_model(model, tol=target[2], xtol=target[3], \ 
     144                                                   maxiter=error[4]): 
    148145      if debug: print "skipping model-invalidity" 
    149146      return inf  #XXX: FORCE TO SATISFY E CONSTRAINTS 
    … …  
    168165if __name__ == '__main__': 
    169166  function_name = model.__name__ 
    170  #dataset_name = data.__name__ 
    171167 
    172168  from mystic.math.measures import mean 
    … …  
    174170  y_mean_error = 0.0 #NOTE: SET THE 'error' HERE! 
    175171  theta = 0.0      #NOTE: SET THE 'failure criteria' HERE! 
    176   Cy = 1.0          #NOTE: SET THE 'cutoff' HERE! 
     172  Cy = 0.1           #NOTE: SET THE 'cutoff' HERE! 
    177173  Cx = (0.0,0.0,0.0) #NOTE: SET THE 'wiggle' HERE! 
    178   valid_tol = 1.e-6  #NOTE: SET THE 'model tolerance' HERE! 
    179   target = (y_mean,theta,Cy,Cx); error = (y_mean_error,None,valid_tol,None) 
     174  valid_tol = 0.0  #NOTE: SET THE 'model tolerance' HERE! 
     175  iter = 200       #NOTE: SET THE 'Cy-validation maxiter' HERE! 
     176  imax = 10        #NOTE: SET THE 'Cx-validation maxiter' HERE! 
     177  target = (y_mean,theta,Cy,Cx,imax) 
     178  error = (y_mean_error,None,valid_tol,None,iter) 
    180179  pars = (target,error) 
    181180 
    … …  
    187186  _n = _npts(npts) 
    188187 
    189   #FIXME: check units versus those in model and in dataset 
     188  #XXX: check units versus those in model 
    190189  w_lower = [0.0]; Y_lower = [theta] 
    191   w_upper = [1.0]; Y_upper = [5.0] 
     190  w_upper = [1.0]; Y_upper = [10.0] 
    192191  h_lower = [0.0]; a_lower = [0.0];  v_lower = [0.0] 
    193192  h_upper = [1.0]; a_upper = [0.0];  v_upper = [0.0] 
    194   w_split = [0.5] 
    195  
    196  #XXX XXX: EDITED TO USE w_split *AND* npts=(2,1,1) *AND* FIX Y1 == Y_lower 
    197  #lower_bounds = (nx * w_lower) + (nx * h_lower) \ 
    198  #             + (ny * w_lower) + (ny * a_lower) \ 
    199  #             + (nz * w_lower) + (nz * v_lower) \ 
    200   lower_bounds = (w_lower) + (w_split) + (nx * h_lower) \ 
    201                + (ny * w_upper) + (ny * a_lower) \ 
    202                + (nz * w_upper) + (nz * v_lower) \ 
    203                + (_n * Y_lower)  
    204  #upper_bounds = (nx * w_upper) + (nx * h_upper) \ 
    205   upper_bounds = (w_split) + (w_upper) + (nx * h_upper) \ 
     193 
     194 #XXX XXX: EDITED TO USE npts=(2,1,1) *AND* FIX Y1 == Y_lower 
     195  lower_bounds = (nx * w_lower) + (nx * h_lower) \ 
     196               + (ny * w_lower) + (ny * a_lower) \ 
     197               + (nz * w_lower) + (nz * v_lower) \ 
     198               + (_n * Y_lower) 
     199  upper_bounds = (nx * w_upper) + (nx * h_upper) \ 
    206200               + (ny * w_upper) + (ny * a_upper) \ 
    207201               + (nz * w_upper) + (nz * v_upper) \ 
    208                + (Y_lower) + (Y_upper)  
    209  #             + (_n * Y_upper)  
     202               + (Y_lower) + (Y_upper) 
    210203  bounds = (lower_bounds,upper_bounds) 
    211204 
    … …  
    220213 
    221214  print " model: f(x) = %s(x)" % function_name 
    222   #print " lipschitz constant: L = %s" % L 
    223215  print " target: %s" % str(target) 
    224216  print " error: %s" % str(error) 
    … …  
    262254 #print "solved: %s" % str( c.flatten(all=True) ) 
    263255 
    264   print "fails valid wrt model:\n%s" % c.valid_wrt_model(model, xtol=Cx, \ 
    265                                                          tol=Cy, blamelist=True) 
    266   print "mean(y): %s >= %s" % (str( c.get_mean_value() ), y_mean) 
     256  print "fails valid wrt model:\n%s" % \ 
     257        c.valid_wrt_model(model, xtol=Cx, tol=Cy, blamelist=True, maxiter=iter) 
     258  print "mean(y): %s >= %s" % (str( c.get_mean_value() ), y_mean - y_mean_error) 
    267259  print "sum_wts: %s == 1.0" % [sum(w) for w in c.wts] 
    268260 
     261  from paramtrans import graphical_distance 
     262  Ry = graphical_distance(model, c, tol=Cy, xtol=Cx, cutoff=0.0, maxiter=0) 
     263  print "vertical_distance: %s <= %s" % (Ry, Cy + max(Cx)) 
     264  Rv = graphical_distance(model, c, tol=Cy, xtol=Cx, cutoff=0.0, maxiter=iter) 
     265  print "graphical_distance: %s <= %s" % (Rv, Cy) 
     266 
    269267# EOF 

  • branches/UQ/math/legacy/TEST_OUQ_StAlData.py

    r546 r556  
    88MINMAX = -1  ## NOTE: sup = maximize = -1; inf = minimize = 1 
    99####################################################################### 
    10 # scaling and mpi info; also optimizer configuration parameters 
    11 # hard-wired: use DE solver, don't use mpi, F-F' calculation 
     10# optimizer configuration parameters 
    1211####################################################################### 
    1312npop = 32  #!!! 
    … …  
    2019 
    2120####################################################################### 
    22 # the model function and the dataset 
    23 ####################################################################### 
    24 #from StAlSurrogate import st_al_surr as model 
     21# the dataset 
     22####################################################################### 
    2523from legacydata import load_dataset 
    2624data = load_dataset('StAlDataset.txt') 
    … …  
    107105    # check mean value, and if necessary use constrain to set mean value 
    108106    y = float(mean(c.values, c.weights)) 
    109     if not (y >= float(target[0])):   #XXX: or target[0]-error[0] ? 
    110      #c.values = impose_mean(target[0]+target[1], c.values, c.weights) 
     107    if not (y >= float(target[0] - error[0])): 
    111108      c.update( constrain( c.flatten(all=True) ) ) 
    112109 
    … …  
    123120      if not [sum(w) for w in c.wts] == [1.0] * len(c.wts): 
    124121        print "norm(w) == 1.0: False" 
    125       if not c.get_mean_value() >= target[0]: 
    126         print "mean(y) >= %s: False"#%s" % (target[0], c.get_mean_value() \ 
    127         #                                                   >= target[0]) 
     122      if not c.get_mean_value() >= float(target[0] - error[0]): 
     123        print "mean(y) >= %s: False" % str(target[0] - error[0]) 
    128124    ###################### end function-specific ###################### 
    129125    # extract weights and positions and values 
    … …  
    142138    ##################### begin function-specific ##################### 
    143139    y = float(mean(c.values, c.weights)) 
    144     if not (y >= float(target[0])): 
     140    if not (y >= float(target[0] - error[0])): 
    145141      if debug: print "skipping mean: %s" % y 
    146142      return inf  #XXX: FORCE TO SATISFY E CONSTRAINTS 
    … …  
    172168####################################################################### 
    173169if __name__ == '__main__': 
    174  #function_name = model.__name__ 
    175170  dataset_name = data.__name__ 
    176  
    177  #model_tol = 2.0  #expect difference is in [-1.677, 1.507], where 2 > |-1.677| 
    178171 
    179172  from mystic.math.measures import mean 
    180173# y_mean = mean(data.values) 
    181174  y_mean = 11.0    #NOTE: SET THE 'mean' HERE! 
    182   y_error = 0.0001 #NOTE: SET THE 'error' HERE! 
     175  y_mean_error = 0.0001 #NOTE: SET THE 'error' HERE! 
    183176  theta = 7.0      #NOTE: SET THE 'failure criteria' HERE! 
    184177  short_tol = 1e-9 #NOTE: SET THE 'short tolerance' HERE! 
    185   target = (y_mean,theta,None); error = (y_error,None,short_tol) 
     178  target = (y_mean,theta,None); error = (y_mean_error,None,short_tol) 
    186179  pars = (target,error) 
    187180 
    … …  
    193186  _n = _npts(npts) 
    194187 
    195   #FIXME: check units versus those in model and in dataset 
     188  #XXX: check units versus those in dataset 
    196189  w_lower = [0.0]; Y_lower = [0.0] 
    197190  w_upper = [1.0]; Y_upper = [100.0] 
    … …  
    201194 
    202195 #XXX XXX: EDITED TO USE w_split 
    203  #lower_bounds = (nx * w_lower) + (nx * h_lower) \ 
    204  #             + (ny * w_lower) + (ny * a_lower) \ 
    205  #             + (nz * w_lower) + (nz * v_lower) \ 
    206196  lower_bounds = (w_lower) + (w_split) + (nx * h_lower) \ 
    207197               + (w_lower) + (w_split) + (ny * a_lower) \ 
    208198               + (w_lower) + (w_split) + (nz * v_lower) \ 
    209199               + (_n * Y_lower)  
    210  #upper_bounds = (nx * w_upper) + (nx * h_upper) \ 
    211  #             + (ny * w_upper) + (ny * a_upper) \ 
    212  #             + (nz * w_upper) + (nz * v_upper) \ 
    213200  upper_bounds = (w_split) + (w_upper) + (nx * h_upper) \ 
    214201               + (w_split) + (w_upper) + (ny * a_upper) \ 
    … …  
    226213  print "..............\n" 
    227214 
    228  #print " model: f(x) = %s(x)" % function_name 
    229215  print " data: z = {%s}" % dataset_name 
    230216  print " target: %s" % str(target) 
    … …  
    273259  print "fails short wrt data:\n%s" % c.short_wrt_data(data, tol=short_tol,\ 
    274260                                                             blamelist=True) 
    275   print "mean(y): %s >= %s" % (str( c.get_mean_value() ), y_mean) 
     261  print "mean(y): %s >= %s" % (str( c.get_mean_value() ), y_mean - y_mean_error) 
    276262  print "sum_wts: %s == 1.0" % [sum(w) for w in c.wts] 
    277263 

  • branches/UQ/math/legacy/TEST_OUQ_StStSurr_Cy.py

    r552 r556  
    22""" 
    33OUQ with model sausage: valid_wrt_model, mean(y), norm(wts) 
    4   with lipschitz and bounds from ExampleDataset 
     4  with 'wiggle room' around the model in Y 
    55""" 
    66 
    7 debug = True#False 
     7debug = False 
    88MINMAX = -1  ## NOTE: sup = maximize = -1; inf = minimize = 1 
    99####################################################################### 
    10 # scaling and mpi info; also optimizer configuration parameters 
    11 # hard-wired: use DE solver, don't use mpi, F-F' calculation 
    12 ####################################################################### 
    13 npop = 20 #20 #40  #!!! 
     10# optimizer configuration parameters 
     11####################################################################### 
     12npop = 40 #50  #!!! 
    1413maxiter = 1000 
    1514maxfun = 1e+6 
    … …  
    2019 
    2120####################################################################### 
    22 # the model function and the dataset 
    23 ####################################################################### 
    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  
    29 def model(x): 
    30   return x[0]; 
     21# the model function 
     22####################################################################### 
     23from StStSurrogate import marc_surr as model 
     24 
    3125 
    3226####################################################################### 
    … …  
    110104    # check mean value, and if necessary use constrain to set mean value 
    111105    y = float(mean(c.values, c.weights)) 
    112     if not (y >= float(target[0])):   #XXX: or target[0]-error[0] ? 
     106    if not (y >= float(target[0] - error[0])): 
    113107      c.update( constrain( c.flatten(all=True) ) ) 
    114108 
    115109    # then test if valid... then impose model validity on product measure 
    116     if not c.valid_wrt_model(model, tol=target[2], xtol=target[3]): 
     110    if not c.valid_wrt_model(model, tol=target[2], xtol=target[3], \ 
     111                                                   maxiter=error[4]): 
    117112      c.set_valid(model, cutoff=target[2], bounds=bounds, tol=error[2], \ 
    118                                  constraints=constrain, xtol=target[3]) 
     113                         constraints=constrain, xtol=target[3], imax=target[4]) 
    119114    ###################### more function-specific ##################### 
    120115    if debug: 
    121       if not c.valid_wrt_model(model, tol=target[2], xtol=target[3]): 
     116      if not c.valid_wrt_model(model, tol=target[2], xtol=target[3], \ 
     117                                                     maxiter=error[4]): 
    122118        print "valid_wrt_model: False" 
    123119      if not [sum(w) for w in c.wts] == [1.0] * len(c.wts): 
    124120        print "norm(w) == 1.0: False" 
    125       if not c.get_mean_value() >= target[0]: 
    126         print "mean(y) >= %s: False" 
     121      if not c.get_mean_value() >= float(target[0] - error[0]): 
     122        print "mean(y) >= %s: False" % str(target[0] - error[0]) 
    127123    ###################### end function-specific ###################### 
    128124    # extract weights and positions and values 
    … …  
    141137    ##################### begin function-specific ##################### 
    142138    y = float(mean(c.values, c.weights)) 
    143     if not (y >= float(target[0])): 
     139    if not (y >= float(target[0] - error[0])): 
    144140      if debug: print "skipping mean: %s" % y 
    145141      return inf  #XXX: FORCE TO SATISFY E CONSTRAINTS 
    146142 
    147     if not c.valid_wrt_model(model, tol=target[2], xtol=target[3]): 
     143    if not c.valid_wrt_model(model, tol=target[2], xtol=target[3], \ 
     144                                                   maxiter=error[4]): 
    148145      if debug: print "skipping model-invalidity" 
    149146      return inf  #XXX: FORCE TO SATISFY E CONSTRAINTS 
    … …  
    168165if __name__ == '__main__': 
    169166  function_name = model.__name__ 
    170  #dataset_name = data.__name__ 
    171167 
    172168  from mystic.math.measures import mean 
    173   y_mean = 0.5     #NOTE: SET THE 'mean' HERE! 
    174   y_mean_error = 0.0 #NOTE: SET THE 'error' HERE! 
    175   theta = 0.0      #NOTE: SET THE 'failure criteria' HERE! 
    176   Cy = 1.0          #NOTE: SET THE 'cutoff' HERE! 
     169  y_mean = 6.5    #NOTE: SET THE 'mean' HERE! 
     170  y_mean_error = 1.0 #NOTE: SET THE 'error' HERE! 
     171  theta = 5.0      #NOTE: SET THE 'failure criteria' HERE! 
     172  Cy = 0.05          #NOTE: SET THE 'cutoff' HERE! 
    177173  Cx = (0.0,0.0,0.0) #NOTE: SET THE 'wiggle' HERE! 
    178   valid_tol = 1.e-6  #NOTE: SET THE 'model tolerance' HERE! 
    179   target = (y_mean,theta,Cy,Cx); error = (y_mean_error,None,valid_tol,None) 
     174  valid_tol = 0.0  #NOTE: SET THE 'model tolerance' HERE! 
     175  iter = 200       #NOTE: SET THE 'Cy-validation maxiter' HERE! 
     176  imax = 10        #NOTE: SET THE 'Cx-validation maxiter' HERE! 
     177  target = (y_mean,theta,Cy,Cx,imax) 
     178  error = (y_mean_error,None,valid_tol,None,iter) 
    180179  pars = (target,error) 
    181180 
    … …  
    187186  _n = _npts(npts) 
    188187 
    189   #FIXME: check units versus those in model and in dataset 
     188  #XXX: check units versus those in model 
    190189  w_lower = [0.0]; Y_lower = [theta] 
    191   w_upper = [1.0]; Y_upper = [5.0] 
    192   h_lower = [0.0]; a_lower = [0.0];  v_lower = [0.0] 
    193   h_upper = [1.0]; a_upper = [0.0];  v_upper = [0.0] 
    194   w_split = [0.5] 
    195  
    196  #XXX XXX: EDITED TO USE w_split *AND* npts=(2,1,1) *AND* FIX Y1 == Y_lower 
    197  #lower_bounds = (nx * w_lower) + (nx * h_lower) \ 
    198  #             + (ny * w_lower) + (ny * a_lower) \ 
    199  #             + (nz * w_lower) + (nz * v_lower) \ 
    200   lower_bounds = (w_lower) + (w_split) + (nx * h_lower) \ 
    201                + (ny * w_upper) + (ny * a_lower) \ 
    202                + (nz * w_upper) + (nz * v_lower) \ 
    203                + (_n * Y_lower)  
    204  #upper_bounds = (nx * w_upper) + (nx * h_upper) \ 
    205   upper_bounds = (w_split) + (w_upper) + (nx * h_upper) \ 
     190  w_upper = [1.0]; Y_upper = [100.0] 
     191  h_lower = [60.0];  a_lower = [0.0];  v_lower = [2.1] 
     192  h_upper = [105.0]; a_upper = [30.0]; v_upper = [2.8] 
     193 
     194 #XXX XXX: EDITED TO USE npts=(2,1,1) *AND* FIX Y1 == Y_lower 
     195  lower_bounds = (nx * w_lower) + (nx * h_lower) \ 
     196               + (ny * w_lower) + (ny * a_lower) \ 
     197               + (nz * w_lower) + (nz * v_lower) \ 
     198               + (_n * Y_lower) 
     199  upper_bounds = (nx * w_upper) + (nx * h_upper) \ 
    206200               + (ny * w_upper) + (ny * a_upper) \ 
    207201               + (nz * w_upper) + (nz * v_upper) \ 
    208                + (Y_lower) + (Y_upper)  
    209  #             + (_n * Y_upper)  
     202               + (Y_lower) + (Y_upper) 
    210203  bounds = (lower_bounds,upper_bounds) 
    211204 
    … …  
    220213 
    221214  print " model: f(x) = %s(x)" % function_name 
    222   #print " lipschitz constant: L = %s" % L 
    223215  print " target: %s" % str(target) 
    224216  print " error: %s" % str(error) 
    … …  
    262254 #print "solved: %s" % str( c.flatten(all=True) ) 
    263255 
    264   print "fails valid wrt model:\n%s" % c.valid_wrt_model(model, xtol=Cx, \ 
    265                                                          tol=Cy, blamelist=True) 
    266   print "mean(y): %s >= %s" % (str( c.get_mean_value() ), y_mean) 
     256  print "fails valid wrt model:\n%s" % \ 
     257        c.valid_wrt_model(model, xtol=Cx, tol=Cy, blamelist=True, maxiter=iter) 
     258  print "mean(y): %s >= %s" % (str( c.get_mean_value() ), y_mean - y_mean_error) 
    267259  print "sum_wts: %s == 1.0" % [sum(w) for w in c.wts] 
    268260 
     261  from paramtrans import graphical_distance 
     262  Ry = graphical_distance(model, c, tol=Cy, xtol=Cx, cutoff=0.0, maxiter=0) 
     263  print "vertical_distance: %s <= %s" % (Ry, Cy + max(Cx)) 
     264  Rv = graphical_distance(model, c, tol=Cy, xtol=Cx, cutoff=0.0, maxiter=iter) 
     265  print "graphical_distance: %s <= %s" % (Rv, Cy) 
     266 
    269267# EOF 

  • branches/UQ/math/legacy/dirac_measure.py

    r547 r556  
    10121012 
    10131013  # construct and configure optimizer 
    1014   debug = False  #!!! 
     1014  debug = True#False  #!!! 
    10151015  maxiter = 1000;  maxfun = 1e+6 
    10161016  crossover = 0.9; percent_change = 0.9 
    … …  
    10281028  if debug: stepmon = VerboseMonitor(10)  #!!! 
    10291029  if npop: # use VTR 
     1030    _i = 2 #XXX: iter returned as results[2] 
    10301031    results = diffev2(cost, guess, npop, ftol=ftol, gtol=gtol, bounds=bounds,\ 
    10311032                      maxiter=maxiter, maxfun=maxfun, constraints=constraints,\ 
    … …  
    10341035                      full_output=1, disp=0, handler=False) 
    10351036  else: # use VTR 
     1037    _i = 3 #XXX: iter returned as results[3]  (results[2] == direc) 
    10361038    results = fmin_powell(cost, guess, ftol=ftol, gtol=gtol, bounds=bounds,\ 
    10371039                      maxiter=maxiter, maxfun=maxfun, constraints=constraints,\ 
    … …  
    10421044  pm.load(results[0], pts)            # params: w,x,y 
    10431045  if debug: print "final cost: %s" % results[1] 
    1044   if debug and results[2] >= maxiter: # iterations 
     1046  if debug and results[_i] >= maxiter: # iterations 
    10451047    print "Warning: constraints solver terminated at maximum iterations" 
    1046  #func_evals = results[3]             # evaluation 
     1048 #func_evals = results[_i+1]           # evaluation 
    10471049  return pm 
    10481050 
    … …  
    11181120  xtol = 0.0 # default 
    11191121  if kwds.has_key('xtol'): xtol = kwds['xtol'] 
    1120   npop = 20 #XXX: tune npop (outer optimization)? 
     1122  npop = 20 #None #XXX: tune npop (outer optimization)? 
    11211123  if kwds.has_key('npop'): npop = kwds.pop('npop') 
     1124  maxiter = 50 #XXX: tune maxiter (outer optimization)? 
     1125  if kwds.has_key('maxiter'): maxiter = kwds.pop('maxiter') 
    11221126  # tolerance for optimization on sum(y) 
    11231127  tol = 0.0 # default 
    11241128  if kwds.has_key('tol'): tol = kwds.pop('tol') 
    1125   ipop = None #XXX: tune npop (inner optimization)? 
     1129  ipop = 10 #XXX: tune npop (inner optimization)? 
    11261130  if kwds.has_key('ipop'): ipop = kwds.pop('ipop') 
     1131  imax = 10  #XXX: tune maxiter (inner optimization)? 
     1132  if kwds.has_key('imax'): imax = kwds.pop('imax') 
    11271133 
    11281134  # if no guess was made, then use bounds constraints 
    … …  
    11411147    points.load(rv, pts) 
    11421148    # calculate infeasibility 
    1143     Rv = graphical_distance(model, points, tol=cutoff, npop=ipop, **kwds) 
     1149    Rv = graphical_distance(model, points, tol=cutoff, npop=ipop, \ 
     1150                                           maxiter=imax, **kwds) 
    11441151    v = infeasibility(Rv, cutoff) 
    11451152    # converting v to E 
    … …  
    11481155  # construct and configure optimizer 
    11491156  debug = False  #!!! 
    1150   maxiter = 1000;  maxfun = 1e+6 
     1157  maxfun = 1e+6 
    11511158  crossover = 0.9; percent_change = 0.8 
    1152   ftol = abs(tol); gtol = None 
     1159  ftol = abs(tol); gtol = None #XXX: optimally, should be VTRCOG... 
    11531160 
    11541161  if debug: 
    … …  
    11611168  from mystic.strategy import Best1Bin, Best1Exp 
    11621169  evalmon = Monitor();  stepmon = Monitor(); strategy = Best1Exp 
    1163   if debug: stepmon = VerboseMonitor(1)  #!!! 
     1170  if debug: stepmon = VerboseMonitor(2)  #!!! 
    11641171  if npop: # use VTR 
     1172    _i = 2 #XXX: iter returned as results[2] 
    11651173    results = diffev2(cost, guess, npop, ftol=ftol, gtol=gtol, bounds=bounds,\ 
    11661174                      maxiter=maxiter, maxfun=maxfun, constraints=constraints,\ 
    … …  
    11691177                      full_output=1, disp=0, handler=False) 
    11701178  else: # use VTR 
     1179    _i = 3 #XXX: iter returned as results[3]  (results[2] == direc) 
    11711180    results = fmin_powell(cost, guess, ftol=ftol, gtol=gtol, bounds=bounds,\ 
    11721181                      maxiter=maxiter, maxfun=maxfun, constraints=constraints,\ 
    … …  
    11771186  pm.load(results[0], pts)            # params: w,x,y 
    11781187 #if debug: print "final cost: %s" % results[1] 
    1179   if debug and results[2] >= maxiter: # iterations 
     1188  if debug and results[_i] >= maxiter: # iterations 
    11801189    print "Warning: constraints solver terminated at maximum iterations" 
    1181  #func_evals = results[3]             # evaluation 
     1190 #func_evals = results[_i+1]           # evaluation 
    11821191  return pm 
    11831192 

  • branches/UQ/math/legacy/paramtrans.py

    r547 r556  
    337337    if cutoff: cutoff = tol 
    338338    else: cutoff = None 
    339   npop = min(20, 4*nxi) #XXX: tune npop? 
     339  npop = min(20, 3*nxi) #XXX: tune npop? 
    340340  if kwds.has_key('npop'): npop = kwds.pop('npop') 
     341  imax = 1000 #XXX: tune maxiter? 
     342  if kwds.has_key('maxiter'): imax = kwds.pop('maxiter') 
    341343 
    342344  ######################################################################### 
    343   def radius(model, point, tol=0.0, xtol=0.0, npop=None): 
     345  def radius(model, point, tol=0.0, xtol=0.0, npop=None, maxiter=None): 
    344346    """graphical distance between a single point x,y and a model F(x')""" 
    345347    # given a single point x,y: find the radius = |y - F(x')| 
    … …  
    359361 
    360362    if debug: 
    361       print "rv: %s" % x 
     363      print "rv: %s" % str(x) 
    362364      print "cost: %s" % cost(x) 
    363365 
    364     if not xtol:  # then radius is difference in x,y and x,F(x) 
    365       return cost(x) 
     366    # if xtol=0, radius is difference in x,y and x,F(x); skip the optimization 
     367    try: 
     368      if not maxiter or not max(xtol): #iterables 
     369        return cost(x) 
     370    except TypeError: 
     371      if not xtol: #non-iterables 
     372        return cost(x) 
    366373 
    367374    # set the range constraints 
    … …  
    373380    if debug: stepmon = VerboseMonitor(1) 
    374381    #XXX: edit settings? 
    375     imax=3000; MINMAX = 1 #XXX: confirm MINMAX=1 is minimization 
     382    MINMAX = 1 #XXX: confirm MINMAX=1 is minimization 
    376383    if npop: # use VTR 
    377384      results = diffev2(cost, x, npop, ftol=tol, gtol=None, \ 
    378                         itermon = stepmon, maxiter=imax, bounds=bounds, \ 
     385                        itermon = stepmon, maxiter=maxiter, bounds=bounds, \ 
    379386                        full_output=1, disp=0, handler=False) 
    380387    else: # use VTR 
    381388      results = fmin_powell(cost, x, ftol=tol, gtol=None, \ 
    382                             itermon = stepmon, bounds=bounds, \ 
     389                            itermon = stepmon, maxiter=maxiter, bounds=bounds, \ 
    383390                            full_output=1, disp=0, handler=False) 
    384391   #solved = results[0]            # x' 
    … …  
    393400 
    394401  #XXX: better to do a single optimization rather than for each point ??? 
    395   d = [radius(model, point, tol, xtol, npop) for point in target] 
     402  d = [radius(model, point, tol, xtol, npop, imax) for point in target] 
    396403  return infeasibility(d, cutoff) 
    397404 

  • branches/UQ/math/legacy/test_ExampleDataset.py

    r546 r556  
    5555  print("testing with some 'wiggle-room'...") 
    5656  Cs = 0.25 
    57   Cx = 1.0; Cy = 2.0 
    58   npop = None # optimizer population per generation 
     57  Cx = 0.5; Cy = 0.25 
    5958 
    6059  # NOTES:   for no 'cutoff', use cutoff=None;  for is_feasible, use raw=False 
    … …  
    6463 
    6564  # Check model validity 
    66   Rv = ex1d_data.valid(F, xtol=Cx, tol=Cy, npop=npop, \ 
     65  Rv = ex1d_data.valid(F, xtol=Cx, tol=Cy, \ 
    6766                       raw=True, blamelist=False, pairs=False, all=False) 
    6867  print("valid(ytol=%s,xtol=%s):\n%s" % (Cy, Cx, Rv)) 

  • branches/UQ/math/legacy/test_ModeledDataset.py

    r546 r556  
    2727  def added(x): 
    2828    return sum(x) 
     29  def first_mass(x): 
     30    return x[0] 
     31 #F = first_mass 
    2932  F = sum_squared 
    3033 #F = added 
    … …  
    3740  # settings 
    3841  print("testing with some 'wiggle-room'...") 
    39   npop = None # optimizer population per generation 
    4042 
    4143  # Check model validity 
    4244  Cx = 0.0; Cy = 1.0 
    43   Rv = data.valid(F, xtol=Cx, tol=Cy, npop=npop, \ 
     45  Rv = data.valid(F, xtol=Cx, tol=Cy, \ 
    4446                       raw=False, blamelist=False, pairs=False, all=True) 
    4547  print("valid(ytol=%s,xtol=%s):\n%s" % (Cy, Cx, Rv)) 
    4648 
    4749  Cx = 0.0; Cy = 2.0 
    48   Rv = data.valid(F, xtol=Cx, tol=Cy, npop=npop, \ 
     50  Rv = data.valid(F, xtol=Cx, tol=Cy, \ 
    4951                       raw=False, blamelist=False, pairs=False, all=True) 
    5052  print("valid(ytol=%s,xtol=%s):\n%s" % (Cy, Cx, Rv)) 
    5153 
    5254  Cx = 1.0; Cy = 2.0 
    53   Rv = data.valid(F, xtol=Cx, tol=Cy, npop=npop, \ 
     55  Rv = data.valid(F, xtol=Cx, tol=Cy, \ 
    5456                       raw=False, blamelist=False, pairs=False, all=True) 
    5557  print("valid(ytol=%s,xtol=%s):\n%s" % (Cy, Cx, Rv)) 
    5658 
    5759  Cx = 2.0; Cy = 2.0 
    58   Rv = data.valid(F, xtol=Cx, tol=Cy, npop=npop, \ 
     60  Rv = data.valid(F, xtol=Cx, tol=Cy, \ 
    5961                       raw=False, blamelist=False, pairs=False, all=True) 
    6062  print("valid(ytol=%s,xtol=%s):\n%s" % (Cy, Cx, Rv)) 

  • branches/UQ/math/legacy/test_StAlDataset.py

    r546 r556  
    4545  print("testing with some 'wiggle-room'...") 
    4646  Cs = 0.25 
    47   Cx = 50.0; Cy = 50.0 
    48   npop = None # optimizer population per generation 
     47  Cx = 0.0; Cy = 1.75 
    4948 
    5049  # NOTES:   for no 'cutoff', use cutoff=None;  for is_feasible, use raw=False 
    … …  
    5453 
    5554  # Check model validity 
    56   Rv = st_al_data.valid(F, xtol=Cx, tol=Cy, npop=npop, \ 
     55  Rv = st_al_data.valid(F, xtol=Cx, tol=Cy, \ 
    5756                        raw=True, blamelist=False, pairs=False, all=False) 
    5857  print("valid(ytol=%s,xtol=%s):\n%s" % (Cy, Cx, Rv)) 
Note: See TracChangeset for help on using the changeset viewer.