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)
voronoi: WIP 23c63ede5bc2ec51d25aed850dc59299db438e34 I am 2022-10-21 12:27:01
IS_stat.sample_alike: use any possible fallback options to sample 80fb6d086e9f17350a1725a319e40bf29591869e I am 2022-10-20 08:56:42
voronoi: WIP ba912d406de12142f2e15a73090118184a0166a4 I am 2022-10-20 02:43:12
IS_stat: add sample_alike() function 2e5394d264b7aaaf5f2486c736ed8cc37afbd030 I am 2022-10-19 13:00:17
voronoi: WIP c3bfa60b6db2042df11528fe95a51d7cc88d3c6b I am 2022-10-19 11:18:13
voronoi: WIP baff3c5ccaf857f60194917fe5d6e33c40b01265 I am 2022-10-19 09:28:33
voronoi.ConvexSpline: little fix of np.sign() for better sleep 968f649e54d97abd596cd60749872fa3ef788c06 I am 2022-10-19 06:28:46
implement even more optimized Convex.. now ConvexSpline class c08511aecb32b880e73c9f10d07feff4dc588e8b I am 2022-10-19 04:11:13
voronoi: WIP cbb61d4544367a67c9a5c3fe6d64a9dc284eb364 I am 2022-10-18 05:05:44
voronoi: WIP e9c247fef30cfb678bad9b955d6f6a94a0ff61e7 I am 2022-10-18 01:02:30
candynodex: hotfix to support numpy masking indexing f359123b72a3d998d5eed4b9e611a5e402a18c8d I am 2022-10-17 09:02:33
add new lightweight candynodes module to replace old heavy ugly candybox ca97509f47c100df90e1e06d82ed7c759627bfd0 I am 2022-10-17 06:31:33
qt_gui: move CandyBox creating to the individual DiceBox widgets 29253625129219b5e550f82cae3da0e53ef5bd47 I am 2022-10-16 03:17:55
voronoi: WIP 4855190869e06df740cc05c64e8b27b2c3c5a88d I am 2022-10-15 22:46:51
voronoi: WIP 0298b3a0563587ea4fc24fa99b51d2e01d517203 I am 2022-10-15 20:00:25
voronoi: WIP 4b80757fa26e3d40dceac888480ec3d814f6abc9 I am 2022-10-15 17:30:11
voronoi: WIP 7cadd1d01d7f4bc3f51e785fe6ebf8cae6ff839a I am 2022-10-15 03:12:54
voronoi: clean up and comment out ContactSolver 543f8482a925ff6cf8485b9be616b5bd3c714d13 I am 2022-10-14 04:24:40
qt_gui.qt_gui: use ConvexSolver instead of ContactSolver (which did not work correctly) c1e49a9b2c3cc13886d317b460d7598f661cf216 I am 2022-10-14 04:01:06
voronoi: add fantastic ConvexSolver 2853826e5ef1abc79d3ac2fb8289b13c45211a31 I am 2022-10-14 04:00:14
Commit 23c63ede5bc2ec51d25aed850dc59299db438e34 - voronoi: WIP
Author: I am
Author date (UTC): 2022-10-21 12:27
Committer name: I am
Committer date (UTC): 2022-10-21 12:27
Parent(s): 80fb6d086e9f17350a1725a319e40bf29591869e
Signer:
Signing key:
Signing status: N
Tree: a633f99bdb5d8bb1332749eabd969b9b33b60431
File Lines added Lines deleted
wellmet/voronoi.py 348 140
File wellmet/voronoi.py changed (mode: 100644) (index 3ae660c..98e7151)
... ... class ContactVoronoi:
478 478 self.base_r = sball.get_Radius_ps(sample_box.nvar, p_base) self.base_r = sball.get_Radius_ps(sample_box.nvar, p_base)
479 479 #č pro funkci IS_stat.sample_like(), víz. diskuse v poznamkách #č pro funkci IS_stat.sample_like(), víz. diskuse v poznamkách
480 480 self.d = (sball.get_Radius_ps(sample_box.nvar - 1, p_base) / self.base_r)**2 self.d = (sball.get_Radius_ps(sample_box.nvar - 1, p_base) / self.base_r)**2
481
481 self.ndim = sample_box.nvar
482 482
483 483 self._nsim = 0 self._nsim = 0
484 484 self._enabled = False # if nsim < 2 and not failsi self._enabled = False # if nsim < 2 and not failsi
 
