4Companion TwoRaySpectrumPropagationLossModel calibration script
6Copyright (c) 2022 SIGNET Lab, Department of Information Engineering,
9SPDX-License-Identifier: GPL-2.0-only
12import argparse
as argp
14from itertools
import product
15from pathlib
import Path
21from matplotlib
import pyplot
as plt
25parser = argp.ArgumentParser(formatter_class=argp.ArgumentDefaultsHelpFormatter)
27 "--num_search_grid_params",
29 help=
"Number of values for each parameter of the search grids",
32 "--num_refinements", default=1, help=
"Number of refinement local search runs to be carried out"
36 default=
"two-ray-to-three-gpp-splm-calibration.csv",
37 help=
"Filename of the fit reference data, obtained from ns-3",
40 "--fit_out_fname", default=
"two-ray-splm-fitted-params.txt", help=
"Filename of the fit results"
43 "--c_plus_plus_out_fname",
44 default=
"two-ray-cplusplus-fitted-params.txt",
45 help=
"Filename of the fit results, encoded as a C++ data structure to be imported in ns-3",
49 default=
"FiguresTwoRayThreeGppChCalibration/",
50 help=
"Output folder for the fit results figures",
52parser.add_argument(
"--epsilon", default=1e-7, help=
"Tolerance value for the preliminary tests")
54 "--preliminary_fit_test",
56 help=
"Whether to run preliminary tests which check the correctness of the script functions",
59 "--fit_ftr_to_threegpp",
61 help=
"Whether to run the calibration with respect to the 3GPP reference channel gains",
66 help=
"Whether to output the code for importing the calibration results in ns-3",
71 help=
"Whether to plot a comparison of the reference data ECDFs vs the fitted FTR distributions",
74args = parser.parse_args()
76num_search_grid_params = int(args.num_search_grid_params)
78num_refinements = int(args.num_refinements)
80ref_data_fname = args.ref_data_fname
82fit_out_fname = args.fit_out_fname
84c_plus_plus_out_fname = args.c_plus_plus_out_fname
86figs_folder = args.figs_folder
88epsilon = float(args.epsilon)
90preliminary_fit_test = bool(args.preliminary_fit_test)
92fit_ftr_to_threegpp = bool(args.fit_ftr_to_threegpp)
94output_ns3_table = bool(args.output_ns3_table)
96plot_fit_results = bool(args.plot_fit_results)
99@contextlib.contextmanager
102 Context manager to patch joblib to report into tqdm progress bar given as argument.
103 Taken from: https://stackoverflow.com/questions/24983493/tracking-progress-of-joblib-parallel-execution
107 class TqdmBatchCompletionCallback(joblib.parallel.BatchCompletionCallBack):
108 def __call__(self, *args, **kwargs):
109 tqdm_object.update(n=self.batch_size)
110 return super().__call__(*args, **kwargs)
112 old_batch_callback = joblib.parallel.BatchCompletionCallBack
113 joblib.parallel.BatchCompletionCallBack = TqdmBatchCompletionCallback
117 joblib.parallel.BatchCompletionCallBack = old_batch_callback
132 def __init__(self, m: float, sigma: float, k: float, delta: float):
133 """! The initializer.
134 @param self: the object pointer
135 @param m: Parameter m for the Gamma variable. Used both as the shape and rate parameters.
136 @param sigma: Parameter sigma. Used as the variance of the amplitudes of the normal diffuse components.
137 @param k: Parameter K. Expresses ratio between dominant specular components and diffuse components.
138 @param delta: Parameter delta [0, 1]. Expresses how similar the amplitudes of the two dominant specular components are.
147 """! The initializer with default values.
148 @param self: the object pointer
157 """! The initializer with default values.
158 @param self: the object pointer
159 @returns A string reporting the value of each of the FTR fading model parameters
162 return f
"m: {self.m}, sigma: {self.sigma}, k: {self.k}, delta: {self.delta}"
166 """! Returns the ECDF for the FTR fading model, for a given parameter grid.
167 @param params: The FTR parameters grid.
168 @param n_samples: The number of samples of the output ECDF
169 @param db: Whether to return the ECDF with the gain expressed in dB
170 @returns The ECDF for the FTR fading model
173 assert params.delta >= 0
and params.delta <= 1.0
176 cmn_sqrt_term = np.sqrt(1 - params.delta**2)
177 v1 = np.sqrt(params.sigma) * np.sqrt(params.k * (1 - cmn_sqrt_term))
178 v2 = np.sqrt(params.sigma) * np.sqrt(params.k * (1 + cmn_sqrt_term))
180 assert abs((v1**2 + v2**2) / (2 * params.sigma) - params.k) < 1e-5
182 assert abs((2 * v1 * v2) / (v1**2 + v2**2) - params.delta) < 1e-4
184 assert v1 == v2 == params.k
186 sqrt_gamma = np.sqrt(np.random.gamma(shape=params.m, scale=1 / params.m, size=n_samples))
189 phi1 = np.random.uniform(low=0, high=1.0, size=n_samples)
190 phi2 = np.random.uniform(low=0, high=1.0, size=n_samples)
193 x = np.random.normal(scale=np.sqrt(params.sigma), size=n_samples)
194 y = np.random.normal(scale=np.sqrt(params.sigma), size=n_samples)
197 compl_phi1 = np.vectorize(complex)(np.cos(phi1), np.sin(phi1))
198 compl_phi2 = np.vectorize(complex)(np.cos(phi2), np.sin(phi2))
199 compl_xy = np.vectorize(complex)(x, y)
201 np.multiply(sqrt_gamma, compl_phi1) * v1
202 + np.multiply(sqrt_gamma, compl_phi2) * v2
207 power = np.square(np.absolute(h))
210 power = 10 * np.log10(power)
212 return np.sort(power)
216 """! Computes the mean of the FTR fading model, given a specific set of parameters.
217 @param params: The FTR fading model parameters.
220 cmn_sqrt_term = np.sqrt(1 - params.delta**2)
221 v1 = np.sqrt(params.sigma) * np.sqrt(params.k * (1 - cmn_sqrt_term))
222 v2 = np.sqrt(params.sigma) * np.sqrt(params.k * (1 + cmn_sqrt_term))
224 mean = v1**2 + v2**2 + 2 * params.sigma
230 """! Computes the mean of the FTR fading model using the formula reported in the corresponding paper,
231 given a specific set of parameters.
232 @param params: The FTR fading model parameters.
235 return 2 * params.sigma * (1 + params.k)
239 """! Computes the Anderson-Darling measure for the specified reference and targets distributions.
240 In particular, the Anderson-Darling measure is defined as:
241 \f$A^2 = -N -S\f$, where \f$S = \sum_{i=1}^N \frac{2i - 1}{N} \left[ ln F(Y_i) + ln F(Y_{N + 1 - i}) \right]\f$.
243 See https://www.itl.nist.gov/div898/handbook/eda/section3/eda35e.htm for further details.
245 @param ref_ecdf: The reference ECDF.
246 @param target_ecdf: The target ECDF we wish to match the reference distribution to.
247 @returns The Anderson-Darling measure for the specified reference and targets distributions.
250 assert len(ref_ecdf) == len(target_ecdf)
253 mult_factors = np.linspace(start=1, stop=n, num=n) * 2 + 1
257 with np.errstate(divide=
"ignore"):
258 log_a_plus_b = np.log(ecdf_values) + np.log(1 - np.flip(ecdf_values))
260 valid_idxs = np.isfinite(log_a_plus_b)
261 A_sq = -np.dot(mult_factors[valid_idxs], log_a_plus_b[valid_idxs])
267 """! Given an ECDF and data points belonging to its domain, returns their associated EDCF value.
268 @param ecdf: The ECDF, represented as a sorted list of samples.
269 @param data_points: A list of data points belonging to the same domain as the samples.
270 @returns The ECDF value of the domain points of the specified ECDF
274 for point
in data_points:
275 idx = np.searchsorted(ecdf, point) / len(ecdf)
276 ecdf_values.append(idx)
278 return np.asarray(ecdf_values)
282 """! Computes the value for the FTR parameter sigma, given k, yielding a unit-mean fading process.
283 @param k: The K parameter of the FTR fading model, which represents the ratio of the average power
284 of the dominant components to the power of the remaining diffuse multipath.
285 @returns The value for the FTR parameter sigma, given k, yielding a unit-mean fading process.
288 return 1 / (2 + 2 * k)
292 ref_data: pd.DataFrame, ref_params_combo: tuple, num_params: int, num_refinements: int
294 """! Estimate the FTR parameters yielding the closest ECDF to the reference one.
296 Uses a global search to estimate the FTR parameters yielding the best fit to the reference ECDF.
297 Then, the search is refined by repeating the procedure in the neighborhood of the parameters
298 identified with the global search. Such a neighborhood is determined as the interval whose center
299 is the previous iteration best value, and the lower and upper bounds are the first lower and upper
300 values which were previously considered, respectively.
302 @param ref_data: The reference data, represented as a DataFrame of samples.
303 @param ref_params_combo: The specific combination of simulation parameters corresponding
304 to the reference ECDF
305 @param num_params: The number of values of each parameter in the global and local search grids.
306 @param num_refinements: The number of local refinement search to be carried out after the global search.
308 @returns An estimate of the FTR parameters yielding the closest ECDF to the reference one.
312 ref_ecdf = ref_data.query(
313 "scen == @ref_params_combo[0] and cond == @ref_params_combo[1] and fc == @ref_params_combo[2]"
317 n_samples = len(ref_ecdf)
324 m_and_k_step = (m_and_k_ub - m_and_k_lb) / n_samples
327 delta_step = 1 / n_samples
330 coarse_search_grid = {
333 np.ones(num_params) * 10,
334 np.linspace(start=m_and_k_lb, stop=m_and_k_ub, endpoint=
True, num=num_params),
338 np.ones(num_params) * 10,
339 np.linspace(start=m_and_k_lb, stop=m_and_k_ub, endpoint=
True, num=num_params),
342 "delta": np.linspace(start=0.0, stop=1.0, endpoint=
True, num=num_params),
346 for element
in product(*coarse_search_grid.values()):
349 params.m = element[0]
350 params.k = element[1]
351 params.delta = element[2]
358 if ad_meas < best_ad:
362 for _
in range(num_refinements):
364 finer_search_grid = {
366 np.ones(num_params) * 10,
368 start=max(0, np.log10(best_params.m) - m_and_k_step),
369 stop=np.log10(best_params.m) + m_and_k_step,
375 np.ones(num_params) * 10,
377 start=max(0, np.log10(best_params.k) - m_and_k_step),
378 stop=np.log10(best_params.k) + m_and_k_step,
383 "delta": np.linspace(
384 start=max(0, best_params.delta - delta_step),
385 stop=min(1, best_params.delta + delta_step),
393 np.log10(best_params.m) + m_and_k_step - max(0, np.log10(best_params.m) - m_and_k_step)
396 min(1, best_params.delta + 1 / num_params) - max(0, best_params.delta - 1 / num_params)
399 for element
in product(*finer_search_grid.values()):
402 params.m = element[0]
403 params.k = element[1]
404 params.delta = element[2]
411 if ad_meas < best_ad:
416 f
"{ref_params_combo[0]}\t{ref_params_combo[1]}\t{ref_params_combo[2]}"
417 + f
" \t{best_params.sigma}\t{best_params.k}\t{best_params.delta}\t{best_params.m}\n"
424 text += f
"TwoRaySpectrumPropagationLossModel::FtrParams({np.format_float_scientific(params.m)}, {np.format_float_scientific(params.sigma)}, \
425 {np.format_float_scientific(params.k)}, {np.format_float_scientific(params.delta)})"
432 Prints to a file the results of the FTR fit, as C++ code.
435 fit (pd.DataFrame): A Pandas Dataframe holding the results of the FTR fit.
436 out_fname (str): The name of the file to print the C++ code to.
441 for scen
in set(fit[
"scen"]):
442 out_str += f
'{{"{scen}",\n{{'
444 for cond
in set(fit[
"cond"]):
445 out_str += f
"{{ChannelCondition::LosConditionValue::{cond}, \n"
448 freqs = np.sort(list(set(fit[
"fc"])))
451 out_str += f
"{float(fc)}, "
452 out_str = out_str[0:-2]
457 fit_line = fit.query(
"scen == @scen and cond == @cond and fc == @fc")
458 assert fit_line.reset_index().shape[0] == 1
461 params.m = fit_line.iloc[0][
"m"]
462 params.k = fit_line.iloc[0][
"k"]
463 params.delta = fit_line.iloc[0][
"delta"]
464 params.sigma = fit_line.iloc[0][
"sigma"]
470 out_str = out_str[0:-2]
474 out_str = out_str[0:-2]
477 out_str = out_str[0:-2]
480 with open(out_fname,
"w", encoding=
"utf-8")
as f:
484if __name__ ==
"__main__":
490 df = pd.read_csv(ref_data_fname, sep=
"\t")
492 df[
"gain"] = 10 * np.log10(df[
"gain"])
495 scenarios = set(df[
"scen"])
496 is_los = set(df[
"cond"])
497 frequencies = np.sort(list(set(df[
"fc"])))
503 if preliminary_fit_test:
510 for delta
in np.linspace(start=0, stop=1, num=100):
514 avg_mean = np.mean(mean_list)
515 assert np.all(np.abs(mean_list - avg_mean) < epsilon)
521 for k
in np.linspace(start=1, stop=500, num=50):
528 assert np.all(np.abs(mean_list - np.float64(1.0)) < epsilon)
529 assert np.all(np.abs(mean_th_list - np.float64(1.0)) < epsilon)
531 if fit_ftr_to_threegpp:
535 desc=
"Fitting FTR to the 3GPP fading model",
536 total=(len(scenarios) * len(is_los) * len(frequencies)),
539 res = joblib.Parallel(n_jobs=10)(
540 joblib.delayed(fit_ftr_to_reference)(
541 df, params_comb, num_search_grid_params, num_refinements
543 for params_comb
in product(scenarios, is_los, frequencies)
546 with open(fit_out_fname,
"w", encoding=
"utf-8")
as f:
547 f.write(
"scen\tcond\tfc\tsigma\tk\tdelta\tm\n")
553 fit = pd.read_csv(fit_out_fname, delimiter=
"\t")
560 sns.set(rc={
"figure.figsize": (7, 5)})
562 sns.set_style(
"darkgrid")
564 fit = pd.read_csv(fit_out_fname, delimiter=
"\t")
566 Path(figs_folder).mkdir(parents=
True, exist_ok=
True)
569 for params_comb
in product(scenarios, is_los, frequencies):
571 "scen == @params_comb[0] and cond == @params_comb[1] and fc == @params_comb[2]"
575 ref_data = df.query(data_query)
578 fit_line = fit.query(data_query)
579 assert fit_line.reset_index().shape[0] == 1
581 params.m = fit_line.iloc[0][
"m"]
582 params.k = fit_line.iloc[0][
"k"]
583 params.delta = fit_line.iloc[0][
"delta"]
584 params.sigma = fit_line.iloc[0][
"sigma"]
591 ad_measures.append(np.sqrt(ad_meas))
593 sns.ecdfplot(data=ref_data, x=
"gain", label=
"38.901 reference model")
594 sns.ecdfplot(ftr_ecdf, label=f
"Fitted FTR, sqrt(AD)={round(np.sqrt(ad_meas), 2)}")
595 plt.xlabel(
"End-to-end channel gain due to small scale fading [dB]")
598 f
"{figs_folder}{params_comb[0]}_{params_comb[1]}_{params_comb[2]/1e9}GHz_fit.png",
605 sns.ecdfplot(ad_measures, label=
"AD measures")
606 plt.xlabel(
"Anderson-Darling goodness-of-fit")
608 plt.savefig(f
"{figs_folder}AD_measures.png", dpi=500, bbox_inches=
"tight")
__init__(self, float m, float sigma, float k, float delta)
The initializer.
__str__(self)
The initializer with default values.
delta
Parameter delta [0, 1].
m
Parameter m for the Gamma variable.
str append_ftr_params_to_cpp_string(str text, FtrParams params)
float compute_anderson_darling_measure(list ref_ecdf, list target_ecdf)
Computes the Anderson-Darling measure for the specified reference and targets distributions.
print_cplusplus_map_from_fit_results(pd.DataFrame fit, str out_fname)
str fit_ftr_to_reference(pd.DataFrame ref_data, tuple ref_params_combo, int num_params, int num_refinements)
Estimate the FTR parameters yielding the closest ECDF to the reference one.
compute_ftr_th_mean(FtrParams params)
Computes the mean of the FTR fading model using the formula reported in the corresponding paper,...
compute_ftr_mean(FtrParams params)
Computes the mean of the FTR fading model, given a specific set of parameters.
get_ftr_ecdf(FtrParams params, int n_samples, db=False)
Returns the ECDF for the FTR fading model, for a given parameter grid.
float get_sigma_from_k(float k)
Computes the value for the FTR parameter sigma, given k, yielding a unit-mean fading process.
np.ndarray compute_ecdf_value(list ecdf, float data_points)
Given an ECDF and data points belonging to its domain, returns their associated EDCF value.