From fa31fb46fdde6304ceb535bf1f2d0fb4bdde27b2 Mon Sep 17 00:00:00 2001 From: HappyZ Date: Wed, 22 May 2019 02:05:58 -0500 Subject: [PATCH] add traditional fitting test --- libs/consts.py | 9 +++- libs/fitting.py | 91 ++++++++++++++++++++++++++++++++++++++++ libs/models.py | 4 +- libs/plotting.py | 98 +++++++++++++++++++++++++++++++++++++++++++ libs/spacemap.py | 84 +++++++++++++++++++++++++------------ tests/dev_map.pickle | Bin 0 -> 32928 bytes tests/fitting.py | 73 ++++++++++++++++++++++++++++++++ tests/spacemap.py | 5 ++- 8 files changed, 333 insertions(+), 31 deletions(-) create mode 100644 libs/fitting.py create mode 100644 libs/plotting.py create mode 100644 tests/dev_map.pickle create mode 100644 tests/fitting.py diff --git a/libs/consts.py b/libs/consts.py index 245e70b..9d537fa 100644 --- a/libs/consts.py +++ b/libs/consts.py @@ -1,11 +1,16 @@ #!/usr/bin/python -RANDOM_SEED = 0 +RANDOM_SEED = 666 GAUSSIAN_NOISE_MEAN = 0.0 # dB GAUSSIAN_NOISE_STD = 3.16 # 10 dB variance -NOISE_FLOOR = -85 #dB +NOISE_FLOOR = -85.0 #dB FLOAT_TOLERANCE = 0.001 + +MAX_RANGE = 10.0 # meter + +COLORMAP_MIN = -60.0 +COLORMAP_MAX = -30.0 diff --git a/libs/fitting.py b/libs/fitting.py new file mode 100644 index 0000000..6d7642b --- /dev/null +++ b/libs/fitting.py @@ -0,0 +1,91 @@ +#!/usr/bin/python + + +import numpy as np +from scipy.optimize import curve_fit +from libs.models import log_gamma_loc + + +def modelfit_log_gamma_func(input_data, *target_data): + # unpack input_data + rx_locs = input_data + + # unpack target_data + if len(target_data) < 4: + print("missing arguments, should have loc_x, loc_y, pwr, gamma at least") + return + tx_loc_x = target_data[0] + tx_loc_y = target_data[1] + tx_power = target_data[2] + env_gamma = target_data[3] + tx_loc_z = None if len(target_data) == 4 else target_data[4] + tx_loc = np.array([tx_loc_x, tx_loc_y]) + if tx_loc_z is not None: + tx_loc = np.array([tx_loc_x, tx_loc_y, tx_loc_z]) + return log_gamma_loc(rx_locs, tx_loc, tx_power, env_gamma) + + +def modelfit_log_gamma( + rx_locs, + rx_rsses, + bounds_pwr=(-60, 0), + bounds_gamma=(2, 6), + bounds_loc_x=(0, 6.2), + bounds_loc_y=(0, 6.2), + bounds_loc_z=None, + monte_carlo_sampling=False, + monte_carlo_sampling_rate=0.8 +): + # initial seeds + seed_pwr = np.random.uniform(bounds_pwr[0], bounds_pwr[1]) + seed_gamma = np.random.uniform(bounds_gamma[0], bounds_gamma[1]) + seed_loc_x = np.random.uniform(bounds_loc_x[0], bounds_loc_x[1]) + seed_loc_y = np.random.uniform(bounds_loc_y[0], bounds_loc_y[1]) + if bounds_loc_z is not None: + seed_loc_z = np.random.uniform(bounds_loc_z[0], bounds_loc_z[1]) + seeds = [ + seed_loc_x, seed_loc_y, + seed_pwr, seed_gamma, + seed_loc_z + ] + bounds = list(zip(*[ + bounds_loc_x, bounds_loc_y, + bounds_pwr, bounds_gamma, + bounds_loc_z + ])) + else: + seeds = [ + seed_loc_x, seed_loc_y, + seed_pwr, seed_gamma + ] + bounds = list(zip(*[ + bounds_loc_x, bounds_loc_y, + bounds_pwr, bounds_gamma + ])) + if monte_carlo_sampling and rx_rsses.shape[0] > 10: + logistics = np.random.choice( + np.arange(rx_rsses.shape[0]), + size=int(monte_carlo_sampling_rate * rx_rsses.shape[0]), + replace=False + ) + rx_locs = rx_locs[logistics, :] + rx_rsses = rx_rsses[logistics, :] + + # fit + popt, pcov = curve_fit( + modelfit_log_gamma_func, rx_locs, rx_rsses, + p0=seeds, bounds=bounds + ) + + # unpack popt + if bounds_loc_z is None: + est_loc_x, est_loc_y, est_tx_pwr, est_env_gamma = popt + est_tx_loc = np.array([est_loc_x, est_loc_y]) + else: + est_loc_x, est_loc_y, est_tx_pwr, est_env_gamma, est_loc_z = popt + est_tx_loc = np.array([est_loc_x, est_loc_y, est_loc_z]) + est_rsses = log_gamma_loc(rx_locs, est_tx_loc, est_tx_pwr, est_env_gamma) + est_errors = est_rsses - rx_rsses + pmse = 1.0 * np.nansum(est_errors * est_errors) / est_errors.shape[0] + + return pmse, est_tx_loc, est_tx_pwr, est_env_gamma, est_rsses diff --git a/libs/models.py b/libs/models.py index 2990f8c..bcc25f9 100644 --- a/libs/models.py +++ b/libs/models.py @@ -12,8 +12,8 @@ np.random.seed(RANDOM_SEED) def log_gamma_loc( - rx_loc: np.ndarray, - tx_loc: np.ndarray, + rx_loc: np.ndarray, # dtype=float, avoid int + tx_loc: np.ndarray, # dtype=float, avoid int pwr: float, gamma: float, loss: float = 0.0, diff --git a/libs/plotting.py b/libs/plotting.py new file mode 100644 index 0000000..766539a --- /dev/null +++ b/libs/plotting.py @@ -0,0 +1,98 @@ +#!/usr/bin/python + + +import numpy as np +import matplotlib.pyplot as plt + +from libs.consts import COLORMAP_MIN +from libs.consts import COLORMAP_MAX +from libs.spacemap import SpaceMap + + +def plotSpace(space_map: SpaceMap): + ''' + ''' + tx_locs = space_map.getTransmitterLocs() + penetration_loss, reflection_loss = space_map.getLosses() + + fig = plt.figure(figsize=(8, 3)) + + ax1 = fig.add_subplot(1, 2, 1) + ax1.set_title("penetration loss") + plt.imshow( + np.transpose(penetration_loss), + cmap='hot', + origin='lower', + interpolation='nearest', + vmin=COLORMAP_MIN, + vmax=COLORMAP_MAX + ) + plt.colorbar() + + ax2 = fig.add_subplot(1, 2, 2) + ax2.set_title("reflection loss") + plt.imshow( + np.transpose(reflection_loss), + cmap='hot', + origin='lower', + interpolation='nearest', + vmin=COLORMAP_MIN, + vmax=COLORMAP_MAX + ) + plt.colorbar() + + # plt.show() + plt.draw() + plt.pause(0.1) + q = input("press enter..") + plt.close() + + +def plotRSS( + ori_rss: np.ndarray, + est_rss: np.ndarray = None, + est_tx_loc: np.ndarray = None, + true_tx_loc: np.ndarray = None +): + ''' + ''' + fig = plt.figure(figsize=(8, 3)) + ax1 = fig.add_subplot(1, 2, 1) + ax1.set_title("original rss map") + plt.imshow( + np.transpose(ori_rss), + cmap='hot', + origin='lower', + interpolation='nearest', + vmin=COLORMAP_MIN, + vmax=COLORMAP_MAX + ) + plt.colorbar() + + ax2 = fig.add_subplot(1, 2, 2) + ax2.set_title("est rss map") + if est_rss is not None: + plt.imshow( + np.transpose(est_rss), + cmap='hot', + origin='lower', + interpolation='nearest', + vmin=COLORMAP_MIN, + vmax=COLORMAP_MAX + ) + plt.colorbar() + if est_tx_loc is not None: + ax2.plot(est_tx_loc[0], est_tx_loc[1], 'bo') + ax2.text(est_tx_loc[0] - 5, est_tx_loc[1] + 2, 'est tx', color='b') + if true_tx_loc is not None: + ax2.plot(true_tx_loc[0], true_tx_loc[1], 'bs') + ax2.text(true_tx_loc[0], true_tx_loc[1], 'true tx', color='b') + + # plt.show() + plt.draw() + plt.pause(0.1) + q = input("press enter..") + plt.close() + + + diff --git a/libs/spacemap.py b/libs/spacemap.py index b3a2cbe..c0ff87c 100644 --- a/libs/spacemap.py +++ b/libs/spacemap.py @@ -3,6 +3,7 @@ import numpy as np +from libs.consts import MAX_RANGE from libs.consts import FLOAT_TOLERANCE from libs.models import log_gamma_dist @@ -17,7 +18,7 @@ class SpaceBlock(): self.mat = material self.has_transmitter = False self.loss_penetration = 0.0 - self.loss_reflection = 0.0 + self.loss_reflection = -100.0 self.related_rays = [] def __mul__(self, val): @@ -95,7 +96,7 @@ class SpaceBlock(): def hasTransmitter(self): return self.has_transmitter - def setLoss(self, penetration=0.0, reflection=0.0): + def setLoss(self, penetration=0.0, reflection=-100.0): self.loss_penetration = penetration self.loss_reflection = reflection @@ -109,7 +110,10 @@ class SpaceRay(): ''' def __init__(self, point1, point2): self.start = point1 - self.end = point2 + if not isinstance(point2, SpaceBlock): + self.end = point1 + SpaceBlock(np.cos(point2), np.sin(point2)) * MAX_RANGE + else: + self.end = point2 # property self.distance = None self.angle_theta = None @@ -127,12 +131,13 @@ class SpaceRay(): self.this_pass_through_loss = None self.this_pass_through_blks = None self.this_gamma = None - self.__prev_factor = None - self.__prev_bounds = None self.loss_total = None self.distance_traveled = None + self.reflection_count = 0 # id self.ray_id = '' + self.prev_ray = None + self.next_rays = [] def getAngle(self, degree=False): if self.angle_theta is None: @@ -172,31 +177,25 @@ class SpaceRay(): def setTravelDistance(self, prior_distance): self.distance_traveled = prior_distance + self.getDistance() - def computeLinePassThroughLoss(self, space_map): - factor = space_map.bs - bounds = space_map.map.shape - - if ( - self.__prev_bounds == bounds and - self.__prev_factor == factor and - self.this_pass_through_blks is not None - ): - return np.sum([ - each.loss_penetration - for each in self.this_pass_through_blks - ]) - - self.__prev_bounds = bounds - self.__prev_factor = factor + def derivePassAndReflectBlocks(self, space_map): self.this_pass_through_blks = [] + self.this_reflection_blks = [] step_blk = SpaceBlock(self.getAngleThetaCos(), self.getAngleThetaSin()) * factor for i in range(1, int(self.getDistance() / factor)): next_blk = self.start + step_blk * i # assume 2D x_idx, y_idx = [int(x) for x in (next_blk / factor).round()] if x_idx < bounds[0] and x_idx > -1 and y_idx < bounds[1] and y_idx > -1: - self.this_pass_through_blks.append(space_map.map[x_idx, y_idx]) - space_map.map[x_idx, y_idx].related_rays.append(self) + if space_map.map[x_idx, y_idx].loss_penetration < 0: + self.this_pass_through_blks.append(space_map.map[x_idx, y_idx]) + if space_map.map[x_idx, y_idx].loss_reflection > -100: + self.this_reflection_blks.append(space_map.map[x_idx, y_idx]) + # space_map.map[x_idx, y_idx].related_rays.append(self) + + def getPassThroughLoss(self): + if self.this_pass_through_blks is None: + print("need to run `derivePassAndReflectBlocks` first") + return self.this_pass_through_loss = np.sum([ each.loss_penetration for each in self.this_pass_through_blks @@ -207,7 +206,7 @@ class SpaceRay(): excluding the end penetration/reflection loss ''' if self.this_pass_through_loss is None: - print("need to run `computeLinePassThroughLoss` first") + print("need to run `getPassThroughLoss` first") return self.loss_total = prior_loss + self.this_pass_through_loss @@ -245,7 +244,9 @@ class SpaceMap(): width: float = 6.4, length: float = 6.4, block_size: float = 0.1 - ): + ): + if MAX_RANGE < width or MAX_RANGE < length: + print("WARNNING! MAX_RANGE {} is less than input".format(MAX_RANGE)) self.width = width self.length = length self.bs = block_size @@ -258,12 +259,21 @@ class SpaceMap(): # initialize the map self.__loss_p = np.zeros(self.map.shape) - self.__loss_r = np.zeros(self.map.shape) + self.__loss_r = np.ones(self.map.shape) * -100 + self.__tx_locs = np.zeros(self.map.shape) for j in range(self.map.shape[1]): y = self.bs * (j + 0.5) for i in range(self.map.shape[0]): self.map[i, j] = SpaceBlock(self.bs * (i + 0.5), y) + # propagation + self.env_gamma = 2.0 + self.ray_trace_deg_step = 1.0 + # according to VTC paper + # "A RAY TRACING METHOD FOR PREDICTING PATH LOSS AND DELAY SPREAD + # IN MICROCELLULAR ENVIRONMENTS" + self.ray_trace_deg_tol = self.ray_trace_deg_step / 180.0 * np.pi / np.sqrt(3) + def getLosses(self): return np.array([self.__loss_p, self.__loss_r]) @@ -279,3 +289,25 @@ class SpaceMap(): self.map[i, j].setLoss(penetration, reflection) self.__loss_p[i, j] = penetration self.__loss_r[i, j] = reflection + + def setHasTransmitter(self, i, j): + self.map[i, j].setTransmitter(flag) + self.__tx_locs[i, j] = 1 if flag else 0 + + def getTransmitterLocs(self): + return self.__tx_locs + + def setEnvGamma(self, val: float): + self.env_gamma = val + + def traceRay(self, tx_power: float, tx_loc: SpaceBlock, rx_loc: SpaceBlock): + rays = [] + for direction in np.arange(0, 359.99, self.ray_trace_deg_step): + ray = SpaceRay(tx_loc, direction / 180.0 * np.pi) + ray.setPower(tx_power, self.env_gamma) + ray.ray_id = "ray_{0:.3f}".format(direction) + rays.append(ray) + while len(rays) > 0: + ray = rays.pop(0) + + diff --git a/tests/dev_map.pickle b/tests/dev_map.pickle new file mode 100644 index 0000000000000000000000000000000000000000..e2f67e5ab22849bbafdcd2f7852a283f1864855b GIT binary patch literal 32928 zcmeI5OOIs55rt#G&++4Tuww0n7ilP`ZQ5diWQm;x1C}hAM~#OCyWH5w|0@2OO?ICz zN}Q|Iz5Os6p{fN(nHi50Co=Qis%eSw=8Mn2`s~ZE-+uAM%TM3_?5lsg`uj&OUw-uU z^S{3Q^u-sS|MS($uU@=*evUh<=bwLaw|ahd6aPBD@$R#a&u@Nn{?P~LAHQ|;?55qc zf1iIM(@$P~{pF|6&p-X}jfnW$Z_YoviJNENoPYko`T1KnZ+`gf+u7fqfARkN@4xro zZ$G#HxrtZjUw-_bcVc|II6v(#dBv*^*L`q32iJSxdJkOhf$KeRy$7!M!1sO+#Cv(W z@8ZvQ*lvV>xWoMYcgI8X@yEM)$ENvYy~|E<#yfXo;)dc4*@4VErtCo02XPFG$DhPQ zyvpA4zHaa3Vf?5Uq>dWPp0W3baiWvi6K@ykuXG*g%NU*Eul@)7!17XO^Lg>y)MN9$ z9~$++$B}jT3164~qZ7Ie3KhWlhH?*F@Kxh1n_a+Z{W#{XI zbksbRF0SHLtkDuXyP0{kDEqKg_=2XZ97J zdHO(o_;GIO8+L|GV=!4)9C%c{#3hgHssFXR|9M^tt4{hx9n?!bLiDIUp{M!PI{FdE z>^vTR;E+AY{83#8buIE!f5nv@^UeH+)#JQ^b@?6T`BZk)3&Mxz|9CfJepvT<()6zS z&<~Vdvwp^6?0Kjce!k9_c*TRf2fyD_XZD<(LcqaWdzrYWDgw*+c zVbupuLwHgzV~}}}dFn^cW^pSoI;`~3JbHLs#^fy?+_OUb_b^X++rH>R@h1KTS$DL0 z;6>dl+~fDu>(t|Xn7=Q2oz)-dldh}gE5Gdr@>4hd!{hAu;eQl2@vXnCJ5n$8z~7N| z$D{n7%FcPM#>G>**BD()f7XL_`i#FJb`PzqcrWr1ALRF+RX+NPJvi)l1ma*{XgrGF zj!sSMD_&1MKj@1geBqPUFY0XO=kZE6*F*6daaq@WO}?Y`Vh=C!I38*6|}AW7Z4N0Ys;sGwdcl@tztt`;fj0 zIsXHG!t6O&Zg`>#>uK?b!##>!&7()p89n$t$8k@Gp*YkBucm!Bdd}$ep5A)SQSFB~ zY3%#7-_5&!m_8x9%{+G0i9Xmn%5JaThpwa32H`1ep!lH|I?%sWy{zNs`jEeQ-|Y7S z)`{o)!pGPZPw{8G;@hl#;$V-zBkP6MXXYDvg79L03Lo0nqkib_7{(jE_1uIn=z1{^ zY90TfG3fbf437n3w?yL=sA&Kn(sP3y!XK6tc`54LlDHb&R6 z(5+cL*ik=YbZPvaWcwWZLUh^V12043ke4yoR{c{P?8@Kl!^UL?57v$UAim}szf;Jb zzJm|$>(Tcn=ch6E=wRJ=qYvnI4IP9vUi2URK;nYK;xVpu^u|vA&L7H$9sAk)*U%q- zA%5@zSN>&>Ki_B2RrAf(iC2C9p6>q!FZQ9~faebG^`AIt=f36_@zEQ*!0wbWKiP9a z)ZcZe?~JSeJ*Qa|PwLd)pETm`AbeWS8E*Uji~ma3p402%Cm!@Yqz}CR-{*fQe~K>{ z8lyXPqc?qEeRE&*oK*%of%W?;`+z*4^rIi8i`L=gy3lXO=s9WBd;TDg#iI_#pr2sR z8Lj$P@q5nMnosQMFZr17IjsY|t4{7m_VMrT?1qqiky>{h?mOQ6T>|FB$g`C@pYrN5Av(7zl`ZW22K|bu*Z^Att_SUP9#9tt~PzU2?{OBW?T5pQW zx#u1r9`->z({Znkfj#Zv-l-iU+ViW7$ACp&Z} z?vnSxdy~ElcnH}Stdkes?EB69feW4yWl(UEx6P2V}k z^b&}mRAck^G0pl|@q145sC)mEx5Mfe_4K>>5Z1m?-+Av1drz)0eJ2ij5zp|b zZb=;K^Rb^7%U`}dBQdU5Yp z|A|XH!=`ZvSGxE-=fLY6#JBqsf7Z#z7)<=(BmPMpii=&y4nE|0G!obN6i?2Tu+b0M z!w>slUFVkm6jC3!s&mEXlrcZ_A4CV^HhjV8{~U? z&q*r=di9*c`W+Ts=}%gJs$V^)+ou1`boM;;zF=|NDR0A}3p_j?bq(dI^?v7Hev5m+ zbtcZR{?+el`hA&mjU9Q*pZ&nO5>`FFFI2DddF&WI=ma)?cT(QL`xj%=ui~(845`bJ z=L7f0P~1Y|c>mpwF}kv^=o9COeE@&qQ5{o#LMQsWLhQnM!NdBGed)>iDPw+)@;q$U zbB5i1&*?d5gT6HC?m4IHo`NUHdFLKr-%vk@J^PIK)%Tt=P1L;I_tA4M*S!?dhxl!; z|L95|4xazY+jFLm%H#dD{RS`ZZ|gb3l%Cmh4ln(Hr|LW0*F`;?H*|&9Kria(Iayqu z^F8OV?(eEYbv6*bc3)5zIC{=#wZ8(to^$$=$Id_X^Y0^#eZl?UajdS?dObfBpSXL7 zUXIjLxaz<7z+=D{>^ZBy^uOmE#W)G6chTQH{&gO^;?Zy8TXo2;t?cEudGD5e)rr3& zbqhndSBI~&9}3wI%2)kjEMy;eyw$kk*uJWH&H=iZuGp2H=wUiE`T}3|zs~Oi`%k?= zT=Kxz>w+KivP+FY`adiV`;C17x(KyO~sh|DESnI3)GoRF{ z_cX@bgTXrSh>xDs>&UwN2Y&Ql2#@8bTboy4OK?r+W0cl3b= z=bpZUd-VbU3ma{Z{U;u7C*sQ)QMAM`DK(D9zntBrojk3QsU z*8h|a317y-4L|mW5M8-%%HI4d|0!dB=sURKZ+)TvU~q3(p3<@BOmlg^kR3ee+g@H@ z$FJfLx97}O`Ow|@FwUOSy5766r!LU*A0DF*NZn|Ge%8R$UVp0oOTuWL5%e(+?T z`gsq0=sv&?geQ8fkbXIqKlM#C=s{hqT@~iL7-f!Km#_$zVALpNYtB^i`>L>f5`UxNG z41+N`vtJqGkB+7bdZN!`F;5xugZ`xtJiwmQEK=j9{nL8R>0lpM9NTA&=c@9RUhKo3 zGfF*gscX~!p3~}*4}Rz<p9t6e$T|O=bWy2&Qs{TH10v)584;QzAx1I z(|YAarSW199Ca^zVh&Q{Jf7D@3ONzynmE__&*V{ z-$u__MZL$POTve9#eVHM&FlB>P=7cNjeW&DeK)NAp}sR_Uk^Qx!RR^DZh1=2)^je` zzxT0!3m3Xqo!0O2LkIRT-{X7EbeH|$`(h{FCO-B(XV}Z%BcU5Q8M06O{Gg-b=sBZb zbm}>Wb+6L@6{0i!PVwPhto}c8pHf#LeLyFWdx3dQix)>in$Ij#E*z1PR zkba?l4*IP5WL^89`oQ_eo_m>ocpVy3_nJKl{onN2^P|6z^Ui)5{QZO01Alw& zG_UWJF+b=`|M9PVNnChvPtf0<9 z?1vyfyy4Yzn%90PJ$laR>POWh|D%0g+40^7)_#DOkbW00_|pd=eNXDt{J@?*OW(mg z>iy_B%iiecIwpR`7oEw=c!lkxe6}y)7yWJ>|4m%u(WqaT=1-i}erec)=!c%@{>b;g z#_w^|ZRf1`S$zY&t?wIMScgCTcT}Cju=1mW@}}`AJzaNrf~;#?I>Sfx=PBd%16qFS zUm|&^d&QS|#rOIdql4q>Tpyk1J)1n-?}dxHiwE@?RzIx&_zmg?@wa?#UvaUQzt8WD z$=}GY{s;L%n@1<-MLdxFHHIfB{uM|0hz~aG4Qu_e=U@KLllX<44|IVKbqL$RzK{tIly{{JtV__{~x`-Hw%;5WEum}g90#ZRjbe#MKvm_K%o@ZZz- zq3gt_PVk}ecKuX;bG+w0<2qG5#bLh{(iigb{{aGb@>l=aC-hN>j*i8XemlZTSoXxV zJgfCI`?2tA%)!I(Zll9xr zebQI;J9+N#dlmZxJ0X6ahq2=jj}$((&Z~dM{9JXo>Tul$*K=^a2d?+P^&Ys_J#hZj KtMjive*1rKJ~__- literal 0 HcmV?d00001 diff --git a/tests/fitting.py b/tests/fitting.py new file mode 100644 index 0000000..1f6b99c --- /dev/null +++ b/tests/fitting.py @@ -0,0 +1,73 @@ +#!/usr/bin/python + +import sys +sys.path.append(".") +sys.path.append("..") +sys.path.append("../..") # Adds higher directory to python modules path + +import os +import time +import pickle +import numpy as np +from libs.consts import NOISE_FLOOR +from libs.fitting import modelfit_log_gamma +from libs.plotting import plotRSS + + +def convert_mat_to_vector(mat, block_size=0.1): + ''' + ''' + rx_locs = [] + rx_rsses = [] + width, length = mat.shape + for i in range(width): + for j in range(length): + if mat[i, j] <= NOISE_FLOOR: + continue + rx_locs.append([(i + 0.5) * block_size, (j + 0.5) * block_size]) + rx_rsses.append(mat[i, j]) + return np.array(rx_locs), np.array(rx_rsses) + + +def convert_vector_to_mat(rx_locs, rx_rsses, shape, block_size=0.1): + ''' + ''' + result = np.ones(shape) * -85.0 + for i in range(rx_locs.shape[0]): + x, y = rx_locs[i, :] / block_size + result[int(x), int(y)] = rx_rsses[i] + return result + + +def test(): + data_ori = pickle.load(open(os.path.join(os.path.dirname(__file__), "dev_map.pickle"), "rb")) + rx_locs, rx_rsses = convert_mat_to_vector(data_ori) + + min_pmse = float('inf') + best_tx_loc = None + best_tx_pwr = None + best_env_gamma = None + best_rsses = None + for i in range(10): + pmse, est_tx_loc, est_tx_pwr, est_env_gamma, est_rsses = modelfit_log_gamma(rx_locs, rx_rsses) + if pmse < min_pmse: + print("found better:", pmse, est_tx_loc, est_tx_pwr, est_env_gamma) + data_est = convert_vector_to_mat(rx_locs, est_rsses, data_ori.shape) + min_pmse = pmse + best_tx_loc = est_tx_loc + best_tx_pwr = est_tx_pwr + best_env_gamma = est_env_gamma + best_rsses = est_rsses + print("final:", min_pmse, best_tx_loc, best_tx_pwr, best_env_gamma) + data_best_est = convert_vector_to_mat(rx_locs, best_rsses, data_ori.shape) + plotRSS( + data_ori, + data_best_est, + est_tx_loc=best_tx_loc / 0.1 + ) + # print(rx_locs[:,1]) + # print(rx_locs[1,:]) + + +if __name__ == "__main__": + test() diff --git a/tests/spacemap.py b/tests/spacemap.py index f21b6a4..875d381 100644 --- a/tests/spacemap.py +++ b/tests/spacemap.py @@ -11,6 +11,7 @@ import numpy as np from libs.spacemap import SpaceBlock from libs.spacemap import SpaceRay from libs.spacemap import SpaceMap +from libs.plotting import plotSpace def test(): @@ -18,7 +19,8 @@ def test(): print(spacemap.getLosses()) print(spacemap.getLoss(np.array([1, 2, 3]), 1)) ray = SpaceRay(SpaceBlock(0.0, 0.0), SpaceBlock(2, 6.38)) - ray.computeLinePassThroughLoss(spacemap) + ray.derivePassAndReflectBlocks(spacemap) + ray.getPassThroughLoss() ray.setTotalLoss(0.0) ray.setInitPower(0, gamma=2.0) ray.setTravelDistance(0.0) @@ -26,6 +28,7 @@ def test(): print("pwr = {}".format(pwr)) print(ray.this_pass_through_blks[0].related_rays) print(ray.this_pass_through_blks[0]) + plotSpace(spacemap) if __name__ == "__main__":