iam-git / WellMet (public) (License: MIT) (since 2021-08-31) (hash sha1)
WellMet is pure Python framework for spatial structural reliability analysis. Or, more specifically, for "failure probability estimation and detection of failure surfaces by adaptive sequential decomposition of the design domain".
Clone URLs: https://rocketgit.com/user/iam-git/WellMet ssh://rocketgit@ssh.rocketgit.com/user/iam-git/WellMet git://git.rocketgit.com/user/iam-git/WellMet
release trunk
List of commits:
Subject Hash Author Date (UTC)
convex_hull.Ghull: integrovaní memikruží je hodně překopano. Breakage. není zatím testováno 2d39dcd67948a99028013b0f7f420adf86eb19e9 I am 2021-07-04 22:02:48
dicebox: replace everywhere "potencial" typo by correct "potential" word. Save English eyes! c4aae036362d93c8a0a5bde9f30537d1b549adaa I am 2021-06-30 14:30:26
welford: add method for sparse data c9dd961c58243ac890b9a518cfff5901dd6a5789 I am 2021-06-28 13:47:09
welford: yet another Welford's algoritm implementation ec31907d6727c408306f4b25c09b7f87632b0ecf I am 2021-06-27 22:27:06
IS_stat.TrueIS: pridal jsem do komentářu vzpomínky 1c5b4596a4da1206b02ed91954cbe8147ec1e208 I am 2021-06-27 15:23:42
IS_stat: add get_IS_estimation() function d048b4cd4e90d8440e7ff2a3ad86bc0f3c160d4f I am 2021-06-26 19:55:14
spring: new brand hustý modul pro vyrovnání odhadů 3fc5eb00de1eb65cff2db71640c340602687f52f I am 2021-06-25 00:52:06
convex_hull.Ghull.rvs: hloupá chyba 7494eef4e2621ac9651037cd4a99e9f570704dee I am 2021-06-21 01:37:49
convex_hull.Ghull.rvs: ensure at least one node outside 1c6d05935c83d11637ae566528cf65eb384215b5 I am 2021-06-19 19:04:40
convex_hull.Ghull: add nonGaussian hulls support, fix MC integration 3e10c2cefff5f3e204c1d4118b195e37023a5a57 I am 2021-06-19 17:44:09
mart: add qhull polygon c3da6bdd1b6412f178c828c1eb1a6368a02eb74c I am 2021-06-18 23:27:48
dicebox.Goal: implement csv export odhadů abdb28f509e664c039a6e5ca4bf1036d00957349 I am 2021-06-15 15:06:15
reader: add export_estimation function. Jednoduchý jako busy ef8fcada4b01be8db29836aa09459d8a7ef382b6 I am 2021-06-15 03:08:53
convex_hull: add meaningless method 96a86b7eabceb15e6c937932c2002409c2df0178 I am 2021-06-14 16:35:39
dicebox: new brand box named Goal is ready 725c2561b62769af6034c690005d1d22e6af41c7 I am 2021-06-05 13:40:15
convex_hull: last node fix in fire() function ca938c11398dfb2e16441b989902dab1558336b3 I am 2021-06-05 13:25:55
convex_hull: small fix for QHull under single simplex integration b0b236cd14ae614ecd1310b8c7cfcafe8d1b0851 I am 2021-05-30 05:57:45
convex_hull: ensure Ghull to always have some nodes for expansion 7110d521d9fece58639a359c83773df904b1177f I am 2021-05-29 07:03:49
dicebox: throw away plot dependency, simplify increment function 1cf0c481b8df1937ff7eda6e0373d5d6a72387db I am 2021-05-28 16:42:19
simplex: add JustCubatureTriangulation class separated from Shull ea291f9deee9f3b5b8615ae1f14a49aecd327461 I am 2021-05-28 16:39:58
Commit 2d39dcd67948a99028013b0f7f420adf86eb19e9 - convex_hull.Ghull: integrovaní memikruží je hodně překopano. Breakage. není zatím testováno
Author: I am
Author date (UTC): 2021-07-04 22:02
Committer name: I am
Committer date (UTC): 2021-07-04 22:02
Parent(s): c4aae036362d93c8a0a5bde9f30537d1b549adaa
Signer:
Signing key:
Signing status: N
Tree: bfad58cbc508a4202ee7f68548b4202ab9825126
File Lines added Lines deleted
convex_hull.py 234 269
File convex_hull.py changed (mode: 100644) (index 4235391..fe06906)
... ... from . import sball
8 8 from scipy import stats from scipy import stats
9 9 from scipy import spatial from scipy import spatial
10 10 from .candybox import CandyBox from .candybox import CandyBox
11 from .IS_stat import PushAndPull # for Shell_IS
11 12
12 13
13 14 #č jako sůl potřebujem statefull třidu, #č jako sůl potřebujem statefull třidu,
 
