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)
estimation: WIP 726058435dcf653f816390e3fc996747fff07908 Олёш 2020-07-14 02:41:20
estimation: simple tesselation added 344f275e586e42c4946f6d4f18a2503ec5f28869 Олёш 2020-07-11 22:25:08
estimation: 2 point Voronoi polishing 18f6563c0b209f4c43c4d19a9dac5f4ff74a790a Олёш 2020-07-11 05:30:10
f_models & whitebox: object generation (creation) fix 3b2e17e86e14f0eaea0e3264e8d6a4f0930505f5 Олёш 2020-07-10 20:54:23
estimation: WIP, Voronoi_2_point worked 364c89027063f0bd3a1af80d7d53db059184dd68 Олёш 2020-07-10 02:51:33
IS_stat: one more IS implementation a23ab3ff5eac9a039f415b9281a02ebbbea778f2 Олёш 2020-07-09 03:09:14
f_models: .pdf() and .sample_pdf() redesigned 40dfc578e14f630c1735a279a3ec41d169cc7606 Олёш 2020-07-08 01:18:29
f_models: linear tranformations and .pdf function added 7ce23a722836eaa3488329870691a332641a92b1 Олёш 2020-07-07 03:36:07
estimation: WIP 602b24888d5b22e7a994a14119af3a12c53a374d Alex 2020-06-29 08:58:54
g_models: ConicSection added e603e051a42dec4694766d536e71300dc62cb3a8 Олёш 2020-06-27 11:11:31
qt_plot: slider and triangulation fix 04a8bf659873bed50b82dc94497f2fc473c883dd Олёш 2020-06-26 19:56:58
qt_plot: show triangulation 0aa36d61be95b6f40f4a01e77a61a701c78789ce Олёш 2020-06-26 10:10:39
f_models: přidany chybějící pdf-ka 1f3c99bcd1eaeed9b8b84175be9189449a7dec5f Олёш 2020-06-25 20:28:33
qt_plot: connection to BlackBox feature added fcff3099dc61e568981e1f1236c29b5326547fd5 Олёш 2020-06-23 09:52:56
qt_plot: estimation graphs are implemented a31c5b3c344893f0f0b7022ee54544493e97c79a Олёш 2020-06-21 22:21:55
qt_plot: WIP 7cf97855991f717873d02ce65058790b324a319d Олёш 2020-06-21 10:17:04
qt_plot: minor changes 931ca9629256f37588e7ed4cb497f71c93954a7b Олёш 2020-06-15 06:01:10
candybox introduced, stores additional data 2d7e227d114a60ab47456bfc62308b16dba4eacf Олёш 2020-06-11 23:54:54
qt_plot: konečně něco interaktivního... 2bcc8842afe87d1e44e066fe168fcfedb8043947 Олёш 2020-06-07 07:35:57
WIP: matplotlib and QtGui ce3d38363636ba60fcddc0ebfe45c1cefadfc7a5 Олёш 2020-06-07 01:57:04
Commit 726058435dcf653f816390e3fc996747fff07908 - estimation: WIP
Author: Олёш
Author date (UTC): 2020-07-14 02:41
Committer name: Олёш
Committer date (UTC): 2020-07-14 02:41
Parent(s): 344f275e586e42c4946f6d4f18a2503ec5f28869
Signer:
Signing key:
Signing status: N
Tree: b2f9655c4670165366a14f5625ed29a20cfab579
File Lines added Lines deleted
estimation.py 521 0
simplex.py 61 0
File estimation.py changed (mode: 100644) (index ce420d6..c2d08d7)
... ... import numpy.ma as ma
10 10 import scipy.stats as stats import scipy.stats as stats
11 11
12 12 from scipy.spatial import cKDTree from scipy.spatial import cKDTree
13 from scipy.spatial import Delaunay
13 14 from scipy import interpolate from scipy import interpolate
14 15
15 16 import collections # for defaultdict import collections # for defaultdict
16 17
17 18 from . import f_models from . import f_models
19 from . import simplex as six
18 20 from .IS_stat import IS from .IS_stat import IS
19 21 from .candybox import CandyBox from .candybox import CandyBox
20 22
23
21 24 # mizí národní slovička :( # mizí národní slovička :(
22 25 # gradient musí být funkcí! # gradient musí být funkcí!
23 26 def progress(estimations, sampled_plan_R, Z, gradient, f_i, nis=500, budget=20000, L2_metric=False): def progress(estimations, sampled_plan_R, Z, gradient, f_i, nis=500, budget=20000, L2_metric=False):
 