... ... class ContactVoronoi:
559 559 def update_contacts(self): def update_contacts(self):
560 560 #č já vím, že sample box pokážde failsi přepočítavá #č já vím, že sample box pokážde failsi přepočítavá
561 561 self.failsi = failsi = self.sample_box.failsi self.failsi = failsi = self.sample_box.failsi
562 self.failure_points = self.sample_box.failure_points # for undiscovered
562 563
563 564 self.CS = ConvexSpline(sampled_plan_model) self.CS = ConvexSpline(sampled_plan_model)
564 565 self.PDF = self.sample_box.pdf(self.model_space) self.PDF = self.sample_box.pdf(self.model_space)
 
... ... class ContactVoronoi:
658 659 self.red_couples[couple] = nodes_idx self.red_couples[couple] = nodes_idx
659 660
660 661
661 def update_couple(self, couple_indices, event_id):
662 pass
663
664 def invalidate_couple(self, couple_indices):
665 pass
666
667 def onboard_couple(self, couple_indices, event_id):
662 def _check_undiscovered(self, gmask, ii, dd, nodes_model, h_pdf=None):
668 663 """ """
669 č Obecně, skoro pro každou dimenzi v rozmězí od 2D do 100D
670 č nejlevnější bylo
671 č 1. strom s pár bodíků
672 č 2. ConvexSolver
673 č 3. strom s nsim bodíků
664 č logika celé třídy je, že máme-li příslušné kontaktu bodíky,
665 č tak stojí za to je uložit.
674 666 """ """
667 couples = self.couples
675 668
676 #č (i, j)? poprvé ta čísla vidíme
677 i, j = couple_indices
669 iim = ii[gmask] #č dovezli nám nedodělané zboží
678 670
679 #č Nejdřív - zda je přímy kontakt, prostě na usečce
680 couple_nodes = self.sampled_plan_model[[i, j]]
681 contact_node_model = np.sum(couple_nodes, axis=0)
682 np.divide(contact_node_model, 2, out=contact_node_model)
671 #č vymaskovat čerstvé víno
672 mask = np.all(iim < self._nsim, axis=1)
673 iim = iim[mask] #č chcu/nechcu globální masku
674 #č vymaskovat zelené kontakty
675 in_failure = np.isin(iim, self.failure_points)
676 has_failure = in_failure.any(axis=1)
677 iim = iim[has_failure]
683 678
684 #č bodík nebude mít vahu - nechť si s tím integrovačka poradí!
685 #č (ten bodík by neměl být omylem doporučen k integraci)
686 idx, dd, ii, mask = self._store_masked(contact_node_model, couple_indices, event_id)
687
688 #č kontakt neexistuje, netřeba nic robiť, j-ko i-čkem není dotčen
689 if not (idx or self.CS.is_couple(couple_indices)):
679 #č zbytečnost
680 if not iim.size: # if empty
690 681 return None return None
691 682
692 #č kontakt existuje, dáme vědět aktualizovačce
693 #č (snad jediné místo kde tuhle funkci voláme)
694 self._add_indices_to_update(j)
683 #č min-max nepočítá se zadarmo
684 #č nejdřív projedeme smyčkou část,
685 #č kde je první index větší a pak ten zbytek
686 mask = iim[:0] > iim[:1]
687 for i, j in iim[mask]:
688 if (i, j) not in couples:
689 #č ok. Ten vzacnej případ nastal.
690 #č kašleme na masky, posíláme všechno na explore
691 self._explore(couple, ii, dd, nodes_model, h_pdf)
692
693 for j, i in iim[~mask]:
694 if (i, j) not in couples:
695 #č ok. Ten vzacnej případ nastal.
696 #č kašleme na masky, posíláme všechno na explore
697 self._explore(couple, ii, dd, nodes_model, h_pdf)
698
699
700
701
702 def _explore(self, couple, ii, dd, nodes_model, h_pdf):
703 print("ContactVoronoi: undiscovered" + str(couple)\
704 + " pair has been found!")
705 i, j = couple
706
707 # similar approach is used in numpy code
708 mask = np.all((ii == i) | (ii == j), axis=1)
709
710 self._check_undiscovered(ii[~mask], dd[~mask], nodes_model[~mask])
711
712 if not np.any(mask):
713 return False, dd, ii, mask
714
715
716 #č přídat, zavolat kólbek
717 f = self.sample_box.f_model
718 #č tady je to super, ťažké datové struktury f-modelu
719 #č vytvaříme na jíž vyfiltrovaných datéch
720 #č Tahle filtrace stromem zaručuje,
721 #č že bodíky (přenejmenším na vstupu f_modelu)
722 #č mají konečné souřadnice bez nan a nekonečno.
723 #č To ale nic neříká o hustotách
724 nodes = f.new_sample(nodes_model[mask], self.model_space)
725 couple_stats = {'event_id': event_id, 'couple': couple_indices}
726 nodes = CandyNodes(nodes, couple_stats)
727 idx = self.nodes.add(nodes)
728 nodes.idx = idx
729
730
731 ddm = dd[mask] #č tohle je kvůli malé pomyšlené
732 nodes.dd1 = ddm[:, 0] #č pytlovince v sample_alike
733
734 #č červené kontakty prostě rovnou integrujeme
735 #č takže netřeba pro ně ukladat zbytečná data
736 # -1 = 'outside', 0=success, 1=failure, 2=mix
737 if event_id == 2:
738 #ddm = dd[mask]
739 iim = ii[mask]
740 #nodes.dd1 = ddm[:, 0] #č teď je pro všechny
741 nodes.dd2 = ddm[:, 1]
742 nodes.ii1 = iim[:, 0]
743 nodes.ii2 = iim[:, 1]
744 if self.on_add_mixed is not None:
745 self.on_add_mixed(idx)
746
747
748 ## checklist
749 self.couples[couple_indices].append(idx)
750 #self.couples_stats #č nevím
751 #self.mixed_couples #č v této fázi nevíme,
752 #self.red_couples #č zda je sada ta nejlepší
753 #self._add_indices_to_update(j) #č vždyť já nevím, jestli se jedná o aktualizaci
754
755 return idx, dd, ii, mask
756
695 757
696 758
697 #č a teď deme vzorkovat
698 ndim = len(contact_node_model)
699 #č vzdalenost tamten strom už má spočítanou (pro zadanou p_norm)
700 r = dd[0,0]
701 # #č ale nechcu kvůli jedinému poloměru komplikovat kód
702 # vec = couple_nodes[0] - couple_nodes[1]
703 # r = np.sqrt(np.inner(vec, vec)) / 2
704 759
705 760 # initially we'll perform IS from multivariate Gaussian distribution # initially we'll perform IS from multivariate Gaussian distribution
706 761 # with center at contact node (i. e. in beetween of two points) # with center at contact node (i. e. in beetween of two points)
707 762 # and single standard deviation (let's say "ball" sampling) # and single standard deviation (let's say "ball" sampling)
708 763 #č r-ko dělíme base_r'kem k získání sigmy. #č r-ko dělíme base_r'kem k získání sigmy.
709 sampling_plan_N, pdf = IS_stat.get_norm_plan(self.ns, ndim, \
764 sampling_plan_N, h_pdf = IS_stat.get_norm_plan(self.ns, ndim, \
710 765 mean=contact_node_model, mean=contact_node_model,
711 766 std=r/self.base_r, design=None) std=r/self.base_r, design=None)
712 767
713 idx, dd, ii, mask = self._store_masked(sampling_plan_N, couple_indices, event_id)
768 idx, score = self.score(sampling_plan_N, h_pdf, couple_indices, event_id)
769
714 770 if not idx: if not idx:
715 771 return None return None
716 772
717 #č uspěch. Dodáme bodíkům vahy
718 nodes = self.nodes[idx]
719 w = nodes.pdf(self.model_space) / pdf[mask] / self.ns
720 #č Tamta filtrace stromem zaručuje,
721 #č že bodíky (přenejmenším na vstupu f_modelu)
722 #č mají konečné souřadnice bez nan a nekonečno.
723 #č To ale nic neříká o hustotách
724 mask = np.isfinite(w)
725 nodes.w = w * mask
726 #č p_couple není finální pravděpodobnost
727 #č bodíky ještě musejí projít konvexní obálkou
728 nodes.p_couple = np.sum(w[mask])
729 nodes.score = nodes.p_couple / np.sqrt(np.var(nodes.w)/len(nodes))
730 #č a doporučíme (dočasně) k integraci
773 #č doporučíme (dočasně) k integraci
731 774 self._recommend_to_integration(idx) self._recommend_to_integration(idx)
732 775
733 if len(nodes) < ndim + 1:
734 return None
735
736 776 #č jdeme do adaptivního IS vzorkování #č jdeme do adaptivního IS vzorkování
737 777 #č zde jenom vzorkujeme; #č zde jenom vzorkujeme;
738 778 #č integrování, vyhodnocování je jinde a později! #č integrování, vyhodnocování je jinde a později!
 
