Changeset 694 for mdsandtrees
 Timestamp:
 09/04/17 16:28:03 (6 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

mdsandtrees/treegenealogy.py
r690 r694 391 391 def xmin_crowd_density(self, x1, x2, y): 392 392 # TODO experimental  requires further work to make it less 'jumpy' and more predictable 393 CONST_LOCAL_AREA_RADIUS = 5 394 CONST_GLOBAL_AREA_RADIUS = 10 395 CONST_WINDOW_SIZE = 20000 #TODO should depend on the maxY ? 393 396 x1_dist_loc = 0 394 397 x2_dist_loc = 0 … … 397 400 x2_dist_glob = 0 398 401 count_glob = 1 399 miny = y 2000400 maxy = y+ 2000402 miny = yCONST_WINDOW_SIZE 403 maxy = y+CONST_WINDOW_SIZE 401 404 i_left = bisect.bisect_left(self.y_sorted, miny) 402 405 i_right = bisect.bisect_right(self.y_sorted, maxy) 403 # print("i " + str(i) + " len " + str(len(self.positions))) 404 # 405 # i = bisect.bisect_left(self.y_sorted, y) 406 # i_left = max(0, i  25) 407 # i_right = min(len(self.y_sorted), i + 25) 406 #TODO test: maxy=y should give the same results, right? 408 407 409 408 def include_pos(pos): 410 409 nonlocal x1_dist_loc, x2_dist_loc, x1_dist_glob, x2_dist_glob, count_loc, count_glob 411 dysq = (pos['y']y)**2 412 dx1 = pos['x']x1 413 dx2 = pos['x']x2 410 411 dysq = (pos['y']y)**2 + 1 #+1 so 1/dysq is at most 1 412 dx1 = math.fabs(pos['x']x1) 413 dx2 = math.fabs(pos['x']x2) 414 414 415 415 d = math.fabs(pos['x']  (x1+x2)/2) 416 416 417 if d < 10:418 x1_dist_loc += math.sqrt(d ysq + dx1**2)419 x2_dist_loc += math.sqrt(d ysq + dx2**2)417 if d < CONST_LOCAL_AREA_RADIUS: 418 x1_dist_loc += math.sqrt(dx1/dysq + dx1**2) 419 x2_dist_loc += math.sqrt(dx2/dysq + dx2**2) 420 420 count_loc += 1 421 elif d > 20:422 x1_dist_glob += math.sqrt(d ysq + dx1**2)423 x2_dist_glob += math.sqrt(d ysq + dx2**2)421 elif d > CONST_GLOBAL_AREA_RADIUS: 422 x1_dist_glob += math.sqrt(dx1/dysq + dx1**2) 423 x2_dist_glob += math.sqrt(dx2/dysq + dx2**2) 424 424 count_glob += 1 425 425 … … 433 433 pos = self.positions_sorted[random.randrange(i_left, i_right)] 434 434 include_pos(pos) 435 436 # return ((x1 if x1_dist > x2_dist else x2)437 # if x1_dist < 10000 else438 # (x1 if x1_dist < x2_dist else x2))439 440 435 441 436 return (x1 if (x1_dist_locx2_dist_loc)/count_loc(x1_dist_globx2_dist_glob)/count_glob > 0 else x2)
Note: See TracChangeset
for help on using the changeset viewer.