... ... def Voronoi_2_point_estimation(sample_box, model_space='Rn', sampling_space=None
474 477 # near_neighbors[ii[k]] = near_neighbors[ii[k]] + 1 # near_neighbors[ii[k]] = near_neighbors[ii[k]] + 1
475 478 # Vor_mask[k] = failsi[ii[k]] # Vor_mask[k] = failsi[ii[k]]
476 479
480
481
482
483
484
485
486 def simple_simplex_estimation(sample_box, model_space='Rn', sampling_space=None, budget=20000, callback=None):
487 """
488 Delaunay triangulation
489 budget=20000
490 """
491
492 if callback is None:
493 callback = lambda *_, **__: None
494
495 nsim = sample_box.nsim
496 nis = max(round(budget/nsim), 100)
497
498 # jsou to informace pro callback
499 estimation={'method': "simple_simplex_estimation", 'nis':nis, 'nsim':nsim}
500 estimation['model_space'] = model_space
501 estimation['sampling_space'] = sampling_space
502
503
504 # vytahneme ze sample_boxu rozdělení
505 f = sample_box.sampled_plan
506
507 # já vím, že sample box pokážde failsi přepočítavá
508 failsi = sample_box.failsi
509
510
511
512 # zde provadím rozdělení na prostor, ve kterém vzorkujem
513 # a prostor "modelu", vô ktôrom, v podstatě, měříme vzdaleností
514 sampled_plan_model = getattr(sample_box, model_space)
515
516 # žádné chyby nechytám
517 tri = spatial.Delaunay(sampled_plan_model)
518 if len(tri.coplanar): # pokud triangulace je v pořadku
519 #print('triangulace v pořádku není')
520 bx._log('triangulation is coplanar')
521 else:
522 #print('triangulace je v pořádku')
523 bx._log('triangulation is OK')
524
525
526
527 if sampling_space is None:
528 sampling_space = model_space
529 # sing like sampling
530 sampled_plan_sing = sampled_plan_model
531 else:
532 sampled_plan_sing = getattr(sample_box, sampling_space)
533
534 sampled_plan_sing_ma = np.ma.asarray(sampled_plan_sing)
535
536 if sampling_space is None:
537 sampling_space = model_space
538 # sing like sampling
539 sampled_plan_sing = sampled_plan_model
540 tree_sampling = tree
541 else:
542 sampled_plan_sing = getattr(sample_box, sampling_space)
543 # narychlo, moc to nepomůže, neboť asi po 16 vzorcích počítá brut forsem
544 tree_sampling = cKDTree(sampled_plan_sing, compact_nodes=False, balanced_tree=False)
545
546 # find distance to the nearest sampling point (from all points)
547 dd2, ii2 = tree_sampling.query(sampled_plan_sing, k=[2], p=p_norm)
548 mindist_sing = dd2.flatten()
549
550 # bacha, rozdíl! Zde hrajeme s hustotou pro zadaní vzorkovacího rozdělení
551 # proto i hustoty bereme z vzorkovacího prostoru
552 PDF_sing = sample_box.pdf(sampling_space)
553
554
555 global_stats = collections.defaultdict(int)
556
557
558
559 # zde postupně v cyklu prochazíme všemi simplexy
560 # tynhlenstím zajišťujeme disjunktnost
561 # a môžeme všechny nasbírané pravděpodobnosti jednoduše sčítat
562 for simplex in tri.simplices:
563
564 vertices = sampled_plan_sing[simplex]
565
566 # hned ošidíme čístou matematiku našimi hustotama (jen tak, pro radost, nikdo nás nenutí...)
567 masses = PDF_sing[simplex] / np.mean(PDF_sing[simplex])
568 barycenter = np.average(vertices, axis=0, weights=masses)
569 I_bc = six.simplex_barycenter_inertia_tensor(vertices-barycenter, masses=masses)
570
571 # matika
572 w, v = np.linalg.eig(I_bc)
573
574 # use IS sampling density with center equal to the simplex's barycenter
575 # set the minimum distance as the standard deviation of IS densisty
576 # u stats.norm zadáváme směrodatnou odchylku, to je asi správné
577 sigmas = np.sqrt(w)
578 h_i = [stats.norm(0, sigmas[j]) for j in range(sample_box.nvar)]
579 # rozdělení ve vlastním prostoru
580 # select nis = 100 points from IS density
581 h_L = f_models.UnCorD(h_i)(nis)
582
583 # здесь уже так легко не отделаемся. Трансформовать кароно.
584 h_plan_bc = (v @ h_L.R.T).T
585 h_plan_sing = h_plan_bc + barycenter
586
587 # sice to má nazev h_plan, ale nese rozdělení a hustoty v f-ku
588 h_plan = f.new_sample(h_plan_sing, sampling_space)
589 # mě příjde, že je to legalní
590 weights_sim = to_sample.pdf(sampling_space) / h.pdf('R') # snad je to správně
591
592 # součet váh nemá cenu kontrolovat, jedná tam nebude, nevyjde
593
594
595 h_plan_model = getattr(h_plan, model_space)
596 dd, ii = tree.query(h_plan_model, k=1, p=p_norm)
597
598 # nechám s velkým písmenkem
599 Vor_mask = ii==i
600 h_plan_model_ma = h_plan_model[Vor_mask]
601 weights_sim = h_plan.w[Vor_mask]
602
603
604 cell_stats = dict()
605 # musí sa (na konci) rovnat jedne
606 # opravdu dělíme nis'em, jako v normálním IS
607 # nikoliv počtem příjatých bodíků,
608 # protože průměrná vaha je o hodně mene významná metrika
609 cell_stats['cell_probability'] = np.sum(weights_sim) / nis
610
611
612 # indexy ii nás moc nezajimajou
613 # vzdalenosti snad byjsme zvladli použit?
614 dd1 = dd[Vor_mask]
615
616 dd2, ii2 = tree.query(h_plan_model_ma, k=[2], p=p_norm)
617 dd2 = dd2.flatten()
618 ii2 = ii2.flatten()
619
620 # tahle hračka s indexy je pro numpy poměrně drahá
621 failsii_2 = failsi[ii2]
622
623 # jeden vzorek (včetně hustoty PDF[i]) je nám vždy znám
624 # porucha
625 if failsi[i]:
626 points_1 = PDF[i] * dd2
627 node_pf_estimations = points_1 / (points_1 + PDF[ii2] * dd1)
628 node_pf_pure_estimations = dd2 / (dd1 + dd2)
629
630 cell_stats['Voronoi_2_point_upper_bound'] = cell_stats['cell_probability']
631 cell_stats['Voronoi_2_point_failure_rate'] = np.sum(weights_sim*np.where(failsii_2,1, node_pf_estimations)) / nis
632 cell_stats['Voronoi_2_point_pure_failure_rate'] = np.sum(weights_sim*np.where(failsii_2,1, node_pf_pure_estimations)) / nis
633 cell_stats['Voronoi_2_point_lower_bound'] = np.sum(weights_sim[failsii_2]) / nis
634 cell_stats['Voronoi_failure_rate'] = cell_stats['cell_probability']
635 nodes=CandyBox(h_plan[Vor_mask], w=h_plan.w[Vor_mask], node_pf_estimations=node_pf_estimations,\
636 node_pf_pure_estimations=node_pf_pure_estimations)
637
638 # neporucha
639 else:
640 dd1 = dd1[failsii_2]
641 dd2 = dd2[failsii_2]
642 points_1 = PDF[i] * dd2
643 points_2 = PDF[ii2[failsii_2]] * dd1
644
645 node_pf_estimations = points_2 / (points_1 + points_2)
646 node_pf_pure_estimations = dd1 / (dd1 + dd2)
647
648 cell_stats['Voronoi_2_point_upper_bound'] = np.sum(weights_sim[failsii_2]) / nis
649 cell_stats['Voronoi_2_point_failure_rate'] = np.sum(weights_sim[failsii_2]*node_pf_estimations) / nis
650 cell_stats['Voronoi_2_point_pure_failure_rate'] = np.sum(weights_sim[failsii_2] * node_pf_pure_estimations) / nis
651 cell_stats['Voronoi_2_point_lower_bound'] = 0
652 cell_stats['Voronoi_failure_rate'] = 0
653 nodes=CandyBox(h_plan[Vor_mask][failsii_2], w=weights_sim[failsii_2], node_pf_estimations=node_pf_estimations,\
654 node_pf_pure_estimations=node_pf_pure_estimations)
655
656
657
658
659
660 for key, value in cell_stats.items():
661 global_stats[key] += value
662
663 # kolbek ↓
664 callback(estimation=estimation, nodes=nodes, cell_stats=cell_stats, out_nodes=h_plan[~Vor_mask])
665
666
667
668
669
670
671
672 def filter(bx, candidates=[]):
673 """
674 supports both sample and sample.R input
675
676 logika metody, nebo, přesněji, implementaci
677 je taková, že bule-li někdo-něco mimo doménu,
678 tak funkce je vrátí a zbytek už neřeší
679 """
680 # co to bylo na vstupu?
681 # když nebyl žádný,
682 # projdeme vlastními kandidaty
683 if len(candidates)==0:
684 candidates = bx.candidates
685
686
687 # tady byl problém. Funkce byla původně navržena tak,
688 # aby ji nezajimalo co je na vstupu
689 # to ale nefunguje
690 # další funkce jako výstup očekavají něco s validním R-kem
691 candidates_to_sample_node = getattr(bx.f.new_sample(candidates), bx.tri_space)
692
693
694 found_simplices = bx.tri.find_simplex(candidates_to_sample_node)
695 # ouside of domain - it's easy
696 outside = candidates[found_simplices < 0]
697 if len(outside) > 0:
698 # my hodnotili svých kandidatov?
699 if bx.candidates == candidates:
700 bx.candidates = outside
701 return outside, len(candidates)/len(outside)
702
703 # tady já chcu vrátit první vhodný vzorek a tím končít
704 for candidate_id in range(len(candidates)):
705 # simplex_id >= 0: # inside of domain
706 simplex_id = found_simplices[candidate_id]
707 simplex = bx.tri.simplices[simplex_id]
708
709 # fp like a failure points. Number of failure points
710 fp = len(np.setdiff1d(simplex, bx.failure_points))
711
712 # pokud je simplex není jednobarevny..
713 if (fp != 0) and (fp != bx.nvar+1):
714 # my hodnotili svých kandidatov?
715 if bx.candidates == candidates:
716 bx.candidates = candidates[candidate_id:]
717 return candidates[candidate_id], candidate_id
718
719 # nepovedlo. nic
720 # mě nenapadá žádný lepší způsob vrátit prázdnou matici
721 return candidates[0:0], len(candidates)
722
723
724
725
726 # pro jistotu pridame
727 bx.simplex_index = {'failure':[], 'success':[], 'mix':[]}
728
729
730 # overall estimations
731 #bx.oiss = IS_stat.ISSI(['failure', 'success', 'out', 'mix'])
732 # -1 = 'out', 0=success, 1=failure, 2=mix
733 bx.oiss = IS_stat.ISSI([-1,0,1,2])
734 # current estimations
735 bx.ciss = IS_stat.ISSI([-1,0,1,2])
736
737
738 # dropneme (pro jistotu) odhady
739 bx.oiss = bx.ciss
740 # -1 = 'out', 0=success, 1=failure, 2=mix
741 bx.ciss = IS_stat.ISSI([-1,0,1,2])
742
743 # drop indexes
744 bx.simplex_index = {'failure':[], 'success':[], 'mix':[]}
745
746 def increment(bx, input_sample):
747 super().increment(input_sample)
748
749 # drop indexes
750 bx.simplex_index = {'failure':[], 'success':[], 'mix':[]}
751
752 # current estimations
753 try: # to čo já vidím v kódu - ISSI slovníky se pokažde generujóu znovu,
754 # není nutně je explicitně kopirovať
755 bx.guessbox.guess('TRI_overall_estimations', bx.nsim-1, bx.oiss.estimations)
756 bx.guessbox.guess('TRI_current_estimations', bx.nsim-1, bx.ciss.estimations)
757
758 # a znovu začneme počítat
759 # -1 = 'out', 0=success, 1=failure, 2=mix
760 bx.ciss = IS_stat.ISSI([-1,0,1,2])
761
762
763
764 # chcu zachytit spadnuti QHull na začatku, kdy ještě není
765 # dostatek teček. Je-li bx.tri fakticky existuje, tj.
766 # triangulace jíž existovala - je třeba nechat QHull spadnout
767 except AttributeError:
768 choose = bx.LHS_like_correction(bx.get_candidates(1))
769 if bx.nsim < bx.nvar + 1: # je to legální
770 bx._log("we have no enough points to build triangulation, so", str(choose), "is our recommendation")
771 return choose
772
773 elif bx.nsim < 2*bx.nvar + 3: # to je ještě budiž
774 bx._log("we have troubles with triangulation, so we offer random sample for now:", str(choose))
775 return choose
776 else: # no to teda ne!
777 raise ValueError("AdaptiveCensoring: s tou triangulací je fakt něco není v pořadku")
778
779
780 def filter(bx, candidates=[]):
781 """
782 za pvré, jako vstup očekávám kandidaty od .get_candidates() funkce,
783 zabalené do полукустарного sample_boxu s zadaným .implicit_multiplicator
784 (je to drobnost pro přesnějši zpracování sad IS IS_statem).
785
786 Metoda musí souřádnicím přiřazovat jev
787 "success", "failure", "mix", "outside"
788
789 TATO metoda jakmile narazí na "mix" nebo "outside"
790 ukladá zjištěné informace do ISSI a nalezeného kandidata vrací
791 """
792 # co to bylo na vstupu?
793 # když nebyl žádný,
794 # projdeme vlastními kandidaty
795 if len(candidates)==0:
796 candidates = bx.candidates
797
798 bx._logi("kandidaty:", candidates)
799
800 # je třeba lokálně zachovat implicit_multiplicator
801 # jinak se ztrací při slajsingu
802 # nechceš přepsat SampleBox, Alexi?
803 try:
804 implicit_multiplicator = candidates.implicit_multiplicator
805 except AttributeError: # kandidaty můžou bejt odkudkoliv
806 # s nekonečným rozptylem nebudou mít váhu "absenční" jevy
807 # moc to odhadům nepomůže, protože je-li kandidaty
808 # nemajú .implicit_multiplicator
809 # asi nebudou mať ani váhy IS v .g_values
810 implicit_multiplicator = float("inf")
811 bx._logi("Dobrý den, kandidaty nemajú .implicit_multiplicator")#. S pozdravem, AdaptiveCensoring")
812
813 # tady byl problém. Funkce byla původně navržena tak,
814 # aby ji nezajimalo co je na vstupu
815 # to ale nefunguje
816 # další funkce jako výstup očekavají něco s validním R-kem
817 # no tj. já zde provádím posouzení transformací z R-ka vstupních souřadnic
818 candidates_to_sample_node = getattr(bx.f.new_sample(candidates), bx.tri_space)
819
820
821 current_simplices = bx.tri.find_simplex(candidates_to_sample_node)
822
823
824 # tak bacha
825 # budeme přepísovat jevy in-place
826 found_simplices = np.ma.array(current_simplices.copy()) #.copy()
827 # nemaskované - obsahuji číslo simplexu
828 # maskované - číslo jevu
829 # -1 = 'out', 0=success, 1=failure, 2=mix
830 # tj. procházíme simplexy z náhodné sady vzorků,
831 # nahrazujeme čislo simplexu odpovidajicím mu jevem
832 # a skryváme ho
833 # pote ty "skryté", "projduté" vzorky využiváme k žískání odhadů
834
835 while len(current_simplices):# > 0:
836 bx._logi("current simplices", current_simplices)
837 # berem hned prvního kandidata
838 # a posuzujeme co je zač
839 simplex_id = current_simplices[0]
840 mask = found_simplices==simplex_id
841
842 if simplex_id < 0: # -1 means ouside
843 # berem kramle
844 break
845 elif simplex_id in bx.simplex_index['success']:
846 found_simplices[mask] = 0
847 elif simplex_id in bx.simplex_index['failure']:
848 found_simplices[mask] = 1
849 elif simplex_id in bx.simplex_index['mix']:
850 found_simplices[mask] = 2
851 # kramle
852 break
853 else: # no index information
854 # tady já chcu vrátit první vhodný vzorek a tím končít
855 # simplex_id >= 0: # inside of domain
856 # asi tady získavam množinu s čísly vrcholů
857 # kteří zakladají simplex
858 simplex = bx.tri.simplices[simplex_id]
859
860 # for debug
861 bx._logi("провал индексу", simplex_id, indent=2)
862
863 # fp like a failure points. Number of failure points
864 # setdiff "Return the unique values in ar1 that are not in ar2."
865 fp = len(np.setdiff1d(simplex, bx.failure_points))
866
867 if fp == bx.nvar+1: #
868 bx.simplex_index['success'].append(simplex_id)
869 found_simplices[mask] = 0
870 elif fp == 0:
871 bx.simplex_index['failure'].append(simplex_id)
872 found_simplices[mask] = 1
873 #print("failure simplex", simplex_id)
874 else:
875 bx.simplex_index['mix'].append(simplex_id)
876 found_simplices[mask] = 2
877 bx._logi("mixed simplex", simplex_id)
878 # bacha! kramle
879 break
880
881 # pridame do seznamu známého
882 found_simplices[mask] = np.ma.masked
883 # eště raz
884 cmask = current_simplices==simplex_id
885 # vyfiltrujeme
886 current_simplices = current_simplices[~cmask]
887
888
889
890 # zde je třeba перехватить ситуацию, куке одӥг но кандидат ӧвӧл
891 # нужно ли?
892 # if len(current_simplices) == 0:
893 # # nepovedlo. nic
894 # bx.candidates = candidates[0:0]
895 # # mě nenapadá žádný lepší způsob vrátit prázdnou matici
896 # return candidates[0:0], len(candidates), -2 # simple_id
897
898 # nemaskované, včetně současného kandidata (nevím proč) - ke kandidatům
899 # землю - крестьянам, фабрики - рабочим
900 # předpokladam, že kandidaty jsou se všim všudy
901 # s vahami (.g_value) a se svým .implicit_multiplicator'em
902 ## zde True hodnoty roušky - to co jíž bylo skryto
903 bx.candidates = candidates[~np.ma.getmaskarray(found_simplices)][1:] # toho prvního prečo nechcem
904 try: # na zacatku je tam prazdný f_model, kterému atribut pripsat nemůžeme
905 bx.candidates.implicit_multiplicator = implicit_multiplicator
906 except:
907 pass
908
909 # vrátíme kandidaty, všechny-ne všechny?
910 # малы одӥг гинэ? уг тодӥськы чик...
911 selected_candidate = candidates[~np.ma.getmaskarray(found_simplices)][:1] # chcem toho prvního
912
913
914 # odešleme ISSI
915 try:
916 # pridame do seznamu známého
917 # rouška musí zůstat z cyklu
918 # proboha, Alexi, co to je za roušku, co se tu děje?
919 # jakmile v tom hlavním cyklu nalezli jsme mix nebo outside
920 # my hned z cyklu vylezli a je neskryli - abychom je vzali jako kandidaty
921 # teď je však skryváme s tou "rouškou", co musela být před opuštěním cyklu nastavena
922 # tak ISSI bude mít možnost odhadovat i pravděpodobnosti mix a outside
923 found_simplices[mask] = np.ma.masked
924 # zde získáme True hodnoty roušek
925 # ukazatel
926 imask = found_simplices.mask
927 found_simplices.mask = ~imask # invertujem, dotkne to i samotnou imask
928 events = found_simplices.compressed()
929 print(candidates)
930 print(candidates.g_values)
931 print("imask", imask)
932 print(candidates.g_values[~imask])
933 print(events)
934 bx.oiss.add_IS_serie(candidates.g_values[~imask], events, implicit_multiplicator)
935 print("global estimation", bx.oiss.estimations)
936 bx.ciss.add_IS_serie(candidates.g_values[~imask], events, implicit_multiplicator)
937 print("current estimation", bx.ciss.estimations)
938 except AttributeError:
939 bx._logi("Это вы мне прислали неваженых кандидатов?")
940 except UnboundLocalError: # čo tu chybu způsobuje?
941 # asi nebyly žádné kandidaty (třeba hned na začátku)
942 assert len(candidates)==0 and len(bx.candidates)==0, "AdaptiveCensoring: Что за бурда с этими кандидатама?"
943 return candidates, 0, -2
944
945
946 # rate = kolík bylo - kolik zůstalo
947 return selected_candidate, len(candidates) - len(bx.candidates), simplex_id
948
949
950
951
952 def get_candidates(bx, Nsim=int(1e4)):
953 # -1 = 'out', 0=success, 1=failure, 2=mix
954
955 # не мудрствуя лукаво
956 user_pf = np.mean(bx.pf_lim)
957 try:
958 low_pf = bx.oiss.estimations[1] # failure
959 upper_pf = 1 - bx.ciss.estimations[0] # sucess
960 self_pf = (low_pf + upper_pf)/2
961 except AttributeError:
962 self_pf = 0.5
963
964 # bereme *mean* od svého a uživatelského odhadu
965 # minimum nejede
966 sampling_r, __ = bx.sball.get_r_iteration(np.mean((self_pf, user_pf)))
967 # asi tam bylo sampling_r/bx.base_r, že?
968 # u stats.norm zadáváme směrodatnou odchylku, je to asi správné
969 h = f_models.UnCorD([stats.norm(0, sampling_r/bx.base_r) for i in range(bx.nvar)])
970
971 # for IS_stats
972 svar = (sampling_r/bx.base_r)**2 # svar like sampling_variance
973 # něco takovýho bych nahrubo placnul
974 implicit_multiplicator = svar**bx.nvar * np.exp(bx.nvar/svar - bx.nvar)
975
976 #
977 # jdeme na to, koťě!
978 #
979
980 # zgenerujeme vzorky
981 # nic zajimavýho
982 h = h(Nsim)
983 # dropneme priliš vzdálené kandidaty
984 distance_from_zero = np.sum(h.R**2, axis=1)
985 mask = distance_from_zero < bx.drop_r
986
987 # a teď bacha!
988 # tady musíme provést jeden trik
989 to_sample = bx.f.new_sample(h.R[mask], 'G') # R-ko smerdžíme ako G-čko
990 w = to_sample.pdf_G / h.pdf_R # snad je to správně
991 # zabalme do boxu
992 candidates = SampleBox(to_sample, w, 'BlackBox internal samples and weights')
993 candidates.implicit_multiplicator = implicit_multiplicator
994 # vahy máme, zbytek už nejsou naši starosti
995 return candidates
996
997 return {**estimation, **global_stats}
File simplex.py added (mode: 100644) (index 0000000..df3ed14)
1 #!/usr/bin/env python
2 # coding: utf-8
3
4 import numpy as np
5
6
7 # unbelivable: I was unable to find a package to calculate the second moments of an simplex
8 def points_inertia_tensor(vertices, masses=1):
9 """
10 coordinates of vertices
11 """
12 nsim, nvar = np.shape(vertices)
13
14 inertia_tensor = np.empty((nvar, nvar))
15 for i in range(nvar):
16 for j in range(i + 1):
17 if i==j:
18 inertia_tensor[i,j] = np.sum( masses * (np.sum(np.square(vertices), axis=1) - np.square(vertices[:, i])))
19 else:
20 inertia_tensor[i,j] = inertia_tensor[j,i] = - np.sum(masses * vertices[:, i]*vertices[:, j])
21 return inertia_tensor
22
23
24
25
26 def simplex_volume(vertices):
27 """
28 coordinates of vertices
29 """
30 nsim, nvar = np.shape(vertices)
31
32 return abs(np.linalg.det(vertices[1:] - vertices[0])) / np.math.factorial(nvar)
33
34
35
36 def simplex_barycenter_inertia_tensor(vertices, masses=1):
37 """
38 Returns the inertia matrix relative to the center of mass
39 coordinates of vertices
40 """
41 nsim, nvar = np.shape(vertices)
42 return points_inertia_tensor(vertices, masses) /(nvar+1)/(nvar+2) * simplex_volume(vertices)
43
44
45
46 def simplex_inertia_tensor(vertices, masses=1):
47 """
48 coordinates of vertices
49 """
50 nsim, nvar = np.shape(vertices)
51 masses = masses * np.append(np.full(nsim, 1/(nvar+1)/(nvar+2)), (nvar+1)/(nvar+2))
52
53 barycenter = np.mean(vertices, axis=0)
54 # barycenter beztak sepri4te ke každemu řádku tensoru
55 #_tensor = points_inertia_tensor(vertices, masses)/(nvar+1) + np.diag(np.square(barycenter))#*(nvar+1))
56 _tensor = points_inertia_tensor(np.vstack((vertices, barycenter)), masses=masses)
57 #return _tensor / (nvar+2) * simplex_volume(vertices)
58 return _tensor * simplex_volume(vertices)
59
60
61
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