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)
ghull: 1DS integration scheme has been implemented. convex_hull has been divided into convex_hull itself and separated ghull module. 29b21b4fc1396384bf8c9bb21b91f82005850574 Aleksei Gerasimov 2021-07-20 13:26:47
IS_stat: add 1D IS-like helpful sampling function 58416707171522e8129333e8bd93f512c7578879 Aleksei Gerasimov 2021-07-20 11:06:27
convex_hull: orth estimation is ready. Uses QR decomposition. 75f20d8771021b037c360df9c7207b8f46f65d91 Aleksei Gerasimov 2021-07-20 08:25:20
convex_hull: 2FORM za mně je Ready. Orth potřebuje Gram-Schmidtovou ortogonalizaci. Oboje jsou otestovany ve 2D. There is some weirdness with candidates inherited from WIP recent commit. It need to be fixed. ffae78d5a15884476760b8ab759963d7ef8ff5fc Aleksei Gerasimov 2021-07-19 14:48:34
convex_hull: 2FORM estimation WIP a0f04ef282d6dacbda52c9734f36f9e7ec9b46da I am 2021-07-19 06:14:54
convex_hull.fire: use sf instead of cdf function to fix precision issue in distant areas. e3b9f036f9d50994357e32827586fb98d85d43a7 I am 2021-07-15 09:13:58
qt_plot: add matplotlib menu 73835a59908b82f08d618a06c4c05a96703715b1 I am 2021-07-14 17:01:54
mart: add candidates plot functions 56b3f5320e83071886c0d8f65a3d9df7c5a231a8 I am 2021-07-13 10:55:13
dicebox.Goal.assess_candidates: put back dd attribute 5da7c67800b8cd6eb3f221273a8270393a571e2b I am 2021-07-12 08:47:00
sball: fix Radial pdf 830fbf3119e2e54e42ea48640dfe553ffa2b3c1f I am 2021-07-08 19:50:42
spring: try to fix infinite loop 0e9a6161996c1b078698b1631d10f1a171332818 I am 2021-07-08 17:54:04
dicebox.Goal: add q_psee potential 4c151252460f6e96f7c1ac5662f81474eae305a0 I am 2021-07-08 00:34:27
dicebox.Goal: Goal is back after the breakage! 19c3c63001ebc19ed5e17892e4b11edb4dbfe828 I am 2021-07-06 16:14:14
IS_stat.PushAndPull: small variance fix, spíše formální 7b5b87655b6978828ed06ba3743ada46a2bfde6b I am 2021-07-06 04:31:43
qt_plot.Ghull: update Ghull related code. A lot of fixes. Вроде фурычит. 29cd8fce14cecc872e1be856d73d9669a6395b08 I am 2021-07-05 18:47:53
IS_stat: new class added (for Ghull, actually) 27b2685694655cb75657c6629d10706669d1f2e2 I am 2021-07-04 22:05:48
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
Commit 29b21b4fc1396384bf8c9bb21b91f82005850574 - ghull: 1DS integration scheme has been implemented. convex_hull has been divided into convex_hull itself and separated ghull module.
Author: Aleksei Gerasimov
Author date (UTC): 2021-07-20 13:26
Committer name: Aleksei Gerasimov
Committer date (UTC): 2021-07-20 13:26
Parent(s): 58416707171522e8129333e8bd93f512c7578879
Signer:
Signing key:
Signing status: N
Tree: e2cb3fc0c33e56b02d133ee2aeb55a480dbec997
File Lines added Lines deleted
convex_hull.py 2 408
dicebox.py 3 2
ghull.py 531 0
qt_plot.py 13 3
File convex_hull.py changed (mode: 100644) (index aee66a8..4b78c5d)
3 3
4 4 import numpy as np import numpy as np
5 5 import mpmath import mpmath
6 from . import sball
7 6 from scipy import stats from scipy import stats
8 from scipy import spatial
9 from .candybox import CandyBox
10 from .IS_stat import PushAndPull # for Shell_IS
7 from scipy import spatial # for QHull
11 8
12 mpmath.mp.dps = 325 # to cover anything double-presigion float can handle
9 mpmath.mp.dps = 325 # to cover anything that double-presigion float can handle
13 10
14 11 #č jako sůl potřebujem statefull třidu, #č jako sůl potřebujem statefull třidu,
15 12 #č která by schovávala vnitřní implementaciju, #č která by schovávala vnitřní implementaciju,
 