... ... class ContactVoronoi:
1035 1075
1036 1076
1037 1077
1038 def score(self, nodes_model, h_pdf, couple_indices, event_id):
1039 "returns: node_index, score"
1040 1078
1041 idx, _dd, _ii, mask = self._store_masked(nodes_model, couple_indices, event_id)
1042 if not idx:
1043 return idx, -1
1079
1044 1080
1045 #č uspěch. Dodáme bodíkům vahy
1046 nodes = self.nodes[idx]
1047 nis = len(h_pdf)
1048 w = nodes.pdf(self.model_space) / h_pdf[mask] / nis
1049 #č Tamta filtrace stromem zaručuje,
1081 #č přídat, zavolat kólbek
1082 f = self.sample_box.f_model
1083 #č tady je to super, ťažké datové struktury f-modelu
1084 #č vytvaříme na jíž vyfiltrovaných datéch
1085 #č Tahle filtrace stromem zaručuje,
1050 1086 #č že bodíky (přenejmenším na vstupu f_modelu) #č že bodíky (přenejmenším na vstupu f_modelu)
1051 1087 #č mají konečné souřadnice bez nan a nekonečno. #č mají konečné souřadnice bez nan a nekonečno.
1052 1088 #č To ale nic neříká o hustotách #č To ale nic neříká o hustotách
1053 mask = np.isfinite(w) #č bacha! Nová, již jiná maska!
1054 nodes.w = w * mask
1055 #č p_couple není finální pravděpodobnost
1056 #č bodíky ještě musejí projít konvexní obálkou
1057 nodes.p_couple = np.sum(w[mask])
1058 nodes.score = nodes.p_couple / np.sqrt(np.var(nodes.w)/len(nodes))
1059 return idx, nodes.score
1089 nodes = f.new_sample(nodes_model[mask], self.model_space)
1090 couple_stats = {'event_id': event_id, 'couple': couple_indices}
1091 nodes = CandyNodes(nodes, couple_stats)
1092 idx = self.nodes.add(nodes)
1093 nodes.idx = idx
1094
1095
1096 ddm = dd[mask] #č tohle je kvůli malé pomyšlené
1097 nodes.dd1 = ddm[:, 0] #č pytlovince v sample_alike
1098
1099 #č červené kontakty prostě rovnou integrujeme
1100 #č takže netřeba pro ně ukladat zbytečná data
1101 # -1 = 'outside', 0=success, 1=failure, 2=mix
1102 if event_id == 2:
1103 #ddm = dd[mask]
1104 iim = ii[mask]
1105 #nodes.dd1 = ddm[:, 0] #č teď je pro všechny
1106 nodes.dd2 = ddm[:, 1]
1107 nodes.ii1 = iim[:, 0]
1108 nodes.ii2 = iim[:, 1]
1109 if self.on_add_mixed is not None:
1110 self.on_add_mixed(idx)
1111
1112
1113 ## checklist
1114 self.couples[couple_indices].append(idx)
1115 #self.couples_stats #č nevím
1116 #self.mixed_couples #č v této fázi nevíme,
1117 #self.red_couples #č zda je sada ta nejlepší
1118 #self._add_indices_to_update(j) #č vždyť já nevím, jestli se jedná o aktualizaci
1119
1120 return idx, dd, ii, mask
1060 1121
1061 1122
1123 def update_couple(self, couple_indices, event_id):
1124 #č při aktualizaci páry kontrolujeme,
1125 #č zda nejsou změny u sady, doporučenou k integraci
1126 #č teček je třeba mít na vědomí,
1127 #č že se můžou zůstávat oblasti, ve kterých jsou
1128 pass
1062 1129
1063 1130
1131 def invalidate_couple(self, couple_indices):
1132 pass
1064 1133
1065 # for simplex: d = nvar+2
1066 # for cell: d = base_r**2
1067 def IS_like(f_plan, sampling_space='G', weights=None, nis=1000, d=1, design=None):
1134 def onboard_couple(self, couple_indices, event_id):
1068 1135 """ """
1069 takes sample and returns sampling plan with the same cov (in case d=1)
1070 covariance matrix we'll divide by d
1136 č Obecně, skoro pro každou dimenzi v rozmězí od 2D do 100D
1137 č nejlevnější bylo
1138 č 1. strom s pár bodíků
1139 č 2. ConvexSolver
1140 č 3. strom s nsim bodíků
1071 1141 """ """
1072 1142
1073 plan = getattr(f_plan, sampling_space)
1143 #č (i, j)? poprvé ta čísla vidíme
1144 i, j = couple_indices
1074 1145
1075
1076 S_bc = np.cov(plan, rowvar=False, bias=True, aweights=weights)
1077 barycenter = np.average(plan, axis=0, weights=weights)
1146 #č Nejdřív - zda je přímy kontakt, prostě na usečce
1147 couple_nodes = self.sampled_plan_model[[i, j]]
1148 contact_node_model = np.sum(couple_nodes, axis=0)
1149 np.divide(contact_node_model, 2, out=contact_node_model)
1150
1151 #č bodík nebude mít vahu - nechť si s tím integrovačka poradí!
1152 #č (ten bodík by neměl být omylem doporučen k integraci)
1153 idx, dd, ii, mask = self._store_masked(contact_node_model, couple_indices, event_id)
1154
1155 #č kontakt neexistuje, netřeba nic robiť, j-ko i-čkem není dotčen
1156 if not (idx or self.CS.is_couple(couple_indices)):
1157 return None
1078 1158
1159 #č kontakt existuje, dáme vědět aktualizovačce
1160 #č (snad jediné místo kde tuhle funkci voláme)
1161 self._add_indices_to_update(j)
1079 1162
1080 #č matika
1081 w, v = np.linalg.eig(S_bc)
1082 1163
1083 # use IS sampling density with center equal to the simplex's barycenter
1084 # set the minimum distance as the standard deviation of IS densisty
1085 #č u stats.norm zadáváme směrodatnou odchylku, to je asi správné
1086 sigmas = np.sqrt(w/d)
1087 mean = 0
1088 h_plan_N, pdf_N = get_norm_plan(nis, f_plan.nvar, mean, sigmas, design)
1164 #č a teď deme vzorkovat
1165 ndim = len(contact_node_model)
1166 #č vzdalenost tamten strom už má spočítanou (pro zadanou p_norm)
1167 r = dd[0,0]
1168 # #č ale nechcu kvůli jedinému poloměru komplikovat kód
1169 # vec = couple_nodes[0] - couple_nodes[1]
1170 # r = np.sqrt(np.inner(vec, vec)) / 2
1089 1171
1090 #ёӵ здесь уже так легко не отделаемся. Трансформовать кароно.
1091 h_plan_bc = (v @ h_plan_N.T).T
1092 h_plan_sing = h_plan_bc + barycenter
1172 # initially we'll perform IS from multivariate Gaussian distribution
1173 # with center at contact node (i. e. in beetween of two points)
1174 # and single standard deviation (let's say "ball" sampling)
1175 #č r-ko dělíme base_r'kem k získání sigmy.
1176 sampling_plan_N, h_pdf = IS_stat.get_norm_plan(self.ns, ndim, \
1177 mean=contact_node_model,
1178 std=r/self.base_r, design=None)
1179
1180 idx, score = self.score(sampling_plan_N, h_pdf, couple_indices, event_id, h_pdf)
1181
1182 if not idx:
1183 return None
1184
1185 #č doporučíme (dočasně) k integraci
1186 self._recommend_to_integration(idx)
1093 1187
1188 if len(nodes) < ndim + 1:
1189 return None
1190
1191 #č jdeme do adaptivního IS vzorkování
1192 #č zde jenom vzorkujeme;
1193 #č integrování, vyhodnocování je jinde a později!
1094 1194
1195 nodes_model = getattr(nodes, self.model_space)
1196 status, nodes_model, h_pdf = IS_stat.sample_alike(nodes_model, \
1197 weights=nodes.w, nis=self.ns, d=self.d)
1198 if not status:
1199 return None
1200
1201 idx, dd, ii, mask = self._store_masked(nodes_model, couple_indices, event_id, h_pdf)
1202 if not idx:
1203 return None
1204
1205 #č uspěch. Dodáme bodíkům vahy
1206 nodes = self.nodes[idx]
1207 w = nodes.pdf(self.model_space) / h_pdf[mask] / self.ns
1208 #č Tamta filtrace stromem zaručuje,
1209 #č že bodíky (přenejmenším na vstupu f_modelu)
1210 #č mají konečné souřadnice bez nan a nekonečno.
1211 #č To ale nic neříká o hustotách
1212 mask = np.isfinite(w)
1213 nodes.w = w * mask
1214 #č p_couple není finální pravděpodobnost
1215 #č bodíky ještě musejí projít konvexní obálkou
1216 nodes.p_couple = np.sum(w[mask])
1217 nodes.score = nodes.p_couple / np.sqrt(np.var(nodes.w)/len(nodes))
1218 #č a doporučíme (dočasně) k integraci
1219 self._recommend_to_integration(idx)
1095 1220
1096 #č sice to má nazev h_plan, ale nese rozdělení a hustoty v f-ku
1097 #
1098 #č jelikož výtvaříme nový vzorek z "čístých" souřadnic
1099 #č máme zaručeno, že h_plan vždy bude f_modelem (podle současného kódu CandyBox)
1100 h_plan = f_plan.new_sample(h_plan_sing, sampling_space)
1101 ## desired: w = h_plan.pdf(sampling_space) / pdf_N # snad je to správně
1102 w = np.divide(h_plan.pdf(sampling_space), pdf_N, out=pdf_N); del(pdf_N)
1103 1221
1104 #č vahy máme
1105 #č zabalme do boxu
1106 #č zbytek už nejsou naši starosti
1107 return CandyBox(h_plan, w=w)
1222 if len(nodes) < nodes.nvar + 1:
1223 # in case of update reset the final score only
1224 nodes.score = -2
1225 return None
1108 1226
1227 nodes_model = getattr(nodes, self.model_space)
1228 nodes_pdf = nodes.pdf(self.model_space)
1109 1229
1110 def sample_alike(self, nodes):
1111 """
1112 č zde jenom vzorkujeme;
1113 č integrování, vyhodnocování je jinde a později!
1114 """
1230 S_bc = np.cov(nodes_model, rowvar=False, bias=True, aweights=nodes_pdf)
1115 1231
1232 #č matika
1233 w, v = eig = np.linalg.eig(S_bc)
1116 1234
1235 if not np.all(w > 0):
1236 nodes.score = -1
1237 return None
1117 1238
1118 # use IS sampling density with center equal to the current point
1119 # set the minimum distance as the standard deviation of IS densisty
1120 # mindist dělíme dvěma - hranice buňky je na půlcesty
1121 # a dělíme base_r'kem k získání sigmy. Přemejmenším empiricky mně to vychází
1122 # (u Norm'u vždy zadaváme směrodatnou odchylku)
1123 h_i = [stats.norm(sampled_plan_sing[i,j], mindist_sing[i]/2/base_r) for j in range(sample_box.nvar)]
1124 h = f_models.UnCorD(h_i)
1239 barycenter = np.average(nodes_model, axis=0, weights=nodes_pdf)
1240 nodes.set_attrs(eig=eig, barycenter=barycenter)
1241 nodes.score = np.trace(S_bc) * np.sqrt(len(nodes))
1125 1242
1126
1127 # select nis points from IS density
1128 # sice to má nazev h_plan, ale nese rozdělení a hustoty v f-ku
1129 h_plan = IS(f, h, space_from_h='R', space_to_f=sampling_space, Nsim=nis)
1130 1243
1131 # součet váh nemá cenu kontrolovat, jedná tam nebude, nevyjde
1132 1244
1133 h_plan_model = getattr(h_plan, model_space)
1134 dd, ii = tree.query(h_plan_model, k=1, p=p_norm)
1135 1245
1136 1246
1137 1247 # nechám s velkým písmenkem # nechám s velkým písmenkem
 