... ... except BaseException as e:
572 573
573 574
574 575 class Ghull: class Ghull:
575 def __init__(self, hull, non_Gaussian_reduction=0.9, design=None):
576 def __init__(self, hull, use_MC=False, non_Gaussian_reduction=0.9):
576 577 self.hull = hull self.hull = hull
577 self.design = design
578 578 self.shell = sball.Shell(hull.sample.nvar) self.shell = sball.Shell(hull.sample.nvar)
579 579 self.outside_dist = sball.Shell(hull.sample.nvar) self.outside_dist = sball.Shell(hull.sample.nvar)
580 580 self.sample = hull.sample self.sample = hull.sample
581 581
582 self._r = 0 # estimation for nonGaussian distributions
583 self.non_Gaussian_reduction = non_Gaussian_reduction
584 self.integration_cutoff = np.inf
582 if use_MC:
583 self.gint = Shell_MC(hull, shell, non_Gaussian_reduction)
584 else:
585 self.gint = Shell_IS(hull, shell, non_Gaussian_reduction)
586
585 587
586 588 def fire(self, ns): def fire(self, ns):
587 589 nodes = self.hull.fire(ns) nodes = self.hull.fire(ns)
 
... ... class Ghull:
607 609 if self.hull.space == 'G': if self.hull.space == 'G':
608 610 return self.hull.get_r() return self.hull.get_r()
609 611 else: else:
610 return self._r * self.non_Gaussian_reduction
612 #č get_design_points() nemusí být. Ale zatím to neřeším
613 hull_points = self.hull.get_design_points()
614 sum_squared = np.sum(np.square(hull_points.G), axis=1)
615 hull_r = np.sqrt(np.nanmin(sum_squared))
616
617 # ask our experimentation team
618 gint_r = self.gint.get_r()
619
620 #č podle mně, nemůže se stat, že by metoda vratila:
621 #č 1. nenůlový poloměr a že by dokonce i centralní bod byl venku.
622 #č 2. poloměr větší jak R
623 #č Odpovědnost za to kladu na gint.
624 return min(hull_r, gint_r)
611 625
612 626 def get_shell_estimation(self): def get_shell_estimation(self):
613 627 shell = self.shell shell = self.shell
 
... ... class Ghull:
626 640 "inner":shell.ps, "shell":shell.p_shell, "outer":shell.pf} "inner":shell.ps, "shell":shell.p_shell, "outer":shell.pf}
627 641 return shell_estimation, global_stats return shell_estimation, global_stats
628 642
629 def integrate(self, nis):
630 candy_nodes, _nis, nf, shell_estimation, global_stats = self.check_it_out(nis)
631 assert _nis == nis
632
633 ns = nis - nf
643 #č nie
644 ##č pro mně je důležité, aby bylo možné rychle přepinat mezi metodama,
645 ##č aby bylo možné výsledky jednoduše porovnavat.
646 ##č Proto explicitně posílám priznak, jakou metodu použit.
647 def integrate(self, nis, callback_all=None, callback_outside=None):
648 #č no to teda disajn třidy je doopravdy hroznej
649 # .get_shell_estimation() will calculate radiuses and will update shell
650 shell_estimation, global_stats = self.get_shell_estimation()
634 651 # shell_estimation -22: inner, -3: shell, -11: outer # shell_estimation -22: inner, -3: shell, -11: outer
635 p_shell = shell_estimation.pop(-3)
636 shell_pf = nf/nis * p_shell
637 shell_ps = ns/nis * p_shell
638
652 shell_estimation.pop(-3)
639 653 # ghull_estimation -22: inner, -21: shell inside, -12: shell outside, -11: outer # ghull_estimation -22: inner, -21: shell inside, -12: shell outside, -11: outer
640 654 ghull_estimation = shell_estimation; del(shell_estimation) ghull_estimation = shell_estimation; del(shell_estimation)
655
656
657 shell_ps, shell_pf, stats = self.gint.integrate(nis, callback_all, callback_outside)
641 658 ghull_estimation[-21] = shell_ps ghull_estimation[-21] = shell_ps
642 659 ghull_estimation[-12] = shell_pf ghull_estimation[-12] = shell_pf
643 global_stats["shell_budget"] = nis
644 global_stats["shell_inside"] = shell_ps
645 global_stats["shell_outside"] = shell_pf
660 global_stats.update(stats)
646 661
647 662 # convex_hull_estimation -2: inside, -1: outside # convex_hull_estimation -2: inside, -1: outside
648 663 inside = shell_ps + self.shell.ps inside = shell_ps + self.shell.ps
 
