File simplex.py changed (mode: 100644) (index d74b680..5f20ce6) |
3 |
3 |
|
|
4 |
4 |
import numpy as np |
import numpy as np |
5 |
5 |
from . import IS_stat |
from . import IS_stat |
|
6 |
|
from . import sball |
|
7 |
|
from . import f_models |
6 |
8 |
from scipy import spatial |
from scipy import spatial |
|
9 |
|
from scipy import stats |
|
10 |
|
from .candybox import CandyBox |
|
11 |
|
|
|
12 |
|
|
|
13 |
|
|
|
14 |
|
|
|
15 |
|
|
|
16 |
|
#č napadlo mě zababáchnuť třidu, která by se sama starala o všem co se tyče |
|
17 |
|
#č vnější domény. Nešlo mě totíž to udělat jednou funkcí, bylo by velmi |
|
18 |
|
#č špatné z hlediska zodpovednosti kódu. Tak to všecko zabalíme to třidy |
|
19 |
|
#č a odebereme z už beztak přetíženého blackboxu část komplexity |
|
20 |
|
# keywords: ISSI, estimation, outside, ConvexHull, Sball, IS kolem středních hodnot |
|
21 |
|
class Shull: # issi_estimate_outside |
|
22 |
|
def __init__(sx, f, nis, model_space, sampling_space=None, powerset_correction=True): |
|
23 |
|
#č tím powerset_corretion je myšlena úplná soustava jevů, |
|
24 |
|
#č tj. vyrovnaní s použitím míry vnější i vnitřní |
|
25 |
|
#č powerset_correction=True přídá -2 (inside) jev do ISSI |
|
26 |
|
|
|
27 |
|
#č zde f-ko musí taky obsahovat vzorky! |
|
28 |
|
sx.f = f |
|
29 |
|
sx.budget = nis |
|
30 |
|
sx.model_space = model_space |
|
31 |
|
if sampling_space is None: |
|
32 |
|
sx.sampling_space = model_space |
|
33 |
|
else: |
|
34 |
|
sx.sampling_space = sampling_space |
|
35 |
|
|
|
36 |
|
sampled_plan_model = getattr(f, model_space) |
|
37 |
|
#č žádná kontrola chyb - nechť to spadné, když má spadnout! |
|
38 |
|
sx.convex_hull = spatial.ConvexHull(sampled_plan_model, incremental=True) |
|
39 |
|
|
|
40 |
|
# current outside probability estimation |
|
41 |
|
sx.p_out = 0.5 # implicit value |
|
42 |
|
sx.sball = sball.Sball(f.nvar) |
|
43 |
|
sx.base_r = sx.sball.get_r(0.5) # we want in average 50/50 ratio |
|
44 |
|
|
|
45 |
|
|
|
46 |
|
# -1 = 'outside' |
|
47 |
|
# -2 = 'inside' |
|
48 |
|
#č založme ISSI |
|
49 |
|
sx.powerset_correction = powerset_correction |
|
50 |
|
#č potřebuji pro korektnost mít před integrací zadané jevy |
|
51 |
|
if powerset_correction: |
|
52 |
|
sx.oiss = IS_stat.ISSI([-1, -2]) |
|
53 |
|
else: |
|
54 |
|
sx.oiss = IS_stat.ISSI([-1]) |
|
55 |
|
|
|
56 |
|
|
|
57 |
|
|
|
58 |
|
def increment(sx, input_sample): |
|
59 |
|
#č sample by měl byt jíž převeden na f (v .add_sample()), |
|
60 |
|
#č zodpovidá za to volajicí kód! |
|
61 |
|
sx.convex_hull.add_points(getattr(input_sample, sx.model_space)) |
|
62 |
|
|
|
63 |
|
|
|
64 |
|
|
|
65 |
|
|
|
66 |
|
def integrate(sx): |
|
67 |
|
# getting rid of the old estimations |
|
68 |
|
sx.oiss.delete_event_data(-1) |
|
69 |
|
sx.oiss.events.append(-1) #č už tak trošku sahám do vnitřku cizí třidy |
|
70 |
|
if sx.powerset_correction: |
|
71 |
|
sx.oiss.delete_event_data(-2) |
|
72 |
|
sx.oiss.events.append(-2) #č a záse |
|
73 |
|
|
|
74 |
|
|
|
75 |
|
# first step |
|
76 |
|
nodes = sx.sample_sball() |
|
77 |
|
mask = nodes.is_outside |
|
78 |
|
|
|
79 |
|
|
|
80 |
|
cut_off = int(sx.budget/3) |
|
81 |
|
#č robím cyklus dokud nesberu dostatečně teček. |
|
82 |
|
#č To je fakt nejrobustnější řešení, co mě napadá |
|
83 |
|
# while (number_of_out_nodes or number_of_nodes_inside is too little) |
|
84 |
|
while (len(mask[mask]) < cut_off) or (len(mask[~mask]) < cut_off): |
|
85 |
|
#č je třeba jenom sehnat dostatečně bodíků a utikat |
|
86 |
|
nodes.add_sample(sx.sample_sball()) |
|
87 |
|
mask = nodes.is_outside |
|
88 |
|
|
|
89 |
|
|
|
90 |
|
return nodes |
|
91 |
|
|
|
92 |
|
|
|
93 |
|
|
|
94 |
|
def sample_sball(sx): |
|
95 |
|
nvar = sx.f.nvar |
|
96 |
|
sampling_r, sx.p_out = sx.sball.get_r_iteration(sx.p_out) |
|
97 |
|
#č asi tam bylo sampling_r/base_r, že? |
|
98 |
|
#č u stats.norm zadáváme směrodatnou odchylku, je to asi správné |
|
99 |
|
h = f_models.UnCorD([stats.norm(0, sampling_r/sx.base_r) for i in range(nvar)]) |
|
100 |
|
nodes = IS_stat.IS(sx.f, h, space_from_h='R', space_to_f=sx.sampling_space, Nsim=sx.budget) |
|
101 |
|
|
|
102 |
|
#č indikatorová funkce |
|
103 |
|
sx.is_outside(nodes) |
|
104 |
|
|
|
105 |
|
# for IS_stats |
|
106 |
|
if sx.powerset_correction: |
|
107 |
|
#č získáme výrovnaný odhad - je to skoro zdarma |
|
108 |
|
#svar = (sampling_r/sx.base_r)**2 # svar like sampling_variance |
|
109 |
|
#č kdysi snažil jsem něco odvést, moc se mi to nepovedlo |
|
110 |
|
#č je to jen tak, jeden z pokusu, hrubej nastřel |
|
111 |
|
#im = svar**nvar * np.exp(nvar/svar - nvar) |
|
112 |
|
|
|
113 |
|
#č radší ne. IM špatně zachycuje nizkou důvěru k tomu, co nemá vlastní tečky |
|
114 |
|
sx.oiss.add_IS_serie(nodes.w, nodes.event, implicit_multiplicator=np.inf) |
|
115 |
|
outside_measure = sx.oiss.estimations[-1] |
|
116 |
|
else: |
|
117 |
|
weights = nodes.w[nodes.is_outside] |
|
118 |
|
#č IM všecko pokazí, jakmile začnu přídávát další jevy |
|
119 |
|
sx.oiss.add_single_event_data(weights, event=-1, nis=sx.budget) |
|
120 |
|
# add_single_event_data() do not calculate estimations itself |
|
121 |
|
weighted_mean, __, events = sx.oiss.get_means() |
|
122 |
|
# outside probability |
|
123 |
|
#č nevyrovnané! |
|
124 |
|
outside_measure = weighted_mean[events==-1][0] |
|
125 |
|
|
|
126 |
|
#č pro přiště |
|
127 |
|
sx.p_out = (sx.p_out + outside_measure) / 2 |
|
128 |
|
|
|
129 |
|
return nodes |
|
130 |
|
|
|
131 |
|
|
|
132 |
|
def is_outside(sx, nodes): |
|
133 |
|
node_coordinates = getattr(nodes, sx.model_space) |
|
134 |
|
mask = is_outside(sx.convex_hull, node_coordinates) |
|
135 |
|
|
|
136 |
|
# -1 = 'outside' |
|
137 |
|
# -2 = 'inside' |
|
138 |
|
nodes.event = mask - 2 |
|
139 |
|
nodes.is_outside = mask |
|
140 |
|
|
|
141 |
|
|
|
142 |
|
|
|
143 |
|
|
|
144 |
|
|
|
145 |
|
|
|
146 |
|
|
|
147 |
|
|
|
148 |
|
|
|
149 |
|
|
7 |
150 |
|
|
8 |
151 |
|
|
9 |
152 |
#č tato metoda je vlastně pro MinEnergyCensoredSampling |
#č tato metoda je vlastně pro MinEnergyCensoredSampling |
|
... |
... |
def sample_simplex(vertices, model_space, sampling_space, nis): |
188 |
331 |
|
|
189 |
332 |
|
|
190 |
333 |
|
|
191 |
|
|
|
192 |
334 |
def is_outside(convex_hull, node_coordinates): |
def is_outside(convex_hull, node_coordinates): |
193 |
335 |
|
|
194 |
336 |
x = node_coordinates |
x = node_coordinates |