... ... class ContactVoronoi:
1141 1251 # dd1 jsou vzdalenosti tečiček do centra Voroneho buňky # dd1 jsou vzdalenosti tečiček do centra Voroneho buňky
1142 1252 dd1 = dd[Vor_mask] dd1 = dd[Vor_mask]
1143 1253
1144 trace = 0
1254
1145 1255 nis_ma = len(weights_sim) nis_ma = len(weights_sim)
1146 1256
1147 1257 h_plan_sing_ma = getattr(h_plan, sampling_space)[Vor_mask] h_plan_sing_ma = getattr(h_plan, sampling_space)[Vor_mask]
1148 S_bc = np.cov(h_plan_sing_ma, rowvar=False, bias=True, aweights=weights_sim)
1149 barycenter = np.average(h_plan_sing_ma, axis=0, weights=weights_sim)
1258
1259 nodes_model = getattr(nodes, self.model_space)
1260 nodes_pdf = nodes.pdf(self.model_space)
1261 S_bc = np.cov(nodes_model, rowvar=False, bias=True, aweights=nodes_pdf)
1262 barycenter = np.average(nodes_model, axis=0, weights=nodes_pdf)
1263 trace = np.trace(S_bc)
1264 nodes.set_attrs(S_bc=S_bc, barycenter=barycenter, trace=trace)
1265
1266 # matika
1267 w, v = np.linalg.eig(S_bc)
1268
1269 if np.any(w <= 0):
1270 return None
1271
1272 goal = trace*len(nodes)
1273
1150 1274 # effective nis # effective nis
1151 1275 nis_eff = nis nis_eff = nis
1152 1276 #print("start", i) #print("start", i)
1153 1277
1154 while trace < np.trace(S_bc):
1155 # number of nodes inside the cell matters too
1156 trace = np.trace(S_bc)# * nis_ma
1157 #print(S_bc)
1158 #print(nis_ma)
1159
1160
1161
1162
1278 max_trace = 0
1279 while max_trace < trace:
1280 trace = np.trace(S_bc)
1281
1163 1282 # matika # matika
1164 1283 w, v = np.linalg.eig(S_bc) w, v = np.linalg.eig(S_bc)
1165 1284
1285 # r_ball/r_base = sigma
1286
1166 1287 # use IS sampling density with center equal to the cell's barycenter # use IS sampling density with center equal to the cell's barycenter
1167 1288 # set the minimum distance as the standard deviation of the IS densisty # set the minimum distance as the standard deviation of the IS densisty
1168 1289 # u stats.norm zadáváme směrodatnou odchylku, to je asi správné # u stats.norm zadáváme směrodatnou odchylku, to je asi správné
 