... ... class Ghull:
651 666 global_stats["inside"] = inside global_stats["inside"] = inside
652 667 global_stats["outside"] = outside global_stats["outside"] = outside
653 668
654 return candy_nodes, ghull_estimation, convex_hull_estimation, global_stats
669 return ghull_estimation, convex_hull_estimation, global_stats
670
671
672
673
674 #č pomocná třída pro integrování
675 #č nejsem tedy jist, jestli je to nejlepší napad -
676 #č dělit Ghull na podtřídy, delegovat funkcionalitu
677 #č ale jinak Ghull se stavá hodně překomplikovaným.
678 #č nic lepšího mně nenapadá, přemyšlel jsem dlouho.
679 class Shell_MC:
680 def __init__(self, shell, hull, non_Gaussian_reduction=0.98):
681 self.shell = shell
682 self.hull = hull
683 self.nvar = hull.sample.nvar
684
685 self.integration_cutoff = np.inf
686 self.non_Gaussian_reduction = non_Gaussian_reduction
687
688 # test if r>0
689 self.r = 0
690 central_G = np.full(hull.sample.nvar, 0, dtype=np.int8)
691 self.DP = self.hull.sample.f_model.new_sample(central_G, space='G')
655 692
693 def reset(self): # clear
694 self.r_needed = (self.hull.space!='G')
656 695
696 self.nsampled = 0
697 self.nfailed = 0
698
699 @property
700 def DP_is_valid(self):
701 #č dalo by se kontrolovat změny obálky
702 #č a provadět kontrolu DP pouze po změně obálky.
703 #č Nejsem liný, jeden bodík prostě za to doopravdy nestoji.
704 mask = self.hull.is_outside(self.DP)
705 return np.any(mask)
706
707 #č na rozdil od rodičovské třídy
708 #č vrací odhad r na základě předchozích integrací
709 #č metoda je navržena tak, aby Shell_IS jú mohl zdědit.
710 def get_r(self):
711 #č pokud návrhový bod z minula je pořad venku,
712 #č tak použijeme redukci.
713 #č Bacha, metoda bude vracet nuly pro obálky v Gaussovskem prostoru!
714 if self.DP_is_valid:
715 return self.r * self.non_Gaussian_reduction
716 else:
717 return self.r
718
719 # stateless
657 720 def rvs(self, nsampled, seats, ns): def rvs(self, nsampled, seats, ns):
658 "Generování­ vzorků (kandidátů a integračních bodů)"
721 "Generování vzorků (kandidátů a integračních bodů)"
659 722 # rand_dir: prepare ns random directions on a unit d-sphere # rand_dir: prepare ns random directions on a unit d-sphere
660 rand_dir = sball.get_random_directions(seats, self.sample.nvar) #random directions
723 rand_dir = sball.get_random_directions(seats, self.nvar) #random directions
661 724
662 725 # generate sampling probabilites # generate sampling probabilites
663 726 left = (ns-nsampled) / ns left = (ns-nsampled) / ns
 
