From 78364194176290280c456c0cb724b2acb1de646a Mon Sep 17 00:00:00 2001 From: meganrm Date: Thu, 21 Sep 2023 16:24:39 -0700 Subject: [PATCH] cleanup --- cellpack/autopack/Environment.py | 52 ++++- cellpack/autopack/ingredient/Ingredient.py | 246 ++++----------------- cellpack/autopack/ingredient/grow.py | 30 ++- 3 files changed, 98 insertions(+), 230 deletions(-) diff --git a/cellpack/autopack/Environment.py b/cellpack/autopack/Environment.py index ccf2df5e8..0b1f09d3e 100644 --- a/cellpack/autopack/Environment.py +++ b/cellpack/autopack/Environment.py @@ -1099,6 +1099,42 @@ def get_ingredient_by_name(self, name, compartment_id=None): else: return None + def get_ingredients_in_tree(self, closest_ingredients): + ingredients = [] + packed_objects = self.packed_objects.get() + if len(packed_objects): + nearby_packed_objects = [ + packed_objects[i] for i in closest_ingredients["indices"] + ] + for obj in nearby_packed_objects: + ingredients.append([obj, closest_ingredients["distances"]]) + return ingredients + + def get_closest_ingredients(self, point, cutoff=10.0): + to_return = {"indices": [], "distances": []} + numpy.zeros(self.totalNbIngr).astype("i") + nb = 0 + number_packed = len(self.packed_objects.get()) + if not number_packed: + return to_return + if self.close_ingr_bhtree is not None: + # request kdtree + nb = [] + self.log.info("finding partners") + + if number_packed >= 1: + distance, nb = self.close_ingr_bhtree.query( + point, number_packed, distance_upper_bound=cutoff + ) # len of ingr posed so far + if number_packed == 1: + distance = [distance] + nb = [nb] + to_return["indices"] = nb + to_return["distances"] = distance # sorted by distance short -> long + return to_return + else: + return to_return + def setExteriorRecipe(self, recipe): """ Set the exterior recipe with the given one. Create the weakref. @@ -1264,14 +1300,13 @@ def onePrevIngredient(self, i, mingrs, distance, nbFreePoints, marray): spacing = self.smallestProteinSize jitter = ingr.getMaxJitter(spacing) dpad = ingr.min_radius + mr + jitter - insidePoints, newDistPoints = ingr.getInsidePoints( - self.grid, - self.grid.masterGridPositions, - dpad, - distance, - centT=centT, + insidePoints, newDistPoints = ingr.get_new_distance_values( jtrans=jtrans, rotMatj=rotMatj, + gridPointsCoords=self.grid.masterGridPositions, + distance=distance, + dpad=dpad, + centT=centT, ) # update free points if len(insidePoints) and self.place_method.find("panda") != -1: @@ -1816,6 +1851,11 @@ def set_result_file_name(self, seed_basename): """ self.result_file = str(self.out_folder / f"results_{seed_basename}") + def update_after_place(self, grid_point_index): + self.order[grid_point_index] = self.lastrank + self.lastrank += 1 + self.nb_ingredient += 1 + def pack_grid( self, seedNum=0, diff --git a/cellpack/autopack/ingredient/Ingredient.py b/cellpack/autopack/ingredient/Ingredient.py index 6bb984bcd..1b45bb80d 100644 --- a/cellpack/autopack/ingredient/Ingredient.py +++ b/cellpack/autopack/ingredient/Ingredient.py @@ -487,19 +487,6 @@ def DecomposeMesh(self, poly, edit=True, copy=False, tri=True, transform=True): ) return faces, vertices, vnormals - def rejectOnce(self, rbnode, moving, afvi): - if rbnode: - self.env.callFunction(self.env.delRB, (rbnode,)) - if afvi is not None and moving is not None: - afvi.vi.deleteObject(moving) - self.haveBeenRejected = True - self.rejectionCounter += 1 - if ( - self.rejectionCounter >= self.rejection_threshold - ): # Graham set this to 6000 for figure 13b (Results Fig 3 Test1) otherwise it fails to fill small guys - self.log.info("PREMATURE ENDING of ingredient rejectOnce", self.name) - self.completion = 1.0 - def addRBsegment(self, pt1, pt2): # ovewrite by grow ingredient pass @@ -930,12 +917,12 @@ def apply_rotation(self, rot, point, origin=[0, 0, 0]): new_pos = r.apply(point) return new_pos + numpy.array(origin) - def alignRotation(self, jtrans): + def alignRotation(self, jtrans, gradients): # for surface points we compute the rotation which # aligns the principal_vector with the surface normal vx, vy, vz = v1 = self.principal_vector # surfacePointsNormals problem here - gradient_center = self.env.gradients[self.gradient].direction + gradient_center = gradients[self.gradient].direction v2 = numpy.array(gradient_center) - numpy.array(jtrans) try: rotMat = numpy.array(rotVectToVect(v1, v2), "f") @@ -987,33 +974,6 @@ def correctBB(self, p1, p2, radc): return numpy.array([numpy.array(mini).flatten(), numpy.array(maxi).flatten()]) # precised: - def checkDistSurface(self, point, cutoff): - if not hasattr(self, "histoVol"): - return False - if self.compartment_id == 0: - compartment = self.env - else: - compartment = self.env.compartments[abs(self.compartment_id) - 1] - compartment_id = self.compartment_id - if compartment_id < 0: - sfpts = compartment.surfacePointsCoords - delta = numpy.array(sfpts) - numpy.array(point) - delta *= delta - distA = numpy.sqrt(delta.sum(1)) - test = distA < cutoff - if True in test: - return True - elif compartment_id == 0: - for o in self.env.compartments: - sfpts = o.surfacePointsCoords - delta = numpy.array(sfpts) - numpy.array(point) - delta *= delta - distA = numpy.sqrt(delta.sum(1)) - test = distA < cutoff - if True in test: - return True - return False - def getListCompFromMask(self, cId, ptsInSphere): # cID ie [-2,-1,-2,0...], ptsinsph = [519,300,etc] current = self.compartment_id @@ -1030,48 +990,6 @@ def getListCompFromMask(self, cId, ptsInSphere): liste = [i for i, x in enumerate(cId) if x == current] return liste - def isInGoodComp(self, pId, nbs=None): - # cID ie [-2,-1,-2,0...], ptsinsph = [519,300,etc] - current = self.compartment_id - cId = self.env.grid.compartment_ids[pId] - if current <= 0: # inside - if current != cId: - return False - return True - if current > 0: # surface - if current != cId and -current != cId: - return False - return True - return False - - def checkCompartment(self, ptsInSphere, nbs=None): - trigger = False - if self.compareCompartment: - cId = numpy.take( - self.env.grid.compartment_ids, ptsInSphere, 0 - ) # shoud be the same ? - if nbs is not None: - if self.compartment_id <= 0 and nbs != 0: - return trigger, True - L = self.getListCompFromMask(cId, ptsInSphere) - - if len(cId) <= 1: - return trigger, True - p = float(len(L)) / float( - len(cId) - ) # ratio accepted compId / totalCompId-> want 1.0 - if p < self.compareCompartmentTolerance: - trigger = True - return trigger, True - # threshold - if ( - self.compareCompartmentThreshold != 0.0 - and p < self.compareCompartmentThreshold - ): - return trigger, True - # reject the ingr - return trigger, False - def get_new_distances_and_inside_points( self, env, @@ -1205,20 +1123,6 @@ def get_new_jitter_location_and_rotation( return self.oneJitter(env, starting_pos, starting_rotation) - def getInsidePoints( - self, - grid, - gridPointsCoords, - dpad, - distance, - centT=None, - jtrans=None, - rotMatj=None, - ): - return self.get_new_distance_values( - jtrans, rotMatj, gridPointsCoords, distance, dpad - ) - def getIngredientsInBox(self, env, jtrans, rotMat, compartment): if env.windowsSize_overwrite: radius = env.windowsSize @@ -1278,28 +1182,16 @@ def getIngredientsInBox(self, env, jtrans, rotMat, compartment): return ingredients - def getIngredientsInTree(self, closest_ingredients): - ingredients = [] - packed_objects = self.env.packed_objects.get() - if len(packed_objects): - nearby_packed_objects = [ - packed_objects[i] for i in closest_ingredients["indices"] - ] - for obj in nearby_packed_objects: - ingredients.append([obj, closest_ingredients["distances"]]) - return ingredients - - def get_partners(self, jtrans, rotMat, organelle, afvi): - env = self.env - closest_ingredients = self.get_closest_ingredients( - jtrans, cutoff=self.env.grid.diag + def get_partners(self, env, jtrans, rotMat, organelle, afvi): + closest_ingredients = env.get_closest_ingredients( + jtrans, cutoff=env.grid.diag ) if not len(closest_ingredients["indices"]): near_by_ingredients = self.getIngredientsInBox( env, jtrans, rotMat, organelle ) else: - near_by_ingredients = self.getIngredientsInTree(closest_ingredients) + near_by_ingredients = env.get_ingredients_in_tree(closest_ingredients) placed_partners = [] if not len(near_by_ingredients): self.log.info("no close ingredient found") @@ -1340,16 +1232,6 @@ def get_partners(self, jtrans, rotMat, organelle, afvi): else: return near_by_ingredients, placed_partners - def getTransform(self): - tTrans = self.vi.ToVec(self.vi.getTranslation(self.moving)) - self.htrans.append(tTrans) - avg = numpy.average(numpy.array(self.htrans)) - d = self.vi.measure_distance(tTrans, avg) - if d < 5.0: - return True - else: - return False - def get_new_pos(self, ingr, pos, rot, positions_to_adjust): """ Takes positions_to_adjust, such as an array of spheres at a level in a @@ -1541,50 +1423,6 @@ def get_rbNodes( self.env.nodes = nodes return nodes - def getClosePairIngredient(self, point, histoVol, cutoff=10.0): - R = {"indices": [], "distances": []} - radius = [ingr.encapsulating_radius for ingr in self.env.rIngr] - radius.append(self.encapsulating_radius) - pos = self.env.rTrans[:] # ).tolist() - pos.append([point[0], point[1], point[2]]) - ind = len(pos) - 1 - bht = spatial.cKDTree(pos, leafsize=10) - # find all pairs for which the distance is less than 1.1 - # times the sum of the radii - pairs = bht.query_ball_point(pos, radius) - for p in pairs: - if p[0] == ind: - R["indices"].append(p[1]) - elif p[1] == ind: - R["indices"].append(p[0]) - return R - - def get_closest_ingredients(self, point, cutoff=10.0): - to_return = {"indices": [], "distances": []} - env = self.env - numpy.zeros(env.totalNbIngr).astype("i") - nb = 0 - number_packed = len(env.packed_objects.get()) - if not number_packed: - return to_return - if env.close_ingr_bhtree is not None: - # request kdtree - nb = [] - self.log.info("finding partners") - - if number_packed >= 1: - distance, nb = env.close_ingr_bhtree.query( - point, number_packed, distance_upper_bound=cutoff - ) # len of ingr posed so far - if number_packed == 1: - distance = [distance] - nb = [nb] - to_return["indices"] = nb - to_return["distances"] = distance # sorted by distance short -> long - return to_return - else: - return to_return - def update_data_tree(self): if len(self.env.packed_objects.get()) >= 1: self.env.close_ingr_bhtree = spatial.cKDTree( @@ -1665,13 +1503,10 @@ def place( dropped_rotation, grid_point_index, new_inside_points, - new_dist_values, ): self.nbPts = self.nbPts + len(new_inside_points) - env.order[grid_point_index] = env.lastrank - env.lastrank += 1 - env.nb_ingredient += 1 + env.update_after_place(grid_point_index) if self.packing_mode[-4:] == "tile": nexthexa = self.tilling.dropTile( @@ -1747,6 +1582,7 @@ def attempt_to_pack_at_grid_location( # I think this should be done ealier, when we're getting the point index if self.packing_mode == "closePartner": target_grid_point_position, rotation_matrix = self.close_partner_check( + env, target_grid_point_position, rotation_matrix, compartment, @@ -1869,7 +1705,7 @@ def attempt_to_pack_at_grid_location( autopack.helper.set_object_static( current_visual_instance, jtrans, rotMatj ) - self.place(env, jtrans, rotMatj, ptInd, insidePoints, newDistPoints) + self.place(env, jtrans, rotMatj, ptInd, insidePoints) else: if is_realtime: self.remove_from_realtime_display(current_visual_instance) @@ -1902,9 +1738,9 @@ def get_rotation(self, pt_ind, env, compartment): elif ( self.use_orient_bias and self.packing_mode == "gradient" ): # you need a gradient here - rot_mat = self.alignRotation(env.grid.masterGridPositions[pt_ind]) + rot_mat = self.alignRotation(env.grid.masterGridPositions[pt_ind], env.gradients) else: - rot_mat = self.env.helper.rotation_matrix( + rot_mat = env.helper.rotation_matrix( random() * self.rotation_range, self.rotation_axis ) # for other points we get a random rotation @@ -1912,7 +1748,7 @@ def get_rotation(self, pt_ind, env, compartment): rot_mat = env.randomRot.get() return rot_mat - def randomize_rotation(self, rotation, histovol): + def randomize_rotation(self, rotation, env): # randomize rotation about axis jitter_rotation = numpy.identity(4) if self.compartment_id > 0: @@ -1927,12 +1763,12 @@ def randomize_rotation(self, rotation, histovol): jitter_rotation = self.getBiasedRotation(rotation, weight=None) else: # should we align to this rotation_axis ? - jitter_rotation = self.env.helper.rotation_matrix( + jitter_rotation = env.helper.rotation_matrix( random() * self.rotation_range, self.rotation_axis ) else: - if histovol is not None: - jitter_rotation = histovol.randomRot.get() + if env is not None: + jitter_rotation = env.randomRot.get() if self.rotation_range != 0.0: return jitter_rotation else: @@ -1980,7 +1816,7 @@ def update_display_rt(self, current_instance, translation, rotation): def rigid_place( self, - histoVol, + env, ptInd, compartment, target_grid_point_position, @@ -1993,14 +1829,14 @@ def rigid_place( """ drop the ingredient on grid point ptInd """ - afvi = histoVol.afviewer - simulationTimes = histoVol.simulationTimes - runTimeDisplay = histoVol.runTimeDisplay - springOptions = histoVol.springOptions + afvi = env.afviewer + simulationTimes = env.simulationTimes + runTimeDisplay = env.runTimeDisplay + springOptions = env.springOptions is_realtime = moving is not None jtrans, rotMatj = self.oneJitter( - histoVol, target_grid_point_position, rotation_matrix + env, target_grid_point_position, rotation_matrix ) # here should go the simulation @@ -2016,7 +1852,7 @@ def rigid_place( self.update_display_rt(moving, jtrans, rotMatj) # 2- get the neighboring object from ptInd near_by_ingredients, placed_partners = self.get_partners( - histoVol, jtrans, rotation_matrix, compartment, afvi + env, jtrans, rotation_matrix, compartment, afvi ) for i, elem in enumerate(near_by_ingredients): ing = elem[2] @@ -2061,7 +1897,7 @@ def rigid_place( targetPoint = jtrans # setup the target position if self.place_method == "spring": - afvi.vi.setRigidBody(afvi.movingMesh, **histoVol.dynamicOptions["spring"]) + afvi.vi.setRigidBody(afvi.movingMesh, **env.dynamicOptions["spring"]) # target can be partner position? target = afvi.vi.getObject("target" + name) if target is None: @@ -2080,9 +1916,9 @@ def rigid_place( ) else: # before assigning should get outside thge object - afvi.vi.setRigidBody(afvi.movingMesh, **histoVol.dynamicOptions["moving"]) + afvi.vi.setRigidBody(afvi.movingMesh, **env.dynamicOptions["moving"]) afvi.vi.setTranslation(self.moving, pos=targetPoint) - afvi.vi.setRigidBody(afvi.staticMesh, **histoVol.dynamicOptions["static"]) + afvi.vi.setRigidBody(afvi.staticMesh, **env.dynamicOptions["static"]) # 4- we run the simulation # c4d.documents.RunAnimation(c4d.documents.GetActiveDocument(), False,True) # if runTimeDisplay : @@ -2114,8 +1950,8 @@ def rigid_place( insidePoints = {} newDistPoints = {} insidePoints, newDistPoints = self.get_new_distance_values( - histoVol.grid, - histoVol.masterGridPositions, + env.grid, + env.masterGridPositions, dpad, distance, centT, @@ -2126,10 +1962,10 @@ def rigid_place( # save dropped ingredient - histoVol.rTrans.append(jtrans) - histoVol.result.append([jtrans, rotMatj, self, ptInd]) - histoVol.rRot.append(rotMatj) - histoVol.rIngr.append(self) + env.rTrans.append(jtrans) + env.result.append([jtrans, rotMatj, self, ptInd]) + env.rRot.append(rotMatj) + env.rIngr.append(self) self.rRot.append(rotMatj) self.tTrans.append(jtrans) @@ -2284,9 +2120,9 @@ def jitter_place( ) return False, packing_location, packing_rotation, {}, {} - def lookForNeighbours(self, trans, rotMat, organelle, afvi): + def lookForNeighbours(self, env, trans, rotMat, organelle, afvi): near_by_ingredients, placed_partners = self.get_partners( - trans, rotMat, organelle, afvi + env, trans, rotMat, organelle, afvi ) targetPoint = trans found = False @@ -2330,7 +2166,7 @@ def pandaBullet_collision(self, pos, rot, rbnode, getnodes=False): collision = False liste_nodes = [] if len(self.env.rTrans) != 0: - closesbody_indice = self.get_closest_ingredients( + closesbody_indice = self.env.get_closest_ingredients( pos, cutoff=self.env.largestProteinSize + self.encapsulating_radius * 2.0, ) # vself.radii[0][0]*2.0 @@ -2362,8 +2198,9 @@ def get_compartment(self, env): else: return env.compartments[abs(self.compartment_id) - 1] - def close_partner_check(self, translation, rotation, compartment, afvi, moving): + def close_partner_check(self, env, translation, rotation, compartment, afvi, moving): target_point, rot_matrix, found = self.lookForNeighbours( + env, translation, rotation, compartment, @@ -2496,8 +2333,6 @@ def pandaBullet_place( inside_points = {} new_dist_points = {} t3 = time() - - # self.update_data_tree(jtrans,rotMatj,ptInd=ptInd)? self.env.static.append(rbnode) self.env.moving = None @@ -2515,11 +2350,6 @@ def pandaBullet_place( ) self.log.info("compute distance loop %d", time() - t3) - # rebuild kdtree - # if len(self.env.rTrans) >= 1: - # self.env.close_ingr_bhtree = spatial.cKDTree( - # self.env.rTrans, leafsize=10 - # ) if self.packing_mode[-4:] == "tile": self.tilling.dropTile( self.tilling.idc, @@ -2637,8 +2467,8 @@ def spheres_SST_place( # ) # collision_results.extend([collision]) if True not in collision_results: - self.env.static.append(rbnode) - self.env.moving = None + env.static.append(rbnode) + env.moving = None for pt in pts_to_check: new_inside_pts, new_dist_points = self.pack_at_grid_pt_location( diff --git a/cellpack/autopack/ingredient/grow.py b/cellpack/autopack/ingredient/grow.py index 37b74560c..08ea6aa29 100644 --- a/cellpack/autopack/ingredient/grow.py +++ b/cellpack/autopack/ingredient/grow.py @@ -905,7 +905,7 @@ def walkSpherePanda( self.vi.update() liste_nodes = [] cutoff = self.env.largestProteinSize + self.uLength - closesbody_indice = self.get_closest_ingredients(pt2, self.env, cutoff=cutoff) + closesbody_indice = self.env.get_closest_ingredients(pt2, self.env, cutoff=cutoff) liste_nodes = self.get_rbNodes( closesbody_indice, pt2, prevpoint=pt1, getInfo=True ) @@ -1138,7 +1138,7 @@ def walkSpherePanda( prev = None if len(self.env.rTrans) > 2: prev = self.env.rTrans[-1] - closesbody_indice = self.get_closest_ingredients( + closesbody_indice = self.env.get_closest_ingredients( newPt, histoVol, cutoff=cutoff ) # vself.radii[0][0]*2.0 if len(closesbody_indice) == 0: @@ -1330,7 +1330,7 @@ def resetLastPoint(self, listePtCurve): self.env.result.pop(len(self.env.result) - 1) # rebuild kdtree if len(self.env.rTrans) > 1: - self.env.close_ingr_bhtree = spatial.cKDTree(self.env.rTrans, leafsize=10) + self.env.close_ingr_bhtree = spatial.cKDTree(self.env.packed_objects.get_positions(), leafsize=10) # also remove from the result ? self.results.pop(len(self.results) - 1) @@ -1527,14 +1527,13 @@ def grow( # self.positions2=[[cent2T],] # rbnode = histoVol.callFunction(histoVol.addRB,(self, numpy.array(jtrans), numpy.array(rotMatj),),{"rtype":self.type},)#cylinder # histoVol.callFunction(histoVol.moveRBnode,(rbnode, jtrans, rotMatj,)) - insidePoints, newDistPoints = self.getInsidePoints( - histoVol.grid, - gridPointsCoords, - dpad, - distance, - centT=cent1T, + insidePoints, newDistPoints = self.get_new_distance_values( jtrans=jtrans, rotMatj=rotMatj, + gridPointsCoords=gridPointsCoords, + distance=distance, + dpad=dpad, + centT=cent1T, ) nbFreePoints = BaseGrid.updateDistances( @@ -1593,14 +1592,13 @@ def updateGrid( for i in range(rg): # len(self.results)): jtrans, rotMatj = self.results[-i] cent1T = self.transformPoints(jtrans, rotMatj, self.positions[-1]) - new_inside_pts, new_dist_points = self.getInsidePoints( - histoVol.grid, - gridPointsCoords, - dpad, - distance, - centT=cent1T, + new_inside_pts, new_dist_points = self.get_new_distance_values( jtrans=jtrans, rotMatj=rotMatj, + gridPointsCoords=gridPointsCoords, + distance=distance, + dpad=dpad, + centT=cent1T, ) insidePoints = self.merge_place_results(new_inside_pts, insidePoints) newDistPoints = self.merge_place_results(new_dist_points, newDistPoints) @@ -1772,7 +1770,7 @@ def grow_place( # rebuild kdtree if len(self.env.rTrans) > 1: - self.env.close_ingr_bhtree = spatial.cKDTree(self.env.rTrans, leafsize=10) + self.env.close_ingr_bhtree = spatial.cKDTree(self.env.packed_objects.get_positions(), leafsize=10) self.currentLength = 0.0 # self.Ptis=[ptInd,histoVol.grid.getPointFrom3D(secondPoint)]