... ... class ContactVoronoi:
1368 1489 mindist_sing = dd2.flatten() mindist_sing = dd2.flatten()
1369 1490
1370 1491
1371 def _store_masked(self, nodes_model, couple_indices, event_id):
1492
1493 def score(self, nodes_model, h_pdf, couple_indices, event_id):
1494 "returns: node_index, score"
1495
1496 idx, _dd, _ii, mask = self._store_masked(nodes_model, couple_indices, event_id, h_pdf)
1497 if not idx:
1498 return idx, -1
1499
1500 #č uspěch. Dodáme bodíkům vahy
1501 nodes = self.nodes[idx]
1502 nis = len(h_pdf)
1503 w = nodes.pdf(self.model_space) / h_pdf[mask] / nis
1504 #č Tamta filtrace stromem zaručuje,
1505 #č že bodíky (přenejmenším na vstupu f_modelu)
1506 #č mají konečné souřadnice bez nan a nekonečno.
1507 #č To ale nic neříká o hustotách
1508 mask = np.isfinite(w) #č bacha! Nová, již jiná maska!
1509 nodes.w = w * mask
1510 #č p_couple není finální pravděpodobnost
1511 #č bodíky ještě musejí projít konvexní obálkou
1512 nodes.p_couple = np.sum(w[mask])
1513 #č rozptyl může být i nulový
1514 #nodes.score = nodes.p_couple / np.sqrt(np.var(nodes.w)/len(nodes))
1515 nodes.score = nodes.p_couple * np.sqrt(len(nodes))
1516 return idx, nodes.score
1517
1518
1519
1520
1521
1522
1523 def sample_alike(self, nodes):
1524 """
1525 č zde jenom vzorkujeme;
1526 č integrování, vyhodnocování je jinde a později!
1527 returns: node_idx, score
1528 """
1529 assert len(nodes) > 0
1530
1531 nodes_model = getattr(nodes, self.model_space)
1532
1533 if len(nodes) < 2:
1534 #č máme jeden bodík. To znaméná, že kontakt existuje
1535 #č a zákazník hodla provést vzorkování.
1536 #č Pro tento vzácný případ naše společnost nabízí řešení -
1537 #č zapojít informace o těch dvou vzorkách, ke kterým patří.
1538 #č Na ty dva vzorky se nemůžeme spoléhat z hlediska směrových
1539 #č preferencí, použijme je jen jako zdroj jakési
1540 #č charakteristické délky. Těžiště necháme v tamtamtom bodíku.
1541
1542 #č Specifika CandyNodes je, že dd1 pro jeden bodík může vrátit
1543 #č jak skalarní hodnotu, tak i pole s jediným prvkem
1544 #č numpy si poradí s obojí
1545 r = nodes.dd1
1546
1547 # we'll perform IS from multivariate Gaussian distribution
1548 # centered at the isolated node
1549 # and some standard deviation
1550 #č r-ko dělíme base_r'kem k získání sigmy.
1551 sampling_plan_N, h_pdf = IS_stat.get_norm_plan(self.ns, self.ndim,\
1552 mean=nodes_model[0],
1553 std=r/self.base_r, design=None)
1554
1555 return self.score(sampling_plan_N, h_pdf, nodes.couple, nodes.event_id)
1556
1557
1558 try:
1559 weights = nodes.w
1560 except AttributeError:
1561 weights = None
1562
1563 sampling_plan, h_pdf = IS_stat.sample_alike(nodes_model, \
1564 weights=weights, nis=self.ns, d=self.d)
1565
1566 return self.score(sampling_plan, h_pdf, nodes.couple, nodes.event_id)
1567
1568
1569
1570 #č h_pdf je jen pro _check_undiscovered()
1571 def _store_masked(self, nodes_model, couple_indices, event_id, h_pdf=None):
1372 1572 """ """
1373 1573 č logika celé třídy je, že máme-li příslušné kontaktu bodíky, č logika celé třídy je, že máme-li příslušné kontaktu bodíky,
1374 1574 č tak stojí za to je uložit. č tak stojí za to je uložit.
 
