Skip to content

Commit 233b3c9

Browse files
inti-abbateinti-abbate
authored
New function to run the optimization with kNN and MLCV methods (#42)
* New function to fit with kNN and MLCV * Modified the Verification.ipynb notebook to include the usage of the new 'run' function
1 parent f955fb2 commit 233b3c9

11 files changed

Lines changed: 443 additions & 209 deletions

File tree

‎docs/examples/Verification.ipynb‎

Lines changed: 230 additions & 95 deletions
Large diffs are not rendered by default.

‎docs/examples/resampled.mcpl.gz‎

-24.1 MB
Binary file not shown.

‎docs/examples/samples.mcpl.gz‎

-2.72 MB
Binary file not shown.

‎docs/examples/samples_bws‎

-391 KB
Binary file not shown.

‎docs/examples/source.xml‎

Lines changed: 0 additions & 35 deletions
This file was deleted.

‎python/kdsource/__init__.py‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -30,10 +30,10 @@
3030
"""
3131

3232
__version__ = "0.1.0"
33-
__author__ = "inti-abbate"
33+
__author__ = "KDSource"
3434

3535
from .geom import Geometry, Metric
36-
from .kde import bw_silv
36+
from .kde import bw_methods, bw_silv
3737
from .kdsource import KDSource, load
3838
from .plist import PList, appendssv, convert2mcpl, join2mcpl, savessv
3939
from .surfsource import SurfaceSourceFile, create_source_file

‎python/kdsource/geom.py‎

Lines changed: 8 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -408,7 +408,7 @@ def transform(self, ts):
408408

409409
def inverse_transform(self, ds):
410410
"""Transform deacde to time."""
411-
return 10**ds
411+
return 10 ** ds
412412

413413
def jac(self, ts):
414414
"""Jacobian of lethargy transformation."""
@@ -549,7 +549,7 @@ def transform(self, poss):
549549
Transform volume position (x,y,z)
550550
to circular flat position (rho,psi).
551551
"""
552-
rhos = np.sqrt(poss[:, 0]**2 + poss[:, 1]**2)
552+
rhos = np.sqrt(poss[:, 0] ** 2 + poss[:, 1] ** 2)
553553
psis = np.arctan2(poss[:, 1], poss[:, 0]) * 180 / np.pi
554554
return np.stack((rhos, psis), axis=1)
555555

@@ -574,7 +574,7 @@ def save(self, mtree):
574574

575575
def jac(self, poss):
576576
"""Jacobian of polar transformation."""
577-
rhos = np.sqrt(poss[:, 0]**2 + poss[:, 1]**2)
577+
rhos = np.sqrt(poss[:, 0] ** 2 + poss[:, 1] ** 2)
578578
return 1 / rhos.reshape(-1)
579579

580580
@staticmethod
@@ -614,7 +614,7 @@ def transform(self, poss):
614614
Transform volume position (x,y,z)
615615
to circular flat position (rho2,psi).
616616
"""
617-
rhos = poss[:, 0]**2 + poss[:, 1]**2
617+
rhos = poss[:, 0] ** 2 + poss[:, 1] ** 2
618618
psis = np.arctan2(poss[:, 1], poss[:, 0]) * 180 / np.pi
619619
return np.stack((rhos, psis), axis=1)
620620

@@ -624,8 +624,8 @@ def inverse_transform(self, poss):
624624
to volume position (x,y,z).
625625
"""
626626
z_col = np.broadcast_to(self.z, (*poss.shape[:-1], 1))
627-
x_col = poss[:, 0]**0.5 * np.cos(poss[:, 1] * np.pi / 180)
628-
y_col = poss[:, 0]**0.5 * np.sin(poss[:, 1] * np.pi / 180)
627+
x_col = poss[:, 0] ** 0.5 * np.cos(poss[:, 1] * np.pi / 180)
628+
y_col = poss[:, 0] ** 0.5 * np.sin(poss[:, 1] * np.pi / 180)
629629
return np.stack((x_col, y_col, z_col), axis=1)
630630

631631
def save(self, mtree):
@@ -1085,7 +1085,7 @@ def GeomFlat(
10851085
[
10861086
Lethargy(E0),
10871087
SurfXY(xmin, xmax, ymin, ymax, z),
1088-
Isotrop(keep_zdir=keep_zdir)
1088+
Isotrop(keep_zdir=keep_zdir),
10891089
],
10901090
trasl=trasl,
10911091
rot=rot,
@@ -1116,7 +1116,7 @@ def GeomFlatTemp(
11161116
Lethargy(E0),
11171117
SurfXY(xmin, xmax, ymin, ymax, z),
11181118
Isotrop(keep_zdir=keep_zdir),
1119-
Decade()
1119+
Decade(),
11201120
],
11211121
trasl=trasl,
11221122
rot=rot,

‎python/kdsource/kde.py‎

Lines changed: 78 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -183,7 +183,7 @@ def bw_mlcv(data, weights=None, n_splits=10, seed=None, grid=None, show=True):
183183
Rule.
184184
grid: array-like
185185
Grid of factors for computing bandwidth grid. Default:
186-
numpy.logspace(-1, 1, 20)
186+
numpy.logspace(-1, 1, 10)
187187
show: bool
188188
Whether to show the scores plot.
189189
@@ -200,7 +200,9 @@ def bw_mlcv(data, weights=None, n_splits=10, seed=None, grid=None, show=True):
200200
)
201201
seed = bw_silv(np.shape(data)[1], N_eff)
202202
if grid is None:
203-
grid = np.logspace(-1, 1, 20)
203+
print("No grid specified. Using np.logspace(-0.1, 0.1, 10).")
204+
print("If fitting fails, change the grid.")
205+
grid = np.logspace(-0.1, 0.1, 10)
204206
bw_grid = np.reshape(grid, (-1, *np.ones(np.ndim(seed), dtype=int))) * seed
205207
if n_splits > len(data):
206208
n_splits = len(data)
@@ -225,6 +227,77 @@ def bw_mlcv(data, weights=None, n_splits=10, seed=None, grid=None, show=True):
225227
return bw_mlcv
226228

227229

230+
def bw_auto(
231+
data,
232+
weights,
233+
grid=None,
234+
Nmlcv=1e4,
235+
Nbatch=1e4,
236+
Nknn=10,
237+
show=True,
238+
n_splits=10,
239+
):
240+
241+
"""
242+
Automatic bandwidth optimization with kNN and MLCV methods combined.
243+
244+
Parameters
245+
----------
246+
data: array-like
247+
Array of samples. Must have shape (obs, dim).
248+
weights: array-like, optional
249+
Array of sample weights. By default all weights are 1.
250+
grid: array-like
251+
Grid to be used for the MLCV bandwidth optimization. Default:
252+
numpy.logspace(-1, 1, 10)
253+
Nmlcv: int
254+
Number of particles to use for the MLCV bandwidth optimization.
255+
Nbatch: int
256+
Number of particles per batch to use for the kNN bandwidth
257+
optimization.
258+
Nknn: int
259+
Number of closest neighbors used for the kNN bandiwdth
260+
optimization.
261+
show: bool
262+
Whether to show the scores plot of cross-validation.
263+
n_splits: int, optional
264+
Number of folds for cross-validation. Must be at least 2.
265+
"""
266+
267+
N, dim = data.shape
268+
269+
# First fit with kNN to generate a seed adaptative bandwidth
270+
print("Fitting first with kNN method:")
271+
if Nbatch == -1:
272+
print("No batch size for kNN selected. Using all the particles list.")
273+
Nbatch = N
274+
BW_knn = bw_knn(data, weights=weights, k=Nknn, batch_size=Nbatch)
275+
print("")
276+
277+
# Then fit with MLCV
278+
print("Fitting now with MLCV using the previous fitting as seed:")
279+
if Nmlcv == -1:
280+
print("No size for MLCV selected. Using all the particles list.")
281+
Nmlcv = N
282+
seed = BW_knn[:Nmlcv]
283+
print("If fitting takes too long, consider reducing Nmlcv.")
284+
BW_mlcv = bw_mlcv(
285+
data[:Nmlcv],
286+
weights=weights,
287+
seed=seed,
288+
grid=grid,
289+
show=show,
290+
n_splits=n_splits,
291+
)
292+
print("")
293+
294+
print("Extending the MLCV optimization to full kNN bandwidth:")
295+
BW_knn_mlcv = BW_knn * BW_mlcv[0] / BW_knn[0]
296+
BW_knn_mlcv *= bw_silv(dim, len(BW_knn)) / bw_silv(dim, len(BW_mlcv))
297+
print("Done.")
298+
return BW_knn_mlcv
299+
300+
228301
def optimize_bw(
229302
bw_method, vecs, ws=None, weightfun=None, maskfun=None, **kwargs
230303
):
@@ -275,9 +348,11 @@ def optimize_bw(
275348
return bw_knn(vecs, weights=ws, **kwargs)
276349
elif bw_method == "mlcv": # Maximum Likelihood Cross-Validation method
277350
return bw_mlcv(vecs, weights=ws, **kwargs)
351+
elif bw_method == "auto": # kNN and MLCV combination
352+
return bw_auto(vecs, weights=ws, **kwargs)
278353
else:
279354
keys = list(bw_methods.keys())
280355
raise Exception("Invalid bw_method. Available: {}".format(keys))
281356

282357

283-
bw_methods = {"silv": bw_silv, "knn": bw_knn, "mlcv": bw_mlcv}
358+
bw_methods = {"auto": bw_auto, "silv": bw_silv, "knn": bw_knn, "mlcv": bw_mlcv}

‎python/kdsource/kdsource.py‎

Lines changed: 52 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -15,13 +15,15 @@
1515
import numpy as np
1616

1717
from .geom import Geometry
18-
from .kde import bw_silv
18+
from .kde import bw_methods, bw_silv
1919
from .kde import optimize_bw
2020
from .plist import PList
2121

22-
R = {'g': 1 / (2 * np.sqrt(np.pi)),
23-
'e': 3 / 5,
24-
'b': 1 / 3} # Roughness of the kernel
22+
R = {
23+
"g": 1 / (2 * np.sqrt(np.pi)),
24+
"e": 3 / 5,
25+
"b": 1 / 3,
26+
} # Roughness of the kernel
2527

2628
STD_DEFAULT = 1
2729

@@ -57,13 +59,13 @@ def load(xmlfilename, N=-1):
5759
kern = kelem.text
5860
else:
5961
print("No kernel specified. Using gaussian as default.")
60-
kern = 'g'
61-
if kern == 'g':
62-
kern = 'gaussian'
63-
elif kern == 'e':
64-
kern = 'epa'
65-
elif kern == 'b':
66-
kern = 'box'
62+
kern = "g"
63+
if kern == "g":
64+
kern = "gaussian"
65+
elif kern == "e":
66+
kern = "epa"
67+
elif kern == "b":
68+
kern = "box"
6769
plist = PList.load(root.find("PList"))
6870
geom = Geometry.load(root.find("Geom"))
6971
scaling = np.array(root.find("scaling").text.split(), dtype="float64")
@@ -78,7 +80,7 @@ def load(xmlfilename, N=-1):
7880

7981

8082
class KDSource:
81-
def __init__(self, plist, geom, bw="silv", J=1.0, kernel="gaussian"):
83+
def __init__(self, plist, geom, bw="auto", J=1.0, kernel="gaussian"):
8284
"""
8385
Object representing Kernel Density Estimation (KDE) sources.
8486
@@ -109,6 +111,8 @@ def __init__(self, plist, geom, bw="silv", J=1.0, kernel="gaussian"):
109111
The source total current, in [1/s]. If set, the density
110112
plots will have the correct units.
111113
kernel: string, optional
114+
The function kernel to fit the variables. Available options
115+
are 'gaussian', 'box', and 'epa'
112116
"""
113117
self.plist = plist
114118
self.geom = geom
@@ -128,7 +132,15 @@ def __init__(self, plist, geom, bw="silv", J=1.0, kernel="gaussian"):
128132
self.R = R[self.kernel[0]]
129133
self.fitted = False
130134

131-
def fit(self, N=-1, skip=0, scaling=None, **kwargs):
135+
def fit(
136+
self,
137+
N=-1,
138+
skip=0,
139+
scaling=None,
140+
importance=None,
141+
bw_method=None,
142+
**kwargs,
143+
):
132144
"""
133145
Fit KDE model to particle list.
134146
@@ -149,8 +161,18 @@ def fit(self, N=-1, skip=0, scaling=None, **kwargs):
149161
scaling: array-like, optional
150162
Scaling to be applied to each variable. This means each
151163
particle variable will be divided by the corresponding
152-
scaling element before applying KDE. By default, the
153-
standard deviation of each variable is used.
164+
scaling element before applying KDE. If set, overrides
165+
importance parameter. By default, the standard deviation of
166+
each variable is used.
167+
importance: array-like, optional
168+
Scaling to be applied to each variable. This means each
169+
particle variable standard deviation will be divided by
170+
the corresponding scaling element before applying KDE.
171+
By default, the importance is set to 1, i.e.
172+
the standard deviation of each variable is used.
173+
bw_method: str
174+
Bandwidth selection method. See bw_methods for available
175+
methods. If None, self.bw_method will be used.
154176
**kwargs: optional
155177
Parameters to be passed to bandwidth selection method. Refer
156178
corresponding method for docs (see bw_methods for method
@@ -165,10 +187,25 @@ def fit(self, N=-1, skip=0, scaling=None, **kwargs):
165187
self.N_eff = np.sum(ws) ** 2 / np.sum(ws ** 2)
166188
if scaling is None:
167189
scaling = self.geom.std(vecs=vecs, weights=ws)
190+
# If an importance is specified, apply to the scaling factors
191+
if importance is None:
192+
importance = [1] * self.geom.dim
193+
elif len(importance) != self.geom.dim:
194+
raise Exception(
195+
"importance must have len = {}.".format(self.geom.dim)
196+
)
197+
scaling /= importance
168198
else:
169199
scaling = np.array(scaling)
170200
scaling[scaling == 0] = STD_DEFAULT
171201
self.scaling = scaling
202+
if bw_method is not None:
203+
keys = list(bw_methods.keys())
204+
if bw_method not in keys:
205+
raise Exception(
206+
"Invalid bw_method. Available: {}".format(keys)
207+
)
208+
self.bw_method = bw_method
172209
if self.bw_method is not None:
173210
print("Calculating bw ... ")
174211
bw = optimize_bw(self.bw_method, vecs / self.scaling, ws, **kwargs)

0 commit comments

Comments
 (0)