File wellmet/shell.py added (mode: 100644) (index 0000000..d4c45fd) |
|
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 .candynodes import CandyNodes |
|
8 |
|
from .IS_stat import get_1DS_sample # for Shell_1DS |
|
9 |
|
|
|
10 |
|
from collections import namedtuple |
|
11 |
|
|
|
12 |
|
try: # try to outthink OS's memory management |
|
13 |
|
import os |
|
14 |
|
mem_bytes = os.sysconf('SC_PAGE_SIZE') * os.sysconf('SC_PHYS_PAGES') |
|
15 |
|
except BaseException as e: |
|
16 |
|
mem_GB = 16 |
|
17 |
|
mem_bytes = mem_GB * 1024**3 # hello, Windows! |
|
18 |
|
#print("ghull failed to get amount of RAM installed. %s GB is supposed."% mem_GB, repr(e)) |
|
19 |
|
#print("BTW, you are using Windows, aren't you?") |
|
20 |
|
|
|
21 |
|
|
|
22 |
|
|
|
23 |
|
ShellStats = namedtuple('ShellStats', ( |
|
24 |
|
'nsim', |
|
25 |
|
'nvar', |
|
26 |
|
'nfacets', |
|
27 |
|
'r', 'R', |
|
28 |
|
'inner', |
|
29 |
|
'shell', |
|
30 |
|
'outer', |
|
31 |
|
'FORM_outside', |
|
32 |
|
'TwoFORM_outside', |
|
33 |
|
'orth_outside' |
|
34 |
|
)) |
|
35 |
|
|
|
36 |
|
ShellEstimation = namedtuple('ShellEstimation', |
|
37 |
|
('shell_budget', 'shell_inside', |
|
38 |
|
'shell_outside', 'inside', 'outside')) |
|
39 |
|
|
|
40 |
|
|
|
41 |
|
#č pomocná třída pro integrování |
|
42 |
|
#č nejsem tedy jist, jestli je to nejlepší napad - |
|
43 |
|
#č dělit Ghull na podtřídy, delegovat funkcionalitu |
|
44 |
|
#č ale jinak Ghull se stavá hodně překomplikovaným. |
|
45 |
|
#č nic lepšího mně nenapadá, přemyšlel jsem dlouho. |
|
46 |
|
class Shell_1DS: |
|
47 |
|
"""1DS stands for 1D sampling. |
|
48 |
|
1DS is, actually, an importance sampling method |
|
49 |
|
that, actually, does not calculate IS weights |
|
50 |
|
as f_density / h_density, but |
|
51 |
|
(arbitrary, nonuniformly) divides interval to subintervals. |
|
52 |
|
By using CDF transformes subintervals to 0-1 measure line. |
|
53 |
|
Then gets one node as middle point of every subinterval, |
|
54 |
|
weights therefore are just interval widths itself. |
|
55 |
|
No sampling imprecisions are introduced, |
|
56 |
|
therefore no spring, no correction are needed. |
|
57 |
|
""" |
|
58 |
|
def __init__(self, hull, shell): |
|
59 |
|
self.shell = shell |
|
60 |
|
self.hull = hull |
|
61 |
|
self.nvar = hull.sample.nvar |
|
62 |
|
|
|
63 |
|
#č objevilo se, že scipy.stats.chi je neunosně nepřesné |
|
64 |
|
#č vzdává se jíž na poloměru 8 |
|
65 |
|
self.sball = sball.Sball(self.nvar) |
|
66 |
|
|
|
67 |
|
self.integration_cutoff = np.inf |
|
68 |
|
|
|
69 |
|
|
|
70 |
|
def integrate(self, nis, callback_all=None, callback_outside=None): |
|
71 |
|
self.reset(nis) |
|
72 |
|
|
|
73 |
|
if self.hull.nsimplex == 0: |
|
74 |
|
bus = self.integration_cutoff |
|
75 |
|
else: |
|
76 |
|
bus = int(mem_bytes / self.hull.nsimplex / 8 / 10) + 1 |
|
77 |
|
while self.nsampled < nis: |
|
78 |
|
|
|
79 |
|
seats = min(nis - self.nsampled, self.integration_cutoff, bus) |
|
80 |
|
try: |
|
81 |
|
self._process_part(seats, nis, callback_all, callback_outside) |
|
82 |
|
except MemoryError as m: |
|
83 |
|
print(self.__class__.__name__ +":", "memory error, %s sedaček" % seats, repr(m)) |
|
84 |
|
self.integration_cutoff = int(np.ceil(seats/2)) |
|
85 |
|
|
|
86 |
|
assert nis == self.nsampled |
|
87 |
|
|
|
88 |
|
return self._get_result() |
|
89 |
|
|
|
90 |
|
|
|
91 |
|
def reset(self, nis): # clear |
|
92 |
|
|
|
93 |
|
self.nsampled = 0 |
|
94 |
|
|
|
95 |
|
#č poloměry bereme ze skořapky |
|
96 |
|
#č za správné nastavení (stejně jako u MC) |
|
97 |
|
#č zodpovidá uživatel třídy |
|
98 |
|
#č třída vlastní odhad r nijak nevyuživá! |
|
99 |
|
r = self.shell.r |
|
100 |
|
R = self.shell.R |
|
101 |
|
# let's predefine 1D sequence at very beginning |
|
102 |
|
if r > np.sqrt(self.nvar - 1): |
|
103 |
|
x_sub = np.geomspace(r, R, nis+1, endpoint=True) |
|
104 |
|
else: |
|
105 |
|
x_sub = np.linspace(r, R, nis+1, endpoint=True) |
|
106 |
|
|
|
107 |
|
|
|
108 |
|
x, weights = get_1DS_sample(self.sball, x_sub) |
|
109 |
|
self.x = x |
|
110 |
|
self.weights = weights |
|
111 |
|
self.mask = np.empty(nis, dtype=np.bool) |
|
112 |
|
|
|
113 |
|
|
|
114 |
|
|
|
115 |
|
# bus analogy |
|
116 |
|
def _process_part(self, seats, nis, callback_all=None, callback_outside=None): |
|
117 |
|
# boarding |
|
118 |
|
left = self.nsampled |
|
119 |
|
right = self.nsampled + seats |
|
120 |
|
rs = self.x[left:right] |
|
121 |
|
|
|
122 |
|
# rand_dir: prepare ns random directions on a unit d-sphere |
|
123 |
|
rand_dir = sball.get_random_directions(seats, self.nvar) #random directions |
|
124 |
|
nodes_G = rand_dir*rs[:,None] |
|
125 |
|
nodes = self.hull.sample.f_model.new_sample(nodes_G, space='G') |
|
126 |
|
|
|
127 |
|
# who is back? |
|
128 |
|
d = self.hull.get_hyperplane_distances(nodes) |
|
129 |
|
mask = d > 0 |
|
130 |
|
assert len(mask) == seats |
|
131 |
|
|
|
132 |
|
# the most important part |
|
133 |
|
self.nsampled += len(mask) |
|
134 |
|
self.mask[left:right] = mask |
|
135 |
|
|
|
136 |
|
if callback_all is not None: |
|
137 |
|
# -2 = 'inside' -1 = 'outside' |
|
138 |
|
candy_nodes = CandyNodes(nodes, event_id=mask-2, is_outside=mask, d=d) |
|
139 |
|
callback_all(candy_nodes) |
|
140 |
|
|
|
141 |
|
if (callback_outside is not None) and np.any(mask): |
|
142 |
|
callback_outside(nodes[mask], d[mask]) |
|
143 |
|
|
|
144 |
|
|
|
145 |
|
|
|
146 |
|
def _get_result(self): |
|
147 |
|
#č mask related to hull.is_outside() |
|
148 |
|
shell_pf = np.sum(self.weights[self.mask]) |
|
149 |
|
shell_ps = np.sum(self.weights[~self.mask]) |
|
150 |
|
|
|
151 |
|
return self.nsampled, shell_ps, shell_pf |
|
152 |
|
|
|
153 |
|
|
|
154 |
|
|
|
155 |
|
|
|
156 |
|
class GaussianAnnulus: |
|
157 |
|
def __init__(self, hull): |
|
158 |
|
assert hull.space == 'G' |
|
159 |
|
|
|
160 |
|
self.hull = hull |
|
161 |
|
self.shell = sball.Shell(hull.sample.nvar) |
|
162 |
|
self.outside_dist = sball.Shell(hull.sample.nvar) |
|
163 |
|
self.sample = hull.sample |
|
164 |
|
|
|
165 |
|
self.gint = Shell_1DS(self.hull, self.shell) |
|
166 |
|
|
|
167 |
|
|
|
168 |
|
def boom(self, ns, use_MC=False): |
|
169 |
|
if use_MC: |
|
170 |
|
self.outside_dist.set_bounds(self.get_R()) |
|
171 |
|
nodes_G = self.outside_dist.rvs(ns) |
|
172 |
|
else: |
|
173 |
|
# rand_dir: prepare ns random directions on a unit d-sphere |
|
174 |
|
rand_dir = sball.get_random_directions(ns, self.sample.nvar) #random directions |
|
175 |
|
|
|
176 |
|
#č deme od vnější .get_R() kružnici směrem ven |
|
177 |
|
r = self.get_R() |
|
178 |
|
|
|
179 |
|
# maximum radius, where norm.pdf() wasn't zero |
|
180 |
|
# -38.575500173381374935388521407730877399444580 |
|
181 |
|
# don't ask me what the magic python use to distinguish |
|
182 |
|
# digits after double precision |
|
183 |
|
max_R_ever = 37 |
|
184 |
|
if r < max_R_ever: |
|
185 |
|
R = max_R_ever |
|
186 |
|
else: |
|
187 |
|
R = r + 10 |
|
188 |
|
r = np.linspace(self.get_R(), max_R_ever, ns, endpoint=True) |
|
189 |
|
nodes_G = rand_dir*r[:,None] |
|
190 |
|
|
|
191 |
|
nodes = self.sample.f_model.new_sample(nodes_G, space='G') |
|
192 |
|
return nodes |
|
193 |
|
|
|
194 |
|
|
|
195 |
|
def get_R(self): |
|
196 |
|
sum_squared = np.sum(np.square(self.sample.G), axis=1) |
|
197 |
|
#index = np.argmax(sum_squared) |
|
198 |
|
return np.sqrt(np.nanmax(sum_squared)) |
|
199 |
|
|
|
200 |
|
def setup_shell(self): |
|
201 |
|
shell = self.shell |
|
202 |
|
r = self.hull.get_r() |
|
203 |
|
R = self.get_R() |
|
204 |
|
|
|
205 |
|
if r<0: |
|
206 |
|
shell.set_bounds(0, R) |
|
207 |
|
else: |
|
208 |
|
shell.set_bounds(r, R) |
|
209 |
|
|
|
210 |
|
return r, R |
|
211 |
|
|
|
212 |
|
def get_shell_estimation(self): |
|
213 |
|
r, R = self.setup_shell() |
|
214 |
|
|
|
215 |
|
return ShellStats( |
|
216 |
|
self.sample.nsim, |
|
217 |
|
self.sample.nvar, |
|
218 |
|
self.hull.nsimplex, |
|
219 |
|
r, R, |
|
220 |
|
self.shell.ps, |
|
221 |
|
self.shell.p_shell, |
|
222 |
|
self.shell.pf, |
|
223 |
|
stats.norm.sf(r), |
|
224 |
|
self.hull.get_2FORM_outside(), |
|
225 |
|
self.hull.get_orth_outside() |
|
226 |
|
) |
|
227 |
|
|
|
228 |
|
def integrate(self, nis, callback_all=None, callback_outside=None): |
|
229 |
|
self.setup_shell() |
|
230 |
|
|
|
231 |
|
nsampled, shell_ps, shell_pf = self.gint.integrate(nis, callback_all, callback_outside) |
|
232 |
|
|
|
233 |
|
inside = shell_ps + self.shell.ps |
|
234 |
|
outside = shell_pf + self.shell.pf |
|
235 |
|
|
|
236 |
|
return ShellEstimation(nsampled, shell_ps, shell_pf, inside, outside) |
|
237 |
|
|
|
238 |
|
|