File convex_hull.py changed (mode: 100644) (index 8e960e2..a83cd27) |
... |
... |
import numpy as np |
6 |
6 |
from . import sball |
from . import sball |
7 |
7 |
#from . import f_models |
#from . import f_models |
8 |
8 |
from scipy import stats |
from scipy import stats |
|
9 |
|
from scipy import spatial |
9 |
10 |
from .candybox import CandyBox |
from .candybox import CandyBox |
10 |
11 |
|
|
11 |
12 |
|
|
|
... |
... |
class GBall: |
49 |
50 |
def npoints(hull): |
def npoints(hull): |
50 |
51 |
return len(hull.sample) |
return len(hull.sample) |
51 |
52 |
|
|
|
53 |
|
@property |
|
54 |
|
def nsimplex(hull): |
|
55 |
|
return 1 |
|
56 |
|
|
52 |
57 |
#def update(hull): pass |
#def update(hull): pass |
53 |
58 |
|
|
54 |
59 |
def is_inside(hull, nodes): |
def is_inside(hull, nodes): |
|
... |
... |
class BrickHull: #č nebo BoundingBrick |
142 |
147 |
def npoints(hull): |
def npoints(hull): |
143 |
148 |
return len(hull.sample) |
return len(hull.sample) |
144 |
149 |
|
|
|
150 |
|
@property |
|
151 |
|
def nsimplex(hull): |
|
152 |
|
return hull.sample.nvar*2 |
|
153 |
|
|
145 |
154 |
@property |
@property |
146 |
155 |
def A(hull): |
def A(hull): |
147 |
156 |
hull._update() |
hull._update() |
|
... |
... |
class DirectHull: |
224 |
233 |
def npoints(hull): |
def npoints(hull): |
225 |
234 |
return len(hull.sample) |
return len(hull.sample) |
226 |
235 |
|
|
|
236 |
|
@property |
|
237 |
|
def nsimplex(hull): |
|
238 |
|
return len(hull.direct_plan) |
|
239 |
|
|
227 |
240 |
@property |
@property |
228 |
241 |
def A(hull): |
def A(hull): |
229 |
242 |
return hull.direct_plan |
return hull.direct_plan |
|
... |
... |
class CompleteHull: |
331 |
344 |
def npoints(hull): |
def npoints(hull): |
332 |
345 |
return len(hull.sample) |
return len(hull.sample) |
333 |
346 |
|
|
|
347 |
|
@property |
|
348 |
|
def nsimplex(hull): |
|
349 |
|
return 2 * (len(hull.direct_plan) + hull.sample.nvar) |
|
350 |
|
|
334 |
351 |
@property |
@property |
335 |
352 |
def A(hull): |
def A(hull): |
336 |
353 |
hull._update() |
hull._update() |
|
... |
... |
class CompleteHull: |
404 |
421 |
return hull.sample.f_model.new_sample(fire_G, space='G') |
return hull.sample.f_model.new_sample(fire_G, space='G') |
405 |
422 |
|
|
406 |
423 |
|
|
|
424 |
|
class QHull: |
|
425 |
|
def __init__(self, sample, space='G', incremental=True): |
|
426 |
|
self.sample = sample |
|
427 |
|
self.incremental = incremental |
|
428 |
|
self.space = space |
|
429 |
|
|
|
430 |
|
def regen(self): |
|
431 |
|
points = getattr(self.sample, self.space) |
|
432 |
|
self.convex_hull = spatial.ConvexHull(points, incremental=self.incremental) |
|
433 |
|
|
|
434 |
|
def enough_points(self): return self.sample.nvar < self.sample.nsim |
|
435 |
|
|
|
436 |
|
|
|
437 |
|
def __getattr__(self, attr): |
|
438 |
|
#č branime se rekurzii |
|
439 |
|
# defend against recursion |
|
440 |
|
#оӵ рекурсилы пезьдэт! |
|
441 |
|
if attr == 'convex_hull': |
|
442 |
|
#č zkusme rychle konvexní obálky sestavit |
|
443 |
|
#č a ihned ji vrátit |
|
444 |
|
self.regen() |
|
445 |
|
return self.convex_hull |
|
446 |
|
|
|
447 |
|
elif attr == 'A': |
|
448 |
|
return self.convex_hull.equations[:,:-1] |
|
449 |
|
elif attr == 'b': |
|
450 |
|
return self.convex_hull.equations[:,-1] |
|
451 |
|
|
|
452 |
|
elif attr == 'points': |
|
453 |
|
if self.enough_points(): |
|
454 |
|
return self.convex_hull.points |
|
455 |
|
else: |
|
456 |
|
return getattr(self.sample, self.space) |
|
457 |
|
|
|
458 |
|
elif attr == 'npoints': |
|
459 |
|
if self.enough_points(): |
|
460 |
|
return self.convex_hull.npoints |
|
461 |
|
else: |
|
462 |
|
return len(self.sample) |
|
463 |
|
|
|
464 |
|
elif attr == 'nsimplex': |
|
465 |
|
if self.enough_points(): |
|
466 |
|
return self.convex_hull.nsimplex |
|
467 |
|
else: |
|
468 |
|
return 0 |
|
469 |
|
|
|
470 |
|
|
|
471 |
|
#ё По всем вопросам обращайтесь |
|
472 |
|
#ё на нашу горячую линию |
|
473 |
|
else: |
|
474 |
|
self._update() #č dycky čerstý chleba! |
|
475 |
|
return getattr(self.convex_hull, attr) |
|
476 |
|
|
|
477 |
|
|
|
478 |
|
#č nejsem jist, jestli ten update vůbec dělat. |
|
479 |
|
#č lze navrhnout třidu aj tak, že sama bude hlídat změny. |
|
480 |
|
#č Jenomže, co kdybychom ten automatický update nechtěli? |
|
481 |
|
def _update(self): |
|
482 |
|
if self.convex_hull.npoints < self.sample.nsim: |
|
483 |
|
if self.incremental: |
|
484 |
|
points = getattr(self.sample, self.space) |
|
485 |
|
self.convex_hull.add_points(points[self.convex_hull.npoints:]) |
|
486 |
|
else: |
|
487 |
|
self.regen() |
|
488 |
|
|
|
489 |
|
|
|
490 |
|
def is_inside(self, nodes): |
|
491 |
|
if self.enough_points(): |
|
492 |
|
self._update() |
|
493 |
|
x = getattr(nodes, self.space) |
|
494 |
|
|
|
495 |
|
#E [normal, offset] forming the hyperplane equation of the facet (see Qhull documentation for more) |
|
496 |
|
A = self.convex_hull.equations[:,:-1] |
|
497 |
|
b = self.convex_hull.equations[:,-1] |
|
498 |
|
|
|
499 |
|
# N=ns, E - number of hyperplane equations |
|
500 |
|
ExN = A @ x.T + np.atleast_2d(b).T |
|
501 |
|
mask = np.all(ExN < 0, axis=0) |
|
502 |
|
return mask |
|
503 |
|
else: |
|
504 |
|
return np.full(len(nodes), False) |
|
505 |
|
|
|
506 |
|
def is_outside(hull, nodes): |
|
507 |
|
return ~hull.is_inside(nodes) |
|
508 |
|
|
|
509 |
|
def get_design_points(hull): |
|
510 |
|
hull._update() |
|
511 |
|
sample_model = -hull.A * hull.b.reshape(-1,1) |
|
512 |
|
return hull.sample.f_model.new_sample(sample_model, space=hull.space) |
|
513 |
|
|
|
514 |
|
|
|
515 |
|
# shortcut for Ghull |
|
516 |
|
# valid only if space==G |
|
517 |
|
def get_r(hull): |
|
518 |
|
if hull.space=='G': |
|
519 |
|
if hull.enough_points(): |
|
520 |
|
hull._update() |
|
521 |
|
b = hull.convex_hull.equations[:,-1] |
|
522 |
|
return -np.nanmax(b) |
|
523 |
|
else: |
|
524 |
|
return 0 |
|
525 |
|
else: |
|
526 |
|
return 0 |
|
527 |
|
|
|
528 |
|
def fire(hull, ns): |
|
529 |
|
if hull.space == 'G': |
|
530 |
|
A = hull.equations[:,:-1] |
|
531 |
|
b = hull.equations[:,-1] |
|
532 |
|
|
|
533 |
|
to_fire = np.argmax(b) |
|
534 |
|
a = A[to_fire] |
|
535 |
|
fire_from = stats.norm.cdf(hull.get_r()) |
|
536 |
|
t = np.linspace(fire_from, 1, ns) |
|
537 |
|
t = stats.norm.ppf(t) |
|
538 |
|
fire_G = t.reshape(-1,1) @ a.reshape(1,-1) |
|
539 |
|
|
|
540 |
|
return hull.sample.f_model.new_sample(fire_G, space='G') |
|
541 |
|
|
|
542 |
|
|
407 |
543 |
#č mým úkolem při návrhu této třidy je pořádně všecko zkomplikovat. |
#č mým úkolem při návrhu této třidy je pořádně všecko zkomplikovat. |
408 |
544 |
#č Dostávám za to peníze. |
#č Dostávám za to peníze. |
409 |
545 |
#č Takže. Udělám 3 druhů estimátorů |
#č Takže. Udělám 3 druhů estimátorů |
|
... |
... |
class Ghull: |
442 |
578 |
|
|
443 |
579 |
# shell_estimation -22: inner, -3: shell, -11: outer |
# shell_estimation -22: inner, -3: shell, -11: outer |
444 |
580 |
shell_estimation = {-22:shell.ps, -3: shell.p_shell, -11: shell.pf} |
shell_estimation = {-22:shell.ps, -3: shell.p_shell, -11: shell.pf} |
445 |
|
global_stats = {"r":r, "R":R, "inner":shell.ps, "shell":shell.p_shell, "outer":shell.pf} |
|
|
581 |
|
global_stats = {"nsim":self.sample.nsim, "nsimplex": self.hull.nsimplex, \ |
|
582 |
|
"r":r, "R":R, "inner":shell.ps, "shell":shell.p_shell, "outer":shell.pf} |
446 |
583 |
return shell_estimation, global_stats |
return shell_estimation, global_stats |
447 |
584 |
|
|
448 |
585 |
def integrate(self, nis): |
def integrate(self, nis): |