File simplex.py changed (mode: 100644) (index b41bdcc..5101b08) |
... |
... |
class FullCubatureTriangulation(_FullTriangulation, _CubatureTriangulation): |
799 |
799 |
|
|
800 |
800 |
|
|
801 |
801 |
|
|
|
802 |
|
#č Triangulation třída byla navržena s těsnou vazbou na Shull in mind. |
|
803 |
|
#č Snahou bylo vyrovnání odhadů, získaných při hojném využití IS |
|
804 |
|
#č Teď ale, chceme-li použit Ghull + quadpy, potřebujem třídu bez vázby na Shull |
|
805 |
|
#č Dávám to do zvláštní třídy jen kvůli logickému rozdělení kódu |
|
806 |
|
#č Jinak by se třída nijak nelišila od FastCubatureTriangulation |
|
807 |
|
class JustCubatureTriangulation(_FastTriangulation, _CubatureTriangulation): |
|
808 |
|
#č Tahle třída stala se komplikovanou jako šroub. |
|
809 |
|
#č zdědené metody, náhled z hlediska vazby na ISSI: |
|
810 |
|
# implicitly inherited from _Triangulation by _FastTriangulation: |
|
811 |
|
# setup - OK (implicitně zadává se issi=None) |
|
812 |
|
# get_events - OK |
|
813 |
|
# update - OK (co vidím, není vazán na Shull, řeší pouze triangulaci) |
|
814 |
|
# is_mixed - OK |
|
815 |
|
|
|
816 |
|
# inherited from _FastTriangulation(_Triangulation) itself: |
|
817 |
|
# integrate - OK |
|
818 |
|
# !get_pf_estimation - neOK |
|
819 |
|
# _invalidate_simplex - OK |
|
820 |
|
# estimate_simplex - OK |
|
821 |
|
|
|
822 |
|
# inherited from _CubatureTriangulation: |
|
823 |
|
# __init__ - OK (implicitně zadává se issi=None) |
|
824 |
|
# integrate_simplex - OK |
|
825 |
|
|
|
826 |
|
#č to je právý pf_estimation |
|
827 |
|
#č jiné třídy ve skutku vracej tri_estimation |
|
828 |
|
def get_pf_estimation(sx): |
|
829 |
|
# TRI-compatible estimation |
|
830 |
|
# -1=outside, 0=success, 1=failure, 2=mix |
|
831 |
|
#č ISSI nepouživáme, outside (-1), ani success (1) nebudou korektní |
|
832 |
|
tri_estimation = {-1:0, 0:0, 1:0, 2:0} |
|
833 |
|
|
|
834 |
|
#č něco konkretnějšího |
|
835 |
|
vertex_estimation = 0 |
|
836 |
|
weighted_vertex_estimation = 0 |
|
837 |
|
|
|
838 |
|
#č nechce se mi počitat přes numpy: np.array(tuple(c.values())) |
|
839 |
|
# let's iterate over all simplices we have |
|
840 |
|
# (calling code is responsible for the .simplex_stats validity) |
|
841 |
|
for event_id, simplex_measure, fm, wfm in sx.simplex_stats.values(): |
|
842 |
|
tri_estimation[event_id] += simplex_measure |
|
843 |
|
vertex_estimation += fm |
|
844 |
|
weighted_vertex_estimation += wfm |
|
845 |
|
|
|
846 |
|
|
|
847 |
|
#ё так, для красоты |
|
848 |
|
global_stats = dict() |
|
849 |
|
global_stats['failure'] = tri_estimation[1] |
|
850 |
|
global_stats['mix'] = tri_estimation[2] |
|
851 |
|
|
|
852 |
|
return {'TRI_estimation': tri_estimation, 'global_stats': global_stats, \ |
|
853 |
|
'vertex_estimation' : vertex_estimation, \ |
|
854 |
|
'weighted_vertex_estimation' : weighted_vertex_estimation, |
|
855 |
|
'coplanar':sx.tri.coplanar} |
802 |
856 |
|
|
803 |
857 |
|
|
804 |
858 |
|
|
|
... |
... |
class FullCubatureTriangulation(_FullTriangulation, _CubatureTriangulation): |
808 |
862 |
|
|
809 |
863 |
|
|
810 |
864 |
|
|
811 |
|
|
|
812 |
|
|
|
813 |
|
# CUbature, one more time |
|
814 |
|
# |
|
815 |
|
#class OneShotTriangulation: |
|
816 |
|
# def setup(sx, sample_box, tri_space='Rn', issi=None, weighting_space=None, \ |
|
817 |
|
# incremental=True, on_add_simplex=None, on_delete_simplex=None): |
|
818 |
|
# |
|
819 |
|
# sx.sample_box = sample_box |
|
820 |
|
# sx.tri_space = tri_space |
|
821 |
|
# sx.issi = issi #č ISSI potřebujem pro tri_estimation |
|
822 |
|
# |
|
823 |
|
# if weighting_space is None: |
|
824 |
|
# sx.weighting_space = tri_space |
|
825 |
|
# else: |
|
826 |
|
# sx.weighting_space = weighting_space |
|
827 |
|
# |
|
828 |
|
# #č kolbeky |
|
829 |
|
# sx.on_add_simplex = on_add_simplex |
|
830 |
|
# sx.on_delete_simplex = on_delete_simplex |
|
831 |
|
# |
|
832 |
|
# #оӵ кылсузъет кылдытом |
|
833 |
|
# sx.simplex_stats = dict() |
|
834 |
|
# #č ty množiny jsou tak trošku overkill, ale budiž |
|
835 |
|
# sx.simplices_set = set() |
|
836 |
|
# |
|
837 |
|
# # create .tri triangulation |
|
838 |
|
# #č tri - Deloneho triangulace |
|
839 |
|
# #č žádné chyby nechytám |
|
840 |
|
# #čs když se tringulace nepovede tak nemáme čo robiť |
|
841 |
|
# # incremental triangulation require one more point |
|
842 |
|
# tri_plan = getattr(sample_box, tri_space) |
|
843 |
|
# sx.tri = spatial.Delaunay(tri_plan, incremental=incremental) |
|
844 |
|
# if len(sx.tri.coplanar): |
|
845 |
|
# #print('triangulace v pořádku není') |
|
846 |
|
# print('Triangulation is coplanar:', sx.tri.coplanar) |
|
847 |
|
# else: |
|
848 |
|
# #print('triangulace je v pořádku') |
|
849 |
|
# print('Triangulation is OK') |
|
850 |
|
# |
|
851 |
|
# |
|
852 |
|
# |
|
853 |
|
# |
|
854 |
|
# def get_events(sx, simplices=None): |
|
855 |
|
# """ |
|
856 |
|
# Metoda musí simplexům přiřazovat jev |
|
857 |
|
# 0=success, 1=failure, 2=mix |
|
858 |
|
# """ |
|
859 |
|
# if simplices is None: |
|
860 |
|
# simplices = sx.tri.simplices |
|
861 |
|
# |
|
862 |
|
# in_failure = np.isin(simplices, sx.sample_box.failure_points) |
|
863 |
|
# has_failure = in_failure.any(axis=1) |
|
864 |
|
# all_failure = in_failure.all(axis=1) |
|
865 |
|
# return np.int8(np.where(has_failure, np.where(all_failure, 1, 2), 0)) |
|
866 |
|
# |
|
867 |
|
# |
|
868 |
|
# |
|
869 |
|
# #č tato funkce běží 91% času |
|
870 |
|
# # bottleneck function |
|
871 |
|
# def update(sx): |
|
872 |
|
# """ |
|
873 |
|
# Triangulace zajistěně existuje |
|
874 |
|
# """ |
|
875 |
|
# #č tato funkce je koncipována jinač |
|
876 |
|
# #č narozdíl od Shull, zde nám taky záleží |
|
877 |
|
# #č na poruchách i neporuchách |
|
878 |
|
# #č f_model proto nám stačit nebude |
|
879 |
|
# #č a bylo by blbě tečky brat zevnějšku, |
|
880 |
|
# #č a failsi - zevnitřku. |
|
881 |
|
# #č Takže - všechno bereme ze sample_box |
|
882 |
|
# #č reference, jenom dejte vědet, |
|
883 |
|
# #č že máme triangulaci aktualizovat |
|
884 |
|
# |
|
885 |
|
# former_simplices = sx.tri.simplices |
|
886 |
|
# former_nsimplex = sx.tri.nsimplex |
|
887 |
|
# |
|
888 |
|
# tri_plan = getattr(sx.sample_box, sx.tri_space) |
|
889 |
|
# #č jako vždy, chyby nechytáme |
|
890 |
|
# #sx.tri.add_points(tri_plan[sx.tri.npoints:]) |
|
891 |
|
# sx._tri_update(tri_plan) |
|
892 |
|
# |
|
893 |
|
# if len(sx.tri.coplanar): # pokud triangulace není v pořadku |
|
894 |
|
# #print('triangulace v pořádku není') |
|
895 |
|
# print('Triangulation has coplanar points:', sx.tri.coplanar) |
|
896 |
|
# |
|
897 |
|
# |
|
898 |
|
# # zkontrolujeme co se změnilo |
|
899 |
|
# if former_nsimplex <= sx.tri.nsimplex: |
|
900 |
|
# equal_mask = former_simplices == sx.tri.simplices[:former_nsimplex] |
|
901 |
|
# mask = ~equal_mask.all(axis=1) # unequal mask |
|
902 |
|
# sx._validate_simplices(former_simplices[mask]) |
|
903 |
|
# |
|
904 |
|
# |
|
905 |
|
# # změněné simplexy přepočítáme |
|
906 |
|
# to_estimate = sx.tri.simplices[:former_nsimplex][mask] |
|
907 |
|
# sx._estimate_simplices(to_estimate) |
|
908 |
|
# |
|
909 |
|
# # teď nové simplexy |
|
910 |
|
# # simplexy свежего разлива |
|
911 |
|
# sx._estimate_simplices(sx.tri.simplices[former_nsimplex:]) |
|
912 |
|
# |
|
913 |
|
# else: #č koplanar. Stavá se. |
|
914 |
|
# equal_mask = former_simplices[:sx.tri.nsimplex] == sx.tri.simplices |
|
915 |
|
# mask = ~equal_mask.all(axis=1) # unequal mask |
|
916 |
|
# sx._validate_simplices(former_simplices[sx.tri.nsimplex:]) |
|
917 |
|
# sx._validate_simplices(former_simplices[:sx.tri.nsimplex][mask]) |
|
918 |
|
# |
|
919 |
|
# sx._estimate_simplices(sx.tri.simplices[mask]) |
|
920 |
|
# |
|
921 |
|
# def _tri_update(sx, tri_plan): |
|
922 |
|
# #č separujeme, abychom vědeli, kolik času žere samotný QHull |
|
923 |
|
# sx.tri.add_points(tri_plan[sx.tri.npoints:]) |
|
924 |
|
# |
|
925 |
|
# def _validate_simplices(sx, to_validate): |
|
926 |
|
# |
|
927 |
|
# #č zde spolehám na to, že pořadí indexů se nikdy nezmění |
|
928 |
|
# #č tj. [12, 13, 26] najednou nestane [26, 12, 13] |
|
929 |
|
# #č (dá se něco takovýho očekavát podle toho co jsem čet v dokumentaci) |
|
930 |
|
# |
|
931 |
|
# #č u těch přečíslovaných zkolntrolujeme, zda fakt v té triangulaci nejsou |
|
932 |
|
# for simplex in to_validate: |
|
933 |
|
# #č když tam je |
|
934 |
|
# #č každopadně tohle je rychlejší než přepočet spousty simplexů |
|
935 |
|
# # few times faster |
|
936 |
|
# isin = (sx.tri.simplices == simplex).all(axis=1).any() |
|
937 |
|
# if not isin: |
|
938 |
|
# sx._invalidate_simplex(simplex) |
|
939 |
|
# |
|
940 |
|
# |
|
941 |
|
# |
|
942 |
|
# |
|
943 |
|
# |
|
944 |
|
# |
|
945 |
|
# def is_mixed(bx, simplices=None): |
|
946 |
|
# |
|
947 |
|
# if simplices is None: |
|
948 |
|
# simplices = bx.tri.simplices |
|
949 |
|
# |
|
950 |
|
# in_failure = np.isin(simplices, bx.failure_points) |
|
951 |
|
# has_failure = in_failure.any(axis=1) |
|
952 |
|
# all_failure = in_failure.all(axis=1) |
|
953 |
|
# return np.logical_xor(has_failure, all_failure) |
|
954 |
|
# |
|
955 |
|
# |
|
956 |
|
# |
|
957 |
|
# def integrate(sx): |
|
958 |
|
# #č Metoda musí simplexům přiřazovat jev |
|
959 |
|
# # 0=success, 1=failure, 2=mix |
|
960 |
|
# #č vyhodil jsem simplex_id'y |
|
961 |
|
# event_ids = sx.get_events() |
|
962 |
|
# |
|
963 |
|
# #č zde chceme ušetřít, a nechat stranou zelené simplexy |
|
964 |
|
# simplices = sx.tri.simplices[event_ids != 0] |
|
965 |
|
# event_ids = event_ids[event_ids != 0] |
|
966 |
|
# |
|
967 |
|
# |
|
968 |
|
# #č zde postupně v cyklu prochazíme simplexy |
|
969 |
|
# #č tynhlenstím zajišťujeme disjunktnost |
|
970 |
|
# #čs a môžeme všechny nasbírané pravděpodobnosti jednoduše sčítat |
|
971 |
|
# for simplex, event_id in zip(simplices, event_ids): |
|
972 |
|
# |
|
973 |
|
# #č ty množiny jsou tak trošku overkill, ale budiž |
|
974 |
|
# sx.simplices_set.add(tuple(simplex)) |
|
975 |
|
# |
|
976 |
|
# # -1 = 'outside', 0=success, 1=failure, 2=mix |
|
977 |
|
# event, fr, wfr = get_failure_ratio(sx.sample_box,\ |
|
978 |
|
# event_id, simplex, sx.weighting_space) |
|
979 |
|
# |
|
980 |
|
# sx.integrate_simplex(simplex, event, event_id, fr, wfr) |
|
981 |
|
# |
|
982 |
|
# |
|
983 |
|
# def get_pf_estimation(sx): |
|
984 |
|
# # TRI-compatible estimation |
|
985 |
|
# # -1=outside, 0=success, 1=failure, 2=mix |
|
986 |
|
# tri_estimation = {-1:0, 0:0, 1:0, 2:0} |
|
987 |
|
# |
|
988 |
|
# # Shull should be inicialized with powerset_correction=True |
|
989 |
|
# #č dostaneme vyrovnané odhady Brna-města (-2) a Brna-venkova (-1) |
|
990 |
|
# #č současný kód Shull zajišťuje, |
|
991 |
|
# #č že v ISSI estimátory budou spočítány |
|
992 |
|
# #sx.issi.get_estimations() |
|
993 |
|
# tri_estimation[-1] = sx.issi.estimations[-1] |
|
994 |
|
# |
|
995 |
|
# #č něco konkretnějšího |
|
996 |
|
# vertex_estimation = 0 |
|
997 |
|
# weighted_vertex_estimation = 0 |
|
998 |
|
# |
|
999 |
|
# pf_inside = sx.issi.estimations[-2] |
|
1000 |
|
# #č nechce se mi počitat přes numpy: np.array(tuple(c.values())) |
|
1001 |
|
# # let's iterate over all simplices we have |
|
1002 |
|
# # (calling code is responsible for the .simplex_stats validity) |
|
1003 |
|
# for event_id, simplex_measure, fm, wfm in sx.simplex_stats.values(): |
|
1004 |
|
# tri_estimation[event_id] += simplex_measure |
|
1005 |
|
# vertex_estimation += fm |
|
1006 |
|
# weighted_vertex_estimation += wfm |
|
1007 |
|
# |
|
1008 |
|
# #č success klidně může být i záporným |
|
1009 |
|
# tri_estimation[0] = pf_inside - tri_estimation[1] - tri_estimation[2] |
|
1010 |
|
# |
|
1011 |
|
# |
|
1012 |
|
# #ё так, для красоты |
|
1013 |
|
# global_stats = dict() |
|
1014 |
|
# global_stats['outside'] = 0 |
|
1015 |
|
# global_stats['success'] = 0 |
|
1016 |
|
# global_stats['failure'] = tri_estimation[1] |
|
1017 |
|
# global_stats['mix'] = tri_estimation[2] |
|
1018 |
|
# |
|
1019 |
|
# return {'TRI_estimation': tri_estimation, 'global_stats': global_stats, \ |
|
1020 |
|
# 'vertex_estimation' : vertex_estimation, \ |
|
1021 |
|
# 'weighted_vertex_estimation' : weighted_vertex_estimation} |
|
1022 |
|
# |
|
1023 |
|
# |
|
1024 |
|
# |
|
1025 |
|
# #č vyhodil jsem simplex_id'y |
|
1026 |
|
# def _invalidate_simplex(sx, indices): |
|
1027 |
|
# simplex = tuple(indices) |
|
1028 |
|
# |
|
1029 |
|
# #č simplex nemusí tam být, |
|
1030 |
|
# if simplex in sx.simplices_set: |
|
1031 |
|
# sx.simplices_set.remove(simplex) |
|
1032 |
|
# |
|
1033 |
|
# if simplex in sx.simplex_stats: |
|
1034 |
|
# sx.simplex_stats.pop(simplex) |
|
1035 |
|
# |
|
1036 |
|
# if sx.on_delete_simplex is not None: |
|
1037 |
|
# sx.on_delete_simplex(indices=indices) |
|
1038 |
|
# |
|
1039 |
|
# |
|
1040 |
|
# |
|
1041 |
|
# |
|
1042 |
|
# |
|
1043 |
|
# |
|
1044 |
|
# def _estimate_simplices(sx, to_estimate): |
|
1045 |
|
# #č zde spolehám na to, že pořadí indexů se nikdy nezmění |
|
1046 |
|
# #č tj. [12, 13, 26] najednou nestane [26, 12, 13] |
|
1047 |
|
# #č (dá se něco takovýho očekavát podle toho co jsem čet v dokumentaci) |
|
1048 |
|
# |
|
1049 |
|
# for simplex in to_estimate: |
|
1050 |
|
# #č ty množiny jsou tak trošku overkill, ale budiž |
|
1051 |
|
# if tuple(simplex) not in sx.simplices_set: |
|
1052 |
|
# sx.estimate_simplex(simplex) |
|
1053 |
|
# |
|
1054 |
|
# #č vyhodil jsem simplex_id'y |
|
1055 |
|
# def estimate_simplex(sx, indices): |
|
1056 |
|
# |
|
1057 |
|
# #č zkusím funkci návrhnout tak, že |
|
1058 |
|
# #ё вызывающая функция запускает estimate_simplex |
|
1059 |
|
# #ё на всём подряд без разбору. |
|
1060 |
|
# #č Našim úkolem je zjistit co je simplex zač |
|
1061 |
|
# #č a ty zelené ignorovat |
|
1062 |
|
# |
|
1063 |
|
# #č ty množiny jsou tak trošku overkill, ale budiž |
|
1064 |
|
# sx.simplices_set.add(tuple(indices)) |
|
1065 |
|
# |
|
1066 |
|
# # -1 = 'outside', 0=success, 1=failure, 2=mix |
|
1067 |
|
# event, event_id, fr, wfr = get_indices_event(sx.sample_box,\ |
|
1068 |
|
# indices, sx.weighting_space) |
|
1069 |
|
# |
|
1070 |
|
# if event_id != 0: |
|
1071 |
|
# sx.integrate_simplex(indices, event, event_id, fr, wfr) |
|
1072 |
|
# |
|
1073 |
|
# |
|
1074 |
|
# |
|
1075 |
|
# def __init__(sx, sample_box, tn_scheme, tri_space='Rn', issi=None, weighting_space=None, \ |
|
1076 |
|
# incremental=True, on_add_simplex=None, on_delete_simplex=None): |
|
1077 |
|
# |
|
1078 |
|
# sx.tn_scheme = tn_scheme |
|
1079 |
|
# |
|
1080 |
|
# sx.setup(sample_box, tri_space=tri_space, issi=issi,\ |
|
1081 |
|
# weighting_space=weighting_space, incremental=incremental,\ |
|
1082 |
|
# on_add_simplex=on_add_simplex,\ |
|
1083 |
|
# on_delete_simplex=on_delete_simplex) |
|
1084 |
|
# |
|
1085 |
|
# |
|
1086 |
|
# #č vyhodil jsem simplex_id'y |
|
1087 |
|
# def integrate_simplex(sx, indices, event_ids, frs, wfrs): |
|
1088 |
|
# print("estimate", indices) |
|
1089 |
|
# |
|
1090 |
|
# #оӵ чылкыт f_model |
|
1091 |
|
# f = sx.sample_box.f_model |
|
1092 |
|
# # quadpy |
|
1093 |
|
# points = [] # container #č vždyť Python nemá pointery |
|
1094 |
|
# def _get_pdf(x): |
|
1095 |
|
# fx = f.sample_pdf(x.T, sx.tri_space) |
|
1096 |
|
# #print('fx', fx) |
|
1097 |
|
# points.append((x.T, fx)) # side effect |
|
1098 |
|
# return fx |
|
1099 |
|
# |
|
1100 |
|
# # premature optimization |
|
1101 |
|
# if sx.on_add_simplex is not None: |
|
1102 |
|
# vertices = f[indices] |
|
1103 |
|
# vertices_model = getattr(vertices, sx.tri_space) |
|
1104 |
|
# else: |
|
1105 |
|
# vertices_model = getattr(f, sx.tri_space)[indices] |
|
1106 |
|
# #č pokud něčemu rotumím, tak |
|
1107 |
|
# # quad_vertices musejí bejt ve tvaru (nvertices=ndim+1, nsimplex, ndim) |
|
1108 |
|
# quad_vertices = |
|
1109 |
|
# simplex_measure = sx.tn_scheme.integrate(_get_pdf, vertices_model) |
|
1110 |
|
# |
|
1111 |
|
# fm = simplex_measure * fr |
|
1112 |
|
# wfm = simplex_measure * wfr |
|
1113 |
|
# |
|
1114 |
|
# #č ISSI tu nemáme, místo toho prostě ukladáme co máme do slovníku |
|
1115 |
|
# sx.simplex_stats[tuple(indices)] = (event_id, simplex_measure, fm, wfm) |
|
1116 |
|
# |
|
1117 |
|
# if sx.on_add_simplex is not None: |
|
1118 |
|
# cell_stats = dict() |
|
1119 |
|
# |
|
1120 |
|
# cell_stats['cell_probability'] = simplex_measure |
|
1121 |
|
# cell_stats['vertex_estimation'] = fm |
|
1122 |
|
# cell_stats['weighted_vertex_estimation'] = wfm |
|
1123 |
|
# cell_stats['event'] = event |
|
1124 |
|
# |
|
1125 |
|
# |
|
1126 |
|
# # kolbek ↓ |
|
1127 |
|
# nodes_model, pdf = points[0] |
|
1128 |
|
# nodes = f.new_sample(nodes_model, sx.tri_space) |
|
1129 |
|
# sx.on_add_simplex(sx, indices=indices, simplex=vertices, nodes=nodes, cell_stats=cell_stats) |
|
1130 |
|
# |
|
1131 |
|
# |
|
1132 |
|
# |
|
1133 |
|
# |
|
1134 |
|
# |
|
1135 |
|
# |
|
1136 |
865 |
# |
# |
1137 |
866 |
## global sample_box function |
## global sample_box function |
1138 |
867 |
##č tím globálným sample_box'em my šetříme čas na poměrně drahých slajséch |
##č tím globálným sample_box'em my šetříme čas na poměrně drahých slajséch |