... ... class QHull:
803 800
804 801
805 802
806 #č mým úkolem při návrhu této třidy je pořádně všecko zkomplikovat.
807 #č Dostávám za to peníze.
808 #č Takže. Udělám 3 druhů estimátorů
809 # convex_hull_estimation -2: inside, -1: outside
810 # shell_estimation -22: inner, -3: shell, -11: outer
811 # ghull_estimation -22: inner, -21: shell inside, -12: shell outside, -11: outer
812
813 try: # try to outthink f****** OS's memory management
814 import os
815 mem_bytes = os.sysconf('SC_PAGE_SIZE') * os.sysconf('SC_PHYS_PAGES')
816 except BaseException as e:
817 mem_GB = 16
818 mem_bytes = mem_GB * 1024**3 # hello, Windows!
819 print("convex_hull failed to get amount of RAM installed. %s GB is supposed."% mem_GB, repr(e))
820 print("BTW, you are using Windows, aren't you?")
821
822
823 class Ghull:
824 def __init__(self, hull, calculate_orth=True, calculate_2FORM=True, use_MC=False, non_Gaussian_reduction=0.9):
825 self.hull = hull
826 self.shell = sball.Shell(hull.sample.nvar)
827 self.outside_dist = sball.Shell(hull.sample.nvar)
828 self.sample = hull.sample
829
830 #č vlastně nevím, ale nemusejí byť úplně zadarmo...
831 self.calculate_orth = calculate_orth
832 self.calculate_2FORM = calculate_2FORM
833
834 self.use_MC = use_MC
835
836 if use_MC:
837 self.gint = Shell_MC(hull, self.shell, non_Gaussian_reduction)
838 else:
839 self.gint = Shell_IS(hull, self.shell, non_Gaussian_reduction)
840
841
842 def fire(self, ns, use_MC=None):
843 nodes = self.hull.fire(ns, use_MC=self.use_MC)
844 if nodes is not None:
845 return nodes
846 else:
847 return self.boom(ns)
848
849 def boom(self, ns):
850 if self.use_MC:
851 self.outside_dist.set_bounds(self.get_R())
852 nodes_G = self.outside_dist.rvs(ns)
853 else:
854 # rand_dir: prepare ns random directions on a unit d-sphere
855 rand_dir = sball.get_random_directions(ns, self.sample.nvar) #random directions
856
857 # maximum radius, where norm.pdf() wasn't zero
858 # -38.575500173381374935388521407730877399444580
859 # don't ask me what the magic python use to distinguish
860 # digits after double precision
861 max_R_ever = -38.5755
862
863 r = np.linspace(self.get_R(), max_R_ever, ns, endpoint=True)
864 nodes_G = rand_dir*r[:,None]
865
866 nodes = self.sample.f_model.new_sample(nodes_G, space='G')
867 return nodes
868
869 def get_orth_outside(self):
870 if self.hull.space == 'G':
871 return self.hull.get_orth_outside()
872 else:
873 #č bude to horší jak s-ball, ale budiž.
874 x = np.full(self.sample.nvar, self.get_R())
875 return calculate_brick_complement_probability(-x, x)
876
877 def get_2FORM_outside(self):
878 if self.hull.space == 'G':
879 return self.hull.get_2FORM_outside()
880 else:
881 #č nebudem do toho michat nic dalšího.
882 #č Když 2FORM, tak 2FORM!
883 return 2 * stats.norm.sf(self.get_R())
884
885 def get_FORM_outside(self):
886 if self.hull.space == 'G':
887 return stats.norm.sf(self.get_r())
888 else:
889 return stats.norm.sf(self.get_R())
890
891 def get_R(self):
892 sum_squared = np.sum(np.square(self.sample.G), axis=1)
893 #index = np.argmax(sum_squared)
894 return np.sqrt(np.nanmax(sum_squared))
895
896 def get_r(self):
897 "calculates minimal signed distance to planes. Can therefore be negative"
898 #return -np.nanmax(self.hull.b)
899 if self.hull.space == 'G':
900 return self.hull.get_r()
901 else:
902 #č get_design_points() nemusí být.
903 #č QHull může nemít dostatek teček
904 try:
905 hull_points = self.hull.get_design_points()
906 sum_squared = np.sum(np.square(hull_points.G), axis=1)
907 hull_r = np.sqrt(np.nanmin(sum_squared))
908 except BaseException as e:
909 msg = "cannot get hull points from the convex hull"
910 print(self.__class__.__name__ +":", msg, repr(e))
911 hull_r = np.inf
912
913
914 # ask our experimentation team
915 gint_r = self.gint.get_r()
916
917 #č podle mně, nemůže se stat, že by metoda vratila:
918 #č 1. nenůlový poloměr a že by dokonce i centralní bod byl venku.
919 #č 2. poloměr větší jak R
920 #č Odpovědnost za to kladu na gint.
921 return min(hull_r, gint_r)
922
923 def get_shell_estimation(self):
924 shell = self.shell
925 r = self.get_r()
926 R = self.get_R()
927
928 if r<0:
929 shell.set_bounds(0, R)
930 else:
931 shell.set_bounds(r, R)
932
933 # shell_estimation -22: inner, -3: shell, -11: outer
934 shell_estimation = {-22:shell.ps, -3: shell.p_shell, -11: shell.pf}
935 global_stats = {"nsim":self.sample.nsim, "ndim":self.sample.nvar, \
936 "nfacets": self.hull.nsimplex, "r":r, "R":R, \
937 "inner":shell.ps, "shell":shell.p_shell, "outer":shell.pf}
938 global_stats['FORM_outside'] = self.get_FORM_outside()
939 if self.calculate_2FORM:
940 global_stats['2FORM_outside'] = self.get_2FORM_outside()
941 if self.calculate_orth:
942 global_stats['orth_outside'] = self.get_orth_outside()
943 return shell_estimation, global_stats
944
945 #č nie
946 ##č pro mně je důležité, aby bylo možné rychle přepinat mezi metodama,
947 ##č aby bylo možné výsledky jednoduše porovnavat.
948 ##č Proto explicitně posílám priznak, jakou metodu použit.
949 def integrate(self, nis, callback_all=None, callback_outside=None):
950 #č no to teda disajn třidy je doopravdy hroznej
951 # .get_shell_estimation() will calculate radiuses and will update shell
952 shell_estimation, global_stats = self.get_shell_estimation()
953 # shell_estimation -22: inner, -3: shell, -11: outer
954 shell_estimation.pop(-3)
955 # ghull_estimation -22: inner, -21: shell inside, -12: shell outside, -11: outer
956 ghull_estimation = shell_estimation; del(shell_estimation)
957
958
959 shell_ps, shell_pf, stats = self.gint.integrate(nis, callback_all, callback_outside)
960 ghull_estimation[-21] = shell_ps
961 ghull_estimation[-12] = shell_pf
962 global_stats.update(stats)
963
964 # convex_hull_estimation -2: inside, -1: outside
965 inside = shell_ps + self.shell.ps
966 outside = shell_pf + self.shell.pf
967 convex_hull_estimation = {-2: inside, -1: outside}
968 global_stats["inside"] = inside
969 global_stats["outside"] = outside
970
971 return ghull_estimation, convex_hull_estimation, global_stats
972
973
974
975
976 #č pomocná třída pro integrování
977 #č nejsem tedy jist, jestli je to nejlepší napad -
978 #č dělit Ghull na podtřídy, delegovat funkcionalitu
979 #č ale jinak Ghull se stavá hodně překomplikovaným.
980 #č nic lepšího mně nenapadá, přemyšlel jsem dlouho.
981 class Shell_MC:
982 def __init__(self, hull, shell, non_Gaussian_reduction=0.98):
983 self.shell = shell
984 self.hull = hull
985 self.nvar = hull.sample.nvar
986
987 self.integration_cutoff = np.inf
988 self.non_Gaussian_reduction = non_Gaussian_reduction
989
990 # test if r>0
991 self.r = 0
992 central_G = np.full(hull.sample.nvar, 0, dtype=np.int8)
993 self.DP = self.hull.sample.f_model.new_sample(central_G, space='G')
994
995 def reset(self): # clear
996 self.r_needed = (self.hull.space!='G')
997
998 self.nsampled = 0
999 self.nfailed = 0
1000
1001 @property
1002 def DP_is_valid(self):
1003 #č dalo by se kontrolovat změny obálky
1004 #č a provadět kontrolu DP pouze po změně obálky.
1005 #č Nejsem liný, jeden bodík prostě za to doopravdy nestoji.
1006 mask = self.hull.is_outside(self.DP)
1007 return np.any(mask)
1008
1009 #č na rozdil od rodičovské třídy
1010 #č vrací odhad r na základě předchozích integrací
1011 #č metoda je navržena tak, aby Shell_IS jú mohl zdědit.
1012 def get_r(self):
1013 #č bojim sa, rerukce bude aplikována dycky
1014 #č Bacha, metoda bude vracet nuly pro obálky v Gaussovskem prostoru!
1015 return self.r * self.non_Gaussian_reduction
1016
1017 # stateless
1018 def rvs(self, nsampled, seats, ns):
1019 "Generování vzorků (kandidátů a integračních bodů)"
1020 # rand_dir: prepare ns random directions on a unit d-sphere
1021 rand_dir = sball.get_random_directions(seats, self.nvar) #random directions
1022
1023 # generate sampling probabilites
1024 left = (ns-nsampled) / ns
1025 right = (ns-nsampled-seats) / ns
1026 #č přidáme trochu zmatku.
1027 #č globálně jdeme r->R, localně ale R_i+ -> R_i-
1028 #č menší poloměry musejí jít dřív - na to zavazano nonGaussian _r
1029 #č převracení lokalně umožní vždycky mít alespoň jeden bodík outside,
1030 #č což je taky velmi vhodné vzhledem k tomu, že se _r bere z outside bodíků
1031 p = np.linspace(right, left, seats, endpoint=False) # probabilities for the radius
1032
1033 # convert probabilitites into random radii
1034 # (distances from origin that are greater than r and less than R)
1035 r = self.shell.isf(p) # actually, it is the same as CDF inverse
1036
1037 #finally a random sample from the optimal IS density:
1038 sample_G = rand_dir*r[:,None]
1039
1040 return sample_G
1041
1042 # bus analogy
1043 def _process_part(self, seats, nis, callback_all=None, callback_outside=None):
1044 # boarding
1045 nodes_G = self.rvs(self.nsampled, seats, nis)
1046 nodes = self.hull.sample.f_model.new_sample(nodes_G, space='G')
1047 # who is back?
1048 mask = self.hull.is_outside(nodes)
1049 if self.r_needed and np.any(mask):
1050 #č rvs má vzorkovat od měnšího poloměru k většímu.
1051 #č to znamená, že první outside bude mít nejměnší r vůbec
1052 self._r_check_out(nodes[mask])
1053 self.r_needed = False
1054
1055 if callback_all is not None:
1056 # -2 = 'inside' -1 = 'outside'
1057 candy_nodes = CandyBox(nodes, event_id=mask-2, is_outside=mask)
1058 callback_all(candy_nodes)
1059
1060 if (callback_outside is not None) and np.any(mask):
1061 callback_outside(nodes[mask])
1062
1063 assert len(mask) == seats
1064 #č nevím proč, ale kdysi mě to vyšlo rychlejc jak obyčejný součet
1065 self.nfailed += len(mask[mask])
1066 self.nsampled += len(mask)
1067
1068 #č metoda je navržena tak, aby Shell_IS jú mohl zdědit.
1069 def _r_check_out(self, outside_nodes):
1070 sum_squares = np.sum(np.square(outside_nodes.G), axis=1)
1071 arg = np.nanargmin(sum_squares)
1072 new_r = np.sqrt(sum_squares[arg])
1073 if (not self.DP_is_valid) or (new_r < self.r):
1074 self.DP = outside_nodes[arg]
1075 self.r = new_r
1076
1077
1078 #č metoda je navržena tak, aby Shell_IS jú mohl zdědit.
1079 def integrate(self, nis, callback_all=None, callback_outside=None):
1080 self.reset()
1081
1082 if self.hull.nsimplex == 0:
1083 bus = self.integration_cutoff
1084 else:
1085 bus = int(mem_bytes / self.hull.nsimplex / 8 / 3) + 1
1086 while self.nsampled < nis:
1087
1088 seats = min(nis - self.nsampled, self.integration_cutoff, bus)
1089 try:
1090 self._process_part(seats, nis, callback_all, callback_outside)
1091 except MemoryError as m:
1092 print(self.__class__.__name__ +":", "memory error, %s sedaček" % seats, repr(m))
1093 self.integration_cutoff = int(np.ceil(seats/2))
1094
1095 assert nis == self.nsampled
1096
1097 return self._get_result()
1098
1099
1100 def _get_result(self):
1101
1102 nis = self.nsampled
1103 nf = self.nfailed
1104 ns = nis - nf
1105 p_shell = self.shell.p_shell
1106 shell_pf = nf/nis * p_shell
1107 shell_ps = ns/nis * p_shell
1108
1109 stats = dict()
1110 stats["shell_budget"] = nis
1111 stats["shell_inside"] = shell_ps
1112 stats["shell_outside"] = shell_pf
1113
1114 return shell_ps, shell_pf, stats
1115
1116
1117
1118 class Shell_IS(Shell_MC):
1119 def reset(self): # clear
1120 self.r_needed = (self.hull.space!='G')
1121
1122 self.nsampled = 0
1123 self.pp = PushAndPull()
1124
1125
1126 # stateless
1127 def rvs(self, nsampled, seats, ns):
1128 "Generování vzorků (kandidátů a integračních bodů)"
1129 # rand_dir: prepare ns random directions on a unit d-sphere
1130 rand_dir = sball.get_random_directions(seats, self.nvar) #random directions
1131
1132 #č poloměry bereme ze skořapky
1133 #č za správné nastavení (stejně jako u MC)
1134 #č zodpovidá uživatel třídy
1135 #č třída vlastní odhad r nijak nevyuživá!
1136 r = self.shell.r
1137 R = self.shell.R
1138 delta_r = R - r
1139 left = nsampled / ns * delta_r + r
1140 right = (nsampled + seats) / ns * delta_r + r
1141 #č přidáme trochu zmatku.
1142 #č globálně jdeme r->R, localně ale R_i+ -> R_i-
1143 #č menší poloměry musejí jít dřív - na to zavazano nonGaussian _r
1144 #č převracení lokalně umožní vždycky mít alespoň jeden bodík outside,
1145 #č což je taky velmi vhodné vzhledem k tomu, že se _r bere z outside bodíků
1146 rs = np.linspace(right, left, seats, endpoint=False)
1147
1148 #finally a random sample from the optimal IS density:
1149 sample_G = rand_dir*rs[:,None]
1150
1151 #č a jsme doma. Platba za špatný design.
1152 #č Potřebujeme hustotu 1D rozdělení, implementovanou v Radial
1153 #č Shell se síce dědí od Radial, ale implementuje .pdf() jako sdruženou
1154 #č hustotu v nD Gaussovskem prostoru
1155 #
1156 #č vahy jsou definovány jako podil původního rozdělení k vzorkovácímu
1157 #č původní - Chi rozdělení. nemá zde být žádný p_shell
1158 f_pdf = stats.chi.pdf(rs, self.nvar)
1159 #č vzorkovácí.
1160 #h_pdf = 1/delta_r
1161 weights = f_pdf * delta_r
1162
1163 return sample_G, weights
1164
1165 # bus analogy
1166 def _process_part(self, seats, nis, callback_all=None, callback_outside=None):
1167 # boarding
1168 nodes_G, weights = self.rvs(self.nsampled, seats, nis)
1169 nodes = self.hull.sample.f_model.new_sample(nodes_G, space='G')
1170 # who is back?
1171 mask = self.hull.is_outside(nodes)
1172 if self.r_needed and np.any(mask):
1173 #č rvs má vzorkovat od měnšího poloměru k většímu.
1174 #č to znamená, že první outside bude mít nejměnší r vůbec
1175 self._r_check_out(nodes[mask])
1176 self.r_needed = False
1177
1178 assert len(mask) == seats
1179 self.pp.add(weights, mask)
1180 self.nsampled += len(mask)
1181
1182 if callback_all is not None:
1183 # -2 = 'inside' -1 = 'outside'
1184 candy_nodes = CandyBox(nodes, event_id=mask-2, is_outside=mask, weights=weights)
1185 callback_all(candy_nodes)
1186
1187 if (callback_outside is not None) and np.any(mask):
1188 callback_outside(nodes[mask])
1189
1190
1191 def _get_result(self):
1192
1193 nis = self.nsampled
1194 #č nejdřív true, pak false. Posilali jsme is_outside
1195 mean_pf, mean_ps = self.pp.mean
1196 var_pf, var_ps = self.pp.var
1197 shell_pf, shell_ps = self.pp.correct_means(p_overall=self.shell.p_shell)
1198
1199 stats = dict()
1200 stats["shell_budget"] = nis
1201 stats["shell_inside_measured"] = mean_ps
1202 stats["shell_outside_measured"] = mean_pf
1203 stats["shell_inside_var"] = var_ps
1204 stats["shell_outside_var"] = var_pf
1205 stats["shell_inside"] = shell_ps
1206 stats["shell_outside"] = shell_pf
1207
1208 return shell_ps, shell_pf, stats
File dicebox.py changed (mode: 100644) (index 25ee32a..9fbad0a)
... ... from scipy import optimize # for BlackSimpleX
20 20
21 21
22 22 from .samplebox import SampleBox # for candidates packing from .samplebox import SampleBox # for candidates packing
23 from .ghull import Ghull
23 24 from . import simplex as sx from . import simplex as sx
24 25 from . import convex_hull as khull from . import convex_hull as khull
25 26 from . import lukiskon as lk from . import lukiskon as lk
26 27 from . import reader # for Goal from . import reader # for Goal
27 28 from . import estimation as stm # for KechatoLukiskon from . import estimation as stm # for KechatoLukiskon
28 29 import collections #E for Counter in MinEnergyCensoredSampling import collections #E for Counter in MinEnergyCensoredSampling
29
30
30 31
31 32
32 33
 