... ... class ContactVoronoi:
1386 1586
1387 1587 # similar approach is used in numpy code # similar approach is used in numpy code
1388 1588 mask = np.all((ii == i) | (ii == j), axis=1) mask = np.all((ii == i) | (ii == j), axis=1)
1589 if not np.all(mask):
1590 #č aniž by se tu trapil if-ámi s h_pdf, pošlu všecko do funkce,
1591 #č neboť ta to všechno potřebuje jen ve velmi vzacném případě
1592 self._check_undiscovered(~mask, ii, dd, nodes_model, h_pdf)
1389 1593
1390 1594 if not np.any(mask): if not np.any(mask):
1391 1595 return False, dd, ii, mask return False, dd, ii, mask
 
... ... class ContactVoronoi:
1405 1609 idx = self.nodes.add(nodes) idx = self.nodes.add(nodes)
1406 1610 nodes.idx = idx nodes.idx = idx
1407 1611
1612
1613 ddm = dd[mask] #č tohle je kvůli malé pomyšlené
1614 nodes.dd1 = ddm[:, 0] #č pytlovince v sample_alike
1615
1408 1616 #č červené kontakty prostě rovnou integrujeme #č červené kontakty prostě rovnou integrujeme
1409 1617 #č takže netřeba pro ně ukladat zbytečná data #č takže netřeba pro ně ukladat zbytečná data
1410 1618 # -1 = 'outside', 0=success, 1=failure, 2=mix # -1 = 'outside', 0=success, 1=failure, 2=mix
1411 1619 if event_id == 2: if event_id == 2:
1412 dd = dd[mask]
1413 ii = ii[mask]
1414 nodes.dd1 = dd[:, 0]
1415 nodes.dd2 = dd[:, 1]
1416 nodes.ii1 = ii[:, 0]
1417 nodes.ii2 = ii[:, 1]
1620 #ddm = dd[mask]
1621 iim = ii[mask]
1622 #nodes.dd1 = ddm[:, 0] #č teď je pro všechny
1623 nodes.dd2 = ddm[:, 1]
1624 nodes.ii1 = iim[:, 0]
1625 nodes.ii2 = iim[:, 1]
1418 1626 if self.on_add_mixed is not None: if self.on_add_mixed is not None:
1419 1627 self.on_add_mixed(idx) self.on_add_mixed(idx)
1420 1628
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