Changeset 569
- Timestamp:
- 09/24/12 15:59:25 (4 years ago)
- Files:
-
- 14 deleted
- 14 edited
- 4 moved
Legend:
- Unmodified
- Added
- Removed
-
branches/UQ/math/legacy/MINMAX_StAlData.py
r556 r569 17 17 # the dataset 18 18 ####################################################################### 19 from legacydata import load_dataset19 from mystic.math.legacydata import load_dataset 20 20 data = load_dataset('StAlDataset.txt') 21 21 … … 67 67 def maximize(params,npts,bounds): 68 68 69 from dirac_measure import scenario69 from mystic.math.dirac_measure import scenario 70 70 from numpy import inf 71 71 target,error = params … … 96 96 ##################### begin function-specific ##################### 97 97 # impose norm on the weights of the discrete measures 98 from dirac_measure import norm_wts_constraintsFactory as factory98 from mystic.math.dirac_measure import norm_wts_constraintsFactory as factory 99 99 constrain = factory(npts) 100 100 … … 158 158 nz = 1 #NOTE: SET THE NUMBER OF 'v' POINTS HERE! 159 159 npts = (nx,ny,nz) 160 from paramtrans import _npts160 from mystic.math.paramtrans import _npts 161 161 _n = _npts(npts) 162 162 … … 219 219 220 220 from numpy import array 221 from dirac_measure import scenario221 from mystic.math.dirac_measure import scenario 222 222 c = scenario() 223 223 c.load(solved,npts) -
branches/UQ/math/legacy/MM_OUQ_StAlData.py
r556 r569 23 23 # the dataset 24 24 ####################################################################### 25 from legacydata import load_dataset25 from mystic.math.legacydata import load_dataset 26 26 data = load_dataset('StAlDataset.txt') 27 27 … … 75 75 def maximize(params,npts,bounds): 76 76 77 from dirac_measure import scenario77 from mystic.math.dirac_measure import scenario 78 78 from numpy import inf 79 79 target,error = params … … 104 104 ##################### begin function-specific ##################### 105 105 # impose mean on the values of the product measure 106 from dirac_measure import mean_y_norm_wts_constraintsFactory as factory106 from mystic.math.dirac_measure import mean_y_norm_wts_constraintsFactory as factory 107 107 constrain = factory((target[0],error[0]), npts) 108 108 … … 187 187 nz = 2 #NOTE: SET THE NUMBER OF 'v' POINTS HERE! 188 188 npts = (nx,ny,nz) 189 from paramtrans import _npts189 from mystic.math.paramtrans import _npts 190 190 _n = _npts(npts) 191 191 … … 250 250 251 251 from numpy import array 252 from dirac_measure import scenario252 from mystic.math.dirac_measure import scenario 253 253 c = scenario() 254 254 c.load(solved,npts) -
branches/UQ/math/legacy/TEST_OUQ_1dData.py
r556 r569 21 21 # the dataset 22 22 ####################################################################### 23 from legacydata import load_dataset23 from mystic.math.legacydata import load_dataset 24 24 data = load_dataset('ExampleDataset.txt', range(0,1)) 25 25 … … 71 71 def maximize(params,npts,bounds): 72 72 73 from dirac_measure import scenario73 from mystic.math.dirac_measure import scenario 74 74 from numpy import inf 75 75 target,error = params … … 100 100 ##################### begin function-specific ##################### 101 101 # impose mean on the values of the product measure 102 from dirac_measure import mean_y_norm_wts_constraintsFactory as factory102 from mystic.math.dirac_measure import mean_y_norm_wts_constraintsFactory as factory 103 103 constrain = factory((target[0],error[0]), npts) 104 104 … … 183 183 nz = 1 #NOTE: SET THE NUMBER OF 'v' POINTS HERE! 184 184 npts = (nx,ny,nz) 185 from paramtrans import _npts185 from mystic.math.paramtrans import _npts 186 186 _n = _npts(npts) 187 187 … … 246 246 247 247 from numpy import array 248 from dirac_measure import scenario248 from mystic.math.dirac_measure import scenario 249 249 c = scenario() 250 250 c.load(solved,npts) -
branches/UQ/math/legacy/TEST_OUQ_1dSurr_CxCy.py
r557 r569 70 70 def maximize(params,npts,bounds): 71 71 72 from dirac_measure import scenario72 from mystic.math.dirac_measure import scenario 73 73 from numpy import inf 74 74 target,error = params … … 99 99 ##################### begin function-specific ##################### 100 100 # impose mean on the values of the product measure 101 from dirac_measure import mean_y_norm_wts_constraintsFactory as factory101 from mystic.math.dirac_measure import mean_y_norm_wts_constraintsFactory as factory 102 102 constrain = factory((target[0],error[0]), npts) 103 103 … … 185 185 nz = 1 #NOTE: SET THE NUMBER OF 'v' POINTS HERE! 186 186 npts = (nx,ny,nz) 187 from paramtrans import _npts187 from mystic.math.paramtrans import _npts 188 188 _n = _npts(npts) 189 189 … … 247 247 248 248 from numpy import array 249 from dirac_measure import scenario249 from mystic.math.dirac_measure import scenario 250 250 c = scenario() 251 251 c.load(solved,npts) … … 261 261 print "sum_wts: %s == 1.0" % [sum(w) for w in c.wts] 262 262 263 from paramtrans import graphical_distance263 from mystic.math.paramtrans import graphical_distance 264 264 Ry = graphical_distance(model, c, ytol=Cy, xtol=Cx, cutoff=0.0, imax=0) 265 265 print "vertical_distance: %s <= %s" % (Ry, Cy + max(Cx)) -
branches/UQ/math/legacy/TEST_OUQ_1dSurr_Cy.py
r557 r569 70 70 def maximize(params,npts,bounds): 71 71 72 from dirac_measure import scenario72 from mystic.math.dirac_measure import scenario 73 73 from numpy import inf 74 74 target,error = params … … 99 99 ##################### begin function-specific ##################### 100 100 # impose mean on the values of the product measure 101 from dirac_measure import mean_y_norm_wts_constraintsFactory as factory101 from mystic.math.dirac_measure import mean_y_norm_wts_constraintsFactory as factory 102 102 constrain = factory((target[0],error[0]), npts) 103 103 … … 185 185 nz = 1 #NOTE: SET THE NUMBER OF 'v' POINTS HERE! 186 186 npts = (nx,ny,nz) 187 from paramtrans import _npts187 from mystic.math.paramtrans import _npts 188 188 _n = _npts(npts) 189 189 … … 247 247 248 248 from numpy import array 249 from dirac_measure import scenario249 from mystic.math.dirac_measure import scenario 250 250 c = scenario() 251 251 c.load(solved,npts) … … 261 261 print "sum_wts: %s == 1.0" % [sum(w) for w in c.wts] 262 262 263 from paramtrans import graphical_distance263 from mystic.math.paramtrans import graphical_distance 264 264 Ry = graphical_distance(model, c, ytol=Cy, xtol=Cx, cutoff=0.0, imax=0) 265 265 print "vertical_distance: %s <= %s" % (Ry, Cy + max(Cx)) -
branches/UQ/math/legacy/TEST_OUQ_StAlData.py
r558 r569 21 21 # the dataset 22 22 ####################################################################### 23 from legacydata import load_dataset23 from mystic.math.legacydata import load_dataset 24 24 data = load_dataset('StAlDataset.txt') 25 25 … … 71 71 def maximize(params,npts,bounds): 72 72 73 from dirac_measure import scenario73 from mystic.math.dirac_measure import scenario 74 74 from numpy import inf 75 75 target,error = params … … 100 100 ##################### begin function-specific ##################### 101 101 # impose mean on the values of the product measure 102 from dirac_measure import mean_y_norm_wts_constraintsFactory as factory102 from mystic.math.dirac_measure import mean_y_norm_wts_constraintsFactory as factory 103 103 constrain = factory((target[0],error[0]), npts) 104 104 … … 184 184 nz = 2 #NOTE: SET THE NUMBER OF 'v' POINTS HERE! 185 185 npts = (nx,ny,nz) 186 from paramtrans import _npts186 from mystic.math.paramtrans import _npts 187 187 _n = _npts(npts) 188 188 … … 247 247 248 248 from numpy import array 249 from dirac_measure import scenario249 from mystic.math.dirac_measure import scenario 250 250 c = scenario() 251 251 c.load(solved,npts) -
branches/UQ/math/legacy/TEST_OUQ_StStSurr_Cy.py
r557 r569 70 70 def maximize(params,npts,bounds): 71 71 72 from dirac_measure import scenario72 from mystic.math.dirac_measure import scenario 73 73 from numpy import inf 74 74 target,error = params … … 99 99 ##################### begin function-specific ##################### 100 100 # impose mean on the values of the product measure 101 from dirac_measure import mean_y_norm_wts_constraintsFactory as factory101 from mystic.math.dirac_measure import mean_y_norm_wts_constraintsFactory as factory 102 102 constrain = factory((target[0],error[0]), npts) 103 103 … … 185 185 nz = 1 #NOTE: SET THE NUMBER OF 'v' POINTS HERE! 186 186 npts = (nx,ny,nz) 187 from paramtrans import _npts187 from mystic.math.paramtrans import _npts 188 188 _n = _npts(npts) 189 189 … … 247 247 248 248 from numpy import array 249 from dirac_measure import scenario249 from mystic.math.dirac_measure import scenario 250 250 c = scenario() 251 251 c.load(solved,npts) … … 261 261 print "sum_wts: %s == 1.0" % [sum(w) for w in c.wts] 262 262 263 from paramtrans import graphical_distance263 from mystic.math.paramtrans import graphical_distance 264 264 Ry = graphical_distance(model, c, ytol=Cy, xtol=Cx, cutoff=0.0, imax=0) 265 265 print "vertical_distance: %s <= %s" % (Ry, Cy + max(Cx)) -
branches/UQ/math/legacy/test_ExampleDataset.py
r557 r569 7 7 Mathematical Modeling and Numerical Analysis (submitted 2012). 8 8 """ 9 from legacydata import load_dataset9 from mystic.math.legacydata import load_dataset 10 10 datafile = 'ExampleDataset.txt' 11 11 … … 31 31 32 32 # Check shortness 33 from paramtrans import lipschitz_distance, graphical_distance33 from mystic.math.paramtrans import lipschitz_distance, graphical_distance 34 34 print("\nshort: %s" % ex1d_data.short()) 35 35 L = ex1d_data.lipschitz; -
branches/UQ/math/legacy/test_ModeledDataset.py
r557 r569 11 11 12 12 # Build a initial dataset 13 from legacydata import datapoint, dataset13 from mystic.math.legacydata import datapoint, dataset 14 14 x = [[0,0,0],[1,1,1],[2,2,2],[3,3,3]] 15 15 y = [0,2,4,9] -
branches/UQ/math/legacy/test_StAlDataset.py
r557 r569 7 7 Mathematical Modeling and Numerical Analysis (submitted 2012). 8 8 """ 9 from legacydata import load_dataset9 from mystic.math.legacydata import load_dataset 10 10 datafile = 'StAlDataset.txt' 11 11 st_al_data = load_dataset(datafile) … … 21 21 22 22 # Check shortness 23 from paramtrans import lipschitz_distance, graphical_distance23 from mystic.math.paramtrans import lipschitz_distance, graphical_distance 24 24 print("\nshort: %s" % st_al_data.short()) 25 25 L = st_al_data.lipschitz; -
branches/UQ/math/sausage/TEST_OUQ_1dSurr_diam.py
r553 r569 22 22 # the model function and the dataset 23 23 ####################################################################### 24 #from StAlSurrogate import st_al_surr as model25 #from legacydata import load_dataset26 #data = load_dataset('ExampleDataset.txt', range(0,1))27 #L = data.lipschitz28 29 24 def model(x): 30 25 return x[0]; 26 31 27 32 28 ####################################################################### … … 76 72 def maximize(params,npts,bounds): 77 73 78 from dirac_measure import scenario74 from mystic.math.dirac_measure import scenario 79 75 from mystic.math import almostEqual 80 76 from numpy import inf … … 189 185 nz = 1 #NOTE: SET THE NUMBER OF 'v' POINTS HERE! 190 186 npts = (nx,ny,nz) 191 from paramtrans import _npts187 from mystic.math.paramtrans import _npts 192 188 _n = _npts(npts) 193 189 … … 252 248 253 249 from numpy import array 254 from dirac_measure import scenario250 from mystic.math.dirac_measure import scenario 255 251 c = scenario() 256 252 c.load(solved,npts) -
branches/UQ/math/sausage/TEST_OUQ_StStSurr.py
r555 r569 23 23 ####################################################################### 24 24 from StStSurrogate import marc_surr as model 25 #from legacydata import load_dataset26 #data = load_dataset('ExampleDataset.txt', range(0,1))27 #L = data.lipschitz28 25 29 26 … … 74 71 def maximize(params,npts,bounds): 75 72 76 from dirac_measure import scenario73 from mystic.math.dirac_measure import scenario 77 74 from mystic.math import almostEqual 78 75 from numpy import inf … … 200 197 nz = 1 #NOTE: SET THE NUMBER OF 'v' POINTS HERE! 201 198 npts = (nx,ny,nz) 202 from paramtrans import _npts199 from mystic.math.paramtrans import _npts 203 200 _n = _npts(npts) 204 201 … … 263 260 264 261 from numpy import array 265 from dirac_measure import scenario262 from mystic.math.dirac_measure import scenario 266 263 c = scenario() 267 264 c.load(solved,npts) -
mystic/_math/dirac_measure.py
r508 r569 2 2 """ 3 3 Classes for Dirac measure data objects. 4 Includes point, dirac_measure, and product_measureclasses.4 Includes point, dirac_measure, product_measure, and scenario classes. 5 5 """ 6 6 # Adapted from seesaw2d.py in branches/UQ/math/examples2/ … … 12 12 13 13 class 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 """ 15 25 16 26 def __init__(self, position, weight): … … 21 31 def __repr__(self): 22 32 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) 23 40 24 41 pass … … 72 89 def __mean(self): 73 90 from mystic.math.measures import mean 74 return mean(self.coords, self.weights) 91 return mean(self.coords, self.weights) 75 92 76 93 def __range(self): … … 174 191 methods: 175 192 c.pof(f) -- calculate the probability of failure 193 c.sampled_pof(f, npts) -- calculate the pof using sampled points 176 194 c.get_expect(f) -- calculate the expectation 177 195 c.set_expect((center,delta), f) -- impose expectation by adjusting positions 178 196 c.flatten() -- convert measure to a flat list of parameters 179 197 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 180 199 181 200 notes: … … 185 204 - weight wxi should be same for each (yj,zk) at xi; similarly for wyi & wzi 186 205 """ 206 def __val(self): 207 raise NotImplementedError, "'value' is undefined in a dirac_measure" 187 208 188 209 def __pts(self): … … 276 297 # u += self.weights[i] 277 298 # return u #XXX: does this need to be normalized? 278 299 279 300 def sampled_pof(self, f, npts=10000): 280 301 """calculate probability of failure over a given function, f, … … 307 328 from numpy import transpose 308 329 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 334 The 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 309 345 310 346 def load(self, params, pts): … … 327 363 pattern followed for each new dimension desired for the product measure. 328 364 """ 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) ) 332 370 return 333 371 … … 349 387 pattern followed for each dimension of the product measure. 350 388 """ 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 395 string 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 ? 352 423 353 424 #__center = None … … 368 439 369 440 441 class 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 515 Inputs: 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 522 Additional 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 527 Notes: 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 549 Inputs: 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 556 Additional Inputs: 557 tol -- maximum acceptable deviation from shortness 558 cutoff -- zero out distances less than cutoff; typically: tol, 0.0, or None 559 560 Notes: 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 583 Inputs: 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 591 Additional Inputs: 592 tol -- maximum acceptable deviation from shortness 593 cutoff -- zero out distances less than cutoff; typically cutoff = tol or 0.0 594 595 Notes: 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 618 Inputs: 619 data -- a collection of data points 620 cutoff -- acceptable deviation from shortness 621 622 Additional 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 639 Inputs: 640 model -- the model function, y' = F(x'), that approximates reality, y = G(x) 641 cutoff -- acceptable model invalidity |y - F(x')| 642 643 Additional 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 650 Notes: 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, 668 where f takes a list of (scenario) values and returns a single value 669 670 Inputs: 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 683 The 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 699 Inputs: 700 params -- a list of parameters (see 'notes') 701 pts -- number of points in each of the underlying discrete measures 702 703 Notes: 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 725 Returns: 726 params -- a list of parameters (see 'notes') 727 728 Notes: 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 370 748 #--------------------------------------------- 371 749 # creators and destructors from parameter list 372 750 373 751 def _mimic(samples, weights): 374 """Generate a product_measure object from a list N product measure752 """Generate a product_measure object from a list of N product measure 375 753 positions and a list of N weights. The resulting product measure will 376 754 mimic the original product measure's statistics, but be larger in size. … … 395 773 396 774 397 def _list_of_measures(samples, weights): 775 def _uniform_weights(samples): 776 """generate a nested list of N x 1D weights from a nested list of N x 1D 777 discrete 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 788 def _list_of_measures(samples, weights=None): 398 789 """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.""" 790 discrete measure positions and a nested list of N x 1D weights. 791 792 Note this function does not return a product measure, it returns a list.""" 400 793 total = [] 794 if not weights: weights = _uniform_weights(samples) 401 795 for i in range(len(samples)): 402 796 next = dirac_measure() … … 407 801 408 802 409 def compose(samples, weights ):803 def compose(samples, weights=None): 410 804 """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.""" 805 discrete measure positions and a nested list of N x 1D weights. If weights 806 are not provided, a uniform distribution with norm = 1.0 will be used.""" 807 if not weights: weights = _uniform_weights(samples) 412 808 total = _list_of_measures(samples, weights) 413 809 c = product_measure(total) … … 423 819 424 820 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 425 830 def unflatten(params, npts): 426 831 """Map a list of random variables to N x 1D discrete measures … … 431 836 432 837 838 from itertools import chain #XXX: faster, but sloppy to have as importable 433 839 def flatten(c): 434 840 """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] 439 842 # 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 ##### 847 def 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) 869 def 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 890 def 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 918 def impose_feasible(cutoff, data, guess=None, **kwds): 919 """impose shortness on a given list of parameters w,x,y. 920 921 Optimization on w,x,y over the given bounds seeks sum(infeasibility) = 0. 922 (this function is not ???-preserving) 923 924 Inputs: 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 930 Additional 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 936 Outputs: 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 1044 def impose_valid(cutoff, model, guess=None, **kwds): 1045 """impose model validity on a given list of parameters w,x,y 1046 1047 Optimization on w,x,y over the given bounds seeks sum(infeasibility) = 0. 1048 (this function is not ???-preserving) 1049 1050 Inputs: 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 1056 Additional 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 1063 Outputs: 1064 pm -- a scenario with desired model validity 1065 1066 Notes: 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 1174 if __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] 443 1235 444 1236 -
mystic/_math/legacydata.py
r558 r569 328 328 if L is None: L = self.lipschitz 329 329 if data is None: data = self 330 from paramtrans import lipschitz_distance, is_feasible330 from mystic.math.paramtrans import lipschitz_distance, is_feasible 331 331 # calculate the shortness 332 332 Rv = lipschitz_distance(L, self, data, **kwds) … … 372 372 elif cutoff is False: cutoff = None 373 373 374 from paramtrans import graphical_distance, is_feasible374 from mystic.math.paramtrans import graphical_distance, is_feasible 375 375 # calculate the model validity 376 376 Rv = graphical_distance(model, self, **kwds) -
mystic/_math/paramtrans.py
r558 r569 305 305 #L = list of lipschitz constants, for use when lipschitz metric is desired 306 306 #constraints = constraints function for finding minimum distance 307 from legacydata import dataset307 from mystic.math.legacydata import dataset 308 308 from numpy import asarray 309 309 from mystic.solvers import diffev2, fmin_powell … … 422 422 print "\nbuilding a scenario from the params..." 423 423 # [store Y as 'values' OR register(F) for Y=F(X) OR points store y as 'val' ?] 424 from dirac_measure import scenario424 from mystic.math.dirac_measure import scenario 425 425 pm = scenario() 426 426 pm.load(param1, pts) … … 437 437 # build a dataset (using X,Y) 438 438 # [store W as 'weights' ?] 439 from legacydata import dataset439 from mystic.math.legacydata import dataset 440 440 d = dataset() 441 441 d.load(X, Y) -
mystic/scripts/support_hypercube_scenario.py
r566 r569 266 266 # would be nice to use meta = ['wx','wx2','x','x2','wy',...] 267 267 268 from dirac_measure import scenario269 from legacydata import dataset268 from mystic.math.dirac_measure import scenario 269 from mystic.math.legacydata import dataset 270 270 try: # select whether to plot the cones 271 271 cones = parsed_opts.cones … … 280 280 try: # get the name of the dataset file 281 281 file = parsed_args[1] 282 from legacydata import load_dataset282 from mystic.math.legacydata import load_dataset 283 283 data = load_dataset(file) 284 284 except: -
mystic/setup.py
r530 r569 316 316 'scripts/support_convergence.py', 317 317 'scripts/support_hypercube.py', 318 'scripts/support_hypercube_measures.py']) 318 'scripts/support_hypercube_measures.py', 319 'scripts/support_hypercube_scenario.py']) 319 320 """ 320 321
Note: See TracChangeset
for help on using the changeset viewer.