... ... class Goal(DiceBox):
1268 1269
1269 1270 def _regen_outside(bx): def _regen_outside(bx):
1270 1271 bx.convex_hull = khull.QHull(bx.f_model, space=bx.tri_space) # for gl_plot bx.convex_hull = khull.QHull(bx.f_model, space=bx.tri_space) # for gl_plot
1271 bx.ghull = khull.Ghull(bx.convex_hull)
1272 bx.ghull = Ghull(bx.convex_hull)
1272 1273 bx._R = -1 # update outer under R>_R condition bx._R = -1 # update outer under R>_R condition
1273 1274 bx._afacet = None bx._afacet = None
1274 1275 bx._bfacet = np.inf bx._bfacet = np.inf
File ghull.py added (mode: 100644) (index 0000000..c820fcd)
1 #!/usr/bin/env python
2 # coding: utf-8
3
4 import numpy as np
5 from . import sball
6 from scipy import stats
7 from .candybox import CandyBox
8 from .IS_stat import PushAndPull # for Shell_IS
9 from .IS_stat import get_1DS_sample # for Shell_1DS
10
11
12
13 try: # try to outthink f****** OS's memory management
14 import os
15 mem_bytes = os.sysconf('SC_PAGE_SIZE') * os.sysconf('SC_PHYS_PAGES')
16 except BaseException as e:
17 mem_GB = 16
18 mem_bytes = mem_GB * 1024**3 # hello, Windows!
19 print("ghull failed to get amount of RAM installed. %s GB is supposed."% mem_GB, repr(e))
20 print("BTW, you are using Windows, aren't you?")
21
22
23
24
25
26 #č pomocná třída pro integrování
27 #č nejsem tedy jist, jestli je to nejlepší napad -
28 #č dělit Ghull na podtřídy, delegovat funkcionalitu
29 #č ale jinak Ghull se stavá hodně překomplikovaným.
30 #č nic lepšího mně nenapadá, přemyšlel jsem dlouho.
31 class _ShellBaseIntegration:
32 def __init__(self, hull, shell, non_Gaussian_reduction=0.98):
33 self.shell = shell
34 self.hull = hull
35 self.nvar = hull.sample.nvar
36
37 self.integration_cutoff = np.inf
38 self.non_Gaussian_reduction = non_Gaussian_reduction
39
40 # test if r>0
41 self.r = 0
42 central_G = np.full(hull.sample.nvar, 0, dtype=np.int8)
43 self.DP = self.hull.sample.f_model.new_sample(central_G, space='G')
44
45
46 @property
47 def DP_is_valid(self):
48 #č dalo by se kontrolovat změny obálky
49 #č a provadět kontrolu DP pouze po změně obálky.
50 #č Nejsem liný, jeden bodík prostě za to doopravdy nestoji.
51 mask = self.hull.is_outside(self.DP)
52 return np.any(mask)
53
54 #č na rozdil od rodičovské třídy
55 #č vrací odhad r na základě předchozích integrací
56 #č metoda je navržena tak, aby Shell_IS jú mohl zdědit.
57 def get_r(self):
58 #č bojim sa, rerukce bude aplikována vždycky
59 #č Bacha, metoda bude vracet nuly pro obálky v Gaussovskem prostoru!
60 return self.r * self.non_Gaussian_reduction
61
62
63 #č metoda je navržena tak, aby Shell_IS jú mohl zdědit.
64 def _r_check_out(self, outside_nodes):
65 sum_squares = np.sum(np.square(outside_nodes.G), axis=1)
66 arg = np.nanargmin(sum_squares)
67 new_r = np.sqrt(sum_squares[arg])
68 if (not self.DP_is_valid) or (new_r < self.r):
69 self.DP = outside_nodes[arg]
70 self.r = new_r
71
72
73 #č metoda je navržena tak, aby Shell_IS jú mohl zdědit.
74 def integrate(self, nis, callback_all=None, callback_outside=None):
75 self.reset(nis)
76
77 if self.hull.nsimplex == 0:
78 bus = self.integration_cutoff
79 else:
80 bus = int(mem_bytes / self.hull.nsimplex / 8 / 3) + 1
81 while self.nsampled < nis:
82
83 seats = min(nis - self.nsampled, self.integration_cutoff, bus)
84 try:
85 self._process_part(seats, nis, callback_all, callback_outside)
86 except MemoryError as m:
87 print(self.__class__.__name__ +":", "memory error, %s sedaček" % seats, repr(m))
88 self.integration_cutoff = int(np.ceil(seats/2))
89
90 assert nis == self.nsampled
91
92 #č finalizovat by mohl i self._get_result(),
93 #č ale nebudeme michat vodku s pivem!
94 self._finalize()
95
96 return self._get_result()
97
98
99
100
101
102 class Shell_MC(_ShellBaseIntegration):
103
104 def reset(self, nis): # clear
105 self.r_needed = (self.hull.space!='G')
106
107 self.nsampled = 0
108 self.nfailed = 0
109
110
111 # stateless
112 def rvs(self, nsampled, seats, ns):
113 "Generování vzorků (kandidátů a integračních bodů)"
114 # rand_dir: prepare ns random directions on a unit d-sphere
115 rand_dir = sball.get_random_directions(seats, self.nvar) #random directions
116
117 # generate sampling probabilites
118 left = (ns-nsampled) / ns
119 right = (ns-nsampled-seats) / ns
120 #č přidáme trochu zmatku.
121 #č globálně jdeme r->R, localně ale R_i+ -> R_i-
122 #č menší poloměry musejí jít dřív - na to zavazano nonGaussian _r
123 #č převracení lokalně umožní vždycky mít alespoň jeden bodík outside,
124 #č což je taky velmi vhodné vzhledem k tomu, že se _r bere z outside bodíků
125 p = np.linspace(right, left, seats, endpoint=False) # probabilities for the radius
126
127 # convert probabilitites into random radii
128 # (distances from origin that are greater than r and less than R)
129 r = self.shell.isf(p) # actually, it is the same as CDF inverse
130
131 #finally a random sample from the optimal IS density:
132 sample_G = rand_dir*r[:,None]
133
134 return sample_G
135
136 # bus analogy
137 def _process_part(self, seats, nis, callback_all=None, callback_outside=None):
138 # boarding
139 nodes_G = self.rvs(self.nsampled, seats, nis)
140 nodes = self.hull.sample.f_model.new_sample(nodes_G, space='G')
141 # who is back?
142 mask = self.hull.is_outside(nodes)
143 if self.r_needed and np.any(mask):
144 #č rvs má vzorkovat od měnšího poloměru k většímu.
145 #č to znamená, že první outside bude mít nejměnší r vůbec
146 self._r_check_out(nodes[mask])
147 self.r_needed = False
148
149 if callback_all is not None:
150 # -2 = 'inside' -1 = 'outside'
151 candy_nodes = CandyBox(nodes, event_id=mask-2, is_outside=mask)
152 callback_all(candy_nodes)
153
154 if (callback_outside is not None) and np.any(mask):
155 callback_outside(nodes[mask])
156
157 assert len(mask) == seats
158 #č nevím proč, ale kdysi mě to vyšlo rychlejc jak obyčejný součet
159 self.nfailed += len(mask[mask])
160 self.nsampled += len(mask)
161
162
163 #č finalizovat by mohl i self._get_result(),
164 #č ale nebudeme michat vodku s pivem!
165 def _finalize(self):
166 #č nebyl žádný outside?
167 #č to přece není v pořádku!
168 #č (jen to nikomu neříkejte,
169 #č ale v současné implementaci třídy
170 #č podmínka dole by neměla nikdy nastat)
171 if self.r_needed:
172 self.r = self.shell.R * self.non_Gaussian_reduction
173
174 def _get_result(self):
175
176 nis = self.nsampled
177 nf = self.nfailed
178 ns = nis - nf
179 p_shell = self.shell.p_shell
180 shell_pf = nf/nis * p_shell
181 shell_ps = ns/nis * p_shell
182
183 stats = dict()
184 stats["shell_budget"] = nis
185 stats["shell_inside"] = shell_ps
186 stats["shell_outside"] = shell_pf
187
188 return shell_ps, shell_pf, stats
189
190
191
192
193 class Shell_IS(Shell_MC):
194 def reset(self, nis): # clear
195 self.r_needed = (self.hull.space!='G')
196
197 self.nsampled = 0
198 self.pp = PushAndPull()
199
200
201 # stateless
202 def rvs(self, nsampled, seats, ns):
203 "Generování vzorků (kandidátů a integračních bodů)"
204 # rand_dir: prepare ns random directions on a unit d-sphere
205 rand_dir = sball.get_random_directions(seats, self.nvar) #random directions
206
207 #č poloměry bereme ze skořapky
208 #č za správné nastavení (stejně jako u MC)
209 #č zodpovidá uživatel třídy
210 #č třída vlastní odhad r nijak nevyuživá!
211 r = self.shell.r
212 R = self.shell.R
213 delta_r = R - r
214 left = nsampled / ns * delta_r + r
215 right = (nsampled + seats) / ns * delta_r + r
216 #č přidáme trochu zmatku.
217 #č globálně jdeme r->R, localně ale R_i+ -> R_i-
218 #č menší poloměry musejí jít dřív - na to zavazano nonGaussian _r
219 #č převracení lokalně umožní vždycky mít alespoň jeden bodík outside,
220 #č což je taky velmi vhodné vzhledem k tomu, že se _r bere z outside bodíků
221 rs = np.linspace(right, left, seats, endpoint=False)
222
223 #finally a random sample from the optimal IS density:
224 sample_G = rand_dir*rs[:,None]
225
226 #č a jsme doma. Platba za špatný design.
227 #č Potřebujeme hustotu 1D rozdělení, implementovanou v Radial
228 #č Shell se síce dědí od Radial, ale implementuje .pdf() jako sdruženou
229 #č hustotu v nD Gaussovskem prostoru
230 #
231 #č vahy jsou definovány jako podil původního rozdělení k vzorkovácímu
232 #č původní - Chi rozdělení. nemá zde být žádný p_shell
233 f_pdf = stats.chi.pdf(rs, self.nvar)
234 #č vzorkovácí.
235 #h_pdf = 1/delta_r
236 weights = f_pdf * delta_r
237
238 return sample_G, weights
239
240 # bus analogy
241 def _process_part(self, seats, nis, callback_all=None, callback_outside=None):
242 # boarding
243 nodes_G, weights = self.rvs(self.nsampled, seats, nis)
244 nodes = self.hull.sample.f_model.new_sample(nodes_G, space='G')
245 # who is back?
246 mask = self.hull.is_outside(nodes)
247 if self.r_needed and np.any(mask):
248 #č rvs má vzorkovat od měnšího poloměru k většímu.
249 #č to znamená, že první outside bude mít nejměnší r vůbec
250 self._r_check_out(nodes[mask])
251 self.r_needed = False
252
253 assert len(mask) == seats
254 self.pp.add(weights, mask)
255 self.nsampled += len(mask)
256
257 if callback_all is not None:
258 # -2 = 'inside' -1 = 'outside'
259 candy_nodes = CandyBox(nodes, event_id=mask-2, is_outside=mask, weights=weights)
260 callback_all(candy_nodes)
261
262 if (callback_outside is not None) and np.any(mask):
263 callback_outside(nodes[mask])
264
265
266 def _get_result(self):
267
268 nis = self.nsampled
269 #č nejdřív true, pak false. Posilali jsme is_outside
270 mean_pf, mean_ps = self.pp.mean
271 var_pf, var_ps = self.pp.var
272 shell_pf, shell_ps = self.pp.correct_means(p_overall=self.shell.p_shell)
273
274 stats = dict()
275 stats["shell_budget"] = nis
276 stats["shell_inside_measured"] = mean_ps
277 stats["shell_outside_measured"] = mean_pf
278 stats["shell_inside_var"] = var_ps
279 stats["shell_outside_var"] = var_pf
280 stats["shell_inside"] = shell_ps
281 stats["shell_outside"] = shell_pf
282
283 return shell_ps, shell_pf, stats
284
285
286
287 #č slovní zásoba došla mně už před pár rokama
288 class Shell_1DS(_ShellBaseIntegration):
289 """1DS stands for 1D sampling.
290 1DS is, actually, an importance sampling method
291 that, actually, does not calculate IS weights
292 as f_density / h_density, but
293 (arbitrary, nonuniformly) divides interval to subintervals.
294 By using CDF transformes subintervals to 0-1 measure line.
295 Then gets one node as middle point of every subinterval,
296 weights therefore are just interval widths itself.
297 No sampling imprecisions are introduced,
298 therefore no spring, no correction are needed.
299 """
300 def reset(self, nis): # clear
301 self.r_needed = (self.hull.space!='G')
302
303 self.nsampled = 0
304
305 #č poloměry bereme ze skořapky
306 #č za správné nastavení (stejně jako u MC)
307 #č zodpovidá uživatel třídy
308 #č třída vlastní odhad r nijak nevyuživá!
309 r = self.shell.r
310 R = self.shell.R
311 # let's predefine 1D sequence at very beginning
312 x_sub = np.linspace(r, R, nis+1, endpoint=True)
313 x, weights = get_1DS_sample(stats.chi(self.nvar), x_sub)
314 self.x = x
315 self.weights = weights
316 self.mask = np.empty(nis, dtype=np.bool)
317
318
319
320 # bus analogy
321 def _process_part(self, seats, nis, callback_all=None, callback_outside=None):
322 # boarding
323 left = self.nsampled
324 right = self.nsampled + seats
325 rs = self.x[left:right]
326
327 # rand_dir: prepare ns random directions on a unit d-sphere
328 rand_dir = sball.get_random_directions(seats, self.nvar) #random directions
329 nodes_G = rand_dir*rs[:,None]
330 nodes = self.hull.sample.f_model.new_sample(nodes_G, space='G')
331
332 # who is back?
333 mask = self.hull.is_outside(nodes)
334 assert len(mask) == seats
335
336 if self.r_needed and np.any(mask):
337 #č má se vzorkovat od měnšího poloměru k většímu.
338 #č to znamená, že první outside bude mít nejměnší r vůbec
339 self._r_check_out(nodes[mask])
340 self.r_needed = False
341
342 # the most important part
343 self.nsampled += len(mask)
344 self.mask[left:right] = mask
345
346 if callback_all is not None:
347 # -2 = 'inside' -1 = 'outside'
348 candy_nodes = CandyBox(nodes, event_id=mask-2, is_outside=mask)
349 callback_all(candy_nodes)
350
351 if (callback_outside is not None) and np.any(mask):
352 callback_outside(nodes[mask])
353
354
355 def _finalize(self):
356 #č nebyl žádný outside?
357 #č to přece není v pořádku!
358 if self.r_needed:
359 self.r = self.shell.R * self.non_Gaussian_reduction
360
361
362 def _get_result(self):
363 #č mask related to hull.is_outside()
364 shell_pf = np.sum(self.weights[self.mask])
365 shell_ps = np.sum(self.weights[~self.mask])
366
367 stats = dict()
368 stats["shell_budget"] = self.nsampled
369 stats["shell_inside"] = shell_ps
370 stats["shell_outside"] = shell_pf
371
372 return shell_ps, shell_pf, stats
373
374
375
376
377 #č mým úkolem při návrhu této třidy je pořádně všecko zkomplikovat.
378 #č Dostávám za to peníze.
379 #č Takže. Udělám 3 druhů estimátorů
380 # convex_hull_estimation -2: inside, -1: outside
381 # shell_estimation -22: inner, -3: shell, -11: outer
382 # ghull_estimation -22: inner, -21: shell inside, -12: shell outside, -11: outer
383 class Ghull:
384 def __init__(self, hull, calculate_orth=True, calculate_2FORM=True, \
385 Integrator=Shell_1DS, non_Gaussian_reduction=0.9):
386 self.hull = hull
387 self.shell = sball.Shell(hull.sample.nvar)
388 self.outside_dist = sball.Shell(hull.sample.nvar)
389 self.sample = hull.sample
390
391 #č vlastně nevím, ale nemusejí byť úplně zadarmo...
392 self.calculate_orth = calculate_orth
393 self.calculate_2FORM = calculate_2FORM
394
395 self.set_integrator(Integrator, non_Gaussian_reduction)
396
397 def set_integrator(self, Integrator, non_Gaussian_reduction=0.9):
398 self.gint = Integrator(self.hull, self.shell, non_Gaussian_reduction)
399
400 def fire(self, ns, use_MC=False):
401 nodes = self.hull.fire(ns, use_MC=use_MC)
402 if nodes is not None:
403 return nodes
404 else:
405 return self.boom(ns)
406
407 def boom(self, ns, use_MC=False):
408 if use_MC:
409 self.outside_dist.set_bounds(self.get_R())
410 nodes_G = self.outside_dist.rvs(ns)
411 else:
412 # rand_dir: prepare ns random directions on a unit d-sphere
413 rand_dir = sball.get_random_directions(ns, self.sample.nvar) #random directions
414
415 # maximum radius, where norm.pdf() wasn't zero
416 # -38.575500173381374935388521407730877399444580
417 # don't ask me what the magic python use to distinguish
418 # digits after double precision
419 max_R_ever = -38.5755
420
421 r = np.linspace(self.get_R(), max_R_ever, ns, endpoint=True)
422 nodes_G = rand_dir*r[:,None]
423
424 nodes = self.sample.f_model.new_sample(nodes_G, space='G')
425 return nodes
426
427 def get_orth_outside(self):
428 if self.hull.space == 'G':
429 return self.hull.get_orth_outside()
430 else:
431 #č bude to horší jak s-ball, ale budiž.
432 x = np.full(self.sample.nvar, self.get_R())
433 return calculate_brick_complement_probability(-x, x)
434
435 def get_2FORM_outside(self):
436 if self.hull.space == 'G':
437 return self.hull.get_2FORM_outside()
438 else:
439 #č nebudem do toho michat nic dalšího.
440 #č Když 2FORM, tak 2FORM!
441 return 2 * stats.norm.sf(self.get_R())
442
443 def get_FORM_outside(self):
444 if self.hull.space == 'G':
445 return stats.norm.sf(self.get_r())
446 else:
447 return stats.norm.sf(self.get_R())
448
449 def get_R(self):
450 sum_squared = np.sum(np.square(self.sample.G), axis=1)
451 #index = np.argmax(sum_squared)
452 return np.sqrt(np.nanmax(sum_squared))
453
454 def get_r(self):
455 "calculates minimal signed distance to planes. Can therefore be negative"
456 #return -np.nanmax(self.hull.b)
457 if self.hull.space == 'G':
458 return self.hull.get_r()
459 else:
460 #č get_design_points() nemusí být.
461 #č QHull může nemít dostatek teček
462 try:
463 hull_points = self.hull.get_design_points()
464 sum_squared = np.sum(np.square(hull_points.G), axis=1)
465 hull_r = np.sqrt(np.nanmin(sum_squared))
466 except BaseException as e:
467 msg = "cannot get hull points from the convex hull"
468 print(self.__class__.__name__ +":", msg, repr(e))
469 hull_r = np.inf
470
471
472 # ask our experimentation team
473 gint_r = self.gint.get_r()
474
475 #č podle mně, nemůže se stat, že by metoda vratila:
476 #č 1. nenůlový poloměr a že by dokonce i centralní bod byl venku.
477 #č 2. poloměr větší jak R
478 #č Odpovědnost za to kladu na gint.
479 return min(hull_r, gint_r)
480
481 def get_shell_estimation(self):
482 shell = self.shell
483 r = self.get_r()
484 R = self.get_R()
485
486 if r<0:
487 shell.set_bounds(0, R)
488 else:
489 shell.set_bounds(r, R)
490
491 # shell_estimation -22: inner, -3: shell, -11: outer
492 shell_estimation = {-22:shell.ps, -3: shell.p_shell, -11: shell.pf}
493 global_stats = {"nsim":self.sample.nsim, "ndim":self.sample.nvar, \
494 "nfacets": self.hull.nsimplex, "r":r, "R":R, \
495 "inner":shell.ps, "shell":shell.p_shell, "outer":shell.pf}
496 global_stats['FORM_outside'] = self.get_FORM_outside()
497 if self.calculate_2FORM:
498 global_stats['2FORM_outside'] = self.get_2FORM_outside()
499 if self.calculate_orth:
500 global_stats['orth_outside'] = self.get_orth_outside()
501 return shell_estimation, global_stats
502
503 #č nie
504 ##č pro mně je důležité, aby bylo možné rychle přepinat mezi metodama,
505 ##č aby bylo možné výsledky jednoduše porovnavat.
506 ##č Proto explicitně posílám priznak, jakou metodu použit.
507 def integrate(self, nis, callback_all=None, callback_outside=None):
508 #č no to teda disajn třidy je doopravdy hroznej
509 # .get_shell_estimation() will calculate radiuses and will update shell
510 shell_estimation, global_stats = self.get_shell_estimation()
511 # shell_estimation -22: inner, -3: shell, -11: outer
512 shell_estimation.pop(-3)
513 # ghull_estimation -22: inner, -21: shell inside, -12: shell outside, -11: outer
514 ghull_estimation = shell_estimation; del(shell_estimation)
515
516
517 shell_ps, shell_pf, stats = self.gint.integrate(nis, callback_all, callback_outside)
518 ghull_estimation[-21] = shell_ps
519 ghull_estimation[-12] = shell_pf
520 global_stats.update(stats)
521
522 # convex_hull_estimation -2: inside, -1: outside
523 inside = shell_ps + self.shell.ps
524 outside = shell_pf + self.shell.pf
525 convex_hull_estimation = {-2: inside, -1: outside}
526 global_stats["inside"] = inside
527 global_stats["outside"] = outside
528
529 return ghull_estimation, convex_hull_estimation, global_stats
530
531
File qt_plot.py changed (mode: 100644) (index eaa8ca2..cd53aca)
... ... from . import stm_df
22 22 from . import sball from . import sball
23 23 from . import schemes from . import schemes
24 24 from . import convex_hull as khull from . import convex_hull as khull
25 #č pro mě je zvykem jako ghull označovat objekt třídy Ghull
26 #č nikoliv čerstvě oddelený modul
27 from .ghull import Ghull, Shell_MC, Shell_IS, Shell_1DS
25 28
26 29
27 30 def qt_gui_plot_2d(sample_box, space='R', *args, **kwargs): def qt_gui_plot_2d(sample_box, space='R', *args, **kwargs):
 