... ... class Ghull:
666 729 #č globálně jdeme r->R, localně ale R_i+ -> R_i- #č globálně jdeme r->R, localně ale R_i+ -> R_i-
667 730 #č menší poloměry musejí jít dřív - na to zavazano nonGaussian _r #č menší poloměry musejí jít dřív - na to zavazano nonGaussian _r
668 731 #č převracení lokalně umožní vždycky mít alespoň jeden bodík outside, #č převracení lokalně umožní vždycky mít alespoň jeden bodík outside,
669 #č což je taky velmi vhodné, vzhledem k tomu, že _r bere z outside bodíků
732 #č což je taky velmi vhodné vzhledem k tomu, že se _r bere z outside bodíků
670 733 p = np.linspace(right, left, seats, endpoint=False) # probabilities for the radius p = np.linspace(right, left, seats, endpoint=False) # probabilities for the radius
671 734
672 735 # convert probabilitites into random radii # convert probabilitites into random radii
 
... ... class Ghull:
677 740 sample_G = rand_dir*r[:,None] sample_G = rand_dir*r[:,None]
678 741
679 742 return sample_G return sample_G
743
744 # bus analogy
745 def _process_part(self, seats, nis, callback_all=None, callback_outside=None):
746 # boarding
747 nodes_G = self.rvs(self.nsampled, seats, nis)
748 nodes = self.sample.f_model.new_sample(nodes_G, space='G')
749 # who is back?
750 mask = self.hull.is_outside(nodes)
751 if self.r_needed and np.any(mask):
752 #č rvs má vzorkovat od měnšího poloměru k většímu.
753 #č to znamená, že první outside bude mít nejměnší r vůbec
754 self._r_check_out(nodes[mask])
755 self.r_needed = False
756
757 if callback_all is not None:
758 # -2 = 'inside' -1 = 'outside'
759 candy_nodes = CandyBox(nodes, event_id=mask-2, is_outside=mask)
760 callback_all(candy_nodes)
761
762 if (callback_outside is not None) and np.any(mask):
763 callback_outside(nodes[mask])
764
765 assert len(mask) == seats
766 #č nevím proč, ale kdysi mě to vyšlo rychlejc jak obyčejný součet
767 self.nfailed += len(mask[mask])
768 self.nsampled += len(mask)
769
770 #č metoda je navržena tak, aby Shell_IS jú mohl zdědit.
771 def _r_check_out(self, outside_nodes):
772 sum_squares = np.sum(np.square(outside_nodes.G), axis=1)
773 arg = np.nanargmax(sum_squares)
774 new_r = np.sqrt(sum_squares[arg])
775 if (not self.DP_is_valid) or (new_r < self.r):
776 self.DP = outside_nodes[arg]
777 self.r = new_r
680 778
681 def check_it_out(self, nis):
682 #č no to teda disajn třidy je doopravdy hroznej
683 # .get_shell_estimation() will calculate radiuses and will update shell
684 shell_estimation, global_stats = self.get_shell_estimation()
685
686 #č rvs má vzorkovat od měnšího poloměru k většímu.
687 #č to znamená, že když první outside bude mít nejměnší r vůbec
688 r_needed = (self.hull.space!='G')
689
779
780 #č metoda je navržena tak, aby Shell_IS jú mohl zdědit.
781 def integrate(self, nis, callback_all=None, callback_outside=None):
782 self.reset()
690 783
691 nsampled = 0
692 nfailed = 0
693 784 if self.hull.nsimplex == 0: if self.hull.nsimplex == 0:
694 785 bus = self.integration_cutoff bus = self.integration_cutoff
695 786 else: else:
696 787 bus = int(mem_bytes / self.hull.nsimplex / 8 / 3) + 1 bus = int(mem_bytes / self.hull.nsimplex / 8 / 3) + 1
697 while nsampled < nis:
788 while self.nsampled < nis:
698 789
699 seats = min(nis - nsampled, self.integration_cutoff, bus)
700 try:
701 nodes_G = self.rvs(nsampled, seats, nis)
702 nodes = self.sample.f_model.new_sample(nodes_G, space='G')
703 mask = self.hull.is_outside(nodes)
704 if r_needed and np.any(mask):
705 self._r = np.sqrt(np.nanmax(np.sum(np.square(nodes_G[mask]), axis=1)))
706 r_needed = False
707 # -2 = 'inside' -1 = 'outside'
708 candy_nodes = CandyBox(nodes, event_id=mask-2, is_outside=mask)
709
710 #č nevím proč, ale kdysi mě to vyšlo rychlejc jak obyčejný součet
711 nfailed += len(mask[mask])
712 nsampled += len(mask)
713 assert len(mask) == seats
790 seats = min(nis - self.nsampled, self.integration_cutoff, bus)
791 try:
792 self._process_part(seats, nis, callback_all, callback_outside)
714 793 except MemoryError as m: except MemoryError as m:
715 print("Ghull: memory error, %s sedaček" % seats, repr(m))
716 self.integration_cutoff = int(np.ceil(seats/3))
794 print(self.__class__.__name__ +":", "memory error, %s sedaček" % seats, repr(m))
795 self.integration_cutoff = int(np.ceil(seats/2))
796
797 assert nis == self.nsampled
798
799 return self._get_result()
800
801
802 def _get_result(self):
803
804 nis = self.nsampled
805 nf = self.nfailed
806 ns = nis - nf
807 p_shell = self.shell.p_shell
808 shell_pf = nf/nis * p_shell
809 shell_ps = ns/nis * p_shell
717 810
718 #č Candy bodů nemusí být nis!
719 return candy_nodes, nsampled, nfailed, shell_estimation, global_stats
811 stats = dict()
812 stats["shell_budget"] = nis
813 stats["shell_inside"] = shell_ps
814 stats["shell_outside"] = shell_pf
815
816 return shell_ps, shell_pf, stats
720 817
721 818
722 ##č napadlo mě zababáchnuť třidu, která by se sama starala o všem co se tyče
723 ##č vnější domény. Nešlo mě totíž to udělat jednou funkcí, bylo by velmi
724 ##č špatné z hlediska zodpovednosti kódu. Tak to všecko zabalíme to třidy
725 ##č a odebereme z už beztak přetíženého blackboxu část komplexity
726 ## keywords: ISSI, estimation, outside, ConvexHull, Sball, IS kolem středních hodnot
727 #class Shull: # issi_estimate_outside
728 # def __init__(sx, f, model_space, sampling_space='G', \
729 # powerset_correction=True, incremental=True, design=None):
730 # #č tím powerset_corretion je myšlena úplná soustava jevů,
731 # #č tj. vyrovnaní s použitím míry vnější i vnitřní
732 # #č powerset_correction=True přídá -2 (inside) jev do ISSI
733 #
734 # #č zde f-ko musí taky obsahovat vzorky!
735 # sx.f = f
736 # sx.model_space = model_space
737 #
738 # #č kašlu na to, pokud uživatel zadal nesmysl, tak
739 # #č sam převedu do nečěho smyslúplnějšího
740 # _dict = {'R':'Rn', 'aR':'aRn', 'P':'GK', 'aP':'aGK', 'U':'G', 'aU':'aG'}
741 # if sampling_space in _dict:
742 # sx.sampling_space = _dict[sampling_space]
743 # else:
744 # sx.sampling_space = sampling_space
745 #
746 # sx.design = design
747 #
748 # sampled_plan_model = getattr(f, model_space)
749 # #č žádná kontrola chyb - nechť to spadné, když má spadnout!
750 # sx.convex_hull = spatial.ConvexHull(sampled_plan_model, incremental=incremental)
751 #
752 # # current outside probability estimation
753 # sx.p_out = 0.5 # implicit value
754 # sx.sball = sball.Sball(f.nvar)
755 # sx.base_r = sx.sball.get_r(0.5) # we want in average 50/50 ratio
756 #
757 #
758 # # -1 = 'outside'
759 # # -2 = 'inside'
760 # #č založme ISSI
761 # sx.powerset_correction = powerset_correction
762 # #č potřebuji pro korektnost mít před integrací zadané jevy
763 # if powerset_correction:
764 # sx.oiss = IS_stat.ISSI([-1, -2])
765 # else:
766 # sx.oiss = IS_stat.ISSI([-1])
767 #
768 #
769 #
770 # def increment(sx, input_sample):
771 # #č sample by měl byt jíž převeden na f (v .add_sample()),
772 # #č zodpovidá za to volajicí kód!
773 # sx.convex_hull.add_points(getattr(input_sample, sx.model_space))
774 #
775 #
776 #
777 # # require
778 # def integrate(sx, nis):
779 # # getting rid of the old estimations
780 # sx.oiss.delete_event_data(-1)
781 # sx.oiss.events.append(-1) #č už tak trošku sahám do vnitřku cizí třidy
782 # if sx.powerset_correction:
783 # sx.oiss.delete_event_data(-2)
784 # sx.oiss.events.append(-2) #č a záse
785 #
786 # #č posunutí středu trošiňku porušuje předpoklady, za kterých
787 # #č Sball volí rozdělení, ale přečo se mi stavá,
788 # #č že ve více dimenzích Shull na začatku prostě nemůže
789 # #č trefit ten blbej... v podstatě simplex
790 # #č Těžiště, přesnějí, prostě nějaký střed můžeme najít dvěma způsoby:
791 # #č 1. mean(vertices.model_space).sampling_space
792 # #č 2. mean(vertices.sampling_space)
793 # #č myslím si, že ten první (a bez váh) je stabilnější
794 # #č (víme, jak střed vrcholů v nějakém prostoru může ani netrefit do
795 # #č simplexu v původním prostoru)
796 ## vertices_model = sx.convex_hull.points[sx.convex_hull.vertices]
797 ## barycenter_model = np.mean(vertices_model, axis=0)
798 ## if sx.model_space == sx.sampling_space:
799 ## sx.barycenter_sampling = barycenter_model
800 ## else:
801 ## barycenter = sx.f.new_sample(barycenter_model, sx.model_space)
802 ## sx.barycenter_sampling = np.squeeze(getattr(barycenter, sx.sampling_space))
803 #
804 # #č rozhodl jsem, že shull musí odhadovat outside, ne motat se centrem Brna.
805 # #č předpokladám, že uživatel zadal buď G, nebo aspoň Rn prostor
806 # #sx.barycenter_sampling = np.full(sx.f.nvar, 0, dtype=np.int8)
807 #
808 # # first step
809 # nodes = sx._sample_sball(nis)
810 # mask = nodes.is_outside
811 #
812 #
813 # cut_off_out = int(nis/3)
814 # cut_off_in = int(nis/3)
815 # #č robím cyklus dokud nesberu dostatečně teček.
816 # #č To je fakt nejrobustnější řešení, co mě napadá
817 # # while (number_of_out_nodes or number_of_nodes_inside is too_small)
818 # while (len(mask[mask]) < cut_off_out): # or (len(mask[~mask]) < cut_off_in):
819 # print(sx.__class__.__name__ + ":", "cut off! Outsides: %s, insides: %s, p_out=%s"\
820 # % (len(mask[mask]), len(mask[~mask]), sx.p_out))
821 # #č je třeba jenom sehnat dostatečně bodíků a utikat
822 # nodes.add_sample(sx._sample_sball(nis))
823 # mask = nodes.is_outside
824 #
825 # #č když neprovadíme výrovnání, tak ten vnitršek nachren nepotřebujem
826 # if (len(mask[~mask]) < cut_off_in) and sx.powerset_correction:
827 # print(sx.__class__.__name__ + ":", \
828 # "cut off inside (%s of %s needed)! Do simplex-based integration of convex hull"\
829 # % (len(mask[~mask]), cut_off_in))
830 #
831 # nodes.add_sample(sx._sample_simplex(nis))
832 #
833 #
834 # return nodes
835 #
836 #
837 #
838 # def _sample_simplex(sx, nis):
839 #
840 # #č f-ko sice musí odkazovat na aktuální f_model
841 # #č ale na druhou stranu normálně ._integrate_simplex()
842 # #č potřebujeme pouze jednou - hned při vytvaření
843 # vertices = sx.f[sx.convex_hull.vertices]
844 # nvar = vertices.nvar
845 #
846 # # IS_like uses .new_sample method, so vertices can not be a SampleBox object
847 # #
848 # #č IS_like těžiště počítá v sampling_space, ale asi mu to až tak nevadí
849 # #
850 # # already divided by nsim in variance formule
851 # # divide by /(nvar+1)/(nvar+2) from simplex inertia tensor solution
852 # # multiply by simplex_volume, but it looks like it shouldn't be here
853 # # for simplex: d = nvar+2
854 # #č sice to má nazev h_plan, ale nese rozdělení a hustoty v f-ku
855 # nodes = IS_stat.IS_like(vertices, sampling_space=sx.sampling_space, \
856 # nis=nis, d=nvar+2, design=sx.design)
857 #
858 # #č indikatorová funkce
859 # sx.is_outside(nodes)
860 #
861 # # for IS_stats
862 # #č zkoušel jsem zadavat celou serii - zhoršuje to odhady
863 # #č nemůžeme důvěrovat tomu, jak ten slepej simplex vidí svět
864 # weights = nodes.w[~nodes.is_outside]
865 # sx.oiss.add_single_event_data(weights, event=-2, nis=nis)
866 # # add_single_event_data() do not calculate estimations itself
867 # sx.oiss.get_estimations()
868 # return nodes
869 #
870 #
871 #
872 # def _sample_sball(sx, nis):
873 # nvar = sx.f.nvar
874 # sampling_r, sx.p_out = sx.sball.get_r_iteration(sx.p_out)
875 # #sampling_r = sx.sball.get_r(sx.p_out)
876 #
877 # #č asi tam bylo sampling_r/base_r, že?
878 # std_ball = sampling_r/sx.base_r
879 # #č chcu std=1, když p_out -> 1
880 # #č a std=sball_solušn, když p_out -> 0
881 # #č surovější nevymyslíš!
882 # std = std_ball + sx.p_out
883 # #č u stats.norm zadáváme směrodatnou odchylku, je to asi správné
884 # #h = f_models.UnCorD([stats.norm(sx.barycenter_sampling[i], std) for i in range(nvar)])
885 # #nodes = IS_stat.IS(sx.f, h, space_from_h='R', space_to_f=sx.sampling_space, Nsim=nis)
886 #
887 # #norm_params = [(sx.barycenter_sampling[i], std) for i in range(nvar)]
888 # nodes = IS_stat.IS_norm(sx.f, std=std, \
889 # sampling_space=sx.sampling_space, nis=nis, design=sx.design)
890 #
891 # outside_measure = sx._apply_nodes(nodes, nis)
892 #
893 # #č pro přiště
894 # sx.p_out = (sx.p_out + outside_measure) / 2
895 #
896 # return nodes
897 #
898 #
899 # def _apply_nodes(sx, nodes, nis):
900 # #č indikatorová funkce
901 # sx.is_outside(nodes)
902 #
903 # # for IS_stats
904 # if sx.powerset_correction:
905 # #č získáme výrovnaný odhad - je to skoro zdarma
906 # #svar = (sampling_r/sx.base_r)**2 # svar like sampling_variance
907 # #č kdysi snažil jsem něco odvést, moc se mi to nepovedlo
908 # #č je to jen tak, jeden z pokusu, hrubej nastřel
909 # #im = svar**nvar * np.exp(nvar/svar - nvar)
910 #
911 # #č radší ne. IM špatně zachycuje nizkou důvěru k tomu, co nemá vlastní tečky
912 # sx.oiss.add_IS_serie(nodes.w, nodes.event_id, implicit_multiplicator=np.inf)
913 # outside_measure = sx.oiss.estimations[-1]
914 # else:
915 # weights = nodes.w[nodes.is_outside]
916 # #č IM všecko pokazí, jakmile začnu přídávát další jevy
917 # sx.oiss.add_single_event_data(weights, event=-1, nis=nis)
918 # # add_single_event_data() do not calculate estimations itself
919 # weighted_mean, __, events = sx.oiss.get_means()
920 # # outside probability
921 # #č nevyrovnané!
922 # outside_measure = weighted_mean[events==-1][0]
923 #
924 # return outside_measure
925 #
926 #
927 #
928 # def is_outside(sx, nodes):
929 # node_coordinates = getattr(nodes, sx.model_space)
930 # mask = is_outside(sx.convex_hull, node_coordinates)
931 #
932 # # -1 = 'outside'
933 # # -2 = 'inside'
934 # nodes.event_id = mask - 2
935 # nodes.is_outside = mask
936 #
937 #
938 819
820 class Shell_IS(Shell_MC):
821 def reset(self): # clear
822 self.r_needed = (self.hull.space!='G')
823
824 self.nsampled = 0
825 self.pp = PushAndPull()
826
827
828 # stateless
829 def rvs(self, nsampled, seats, ns):
830 "Generování vzorků (kandidátů a integračních bodů)"
831 # rand_dir: prepare ns random directions on a unit d-sphere
832 rand_dir = sball.get_random_directions(seats, self.nvar) #random directions
833
834 #č poloměry bereme ze skořapky
835 #č za správné nastavení (stejně jako u MC)
836 #č zodpovidá uživatel třídy
837 #č třída vlastní odhad r nijak nevyuživá!
838 r = self.shell.r
839 R = self.shell.R
840 delta_r = R - r
841 left = nsampled / ns * delta_r + r
842 right = (nsampled + seats) / ns * delta_r + r
843 #č přidáme trochu zmatku.
844 #č globálně jdeme r->R, localně ale R_i+ -> R_i-
845 #č menší poloměry musejí jít dřív - na to zavazano nonGaussian _r
846 #č převracení lokalně umožní vždycky mít alespoň jeden bodík outside,
847 #č což je taky velmi vhodné vzhledem k tomu, že se _r bere z outside bodíků
848 rs = np.linspace(right, left, seats, endpoint=False)
849
850 #finally a random sample from the optimal IS density:
851 sample_G = rand_dir*rs[:,None]
852
853 #č a jsme doma. Platba za špatný design.
854 #č Potřebujeme hustotu 1D rozdělení, implementovanou v Radial
855 #č Shell se síce dědí od Radial, ale implementuje .pdf() jako sdruženou
856 #č hustotu v nD Gaussovskem prostoru
857 weights = sball.Radial.pdf(self.shell, rs) * delta_r
858
859 return sample_G, weights
939 860
861 # bus analogy
862 def _process_part(self, seats, nis, callback_all=None, callback_outside=None):
863 # boarding
864 nodes_G, weights = self.rvs(self.nsampled, seats, nis)
865 nodes = self.sample.f_model.new_sample(nodes_G, space='G')
866 # who is back?
867 mask = self.hull.is_outside(nodes)
868 if self.r_needed and np.any(mask):
869 #č rvs má vzorkovat od měnšího poloměru k většímu.
870 #č to znamená, že první outside bude mít nejměnší r vůbec
871 self._r_check_out(nodes[mask])
872 self.r_needed = False
873
874 assert len(mask) == seats
875 self.pp.add(weights, mask)
876 self.nsampled += len(mask)
877
878 if callback_all is not None:
879 # -2 = 'inside' -1 = 'outside'
880 candy_nodes = CandyBox(nodes, event_id=mask-2, is_outside=mask, weights=weights)
881 callback_all(candy_nodes)
882
883 if (callback_outside is not None) and np.any(mask):
884 callback_outside(nodes[mask])
885
886
887 def _get_result(self):
888
889 nis = self.nsampled
890 #č nejdřív true, pak false. Posilali jsme is_outside
891 mean_pf, mean_ps = self.pp.mean
892 var_pf, var_ps = self.pp.var
893 shell_pf, shell_ps = self.pp.corrected_mean(p_overall=self.shell.p_shell)
894
895 stats = dict()
896 stats["shell_budget"] = nis
897 stats["shell_inside_measured"] = mean_ps
898 stats["shell_outside_measured"] = mean_pf
899 stats["shell_inside_var"] = var_ps
900 stats["shell_outside_var"] = var_pf
901 stats["shell_inside"] = shell_ps
902 stats["shell_outside"] = shell_pf
903
904 return shell_ps, shell_pf, stats
Hints:
Before first commit, do not forget to setup your git environment:
git config --global user.name "your_name_here"
git config --global user.email "your@email_here"

Clone this repository using HTTP(S):
git clone https://rocketgit.com/user/iam-git/WellMet

Clone this repository using ssh (do not forget to upload a key first):
git clone ssh://rocketgit@ssh.rocketgit.com/user/iam-git/WellMet

Clone this repository using git:
git clone git://git.rocketgit.com/user/iam-git/WellMet

You are allowed to anonymously push to this repository.
This means that your pushed commits will automatically be transformed into a merge request:
... clone the repository ...
... make some changes and some commits ...
git push origin main