4Companion TwoRaySpectrumPropagationLossModel calibration script
6Copyright (c) 2022 SIGNET Lab, Department of Information Engineering,
9This program is free software: you can redistribute it and/or modify it under
10the terms of the GNU General Public License as published by the Free Software
11Foundation, either version 3 of the License, or (at your option) any later
14This program is distributed in the hope that it will be useful, but WITHOUT
15ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
16FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
18You should have received a copy of the GNU General Public License along with
19this program. If not, see <http://www.gnu.org/licenses/>.
23from matplotlib
import pyplot
as plt
26from pathlib
import Path
27from itertools
import product
31import argparse
as argp
35parser = argp.ArgumentParser(formatter_class=argp.ArgumentDefaultsHelpFormatter)
36parser.add_argument(
"--num_search_grid_params", default=30,
37 help=
"Number of values for each parameter of the search grids")
38parser.add_argument(
"--num_refinements", default=1,
39 help=
"Number of refinement local search runs to be carried out")
40parser.add_argument(
"--ref_data_fname", default=
"two-ray-to-three-gpp-splm-calibration.csv",
41 help=
"Filename of the fit reference data, obtained from ns-3")
42parser.add_argument(
"--fit_out_fname", default=
"two-ray-splm-fitted-params.txt",
43 help=
"Filename of the fit results")
44parser.add_argument(
"--c_plus_plus_out_fname", 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")
46parser.add_argument(
"--figs_folder", default=
"FiguresTwoRayThreeGppChCalibration/",
47 help=
"Output folder for the fit results figures")
48parser.add_argument(
"--epsilon", default=1e-7,
49 help=
"Tolerance value for the preliminary tests")
50parser.add_argument(
"--preliminary_fit_test", default=
True,
51 help=
"Whether to run preliminary tests which check the correctness of the script functions")
52parser.add_argument(
"--fit_ftr_to_threegpp", default=
True,
53 help=
"Whether to run the calibration with respect to the 3GPP reference channel gains")
54parser.add_argument(
"--output_ns3_table", default=
True,
55 help=
"Whether to output the code for importing the calibration results in ns-3")
56parser.add_argument(
"--plot_fit_results", default=
False,
57 help=
"Whether to plot a comparison of the reference data ECDFs vs the fitted FTR distributions")
59args = parser.parse_args()
61num_search_grid_params = int(args.num_search_grid_params)
63num_refinements = int(args.num_refinements)
65ref_data_fname = args.ref_data_fname
67fit_out_fname = args.fit_out_fname
69c_plus_plus_out_fname = args.c_plus_plus_out_fname
71figs_folder = args.figs_folder
73epsilon = float(args.epsilon)
75preliminary_fit_test = bool(args.preliminary_fit_test)
77fit_ftr_to_threegpp = bool(args.fit_ftr_to_threegpp)
79output_ns3_table = bool(args.output_ns3_table)
81plot_fit_results = bool(args.plot_fit_results)
84@contextlib.contextmanager
87 Context manager to patch joblib to report into tqdm progress bar given as argument.
88 Taken
from: https://stackoverflow.com/questions/24983493/tracking-progress-of-joblib-parallel-execution
91 class TqdmBatchCompletionCallback(joblib.parallel.BatchCompletionCallBack):
92 def __call__(self, *args, **kwargs):
93 tqdm_object.update(n=self.batch_size)
94 return super().__call__(*args, **kwargs)
96 old_batch_callback = joblib.parallel.BatchCompletionCallBack
97 joblib.parallel.BatchCompletionCallBack = TqdmBatchCompletionCallback
101 joblib.parallel.BatchCompletionCallBack = old_batch_callback
116 def __init__(self, m: float, sigma: float, k: float, delta: float):
117 '''! The initializer.
118 @param self: the object pointer
119 @param m: Parameter m
for the Gamma variable. Used both
as the shape
and rate parameters.
120 @param sigma: Parameter sigma. Used
as the variance of the amplitudes of the normal diffuse components.
121 @param k: Parameter K. Expresses ratio between dominant specular components
and diffuse components.
122 @param delta: Parameter delta [0, 1]. Expresses how similar the amplitudes of the two dominant specular components are.
131 '''! The initializer with default values.
132 @param self: the object pointer
141 '''! The initializer with default values.
142 @param self: the object pointer
143 @returns A string reporting the value of each of the FTR fading model parameters
146 return f
'm: {self.m}, sigma: {self.sigma}, k: {self.k}, delta: {self.delta}'
150 '''! Returns the ECDF for the FTR fading model,
for a given parameter grid.
151 @param params: The FTR parameters grid.
152 @param n_samples: The number of samples of the output ECDF
153 @param db: Whether to
return the ECDF
with the gain expressed
in dB
154 @returns The ECDF
for the FTR fading model
157 assert (params.delta >= 0
and params.delta <= 1.0)
160 cmn_sqrt_term = np.sqrt(1 - params.delta**2)
161 v1 = np.sqrt(params.sigma) * np.sqrt(params.k * (1 - cmn_sqrt_term))
162 v2 = np.sqrt(params.sigma) * np.sqrt(params.k * (1 + cmn_sqrt_term))
164 assert (abs((v1**2 + v2**2)/(2*params.sigma) - params.k) < 1e-5)
166 assert (abs((2*v1*v2)/(v1**2 + v2**2) - params.delta) < 1e-4)
168 assert (v1 == v2 == params.k)
170 sqrt_gamma = np.sqrt(np.random.gamma(
171 shape=params.m, scale=1/params.m, size=n_samples))
174 phi1 = np.random.uniform(low=0, high=1.0, size=n_samples)
175 phi2 = np.random.uniform(low=0, high=1.0, size=n_samples)
178 x = np.random.normal(scale=np.sqrt(params.sigma), size=n_samples)
179 y = np.random.normal(scale=np.sqrt(params.sigma), size=n_samples)
182 compl_phi1 = np.vectorize(complex)(np.cos(phi1), np.sin(phi1))
183 compl_phi2 = np.vectorize(complex)(np.cos(phi2), np.sin(phi2))
184 compl_xy = np.vectorize(complex)(x, y)
185 h = np.multiply(sqrt_gamma, compl_phi1) * v1 + \
186 np.multiply(sqrt_gamma, compl_phi2) * v2 + compl_xy
189 power = np.square(np.absolute(h))
192 power = 10*np.log10(power)
194 return np.sort(power)
198 '''! Computes the mean of the FTR fading model, given a specific set of parameters.
199 @param params: The FTR fading model parameters.
202 cmn_sqrt_term = np.sqrt(1 - params.delta**2)
203 v1 = np.sqrt(params.sigma) * np.sqrt(params.k * (1 - cmn_sqrt_term))
204 v2 = np.sqrt(params.sigma) * np.sqrt(params.k * (1 + cmn_sqrt_term))
206 mean = v1**2 + v2**2 + 2*params.sigma
212 '''! Computes the mean of the FTR fading model using the formula reported in the corresponding paper,
213 given a specific set of parameters.
214 @param params: The FTR fading model parameters.
217 return 2 * params.sigma * (1 + params.k)
221 '''! Computes the Anderson-Darling measure for the specified reference
and targets distributions.
222 In particular, the Anderson-Darling measure
is defined
as:
223 \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$.
225 See https://www.itl.nist.gov/div898/handbook/eda/section3/eda35e.htm
for further details.
227 @param ref_ecdf: The reference ECDF.
228 @param target_ecdf: The target ECDF we wish to match the reference distribution to.
229 @returns The Anderson-Darling measure
for the specified reference
and targets distributions.
232 assert (len(ref_ecdf) == len(target_ecdf))
235 mult_factors = np.linspace(start=1, stop=n, num=n)*2 + 1
239 with np.errstate(divide=
'ignore'):
240 log_a_plus_b = np.log(ecdf_values) + np.log(1 - np.flip(ecdf_values))
242 valid_idxs = np.isfinite(log_a_plus_b)
243 A_sq = - np.dot(mult_factors[valid_idxs], log_a_plus_b[valid_idxs])
249 '''! Given an ECDF and data points belonging to its domain, returns their associated EDCF value.
250 @param ecdf: The ECDF, represented
as a sorted list of samples.
251 @param data_points: A list of data points belonging to the same domain
as the samples.
252 @returns The ECDF value of the domain points of the specified ECDF
256 for point
in data_points:
257 idx = np.searchsorted(ecdf, point) / len(ecdf)
258 ecdf_values.append(idx)
260 return np.asarray(ecdf_values)
264 '''! Computes the value for the FTR parameter sigma, given k, yielding a unit-mean fading process.
265 @param k: The K parameter of the FTR fading model, which represents the ratio of the average power
266 of the dominant components to the power of the remaining diffuse multipath.
267 @returns The value
for the FTR parameter sigma, given k, yielding a unit-mean fading process.
270 return 1 / (2 + 2 * k)
273def fit_ftr_to_reference(ref_data: pd.DataFrame, ref_params_combo: tuple, num_params: int, num_refinements: int) -> str:
274 '''! Estimate the FTR parameters yielding the closest ECDF to the reference one.
276 Uses a global search to estimate the FTR parameters yielding the best fit to the reference ECDF.
277 Then, the search
is refined by repeating the procedure
in the neighborhood of the parameters
278 identified
with the
global search. Such a neighborhood
is determined
as the interval whose center
279 is the previous iteration best value,
and the lower
and upper bounds are the first lower
and upper
280 values which were previously considered, respectively.
282 @param ref_data: The reference data, represented
as a DataFrame of samples.
283 @param ref_params_combo: The specific combination of simulation parameters corresponding
284 to the reference ECDF
285 @param num_params: The number of values of each parameter
in the
global and local search grids.
286 @param num_refinements: The number of local refinement search to be carried out after the
global search.
288 @returns An estimate of the FTR parameters yielding the closest ECDF to the reference one.
292 ref_ecdf = ref_data.query(
293 'scen == @ref_params_combo[0] and cond == @ref_params_combo[1] and fc == @ref_params_combo[2]')
296 n_samples = len(ref_ecdf)
303 m_and_k_step = (m_and_k_ub - m_and_k_lb) / n_samples
306 delta_step = 1/n_samples
309 coarse_search_grid = {
311 'm': np.power(np.ones(num_params)*10, np.linspace(start=m_and_k_lb, stop=m_and_k_ub, endpoint=
True, num=num_params)),
313 'k': np.power(np.ones(num_params)*10, np.linspace(start=m_and_k_lb, stop=m_and_k_ub, endpoint=
True, num=num_params)),
315 'delta': np.linspace(start=0.0, stop=1.0, endpoint=
True, num=num_params)
319 for element
in product(*coarse_search_grid.values()):
323 params.m = element[0]
324 params.k = element[1]
325 params.delta = element[2]
332 if (ad_meas < best_ad):
336 for _
in range(num_refinements):
338 finer_search_grid = {
339 'm': np.power(np.ones(num_params)*10,
340 np.linspace(start=
max(0, np.log10(best_params.m) - m_and_k_step),
341 stop=np.log10(best_params.m) +
343 endpoint=
True, num=num_params)),
344 'k': np.power(np.ones(num_params)*10,
345 np.linspace(start=
max(0, np.log10(best_params.k) - m_and_k_step),
346 stop=np.log10(best_params.k) +
348 endpoint=
True, num=num_params)),
349 'delta': np.linspace(start=
max(0, best_params.delta - delta_step),
350 stop=
min(1, best_params.delta + delta_step),
351 endpoint=
True, num=num_params)
355 m_and_k_step = (np.log10(best_params.m) + m_and_k_step -
356 max(0, np.log10(best_params.m) - m_and_k_step)) / n_samples
357 delta_step = (
min(1, best_params.delta + 1/num_params) -
358 max(0, best_params.delta - 1/num_params)) / n_samples
360 for element
in product(*finer_search_grid.values()):
364 params.m = element[0]
365 params.k = element[1]
366 params.delta = element[2]
373 if (ad_meas < best_ad):
377 out_str = f
"{ref_params_combo[0]}\t{ref_params_combo[1]}\t{ref_params_combo[2]}" + \
378 f
" \t{best_params.sigma}\t{best_params.k}\t{best_params.delta}\t{best_params.m}\n"
384 text += f
'TwoRaySpectrumPropagationLossModel::FtrParams({np.format_float_scientific(params.m)}, {np.format_float_scientific(params.sigma)}, \
385 {np.format_float_scientific(params.k)}, {np.format_float_scientific(params.delta)})'
392 Prints to a file the results of the FTR fit, as C++ code.
395 fit (pd.DataFrame): A Pandas Dataframe holding the results of the FTR fit.
396 out_fname (str): The name of the file to
print the C++ code to.
401 for scen
in set(fit[
'scen']):
402 out_str += f
'{{\"{scen}\",\n{{'
404 for cond
in set(fit[
'cond']):
405 out_str += f
'{{ChannelCondition::LosConditionValue::{cond}, \n'
408 freqs = np.sort(
list(set(fit[
'fc'])))
411 out_str += f
'{float(fc)}, '
412 out_str = out_str[0:-2]
417 fit_line = fit.query(
418 'scen == @scen and cond == @cond and fc == @fc')
419 assert(fit_line.reset_index().shape[0] == 1)
422 params.m = fit_line.iloc[0][
'm']
423 params.k = fit_line.iloc[0][
'k']
424 params.delta = fit_line.iloc[0][
'delta']
425 params.sigma = fit_line.iloc[0][
'sigma']
431 out_str = out_str[0:-2]
435 out_str = out_str[0:-2]
438 out_str = out_str[0:-2]
441 with open(out_fname,
"w", encoding=
"utf-8")
as f:
445if __name__ ==
'__main__':
452 df = pd.read_csv(ref_data_fname, sep=
'\t')
454 df[
'gain'] = 10*np.log10(df[
'gain'])
457 scenarios = set(df[
'scen'])
458 is_los = set(df[
'cond'])
459 frequencies = np.sort(
list(set(df[
'fc'])))
465 if preliminary_fit_test:
473 for delta
in np.linspace(start=0, stop=1, num=100):
477 avg_mean = np.mean(mean_list)
478 assert np.all(np.abs(mean_list - avg_mean) < epsilon)
484 for k
in np.linspace(start=1, stop=500, num=50):
491 assert np.all(np.abs(mean_list - np.float64(1.0)) < epsilon)
492 assert np.all(np.abs(mean_th_list - np.float64(1.0)) < epsilon)
494 if fit_ftr_to_threegpp:
497 with tqdm_joblib(tqdm(desc=
"Fitting FTR to the 3GPP fading model", total=(len(scenarios) * len(is_los) * len(frequencies))))
as progress_bar:
498 res = joblib.Parallel(n_jobs=10)(
499 joblib.delayed(fit_ftr_to_reference)(df, params_comb, num_search_grid_params, num_refinements)
for params_comb
in product(scenarios, is_los, frequencies))
501 with open(fit_out_fname,
"w", encoding=
"utf-8")
as f:
502 f.write(
"scen\tcond\tfc\tsigma\tk\tdelta\tm\n")
509 fit = pd.read_csv(fit_out_fname, delimiter=
'\t')
517 sns.set(rc={
'figure.figsize': (7, 5)})
519 sns.set_style(
'darkgrid')
521 fit = pd.read_csv(fit_out_fname, delimiter=
'\t')
523 Path(figs_folder).mkdir(parents=
True, exist_ok=
True)
526 for params_comb
in product(scenarios, is_los, frequencies):
528 data_query =
'scen == @params_comb[0] and cond == @params_comb[1] and fc == @params_comb[2]'
531 ref_data = df.query(data_query)
534 fit_line = fit.query(data_query)
535 assert(fit_line.reset_index().shape[0] == 1)
537 params.m = fit_line.iloc[0][
'm']
538 params.k = fit_line.iloc[0][
'k']
539 params.delta = fit_line.iloc[0][
'delta']
540 params.sigma = fit_line.iloc[0][
'sigma']
547 np.sort(ref_data[
'gain']), ftr_ecdf)
548 ad_measures.append(np.sqrt(ad_meas))
550 sns.ecdfplot(data=ref_data, x=
'gain',
551 label=
'38.901 reference model')
553 ftr_ecdf, label=f
'Fitted FTR, sqrt(AD)={round(np.sqrt(ad_meas), 2)}')
555 'End-to-end channel gain due to small scale fading [dB]')
558 f
'{figs_folder}{params_comb[0]}_{params_comb[1]}_{params_comb[2]/1e9}GHz_fit.png', dpi=500, bbox_inches=
'tight')
562 sns.ecdfplot(ad_measures, label=
'AD measures')
563 plt.xlabel(
'Anderson-Darling goodness-of-fit')
565 plt.savefig(f
'{figs_folder}AD_measures.png',
566 dpi=500, bbox_inches=
'tight')
def __str__(self)
The initializer with default values.
def __init__(self, float m, float sigma, float k, float delta)
The initializer.
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.
def compute_ftr_th_mean(FtrParams params)
Computes the mean of the FTR fading model using the formula reported in the corresponding paper,...
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.
def get_ftr_ecdf(FtrParams params, int n_samples, db=False)
Returns the ECDF for the FTR fading model, for a given parameter grid.
def print_cplusplus_map_from_fit_results(pd.DataFrame fit, str out_fname)
float get_sigma_from_k(float k)
Computes the value for the FTR parameter sigma, given k, yielding a unit-mean fading process.
def compute_ftr_mean(FtrParams params)
Computes the mean of the FTR fading model, given a specific set of parameters.
def tqdm_joblib(tqdm_object)
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.