... ... class HullEstimationWidget(pg.LayoutWidget):
3033 3036 'title': "number of random directions", \ 'title': "number of random directions", \
3034 3037 'tip': "Used only for random scheme in DirectHull (or CompleteHull)"}) 'tip': "Used only for random scheme in DirectHull (or CompleteHull)"})
3035 3038
3036 params.append({'name': 'use_MC', 'title': "Use MC", 'type': 'bool', 'value': False })
3039 params.append({'name': 'integrator', 'title': "integrator", 'type': 'list', \
3040 'values': ['MC', 'IS', '1DS'], 'value': '1DS' })
3037 3041 params.append({'name': 'nonG_reduction', 'type': 'float', \ params.append({'name': 'nonG_reduction', 'type': 'float', \
3038 3042 'title': "non Gaussian reduction", \ 'title': "non Gaussian reduction", \
3039 3043 'limits': (0, 1), 'value': 0.9, 'default': 0.9,\ 'limits': (0, 1), 'value': 0.9, 'default': 0.9,\
 
... ... class HullEstimationWidget(pg.LayoutWidget):
3065 3069 self.param = pg.parametertree.Parameter.create(name='params', type='group', children=params) self.param = pg.parametertree.Parameter.create(name='params', type='group', children=params)
3066 3070
3067 3071 def get_ghull(self): def get_ghull(self):
3068 use_MC = self.param.getValues()['use_MC'][0]
3072 integrator = self.param.getValues()['integrator'][0]
3073 if integrator == 'MC':
3074 Integrator = Shell_MC
3075 elif integrator == 'IS':
3076 Integrator = Shell_IS
3077 elif integrator == '1DS':
3078 Integrator = Shell_1DS
3069 3079 nonG_reduction = self.param.getValues()['nonG_reduction'][0] nonG_reduction = self.param.getValues()['nonG_reduction'][0]
3070 3080 hull = self.get_hull() hull = self.get_hull()
3071 self.ghull = khull.Ghull(hull, use_MC=use_MC, non_Gaussian_reduction=nonG_reduction)
3081 self.ghull = Ghull(hull, Integrator=Integrator, non_Gaussian_reduction=nonG_reduction)
3072 3082 return self.ghull return self.ghull
3073 3083
3074 3084
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