diff --git a/cellpack/autopack/Analysis.py b/cellpack/autopack/Analysis.py index d0b67b95e..6105a2b06 100644 --- a/cellpack/autopack/Analysis.py +++ b/cellpack/autopack/Analysis.py @@ -170,7 +170,7 @@ def getDistanceFrom(self, target, parents=None, **options): objects : list of object or list of points """ # get distance from object to the target. - # all object are in h.molecules and orga.molecules + # all object are in env.packed_objects # get options if type(target) == list or type(target) == tuple: @@ -447,23 +447,6 @@ def rectangle_circle_area(self, bbox, center, radius): area = bbox[0][1] ** 2 return area - def getAxeValue(self, ingrname, axe=0): - ingredient_positions = [ - self.env.molecules[i][0][axe] - for i in range(len(self.env.molecules)) - if self.env.molecules[i][2].name == ingrname - ] - return ingredient_positions - - def getAxesValues(self, positions): - pp = numpy.array(positions).transpose() - if len(positions) == 0: - return 1, 1, 1 - px = pp[0] - py = pp[1] - pz = pp[2] - return px, py, pz - def getVolumeShell(self, bbox, radii, center): # rectangle_circle_area volumes = [] @@ -473,9 +456,6 @@ def getVolumeShell(self, bbox, radii, center): r2 = radii[i + 1] v1 = self.g.calc_volume(r1, box_size0 / 2.0) v2 = self.g.calc_volume(r2, box_size0 / 2.0) - # if v1 == 0 or v2 == 0 : - # volumes.append((4./3.)*numpy.pi*(numpy.power(r2,3)-numpy.power(r1, 3))) - # else : volumes.append(v2 - v1) return volumes @@ -1774,7 +1754,7 @@ def pack( t2 = time() run_time = t2 - t1 print(f"time to run pack_grid for {self.env.place_method}: {run_time:0.2f}") - print(f"num placed: {len(self.env.molecules)}") + print(f"num placed: {len(self.env.packed_objects.get_ingredients())}") if show_plotly_plot: min_bound, max_bound = self.env.get_bounding_box_limits() width = max_bound - min_bound @@ -1787,7 +1767,7 @@ def pack( range=[min_bound[1] - 0.2 * width[1], max_bound[1] + 0.2 * width[1]] ) self.plotly.update_title( - f"{self.env.place_method} took {str(round(t2 - t1, 2))}s, packed {len(self.env.molecules)}" + f"{self.env.place_method} took {str(round(t2 - t1, 2))}s, packed {len(self.env.packed_objects.get_ingredients())}" ) self.plotly.make_grid_heatmap(self.env) self.plotly.add_ingredient_positions(self.env) @@ -1875,29 +1855,6 @@ def plot_one_result_2d( plt.close() # closes the current figure return ingrpos - # res=plotOneResult(None,filename="results_seed_8.json") - - def plot_one_result_3D(self, filename, width=1000.0): - plt.close("all") # closes the current figure - pos = [] - s = [] - c = [] - for i in range(len(self.env.molecules)): - m = self.env.molecules[i] - pos.append(numpy.array(m[0]).tolist()) - s.append(m[2].encapsulating_radius ** 2) - c.append(m[2].color) - fig = plt.figure() - ax = fig.gca(projection="3d") - x, y, z = numpy.array(pos).transpose() - ax.scatter(x, y, z, s=s, c=c) - ax.legend() - ax.set_xlim3d(0, width) - ax.set_ylim3d(0, width) - ax.set_zlim3d(0, width) - plt.savefig(filename) - return x, y, z, s, c - def set_ingredient_color(self, ingr): """ Sets the color of an ingredient diff --git a/cellpack/autopack/Compartment.py b/cellpack/autopack/Compartment.py index 693172ad6..6d93c80f7 100644 --- a/cellpack/autopack/Compartment.py +++ b/cellpack/autopack/Compartment.py @@ -68,6 +68,10 @@ import cellpack.autopack as autopack from cellpack.autopack import transformation as tr, binvox_rw from cellpack.autopack.BaseGrid import gridPoint +from cellpack.autopack.interface_objects.packed_objects import ( + PackedObject, + PackedObjects, +) from .Recipe import Recipe from .ray import ( makeMarchingCube, @@ -138,6 +142,7 @@ def __init__(self, name, object_info): self.ref_obj = None self.meshType = None self.representations = object_info["representations"] + self.packed_objects = PackedObjects() if self.representations.has_mesh(): self.gname = self.representations.get_mesh_name() self.meshType = self.representations.get_mesh_format() @@ -172,7 +177,6 @@ def __init__(self, name, object_info): # is added to a Environment. Positivefor surface pts # negative for interior points # self.parent = None - self.molecules = [] # list of ( (x,y,z), rotation, ingredient) triplet generated by fill self.overwriteSurfacePts = True # do we discretize surface point per edges @@ -207,7 +211,6 @@ def reset(self): # list of grid point indices on compartment surface self.surfacePoints = None self.surfacePointsNormals = {} # will be point index:normal - # self.molecules = [] def transformMesh(self, pos, rotation): rot = tr.matrix_from_quaternion(rotation).transpose() @@ -266,6 +269,17 @@ def initialize_shape(self, mesh_store): self.encapsulating_radius = radius self.radius = mesh_store.get_smallest_radius(self.gname, center) + def store_packed_object(self, env): + packed_object = PackedObject( + position=[0, 0, 0], + rotation=numpy.identity(4), + radius=self.radius, + pt_index=-1, + ingredient=self, + is_compartment=True, + ) + env.packed_objects.add(packed_object) + def addShapeRB(self, env): # in case our shape is a regular primitive if self.stype == "capsule": @@ -503,10 +517,10 @@ def setInnerRecipe(self, recipe): self.innerRecipe.number = self.number recipe.compartment = self # weakref.ref(self) for ingr in recipe.ingredients: - ingr.compNum = -self.number + ingr.compartment_id = -self.number if hasattr(ingr, "compMask"): if not ingr.compMask: - ingr.compMask = [ingr.compNum] + ingr.compMask = [ingr.compartment_id] def setSurfaceRecipe(self, recipe): """set the inner recipe that define the ingredient to pack at the surface""" @@ -516,7 +530,7 @@ def setSurfaceRecipe(self, recipe): self.surfaceRecipe.number = self.number recipe.compartment = self # weakref.ref(self) for ingr in recipe.ingredients: - ingr.compNum = self.number + ingr.compartment_id = self.number def getCenter(self): """get the center of the mesh (vertices barycenter)""" @@ -3200,12 +3214,19 @@ def create_voxelization( voxel_size, image_size, position, - mesh_store, - hollow=False, + **kwargs, ): """ Creates a voxelization mask at the position of the compartment """ + if "mesh_store" not in kwargs: + raise RuntimeError( + "mesh_store must be provided for compartment voxelization" + ) + + mesh_store = kwargs["mesh_store"] + hollow = kwargs.get("hollow", False) + relative_position = position - bounding_box[0] mask = self.create_voxelized_mask( *image_size, diff --git a/cellpack/autopack/Environment.py b/cellpack/autopack/Environment.py index c26b0918c..9545e6148 100644 --- a/cellpack/autopack/Environment.py +++ b/cellpack/autopack/Environment.py @@ -49,6 +49,7 @@ import os from time import time from random import random, uniform, seed +from cellpack.autopack.interface_objects.packed_objects import PackedObjects from scipy import spatial import numpy import pickle @@ -70,6 +71,9 @@ from cellpack.autopack.utils import ( cmp_to_key, expand_object_using_key, + get_value_from_distribution, + get_max_value_from_distribution, + get_min_value_from_distribution, ingredient_compare0, ingredient_compare1, ingredient_compare2, @@ -86,7 +90,6 @@ # backward compatibility with kevin method from cellpack.autopack.BaseGrid import BaseGrid as BaseGrid -from .trajectory import dcdTrajectory, molbTrajectory from .randomRot import RandomRot from tqdm import tqdm @@ -197,10 +200,6 @@ def __init__(self, config=None, recipe=None): None # Oct 20, 2012 Graham wonders if this is part of the problem ) self.freePointMask = None - self.molecules = ( - [] - ) # list of ( (x,y,z), rotation, ingredient) triplet generated by packing - self.ingr_result = {} self.randomRot = RandomRot() # the class used to generate random rotation self.activeIngr = [] @@ -253,10 +252,9 @@ def __init__(self, config=None, recipe=None): self.close_ingr_bhtree = ( None # RBHTree(a.tolist(),range(len(a)),10,10,10,10,10,9999) ) - self.rTrans = [] - self.result = [] - self.rIngr = [] - self.rRot = [] + + self.packed_objects = PackedObjects() + # should be part of an independent module self.panda_solver = "bullet" # or bullet # could be a problem here for pp @@ -387,36 +385,15 @@ def reportprogress(self, label=None, progress=None): def set_partners_ingredient(self, ingr): if ingr.partners is not None: for partner in ingr.partners.all_partners: - partner_ingr = self.getIngrFromName(partner.name) + partner_ingr = self.get_ingredient_by_name(partner.name) partner.set_ingredient(partner_ingr) if ingr.type == "Grow": # TODO: I don't think this code is needed, # but I haven't dug into it enough to delete it all yet ingr.prepare_alternates() - def get_positions_for_ingredient(self, ingredient_name): - return numpy.array( - [ - self.molecules[i][0] - for i in range(len(self.molecules)) - if self.molecules[i][2].name == ingredient_name - ] - ) - - def get_rotations_for_ingredient(self, ingredient_name): - return numpy.array( - [ - self.molecules[i][1] - for i in range(len(self.molecules)) - if self.molecules[i][2].name == ingredient_name - ] - ) - - def get_all_positions(self): - return numpy.array([self.molecules[i][0] for i in range(len(self.molecules))]) - def get_all_distances(self, position=None): - positions = self.get_all_positions() + positions = self.packed_objects.get_positions() if len(positions) == 0: return numpy.array([]) elif position is not None: @@ -426,7 +403,9 @@ def get_all_distances(self, position=None): def get_distances(self, ingredient_name, center): - ingredient_positions = self.get_positions_for_ingredient(ingredient_name) + ingredient_positions = self.packed_objects.get_positions_for_ingredient( + ingredient_name + ) if len(ingredient_positions): distances_between_ingredients = spatial.distance.pdist(ingredient_positions) @@ -444,7 +423,7 @@ def get_distances(self, ingredient_name, center): ) def get_ingredient_angles(self, ingredient_name, center, ingredient_positions): - ingredient_rotation = self.get_rotations_for_ingredient( + ingredient_rotation = self.packed_objects.get_rotations_for_ingredient( ingredient_name=ingredient_name, ) ingredient_position_vector = numpy.array(ingredient_positions) - numpy.array( @@ -504,16 +483,18 @@ def calc_pairwise_distances(self, ingr1name, ingr2name): """ Returns pairwise distances between ingredients of different types """ - ingr_pos_1 = self.get_positions_for_ingredient(ingredient_name=ingr1name) - ingr_pos_2 = self.get_positions_for_ingredient(ingredient_name=ingr2name) + ingr_pos_1 = self.packed_objects.get_positions_for_ingredient( + ingredient_name=ingr1name + ) + ingr_pos_2 = self.packed_objects.get_positions_for_ingredient( + ingredient_name=ingr2name + ) return numpy.ravel(spatial.distance.cdist(ingr_pos_1, ingr_pos_2)) def save_result( self, free_points, distances, - t0, - vTestid, all_ingr_as_array, save_grid_logs=False, save_result_as_file=False, @@ -530,57 +511,15 @@ def save_result( self.collectResultPerIngredient() if save_result_as_file: self.store() - self.store_asTxt() Writer(format=self.format_output).save( self, - kwds=["compNum"], + kwds=["compartment_id"], result=True, quaternion=True, all_ingr_as_array=all_ingr_as_array, compartments=self.compartments, ) - self.log.info("time to save result file %d", time() - t0) - self.log.info("self.compartments In Environment = %d", len(self.compartments)) - if self.compartments == []: - unitVol = self.grid.gridSpacing**3 - innerPointNum = len(free_points) - self.log.info(" . . . . ") - self.log.info("inner Point Count = %d", innerPointNum) - self.log.info("innerVolume temp Confirm = %d", innerPointNum * unitVol) - usedPts = 0 - unUsedPts = 0 - # fpts = self.freePointsAfterFill - vDistanceString = "" - for i in range(innerPointNum): - pt = free_points[i] # fpts[i] - # for pt in self.histo.freePointsAfterFill:#[:self.histo.nbFreePointsAfterFill]: - d = self.distancesAfterFill[pt] - vDistanceString += str(d) + "\n" - if d <= 0: # >self.smallestProteinSize-0.001: - usedPts += 1 - else: - unUsedPts += 1 - # Graham Note: There is overused disk space- the rotation matrix is 4x4 with an offset of 0,0,0 - # and we have a separate translation vector in the results and molecules arrays. - # Get rid of the translation vector and move it to the rotation matrix to save space... - # will that slow the time it takes to extract the vector from the matrix when we need to call it? - - self.log.info("unitVolume2 = %d", unitVol) - self.log.info("Number of Points Unused = %d", unUsedPts) - self.log.info("Number of Points Used = %d", usedPts) - self.log.info("Volume Used = %d", usedPts * unitVol) - self.log.info("Volume Unused = %d", unUsedPts * unitVol) - self.log.info("vTestid = %s", vTestid) - self.log.info("self.nbGridPoints = %r", self.grid.nbGridPoints) - self.log.info("self.gridVolume = %d", self.grid.gridVolume) - self.log.info("histoVol.timeUpDistLoopTotal = %d", self.timeUpDistLoopTotal) - - # END Analysis Tools: Graham added back this big chunk of code - # for analysis tools and graphic on 5/16/12 - # Needs to be cleaned up into a function and proper uPy code - self.log.info("time to save end %d", time() - t0) - def loadResult( self, resultfilename=None, restore_grid=True, backward=False, transpose=True ): @@ -1064,7 +1003,6 @@ def _add_compartment(self, compartment, parent): compartment.setNumber(self.nbCompartments) self.nbCompartments += 1 self.compartments.append(compartment) - CompartmentList.add_compartment(parent, compartment) def compartment_id_for_nearest_grid_point(self, point): @@ -1122,11 +1060,11 @@ def getIngrFromNameInRecipe(self, name, r): # return ingr return None - def getIngrFromName(self, name, compNum=None): + def get_ingredient_by_name(self, name, compartment_id=None): """ Given an ingredient name and optionally the compartment number, retrieve the ingredient object instance """ - if compNum is None: + if compartment_id is None: r = self.exteriorRecipe ingr = self.getIngrFromNameInRecipe(name, r) if ingr is not None: @@ -1140,15 +1078,15 @@ def getIngrFromName(self, name, compNum=None): ingr = self.getIngrFromNameInRecipe(name, ri) if ingr is not None: return ingr - elif compNum == 0: + elif compartment_id == 0: r = self.exteriorRecipe ingr = self.getIngrFromNameInRecipe(name, r) if ingr is not None: return ingr else: return None - elif compNum > 0: - o = self.compartments[compNum - 1] + elif compartment_id > 0: + o = self.compartments[compartment_id - 1] rs = o.surfaceRecipe ingr = self.getIngrFromNameInRecipe(name, rs) if ingr is not None: @@ -1156,7 +1094,7 @@ def getIngrFromName(self, name, compNum=None): else: return None else: # <0 - o = self.compartments[(compNum * -1) - 1] + o = self.compartments[(compartment_id * -1) - 1] ri = o.innerRecipe ingr = self.getIngrFromNameInRecipe(name, ri) if ingr is not None: @@ -1164,6 +1102,42 @@ def getIngrFromName(self, name, compNum=None): else: return None + def get_ingredients_in_tree(self, closest_ingredients): + ingredients = [] + packed_objects = self.packed_objects.get_ingredients() + 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_ingredients()) + 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. @@ -1172,7 +1146,7 @@ def setExteriorRecipe(self, recipe): self.exteriorRecipe = recipe recipe.compartment = self # weakref.ref(self) for ingr in recipe.ingredients: - ingr.compNum = 0 + ingr.compartment_id = 0 def BuildCompartmentsGrids(self): """ @@ -1286,19 +1260,20 @@ def buildGrid(self, rebuild=True): else: self.grid.distToClosestSurf_store = self.grid.distToClosestSurf[:] - distance = self.grid.distToClosestSurf # [:] + # distance = self.grid.distToClosestSurf # [:] nbFreePoints = nbPoints # -1 - for i, mingrs in enumerate(self.molecules): # ( jtrans, rotMatj, self, ptInd ) - nbFreePoints = self.onePrevIngredient( - i, mingrs, distance, nbFreePoints, self.molecules - ) - for organelle in self.compartments: - for i, mingrs in enumerate( - organelle.molecules - ): # ( jtrans, rotMatj, self, ptInd ) - nbFreePoints = self.onePrevIngredient( - i, mingrs, distance, nbFreePoints, organelle.molecules - ) + # TODO: refactor this to work with new placed_objects data structure + # for i, mingrs in enumerate(self.molecules): # ( jtrans, rotMatj, self, ptInd ) + # nbFreePoints = self.onePrevIngredient( + # i, mingrs, distance, nbFreePoints, self.molecules + # ) + # for organelle in self.compartments: + # for i, mingrs in enumerate( + # organelle.molecules + # ): # ( jtrans, rotMatj, self, ptInd ) + # nbFreePoints = self.onePrevIngredient( + # i, mingrs, distance, nbFreePoints, organelle.molecules + # ) self.grid.nbFreePoints = nbFreePoints if self.use_gradient and len(self.gradients) and rebuild: @@ -1321,22 +1296,21 @@ def onePrevIngredient(self, i, mingrs, distance, nbFreePoints, marray): """ Unused """ - jtrans, rotMatj, ingr, ptInd = mingrs + jtrans, rotMatj, ingr, ptInd, _ = mingrs centT = ingr.transformPoints(jtrans, rotMatj, ingr.positions[-1]) insidePoints = {} newDistPoints = {} - mr = self.get_dpad(ingr.compNum) + mr = self.get_dpad(ingr.compartment_id) 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: @@ -1507,15 +1481,13 @@ def reset(self): self.jitterLength = 0.0 r = self.exteriorRecipe self.resetIngrRecip(r) - self.molecules = [] + self.packed_objects = PackedObjects() for orga in self.compartments: orga.reset() rs = orga.surfaceRecipe self.resetIngrRecip(rs) ri = orga.innerRecipe self.resetIngrRecip(ri) - orga.molecules = [] - self.ingr_result = {} if self.world is not None: # need to clear all node # nodes = self.rb_panda[:] @@ -1527,11 +1499,6 @@ def reset(self): del self.octree self.octree = None # the reset doesnt touch the grid... - - self.rTrans = [] - self.rIngr = [] - self.rRot = [] - self.result = [] # rapid node ? def resetIngrRecip(self, recip): @@ -1572,9 +1539,6 @@ def getActiveIng(self): """Return all remaining active ingredients""" allIngredients = [] recipe = self.exteriorRecipe - if recipe is not None: - if not hasattr(recipe, "molecules"): - recipe.molecules = [] if recipe: for ingr in recipe.ingredients: ingr.counter = 0 # counter of placed molecules @@ -1585,8 +1549,6 @@ def getActiveIng(self): ingr.completion = 1.0 for compartment in self.compartments: - if not hasattr(compartment, "molecules"): - compartment.molecules = [] recipe = compartment.surfaceRecipe if recipe: for ingr in recipe.ingredients: @@ -1643,19 +1605,28 @@ def pickIngredient(self, vThreshStart, verbose=0): # print (r,n,ingr.name,len(self.activeIngr)) #Graham turned back on 5/16/12, but may be costly return ingr - def get_dpad(self, compNum): + def get_dpad(self, compartment_id): """Return the largest encapsulating_radius and use it for padding""" mr = 0.0 - if compNum == 0: # cytoplasm -> use cyto and all surfaces + if compartment_id == 0: # cytoplasm -> use cyto and all surfaces for ingr1 in self.activeIngr: - if ingr1.compNum >= 0: - r = ingr1.encapsulating_radius + if ingr1.compartment_id >= 0: + if hasattr(ingr1, "max_radius"): + r = ingr1.max_radius + else: + r = ingr1.encapsulating_radius if r > mr: mr = r else: for ingr1 in self.activeIngr: - if ingr1.compNum == compNum or ingr1.compNum == -compNum: - r = ingr1.encapsulating_radius + if ( + ingr1.compartment_id == compartment_id + or ingr1.compartment_id == -compartment_id + ): + if hasattr(ingr1, "max_radius"): + r = ingr1.max_radius + else: + r = ingr1.encapsulating_radius if r > mr: mr = r return mr @@ -1829,32 +1800,7 @@ def prep_molecules_for_save(self, distances, free_points, nbFreePoints): if self.runTimeDisplay and autopack.helper.host == "simularium": autopack.helper.writeToFile("./realtime", self.boundingBox) - - if self.afviewer is not None and hasattr(self.afviewer, "vi"): - self.afviewer.vi.progressBar(label="Filling Complete") - self.afviewer.vi.resetProgressBar() - ingredients = {} - all_ingr_as_array = self.molecules - for pos, rot, ingr, ptInd in self.molecules: - if ingr.name not in ingredients: - ingredients[ingr.name] = [ingr, [], [], []] - mat = rot.copy() - mat[:3, 3] = pos - ingredients[ingr.name][1].append(pos) - ingredients[ingr.name][2].append(rot) - ingredients[ingr.name][3].append(numpy.array(mat)) - for compartment in self.compartments: - for pos, rot, ingr, ptInd in compartment.molecules: - if ingr.name not in ingredients: - ingredients[ingr.name] = [ingr, [], [], []] - mat = rot.copy() - mat[:3, 3] = pos - ingredients[ingr.name][1].append(pos) - ingredients[ingr.name][2].append(rot) - ingredients[ingr.name][3].append(numpy.array(mat)) - all_ingr_as_array.append([pos, rot, ingr, ptInd]) - self.ingr_result = ingredients - return all_ingr_as_array + return self.packed_objects.get_ingredients() def check_new_placement(self, new_position): distances = self.get_all_distances(new_position) @@ -1876,44 +1822,33 @@ def distance_check_failed(self): expected_min_distance = self.smallestProteinSize * 2 return min_distance < expected_min_distance + 0.001 - @staticmethod - def get_count_from_options(count_options): - """ - Returns a count from the options - """ - count = None - if count_options.get("distribution") == "uniform": - count = int( - numpy.random.randint( - count_options.get("min", 0), count_options.get("max", 1) - ) - ) - elif count_options.get("distribution") == "normal": - count = int( - numpy.rint( - numpy.random.normal( - count_options.get("mean", 0), count_options.get("std", 1) - ) - ) - ) - elif count_options.get("distribution") == "list": - count = int( - numpy.rint(numpy.random.choice(count_options.get("list_values", None))) - ) - - return count - - def update_count(self, allIngredients): + def update_variable_ingredient_attributes(self, allIngredients): """ - updates the count for all ingredients based on input options + updates variable attributes for all ingredients based on input options """ for ingr in allIngredients: if hasattr(ingr, "count_options") and ingr.count_options is not None: - count = self.get_count_from_options(count_options=ingr.count_options) + + count = get_value_from_distribution( + distribution_options=ingr.count_options + ) if count is not None: ingr.count = count ingr.left_to_place = count + if hasattr(ingr, "size_options") and ingr.size_options is not None: + max_radius = get_max_value_from_distribution( + distribution_options=ingr.size_options + ) + if max_radius is not None: + ingr.max_radius = max_radius + + min_radius = get_min_value_from_distribution( + distribution_options=ingr.size_options + ) + if min_radius is not None: + ingr.min_radius = min_radius + def add_seed_number_to_base_name(self, seed_number): return f"{self.base_name}_seed_{seed_number}" @@ -1923,6 +1858,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, @@ -1949,7 +1889,7 @@ def pack_grid( allIngredients = self.callFunction(self.getActiveIng) # set the number of ingredients to pack - self.update_count(allIngredients) + self.update_variable_ingredient_attributes(allIngredients) # verify partner usePP = False @@ -1977,6 +1917,10 @@ def pack_grid( bb_outside = numpy.nonzero(self.freePointMask) self.grid.compartment_ids[bb_outside] = 99999 compartment_ids = self.grid.compartment_ids + + for compartment in self.compartments: + compartment.store_packed_object(self) + # why a copy? --> can we split ? distances = self.grid.distToClosestSurf[:] spacing = self.spacing or self.smallestProteinSize @@ -2110,20 +2054,20 @@ def pack_grid( ) if self.afviewer.renderDistance: self.afviewer.vi.displayParticleVolumeDistance(distances, self) - current_ingr_compartment = ingr.compNum + current_ingr_compartment = ingr.compartment_id # compute dpad which is the distance at which we need to update # distances after the drop is successfull max_radius = self.get_dpad(current_ingr_compartment) self.log.info( - f"picked Ingr radius {ingr.min_radius}, compNum {current_ingr_compartment}" + f"picked Ingr radius {ingr.min_radius}, compartment_id {current_ingr_compartment}" ) # find the points that can be used for this ingredient ## - if ingr.compNum > 0: - compartment = self.compartments[ingr.compNum - 1] + if ingr.compartment_id > 0: + compartment = self.compartments[ingr.compartment_id - 1] surface_points = compartment.surfacePoints res = [True, surface_points[int(random() * len(surface_points))]] else: @@ -2260,20 +2204,17 @@ def pack_grid( self.activeIngr = self.activeIngr0 + self.activeIngr12 if dump and ((time() - stime) > dump_freq): - print("SAVING", self.result_file) all_ingr_as_array = self.prep_molecules_for_save( distances, free_points, nbFreePoints, ) stime = time() - print(f"placed {len(self.molecules)}") + self.log.info(f"placed {len(self.packed_objects.get_ingredients())}") if self.saveResult: self.save_result( free_points, distances=distances, - t0=stime, - vTestid=vTestid, all_ingr_as_array=all_ingr_as_array, ) @@ -2289,33 +2230,37 @@ def pack_grid( self.save_result( free_points, distances=distances, - t0=time(), - vTestid=vTestid, all_ingr_as_array=all_ingr_as_array, ) def restore_molecules_array(self, ingr): - if len(ingr.results): - for elem in ingr.results: - if ingr.compNum == 0: - self.molecules.append([elem[0], numpy.array(elem[1]), ingr, 0]) - else: - ingr.recipe.compartment.molecules.append( - [elem[0], numpy.array(elem[1]), ingr, 0] - ) + pass + # if len(ingr.results): + # for elem in ingr.results: + # TODO: fix this to reset ingredients from results + # if ingr.compartment_id == 0: + # self.molecules.append( + # [elem[0], numpy.array(elem[1]), ingr, 0, ingr.radius] + # ) + # else: + # ingr.recipe.compartment.molecules.append( + # [elem[0], numpy.array(elem[1]), ingr, 0, ingr.radius] + # ) def restore(self, result, orgaresult, freePoint, tree=False): # should we used the grid ? the freePoint can be computed - # result is [pos,rot,ingr.name,ingr.compNum,ptInd] - # orgaresult is [[pos,rot,ingr.name,ingr.compNum,ptInd],[pos,rot,ingr.name,ingr.compNum,ptInd]...] + # result is [pos,rot,ingr.name,ingr.compartment_id,ptInd] + # orgaresult is [[pos,rot,ingr.name,ingr.compartment_id,ptInd],[pos,rot,ingr.name,ingr.compartment_id,ptInd]...] # after restore we can build the grid and fill! # ingredient based dictionary + # TODO: refactor with new packed_objects + ingredients = {} molecules = [] for elem in result: - pos, rot, name, compNum, ptInd = elem + pos, rot, name, compartment_id, ptInd = elem # needto check the name if it got the comp rule - ingr = self.getIngrFromName(name, compNum) + ingr = self.get_ingredient_by_name(name, compartment_id) if ingr is not None: molecules.append([pos, numpy.array(rot), ingr, ptInd]) if name not in ingredients: @@ -2325,19 +2270,15 @@ def restore(self, result, orgaresult, freePoint, tree=False): ingredients[name][1].append(pos) ingredients[name][2].append(numpy.array(rot)) ingredients[name][3].append(numpy.array(mat)) - self.rTrans.append(numpy.array(pos).flatten()) - self.rRot.append(numpy.array(rot)) # rotMatj - self.rIngr.append(ingr) ingr.results.append([pos, rot]) - self.molecules = molecules if self.exteriorRecipe: self.exteriorRecipe.molecules = molecules if len(orgaresult) == len(self.compartments): for i, o in enumerate(self.compartments): molecules = [] for elem in orgaresult[i]: - pos, rot, name, compNum, ptInd = elem - ingr = self.getIngrFromName(name, compNum) + pos, rot, name, compartment_id, ptInd = elem + ingr = self.get_ingredient_by_name(name, compartment_id) if ingr is not None: molecules.append([pos, numpy.array(rot), ingr, ptInd]) if name not in ingredients: @@ -2347,14 +2288,13 @@ def restore(self, result, orgaresult, freePoint, tree=False): ingredients[name][1].append(pos) ingredients[name][2].append(numpy.array(rot)) ingredients[name][3].append(numpy.array(mat)) - self.rTrans.append(numpy.array(pos).flatten()) - self.rRot.append(numpy.array(rot)) # rotMatj - self.rIngr.append(ingr) ingr.results.append([pos, rot]) o.molecules = molecules # consider that one filling have occured - if len(self.rTrans) and tree: - self.close_ingr_bhtree = spatial.cKDTree(self.rTrans, leafsize=10) + if len(self.packed_objects.get_ingredients()) and tree: + self.close_ingr_bhtree = spatial.cKDTree( + self.packed_objects.get_positions(), leafsize=10 + ) self.cFill = self.nFill self.ingr_result = ingredients if len(freePoint): @@ -2377,25 +2317,10 @@ def store(self, resultfilename=None): if resultfilename is None: resultfilename = self.result_file resultfilename = autopack.fixOnePath(resultfilename) - rfile = open(resultfilename, "wb") - # pickle.dump(self.molecules, rfile) - # OR - result = [] - for pos, rot, ingr, ptInd in self.molecules: - result.append([pos, rot, ingr.name, ingr.compNum, ptInd]) - pickle.dump(result, rfile) - rfile.close() - for i, orga in enumerate(self.compartments): - orfile = open(resultfilename + "_organelle_" + str(i), "wb") - result = [] - for pos, rot, ingr, ptInd in orga.molecules: - result.append([pos, rot, ingr.name, ingr.compNum, ptInd]) - pickle.dump(result, orfile) - # pickle.dump(orga.molecules, orfile) - orfile.close() - rfile = open(resultfilename + "_free_points", "wb") - pickle.dump(self.grid.free_points, rfile) - rfile.close() + with open(resultfilename, "wb") as rfile: + pickle.dump(self.packed_objects.get_ingredients(), rfile) + with open(resultfilename + "_free_points", "wb") as rfile: + pickle.dump(self.grid.free_points, rfile) @classmethod def dropOneIngr(self, pos, rot, ingrname, ingrcompNum, ptInd, rad=1.0): @@ -2437,10 +2362,10 @@ def getOneIngrJson(self, ingr, ingrdic): return ( ingrdic["results"], ingr.composition_name, - ingr.compNum, + ingr.compartment_id, 1, ingr.encapsulating_radius, - ) # ingrdic["compNum"],1,ingrdic["encapsulating_radius"] + ) # ingrdic["compartment_id"],1,ingrdic["encapsulating_radius"] def load_asTxt(self, resultfilename=None): if resultfilename is None: @@ -2493,19 +2418,12 @@ def cb(ingr): ingr.results = [] self.loopThroughIngr(cb) - for pos, rot, ingr, ptInd in self.molecules: + for obj in self.packed_objects.get_ingredients(): + ingr = obj.ingredient if isinstance(ingr, GrowIngredient) or isinstance(ingr, ActinIngredient): pass # already store else: - ingr.results.append([pos, rot]) - for i, orga in enumerate(self.compartments): - for pos, rot, ingr, ptInd in orga.molecules: - if isinstance(ingr, GrowIngredient) or isinstance( - ingr, ActinIngredient - ): - pass # already store - else: - ingr.results.append([pos, rot]) + ingr.results.append([obj.position, obj.rotation]) def load_asJson(self, resultfilename=None): if resultfilename is None: @@ -2645,7 +2563,7 @@ def load_asJson(self, resultfilename=None): def dropOneIngrJson(self, ingr, rdic): adic = OrderedDict() # [ingr.name] - adic["compNum"] = ingr.compNum + adic["compartment_id"] = ingr.compartment_id adic["encapsulating_radius"] = float(ingr.encapsulating_radius) adic["results"] = [] for r in ingr.results: @@ -2717,43 +2635,6 @@ def store_asJson(self, resultfilename=None, indent=True): ) # ,indent=4, separators=(',', ': ') print("ok dump", resultfilename) - def store_asTxt(self, resultfilename=None): - if resultfilename is None: - resultfilename = self.result_file - resultfilename = autopack.fixOnePath(resultfilename) - rfile = open(resultfilename + ".txt", "w") # doesnt work with symbol link ? - # pickle.dump(self.molecules, rfile) - # OR - line = "" - line += "\n" - for pos, rot, ingr, ptInd in self.molecules: - line += self.dropOneIngr( - pos, rot, ingr.name, ingr.compNum, ptInd, rad=ingr.encapsulating_radius - ) - # result.append([pos,rot,ingr.name,ingr.compNum,ptInd]) - rfile.write(line) - # write the curve point - - rfile.close() - for i, orga in enumerate(self.compartments): - orfile = open(resultfilename + "_organelle_" + str(i) + ".txt", "w") - line = "" - for pos, rot, ingr, ptInd in orga.molecules: - line += self.dropOneIngr( - pos, - rot, - ingr.name, - ingr.compNum, - ptInd, - rad=ingr.encapsulating_radius, - ) - orfile.write(line) - # pickle.dump(orga.molecules, orfile) - orfile.close() - # rfile = open(resultfilename+"free_points", 'w') - # pickle.dump(self.free_points, rfile) - # rfile.close() - @classmethod def convertPickleToText(self, resultfilename=None, norga=0): if resultfilename is None: @@ -2770,19 +2651,18 @@ def convertPickleToText(self, resultfilename=None, norga=0): rfile.close() rfile = open(resultfilename + ".txt", "w") line = "" - for pos, rot, ingrName, compNum, ptInd in result: - line += self.dropOneIngr(pos, rot, ingrName, compNum, ptInd) - # result.append([pos,rot,ingr.name,ingr.compNum,ptInd]) + for pos, rot, ingrName, compartment_id, ptInd in result: + line += self.dropOneIngr(pos, rot, ingrName, compartment_id, ptInd) + # result.append([pos,rot,ingr.name,ingr.compartment_id,ptInd]) rfile.write(line) rfile.close() for i in range(norga): orfile = open(resultfilename + "_organelle_" + str(i) + ".txt", "w") result = [] line = "" - for pos, rot, ingrName, compNum, ptInd in orgaresult[i]: - line += self.dropOneIngr(pos, rot, ingrName, compNum, ptInd) + for pos, rot, ingrName, compartment_id, ptInd in orgaresult[i]: + line += self.dropOneIngr(pos, rot, ingrName, compartment_id, ptInd) orfile.write(line) - # pickle.dump(orga.molecules, orfile) orfile.close() # freepoint @@ -3226,17 +3106,6 @@ def exportToReaDDy(self): # Animate # ============================================================================== - def readTraj(self, filename): - self.collectResultPerIngredient() - lenIngrInstance = len(self.molecules) - for orga in self.compartments: - lenIngrInstance += len(orga.molecules) - fileName, fileExtension = os.path.splitext(filename) - if fileExtension == ".dcd": - self.traj = dcdTrajectory(filename, lenIngrInstance) - elif fileExtension == ".molb": - self.traj = molbTrajectory(filename, lenIngrInstance) - self.traj.completeMapping(self) def linkTraj(self): # link the traj usin upy for creating a new synchronized calleback? @@ -3281,45 +3150,28 @@ def create_voxelization(self, image_data, image_size, voxel_size, hollow=False): The updated image data. """ channel_colors = [] - for pos, rot, ingr, _ in self.molecules: - if ingr.name not in image_data: - image_data[ingr.name] = numpy.zeros(image_size, dtype=numpy.uint8) - if ingr.color is not None: - color = ingr.color - if all([x <= 1 for x in ingr.color]): - color = [int(col * 255) for col in ingr.color] - channel_colors.append(color) - - image_data[ingr.name] = ingr.create_voxelization( - image_data=image_data[ingr.name], - bounding_box=self.boundingBox, - voxel_size=voxel_size, - image_size=image_size, - position=pos, - rotation=rot, - ) - for compartment in self.compartments: - if compartment.name not in image_data: - image_data[compartment.name] = numpy.zeros( - image_size, dtype=numpy.uint8 - ) - if hasattr(compartment, "color") and compartment.color is not None: - color = compartment.color - if all([x <= 1 for x in compartment.color]): - color = [int(col * 255) for col in compartment.color] + for obj in self.packed_objects.get_all(): + mesh_store = None + if obj.name not in image_data: + image_data[obj.name] = numpy.zeros(image_size, dtype=numpy.uint8) + if obj.color is not None: + color = obj.color + if all([x <= 1 for x in obj.color]): + color = [int(col * 255) for col in obj.color] channel_colors.append(color) - else: - channel_colors.append([0, 255, 0]) + if obj.is_compartment: + mesh_store = self.mesh_store - image_data[compartment.name] = compartment.create_voxelization( - image_data=image_data[compartment.name], + image_data[obj.name] = obj.ingredient.create_voxelization( + image_data=image_data[obj.name], bounding_box=self.boundingBox, voxel_size=voxel_size, image_size=image_size, - position=compartment.position, - mesh_store=self.mesh_store, + position=obj.position, + rotation=obj.rotation, hollow=hollow, + mesh_store=mesh_store, ) return image_data, channel_colors diff --git a/cellpack/autopack/Graphics.py b/cellpack/autopack/Graphics.py index 25e317b0e..abec13189 100644 --- a/cellpack/autopack/Graphics.py +++ b/cellpack/autopack/Graphics.py @@ -241,10 +241,10 @@ def addMasterIngr(self, ingr, parent=None): @type parent: hostObject @param parent: specify a parent to insert under """ - if ingr.compNum == 0: + if ingr.compartment_id == 0: compartment = self.histo else: - compartment = self.histo.compartments[abs(ingr.compNum) - 1] + compartment = self.histo.compartments[abs(ingr.compartment_id) - 1] name = "%s_%s" % (ingr.name, compartment.name) gi = self.checkCreateEmpty(name, parent=parent) self.orgaToMasterGeom[ingr] = gi @@ -1205,14 +1205,14 @@ def displayIngrResults(self, ingredient, doSphere=False, doMesh=True): print("no ingredient provided") return o = ingredient.recipe.compartment - # comp = ingredient.compNum + # comp = ingredient.compartment_id verts = [] radii = [] matrices = [] # retrieve the result from the molecules array of the recipes res = [ self.collectResult(ingr, pos, rot) - for pos, rot, ingr, ptInd in o.molecules + for pos, rot, ingr, ptInd, _ in o.molecules if ingr == ingredient ] @@ -1302,7 +1302,7 @@ def displayCytoplasmIngredients(self): else: return if self.doSpheres: - for pos, rot, ingr, ptInd in self.histo.molecules: + for pos, rot, ingr, ptInd, _ in self.histo.molecules: level = ingr.deepest_level px = ingr.transformPoints(pos, rot, ingr.positions[level]) if ingr.model_type == "Spheres": @@ -1333,7 +1333,7 @@ def displayCytoplasmIngredients(self): if r: meshGeoms = {} # self.meshGeoms#{} inds = {} - for pos, rot, ingr, ptInd in self.histo.molecules: + for pos, rot, ingr, ptInd, _ in self.histo.molecules: if ingr.mesh: # display mesh mat = rot.copy() mat[:3, 3] = pos @@ -1416,7 +1416,7 @@ def displayCompartmentsIngredients(self): for ingr in ri.ingredients: verts[ingr] = [] radii[ingr] = [] - for pos, rot, ingr, ptInd in orga.molecules: + for pos, rot, ingr, ptInd, _ in orga.molecules: level = ingr.deepest_level px = ingr.transformPoints(pos, rot, ingr.positions[level]) if ingr.model_type == "Spheres": @@ -1460,7 +1460,7 @@ def displayCompartmentsIngredients(self): # for ingr in ri.ingredients: # if ingr.mesh: # display mesh # self.displayIngrMesh(matrices,ingr) - for pos, rot, ingr, ptInd in orga.molecules: + for pos, rot, ingr, ptInd, _ in orga.molecules: if ingr.mesh: # display mesh mat = rot.copy() mat[:3, 3] = pos @@ -1599,13 +1599,13 @@ def appendIngrInstance(self, ingr, sel=None, bb=None): rotMatj = self.helper.getMatRotation(ch) jtrans = self.helper.ToVec(self.helper.getTranslation(ch)) ptId = self.histo.getPointFrom3D(jtrans, bb, spacing, nb) - res.append([jtrans, rotMatj, ingr, ptId]) + res.append([jtrans, rotMatj, ingr, ptId, ingr.radius]) # need to handle the PtID none else: rotMatj = self.helper.getMatRotation(o) jtrans = self.helper.ToVec(self.helper.getTranslation(o)) ptId = self.histo.getPointFrom3D(jtrans, bb, spacing, nb) - res.append([jtrans, rotMatj, ingr, ptId]) + res.append([jtrans, rotMatj, ingr, ptId, ingr.radius]) ingr.counter = 0 orga.molecules.extend(res) @@ -1986,13 +1986,13 @@ def addIngredientFromGeom(self, name, ingrobj, recipe=None, **kw): if ingr is not None: if recipe is None: self.histo.exteriorRecipe.addIngredient(ingr) - ingr.compNum = 0 + ingr.compartment_id = 0 g = self.vi.getObject(self.name + "_cytoplasm") self.addMasterIngr(ingr, parent=g) ingr.env = self.histo else: recipe.addIngredient(ingr) - ingr.compNum = recipe.number + ingr.compartment_id = recipe.number # g = self.vi.getObject("O" + o.name) ingr.env = self.histo rep = self.vi.getObject(ingr.composition_name + "_mesh") @@ -2002,7 +2002,7 @@ def addIngredientFromGeom(self, name, ingrobj, recipe=None, **kw): ingr.meshFile = autopack.cache_geoms + os.sep + ingr.composition_name ingr.meshName = ingr.composition_name + "_mesh" ingr.saveDejaVuMesh(ingr.meshFile) - self.addMasterIngr(ingr, parent=self.orgaToMasterGeom[ingr.compNum]) + self.addMasterIngr(ingr, parent=self.orgaToMasterGeom[ingr.compartment_id]) return ingr def addCompartmentFromGeom(self, name, obj, **kw): @@ -2301,7 +2301,7 @@ def colorByOrder( else: # for each mol get the order and color the poly inds = {} - for pos, rot, ingr, ptInd in self.histo.molecules: + for pos, rot, ingr, ptInd, _ in self.histo.molecules: if ingr not in inds: inds[ingr] = [ptInd] else: diff --git a/cellpack/autopack/ingredient/Ingredient.py b/cellpack/autopack/ingredient/Ingredient.py index f454b8b15..e37d420f8 100644 --- a/cellpack/autopack/ingredient/Ingredient.py +++ b/cellpack/autopack/ingredient/Ingredient.py @@ -54,6 +54,10 @@ from time import time import math +from cellpack.autopack.interface_objects.ingredient_types import INGREDIENT_TYPE +from cellpack.autopack.interface_objects.packed_objects import PackedObject +from cellpack.autopack.utils import get_distance, get_value_from_distribution + from .utils import ( ApplyMatrix, getNormedVectorOnes, @@ -72,15 +76,15 @@ reporthook = helper.reporthook -class CountDistributions(MetaEnum): - "All available count distributions" +class DistributionTypes(MetaEnum): + # All available distribution types UNIFORM = "uniform" NORMAL = "normal" LIST = "list" -class CountOptions(MetaEnum): - "All available count options" +class DistributionOptions(MetaEnum): + # All available distribution options MIN = "min" MAX = "max" MEAN = "mean" @@ -88,10 +92,10 @@ class CountOptions(MetaEnum): LIST_VALUES = "list_values" -REQUIRED_COUNT_OPTIONS = { - CountDistributions.UNIFORM: [CountOptions.MIN, CountOptions.MAX], - CountDistributions.NORMAL: [CountOptions.MEAN, CountOptions.STD], - CountDistributions.LIST: [CountOptions.LIST_VALUES], +REQUIRED_DISTRIBUTION_OPTIONS = { + DistributionTypes.UNIFORM: [DistributionOptions.MIN, DistributionOptions.MAX], + DistributionTypes.NORMAL: [DistributionOptions.MEAN, DistributionOptions.STD], + DistributionTypes.LIST: [DistributionOptions.LIST_VALUES], } @@ -137,7 +141,7 @@ class Ingredient(Agent): of the (+) values - an optional principal vector used to align the ingredient - recipe will be a weakref to the Recipe this Ingredient belongs to - - compNum is the compartment number (0 for cytoplasm, positive for compartment + - compartment_id is the compartment number (0 for cytoplasm, positive for compartment surface and negative compartment interior - Attributes used by the filling algorithm: - count counts the number of placed ingredients during a fill @@ -177,6 +181,7 @@ class Ingredient(Agent): "resolution_dictionary", "rotation_axis", "rotation_range", + "size_options", "type", "use_orient_bias", "use_rotation_axis", @@ -215,6 +220,7 @@ def __init__( resolution_dictionary=None, rotation_axis=None, rotation_range=6.2831, + size_options=None, use_orient_bias=False, use_rotation_axis=False, weight=0.2, @@ -239,6 +245,7 @@ def __init__( self.molarity = molarity self.count = count self.count_options = count_options + self.size_options = size_options self.priority = priority self.log.info( "priority %d, self.priority %r", @@ -290,7 +297,7 @@ def __init__( self.principal_vector = principal_vector self.recipe = None # will be set when added to a recipe - self.compNum = None + self.compartment_id = None self.compId_accepted = ( [] ) # if this list is defined, point picked outise the list are rejected @@ -360,6 +367,26 @@ def __init__( self.organism = "" # add tiling property ? as any ingredient coud tile as hexagon. It is just the packing type + @staticmethod + def validate_distribution_options(distribution_options): + """ + Validates distribution options and returns validated distribution options + """ + if "distribution" not in distribution_options: + raise Exception("Ingredient count options must contain a distribution") + if not DistributionTypes.is_member(distribution_options["distribution"]): + raise Exception( + f"{distribution_options['distribution']} is not a valid distribution" + ) + for required_option in REQUIRED_DISTRIBUTION_OPTIONS.get( + distribution_options["distribution"], [] + ): + if required_option not in distribution_options: + raise Exception( + f"Missing option '{required_option}' for {distribution_options['distribution']} distribution" + ) + return distribution_options + @staticmethod def validate_ingredient_info(ingredient_info): """ @@ -372,20 +399,14 @@ def validate_ingredient_info(ingredient_info): raise Exception("Ingredient count must be greater than or equal to 0") if "count_options" in ingredient_info: - count_options = ingredient_info["count_options"] - if "distribution" not in count_options: - raise Exception("Ingredient count options must contain a distribution") - if not CountDistributions.is_member(count_options["distribution"]): - raise Exception( - f"{count_options['distribution']} is not a valid count distribution" - ) - for required_option in REQUIRED_COUNT_OPTIONS.get( - count_options["distribution"], [] - ): - if required_option not in count_options: - raise Exception( - f"Missing option '{required_option}' for {count_options['distribution']} distribution" - ) + ingredient_info["count_options"] = Ingredient.validate_distribution_options( + ingredient_info["count_options"] + ) + + if "size_options" in ingredient_info: + ingredient_info["size_options"] = Ingredient.validate_distribution_options( + ingredient_info["size_options"] + ) return ingredient_info @@ -466,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 @@ -724,7 +732,7 @@ def jitterPosition(self, position, spacing, normal=None): spacing is the grid spacing this will jitter gauss(0., 0.3) * Ingredient.max_jitter """ - if self.compNum > 0: + if self.compartment_id > 0: vx, vy, vz = v1 = self.principal_vector # surfacePointsNormals problem here v2 = normal @@ -742,7 +750,7 @@ def jitterPosition(self, position, spacing, normal=None): dz = jz * spacing * uniform(-1.0, 1.0) # d2 = dx*dx + dy*dy + dz*dz # if d2 < jitter2: - if self.compNum > 0: # jitter less among normal + if self.compartment_id > 0: # jitter less among normal dx, dy, dz, dum = numpy.dot(rotMat, (dx, dy, dz, 0)) position[0] += dx position[1] += dy @@ -814,7 +822,7 @@ def get_list_of_free_indices( ): allIngrPts = [] allIngrDist = [] - current_comp_id = self.compNum + current_comp_id = self.compartment_id # gets min distance an object has to be away to allow packing for this object cuttoff = self.get_cuttoff_value(spacing) if self.packing_mode == "close": @@ -909,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") @@ -966,36 +974,9 @@ 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.compNum == 0: - compartment = self.env - else: - compartment = self.env.compartments[abs(self.compNum) - 1] - compNum = self.compNum - if compNum < 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 compNum == 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.compNum + current = self.compartment_id if current < 0: # inside ins = [i for i, x in enumerate(cId) if x == current] # surf=[i for i,x in enumerate(cId) if x == -current] @@ -1009,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.compNum - 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.compNum <= 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, @@ -1092,7 +1031,7 @@ def is_point_in_correct_region(self, point): nearest_grid_point_compartment_id = ( self.env.compartment_id_for_nearest_grid_point(point) ) # offset ? - compartment_ingr_belongs_in = self.compNum + compartment_ingr_belongs_in = self.compartment_id if compartment_ingr_belongs_in == 0: compartment = self.env else: @@ -1101,7 +1040,7 @@ def is_point_in_correct_region(self, point): compartment = self.env.compartments[abs(compartment_ingr_belongs_in) - 1] if compartment_ingr_belongs_in > 0: # surface ingredient if self.type == "Grow": - # need a list of accepted compNum + # need a list of accepted compartment_id check = False if len(self.compMask): check = nearest_grid_point_compartment_id in self.compMask @@ -1130,7 +1069,7 @@ def is_point_in_correct_region(self, point): def far_enough_from_surfaces(self, point, cutoff): # check if clear of all other compartment surfaces ingredient_compartment = self.get_compartment(self.env) - ingredient_compartment_id = self.compNum + ingredient_compartment_id = self.compartment_id for compartment in self.env.compartments: if ( ingredient_compartment_id > 0 @@ -1184,37 +1123,21 @@ 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 - ) - # return self.get_new_distance_values( - # grid, gridPointsCoords, dpad, distance, centT, jtrans, rotMatj, dpad - # ) - - def getIngredientsInBox(self, histoVol, jtrans, rotMat, compartment, afvi): - if histoVol.windowsSize_overwrite: - rad = histoVol.windowsSize + def getIngredientsInBox(self, env, jtrans, rotMat, compartment): + if env.windowsSize_overwrite: + radius = env.windowsSize else: - # rad = self.min_radius*2.0# + histoVol.largestProteinSize + \ - # histoVol.smallestProteinSize + histoVol.windowsSize - rad = ( + radius = ( self.min_radius - + histoVol.largestProteinSize - + histoVol.smallestProteinSize - + histoVol.windowsSize + + env.largestProteinSize + + env.smallestProteinSize + + env.windowsSize ) x, y, z = jtrans - bb = ([x - rad, y - rad, z - rad], [x + rad, y + rad, z + rad]) + bb = ( + [x - radius, y - radius, z - radius], + [x + radius, y + radius, z + radius], + ) if self.model_type == "Cylinders": cent1T = self.transformPoints( jtrans, rotMat, self.positions[self.deepest_level] @@ -1240,7 +1163,7 @@ def getIngredientsInBox(self, histoVol, jtrans, rotMat, compartment, afvi): if bb[0][i] > maxBB[i]: maxBB[i] = bb[0][i] bb = [minBB, maxBB] - if histoVol.runTimeDisplay > 1: + if env.runTimeDisplay > 1: box = self.vi.getObject("partBox") if box is None: box = self.vi.Box("partBox", cornerPoints=bb, visible=1) @@ -1249,45 +1172,33 @@ def getIngredientsInBox(self, histoVol, jtrans, rotMat, compartment, afvi): self.vi.updateBox(box, cornerPoints=bb) self.vi.update() # sleep(1.0) - pointsInCube = histoVol.grid.getPointsInCube(bb, jtrans, rad) + # pointsInCube = env.grid.getPointsInCube(bb, jtrans, radius) # should we got all ingre from all recipes? # can use the kdtree for it... # maybe just add the surface if its not already the surface - mingrs = [m for m in compartment.molecules if m[3] in pointsInCube] - return mingrs - - def getIngredientsInTree(self, closest_ingredients): - if len(self.env.rIngr): - ingrs = [self.env.rIngr[i] for i in closest_ingredients["indices"]] - return [ - numpy.array(self.env.rTrans)[closest_ingredients["indices"]], - numpy.array(self.env.rRot)[closest_ingredients["indices"]], - ingrs, - closest_ingredients["distances"], - ] - else: - return [] + ingredients = [] + for obj in compartment.packed_objects.get(): + ingredients.append([obj, get_distance(obj.position, jtrans)]) - def get_partners(self, jtrans, rotMat, organelle, afvi): - env = self.env - closest_ingredients = self.get_closest_ingredients( - jtrans, cutoff=self.env.grid.diag - ) + return ingredients + + def get_partners(self, env, jtrans, rotMat, organelle): + 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, afvi + 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) or not len(near_by_ingredients[2]): + if not len(near_by_ingredients): self.log.info("no close ingredient found") return [], [] else: self.log.info("nb close ingredient %s", self.name) - for i in range(len(near_by_ingredients[2])): - packed_ingredient = near_by_ingredients[2][i] - distance = (near_by_ingredients[3][i],) + for i in range(len(near_by_ingredients)): + packed_ingredient = near_by_ingredients[i][0].ingredient + distance = near_by_ingredients[i][1] if self.packing_mode == "closePartner": if self.partners.is_partner(packed_ingredient.name): placed_partners.append( @@ -1319,16 +1230,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 @@ -1339,12 +1240,13 @@ def get_new_pos(self, ingr, pos, rot, positions_to_adjust): return self.transformPoints(pos, rot, positions_to_adjust) def check_against_one_packed_ingr(self, index, level, search_tree): - overlapped_ingr = self.env.rIngr[index] + ingredient_instance = self.env.packed_objects.get_ingredients()[index] + ingredient_class = ingredient_instance.ingredient positions_of_packed_ingr_spheres = self.get_new_pos( - self.env.rIngr[index], - self.env.rTrans[index], - self.env.rRot[index], - overlapped_ingr.positions[level], + ingredient_class, + ingredient_instance.position, + ingredient_instance.rotation, + ingredient_class.positions[level], ) # check distances between the spheres at this level in the ingr we are packing # to the spheres at this level in the ingr already placed @@ -1354,7 +1256,9 @@ def check_against_one_packed_ingr(self, index, level, search_tree): ) # return index of sph1 closest to pos of packed ingr cradii = numpy.array(self.radii[level])[ind] - oradii = numpy.array(self.env.rIngr[index].radii[level]) + oradii = numpy.array( + self.env.packed_objects.get_ingredients()[index].ingredient.radii[level] + ) sumradii = numpy.add(cradii.transpose(), oradii).transpose() sD = dist_from_packed_spheres_to_new_spheres - sumradii return len(numpy.nonzero(sD < 0.0)[0]) != 0 @@ -1362,12 +1266,13 @@ def check_against_one_packed_ingr(self, index, level, search_tree): def np_check_collision(self, packing_location, rotation): has_collision = False # no ingredients packed yet - if not len(self.env.rTrans): + packed_objects = self.env.packed_objects.get_ingredients() + if not len(packed_objects): return has_collision else: if self.env.close_ingr_bhtree is None: self.env.close_ingr_bhtree = spatial.cKDTree( - self.env.rTrans, leafsize=10 + self.env.packed_objects.get_positions(), leafsize=10 ) # starting at level 0, check encapsulating radii level = 0 @@ -1375,9 +1280,9 @@ def np_check_collision(self, packing_location, rotation): ( distances_from_packing_location_to_all_ingr, ingr_indexes, - ) = self.env.close_ingr_bhtree.query(packing_location, len(self.env.rTrans)) + ) = self.env.close_ingr_bhtree.query(packing_location, len(packed_objects)) radii_of_placed_ingr = numpy.array( - [ing.encapsulating_radius for ing in self.env.rIngr] + self.env.packed_objects.get_encapsulating_radii() )[ingr_indexes] overlap_distance = distances_from_packing_location_to_all_ingr - ( self.encapsulating_radius + radii_of_placed_ingr @@ -1427,10 +1332,10 @@ def get_rbNodes( ): # move around the rbnode and return it # self.env.loopThroughIngr( self.env.reset_rbnode ) - if self.compNum == 0: + if self.compartment_id == 0: organelle = self.env else: - organelle = self.env.compartments[abs(self.compNum) - 1] + organelle = self.env.compartments[abs(self.compartment_id) - 1] nodes = [] # a=numpy.asarray(self.env.rTrans)[close_indice["indices"]] # b=numpy.array([currentpt,]) @@ -1488,7 +1393,7 @@ def get_rbNodes( nodes.append(rbnode) # append organelle rb nodes for o in self.env.compartments: - if self.compNum > 0 and o.name == organelle.name: + if self.compartment_id > 0 and o.name == organelle.name: # this i notworking for growing ingredient like hair. # should had after second segments if self.type != "Grow": @@ -1508,7 +1413,7 @@ def get_rbNodes( nodes.append(orbnode) else: nodes.append([orbnode, [0, 0, 0], numpy.identity(4), o]) - # if self.compNum < 0 or self.compNum == 0 : + # if self.compartment_id < 0 or self.compartment_id == 0 : # for o in self.env.compartments: # if o.rbnode is not None : # if not getInfo : @@ -1516,85 +1421,11 @@ 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 - if not len(env.rTrans): - return to_return - if env.close_ingr_bhtree is not None: - # request kdtree - nb = [] - self.log.info("finding partners") - if len(env.rTrans) >= 1: - # nb = histoVol.close_ingr_bhtree.query_ball_point(point,cutoff) - # else :#use the general query, how many we want - distance, nb = env.close_ingr_bhtree.query( - point, len(env.rTrans), distance_upper_bound=cutoff - ) # len of ingr posed so far - if len(env.rTrans) == 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, jtrans, rotMatj, ptInd=0, pt1=None, pt2=None, updateTree=True - ): - # self.env.static.append(rbnode) - # self.env.moving = None - self.env.nb_ingredient += 1 - self.env.rTrans.append(numpy.array(jtrans).flatten().tolist()) - self.env.rRot.append(numpy.array(rotMatj)) # rotMatj - self.env.rIngr.append(self) - if pt1 is not None: - self.env.result.append( - [ - [ - numpy.array(pt1).flatten().tolist(), - numpy.array(pt2).flatten().tolist(), - ], - rotMatj.tolist(), - self, - ptInd, - ] - ) - else: - self.env.result.append( - [ - numpy.array(jtrans).flatten().tolist(), - numpy.array(rotMatj), - self, - ptInd, - ] + def update_data_tree(self): + if len(self.env.packed_objects.get_ingredients()) >= 1: + self.env.close_ingr_bhtree = spatial.cKDTree( + self.env.packed_objects.get_positions(), leafsize=10 ) - if updateTree: - if len(self.env.rTrans) >= 1: - self.env.close_ingr_bhtree = spatial.cKDTree( - self.env.rTrans, leafsize=10 - ) def pack_at_grid_pt_location( self, @@ -1605,9 +1436,12 @@ def pack_at_grid_pt_location( grid_point_distances, inside_points, new_dist_points, + pt_index, ): + packing_location = jtrans radius_of_area_to_check = self.encapsulating_radius + dpad + self.store_packed_object(packing_location, rotation_matrix, pt_index) bounding_points_to_check = self.get_all_positions_to_check(packing_location) @@ -1647,24 +1481,30 @@ def reject( self.log.info("PREMATURE ENDING of ingredient %s", self.name) self.completion = 1.0 + def store_packed_object(self, position, rotation, index): + packed_object = PackedObject( + position=position, + rotation=rotation, + radius=self.get_radius(), + pt_index=index, + ingredient=self, + ) + self.env.packed_objects.add(packed_object) + if self.compartment_id != 0: + compartment = self.get_compartment(self.env) + compartment.packed_objects.add(packed_object) + def place( self, env, - compartment, dropped_position, dropped_rotation, grid_point_index, new_inside_points, - new_dist_values, ): self.nbPts = self.nbPts + len(new_inside_points) - # self.update_distances(new_inside_points, new_dist_values) - compartment.molecules.append( - [dropped_position, dropped_rotation, self, grid_point_index] - ) - 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( @@ -1678,7 +1518,18 @@ def place( self.counter += 1 self.completion = float(self.counter) / float(self.left_to_place) self.rejectionCounter = 0 - self.update_data_tree(dropped_position, dropped_rotation, grid_point_index) + self.update_data_tree() + + def update_ingredient_size(self): + # update the size of the ingredient based on input options + if hasattr(self, "size_options") and self.size_options is not None: + if self.type == INGREDIENT_TYPE.SINGLE_SPHERE: + radius = get_value_from_distribution( + distribution_options=self.size_options + ) + if radius is not None: + self.radius = radius + self.encapsulating_radius = radius def attempt_to_pack_at_grid_location( self, @@ -1692,6 +1543,8 @@ def attempt_to_pack_at_grid_location( ): success = False jitter = self.getMaxJitter(spacing) + self.update_ingredient_size() + dpad = self.min_radius + max_radius + jitter self.vi = autopack.helper self.env = env # NOTE: do we need to store the env on the ingredient? @@ -1727,6 +1580,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, @@ -1757,13 +1611,12 @@ def attempt_to_pack_at_grid_location( newDistPoints, ) = self.jitter_place( env, - compartment, target_grid_point_position, rotation_matrix, current_visual_instance, grid_point_distances, dpad, - env.afviewer, + ptInd, ) elif self.place_method == "spheresSST": ( @@ -1843,15 +1696,14 @@ def attempt_to_pack_at_grid_location( grid_point_distances, insidePoints, newDistPoints, + ptInd, ) if success: if is_realtime: autopack.helper.set_object_static( current_visual_instance, jtrans, rotMatj ) - self.place( - env, compartment, jtrans, rotMatj, ptInd, insidePoints, newDistPoints - ) + self.place(env, jtrans, rotMatj, ptInd, insidePoints) else: if is_realtime: self.remove_from_realtime_display(current_visual_instance) @@ -1861,7 +1713,7 @@ def attempt_to_pack_at_grid_location( def get_rotation(self, pt_ind, env, compartment): # compute rotation matrix rotMat - comp_num = self.compNum + comp_num = self.compartment_id rot_mat = numpy.identity(4) if comp_num > 0: @@ -1884,9 +1736,11 @@ 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 @@ -1894,10 +1748,10 @@ 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.compNum > 0: + if self.compartment_id > 0: jitter_rotation = self.getAxisRotation(rotation) else: if self.use_rotation_axis: @@ -1907,15 +1761,14 @@ def randomize_rotation(self, rotation, histovol): # set use_rotation_axis = 1 and set rotation_axis = 0, 0, 0 for that ingredient elif self.use_orient_bias and self.packing_mode == "gradient": jitter_rotation = self.getBiasedRotation(rotation, weight=None) - # weight = 1.0 - self.env.gradients[self.gradient].weight[ptInd]) 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: @@ -1946,7 +1799,7 @@ def randomize_translation(self, env, translation, rotation): dz = jz * jitter * uniform(-1.0, 1.0) d2 = dx * dx + dy * dy + dz * dz if d2 < jitter_sq: - if self.compNum > 0: # jitter less among normal + if self.compartment_id > 0: # jitter less among normal dx, dy, dz, _ = numpy.dot(rotation, (dx, dy, dz, 0)) jitter_trans = (tx + dx, ty + dy, tz + dz) found = True @@ -1963,7 +1816,7 @@ def update_display_rt(self, current_instance, translation, rotation): def rigid_place( self, - histoVol, + env, ptInd, compartment, target_grid_point_position, @@ -1976,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 @@ -1999,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 ) for i, elem in enumerate(near_by_ingredients): ing = elem[2] @@ -2044,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: @@ -2063,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 : @@ -2097,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, @@ -2109,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) @@ -2168,13 +2021,12 @@ def get_all_positions_to_check(self, packing_location): def jitter_place( self, env, - compartment, targeted_master_grid_point, rot_mat, moving, distance, dpad, - afvi, + pt_index, ): """ Check if the given grid point is available for packing using the jitter collision detection @@ -2257,6 +2109,8 @@ def jitter_place( len(insidePoints), len(newDistPoints), ) + for pt in points_to_check: + self.store_packed_object(pt, packing_rotation, pt_index) return ( True, packing_location, @@ -2266,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): near_by_ingredients, placed_partners = self.get_partners( - trans, rotMat, organelle, afvi + env, trans, rotMat, organelle ) targetPoint = trans found = False @@ -2288,7 +2142,7 @@ def lookForNeighbours(self, trans, rotMat, organelle, afvi): return targetPoint, rotMat, found else: # maybe get the ptid that can have it found = True - if self.compNum > 0: + if self.compartment_id > 0: # surface d, i = organelle.OGsrfPtsBht.query(targetPoint) vx, vy, vz = v1 = self.principal_vector @@ -2312,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 @@ -2339,17 +2193,19 @@ def pandaBullet_collision(self, pos, rot, rbnode, getnodes=False): return collision def get_compartment(self, env): - if self.compNum == 0: + if self.compartment_id == 0: return env else: - return env.compartments[abs(self.compNum) - 1] + 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, - afvi, ) if not found and self.counter != 0: self.reject() @@ -2478,16 +2334,10 @@ 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 for pt in points_to_check: - self.env.rTrans.append(pt) - self.env.rRot.append(packing_rotation) - self.env.rIngr.append(self) - self.env.result.append([pt, packing_rotation, self, ptInd]) self.pack_at_grid_pt_location( env, @@ -2497,14 +2347,10 @@ def pandaBullet_place( distance, inside_points, new_dist_points, + ptInd, ) 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, @@ -2547,8 +2393,6 @@ def spheres_SST_place( env.setupPanda() is_realtime = moving is not None - gridPointsCoords = env.grid.masterGridPositions - targetPoint = target_grid_point_position if is_realtime: @@ -2624,22 +2468,19 @@ def spheres_SST_place( # ) # collision_results.extend([collision]) if True not in collision_results: - # self.update_data_tree(jtrans,rotMatj,ptInd=ptInd)? - self.env.static.append(rbnode) - self.env.moving = None + env.static.append(rbnode) + env.moving = None for pt in pts_to_check: - self.env.rTrans.append(pt) - self.env.rRot.append(packing_rotation) - self.env.rIngr.append(self) - self.env.result.append([pt, packing_rotation, self, ptInd]) - new_inside_pts, new_dist_points = self.get_new_distance_values( + new_inside_pts, new_dist_points = self.pack_at_grid_pt_location( + env, pt, packing_rotation, - gridPointsCoords, - distance, dpad, - self.deepest_level, + distance, + insidePoints, + newDistPoints, + ptInd, ) insidePoints = self.merge_place_results( new_inside_pts, insidePoints @@ -2647,12 +2488,6 @@ def spheres_SST_place( newDistPoints = self.merge_place_results( new_dist_points, newDistPoints ) - # rebuild kdtree - if len(self.env.rTrans) >= 1: - del self.env.close_ingr_bhtree - self.env.close_ingr_bhtree = spatial.cKDTree( - self.env.rTrans, leafsize=10 - ) success = True return ( diff --git a/cellpack/autopack/ingredient/agent.py b/cellpack/autopack/ingredient/agent.py index 152b94b87..8ae82dd2c 100644 --- a/cellpack/autopack/ingredient/agent.py +++ b/cellpack/autopack/ingredient/agent.py @@ -86,10 +86,8 @@ def pick_partner_grid_index( self, near_by_ingredients, placed_partners, current_packing_position=[0, 0, 0] ): # near_by_ingredient is [ - # ingredient pos[], - # ingredient rot[], - # ingredient[], - # distance[]? + # PackedObject + # distance[] # ] # placed_partners is (index,placed_partner_ingredient,distance_from_current_point) # weight using the distance function @@ -100,13 +98,13 @@ def pick_partner_grid_index( return None partner_index = placed_partners[i][0] # i,part,dist partner = placed_partners[i][1] - partner_ingredient = near_by_ingredients[2][partner_index] + partner_ingredient = near_by_ingredients[partner_index][0].ingredient self.log.info(f"binding to {partner_ingredient.name}") - if self.compNum > 0: + if self.compartment_id > 0: packing_position = self.env.grid.getClosestFreeGridPoint( packing_position, - compId=self.compNum, + compId=self.compartment_id, ball=( partner_ingredient.encapsulating_radius + self.encapsulating_radius ), @@ -120,7 +118,7 @@ def pick_partner_grid_index( if binding_probability > 0: bind = chance <= binding_probability if bind: - partner_position = near_by_ingredients[0][partner_index] + partner_position = near_by_ingredients[partner_index][0].position # get closestFreePoint using freePoint and masterGridPosition # if self.place_method == "rigid-body" or self.place_method == "jitter": diff --git a/cellpack/autopack/ingredient/grow.py b/cellpack/autopack/ingredient/grow.py index 0a8803646..e993220cc 100644 --- a/cellpack/autopack/ingredient/grow.py +++ b/cellpack/autopack/ingredient/grow.py @@ -160,7 +160,7 @@ def __init__( self.orientation = orientation self.seedOnPlus = True # The filament should continue to grow on its (+) end self.seedOnMinus = False # The filamen should continue to grow on its (-) end. - # if self.compNum > 0 : + # if self.compartment_id > 0 : # self.seedOnMinus = False self.vector = [0.0, 0.0, 0.0] self.biased = biased @@ -601,7 +601,7 @@ def walkLatticeSurface( step = histoVol.grid.gridSpacing * size bb = ([cx - step, cy - step, cz - step], [cx + step, cy + step, cz + step]) pointsInCube = histoVol.grid.getPointsInCube(bb, posc, step, addSP=False) - o = self.env.compartments[abs(self.compNum) - 1] + o = self.env.compartments[abs(self.compartment_id) - 1] sfpts = o.surfacePointsCoords found = False @@ -905,7 +905,9 @@ 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 +1140,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: @@ -1211,7 +1213,13 @@ def walkSpherePanda( # self.update_data_tree(jtrans1,rotMatj1,pt1=newPt,pt2=newPts[1]) self.partners[alternate].ingr.update_data_tree(jtrans, rotMatj) self.compartment.molecules.append( - [jtrans, rotMatj, self.partners[alternate].ingr, 0] + [ + jtrans, + rotMatj, + self.partners[alternate].ingr, + 0, + self.partners[alternate].ingr.radius, + ] ) # transpose ? newv, d1 = self.vi.measure_distance(pt2, newPt, vec=True) # v,d2 = self.vi.measure_distance(newPt,newPts[1],vec=True) @@ -1324,7 +1332,9 @@ 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) @@ -1361,8 +1371,8 @@ def grow( counter = 0 mask = None - if self.walkingMode == "lattice" and self.compNum > 0: - o = self.env.compartments[abs(self.compNum) - 1] + if self.walkingMode == "lattice" and self.compartment_id > 0: + o = self.env.compartments[abs(self.compartment_id) - 1] v = o.surfacePointsCoords mask = numpy.ones(len(v), int) alternate = False @@ -1413,7 +1423,7 @@ def grow( marge=self.marge, checkcollision=True, ) - if self.walkingMode == "lattice" and self.compNum > 0: + if self.walkingMode == "lattice" and self.compartment_id > 0: secondPoint, success, mask = self.walkLatticeSurface( startingPoint, distance, @@ -1521,14 +1531,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( @@ -1587,14 +1596,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) @@ -1605,10 +1613,12 @@ def updateGrid( return insidePoints, newDistPoints, nbFreePoints, free_points def getFirstPoint(self, ptInd, seed=0): - if self.compNum > 0: # surfacegrowing: first point is aling to the normal: - v2 = self.env.compartments[abs(self.compNum) - 1].surfacePointsNormals[ - ptInd - ] + if ( + self.compartment_id > 0 + ): # surfacegrowing: first point is aling to the normal: + v2 = self.env.compartments[ + abs(self.compartment_id) - 1 + ].surfacePointsNormals[ptInd] secondPoint = ( numpy.array(self.startingpoint) + numpy.array(v2) * self.uLength ) @@ -1634,7 +1644,7 @@ def getFirstPoint(self, ptInd, seed=0): secondPoint, dist=self.cutoff_boundary, jitter=self.max_jitter ) closeS = False - if inside and self.compNum <= 0: + if inside and self.compartment_id <= 0: # only if not surface ingredient closeS = self.far_enough_from_surfaces( secondPoint, cutoff=self.cutoff_surface @@ -1656,7 +1666,7 @@ def getFirstPoint(self, ptInd, seed=0): inside = self.env.grid.checkPointInside( secondPoint, dist=self.cutoff_boundary, jitter=self.max_jitter ) - if self.compNum <= 0: + if self.compartment_id <= 0: closeS = self.far_enough_from_surfaces( secondPoint, cutoff=self.cutoff_surface ) @@ -1716,8 +1726,10 @@ def grow_place( normal = None # jitter the first point - if self.compNum > 0: - normal = env.compartments[abs(self.compNum) - 1].surfacePointsNormals[ptInd] + if self.compartment_id > 0: + normal = env.compartments[ + abs(self.compartment_id) - 1 + ].surfacePointsNormals[ptInd] self.startingpoint = previousPoint = startingPoint = self.jitterPosition( numpy.array(env.grid.masterGridPositions[ptInd]), env.smallestProteinSize, @@ -1728,10 +1740,10 @@ def grow_place( # if u != self.uLength: # self.positions2 = [[self.vector]] - if self.compNum == 0: + if self.compartment_id == 0: compartment = self.env else: - compartment = self.env.compartments[abs(self.compNum) - 1] + compartment = self.env.compartments[abs(self.compartment_id) - 1] self.compartment = compartment secondPoint = self.getFirstPoint(ptInd) @@ -1766,7 +1778,9 @@ 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)] @@ -1830,7 +1844,7 @@ def grow_place( for i in range(len(self.results)): jtrans, rotMatj = self.results[-i] dist, ptInd = env.grid.getClosestGridPoint(jtrans) - compartment.molecules.append([jtrans, rotMatj, self, ptInd]) + compartment.molecules.append([jtrans, rotMatj, self, ptInd, self.radius]) # reset the result ? self.results = [] # print ("After :",listePtLinear) diff --git a/cellpack/autopack/ingredient/multi_cylinder.py b/cellpack/autopack/ingredient/multi_cylinder.py index 2162e438d..55086ad65 100644 --- a/cellpack/autopack/ingredient/multi_cylinder.py +++ b/cellpack/autopack/ingredient/multi_cylinder.py @@ -205,8 +205,8 @@ def collides_with_compartment( ptsInSphereId = numpy.take(pointsInCube, ptsWithinCaps[0], 0) compIdsSphere = numpy.take(env.grid.compartment_ids, ptsInSphereId, 0) - if self.compNum <= 0: - wrongPt = [cid for cid in compIdsSphere if cid != self.compNum] + if self.compartment_id <= 0: + wrongPt = [cid for cid in compIdsSphere if cid != self.compartment_id] if len(wrongPt): # print wrongPt return True @@ -354,8 +354,10 @@ def checkCylCollisions( ptsInSphereId = numpy.take(pointsInCube, ptsWithinCaps[0], 0) compIdsSphere = numpy.take(env.grid.compartment_ids, ptsInSphereId, 0) # print "compId",compIdsSphere - if self.compNum <= 0: - wrongPt = [cid for cid in compIdsSphere if cid != self.compNum] + if self.compartment_id <= 0: + wrongPt = [ + cid for cid in compIdsSphere if cid != self.compartment_id + ] if len(wrongPt): return True, insidePoints, newDistPoints diff --git a/cellpack/autopack/ingredient/multi_sphere.py b/cellpack/autopack/ingredient/multi_sphere.py index 068cf83c8..7a8a00722 100644 --- a/cellpack/autopack/ingredient/multi_sphere.py +++ b/cellpack/autopack/ingredient/multi_sphere.py @@ -33,6 +33,7 @@ def __init__( offset=[0, 0, 0], orient_bias_range=[-pi, pi], overwrite_distance_function=True, # overWrite + object_name=None, packing_mode="random", packing=0, partners=None, @@ -60,6 +61,7 @@ def __init__( molarity=molarity, name=name, jitter_attempts=jitter_attempts, + object_name=object_name, offset=offset, orient_bias_range=orient_bias_range, packing_mode=packing_mode, @@ -88,6 +90,43 @@ def __init__( if name is None: name = "%s_%f" % (str(self.radii), molarity) + def get_radius(self): + return self.encapsulating_radius + + def get_new_distance_values( + self, + packed_position, + packed_rotation, + gridPointsCoords, + distance, + dpad, + level, + ): + inside_points = {} + new_dist_points = {} + padded_sphere = self.radius + dpad + ptsInSphere = self.env.grid.getPointsInSphere(packed_position, padded_sphere) + delta = numpy.take(gridPointsCoords, ptsInSphere, 0) - packed_position + delta *= delta + distA = numpy.sqrt(delta.sum(1)) + for pti in range(len(ptsInSphere)): + pt = ptsInSphere[pti] + dist = distA[pti] + d = dist - self.radius + if d <= 0: # point is inside dropped sphere + if pt in inside_points: + if abs(d) < abs(inside_points[pt]): + inside_points[pt] = d + else: + inside_points[pt] = d + elif d < distance[pt]: # point in region of influence + if pt in new_dist_points: + if d < new_dist_points[pt]: + new_dist_points[pt] = d + else: + new_dist_points[pt] = d + return inside_points, new_dist_points + def collision_jitter( self, jtrans, @@ -255,6 +294,7 @@ def pack_at_grid_pt_location( grid_point_distances, inside_points, new_dist_points, + pt_index, ): level = self.deepest_level centers = self.positions[level] @@ -263,6 +303,7 @@ def pack_at_grid_pt_location( jtrans, rotation_matrix, centers ) # centers) grid_points_coords = env.grid.masterGridPositions + self.store_packed_object(jtrans, rotation_matrix, pt_index) for radius_of_sphere_in_tree, pos_of_sphere in zip(radii, transformed_centers): radius_of_area_to_check = ( radius_of_sphere_in_tree + dpad @@ -321,7 +362,7 @@ def collides_with_compartment( for radc, posc in zip(radii, transformed_centers): ptsInSphere = env.grid.getPointsInSphere(posc, radc[0]) # indices compIdsSphere = numpy.take(env.grid.compartment_ids, ptsInSphere, 0) - wrongPt = [cid for cid in compIdsSphere if cid != self.compNum] + wrongPt = [cid for cid in compIdsSphere if cid != self.compartment_id] if len(wrongPt): return True return False diff --git a/cellpack/autopack/ingredient/single_cube.py b/cellpack/autopack/ingredient/single_cube.py index 0f83cf506..6dd639df6 100644 --- a/cellpack/autopack/ingredient/single_cube.py +++ b/cellpack/autopack/ingredient/single_cube.py @@ -251,8 +251,8 @@ def collides_with_compartment( ptinsideId = numpy.take(pointsInCube, ptinside, 0) compIdsSphere = numpy.take(env.grid.compartment_ids, ptinsideId, 0) # print "compId",compIdsSphere - if self.compNum <= 0: - wrongPt = [cid for cid in compIdsSphere if cid != self.compNum] + if self.compartment_id <= 0: + wrongPt = [cid for cid in compIdsSphere if cid != self.compartment_id] if len(wrongPt): # print wrongPt return True diff --git a/cellpack/autopack/ingredient/single_cylinder.py b/cellpack/autopack/ingredient/single_cylinder.py index 1b3ab3a61..b1e291f76 100644 --- a/cellpack/autopack/ingredient/single_cylinder.py +++ b/cellpack/autopack/ingredient/single_cylinder.py @@ -206,8 +206,8 @@ def collides_with_compartment( ptsInSphereId = numpy.take(pointsInCube, ptsWithinCaps[0], 0) compIdsSphere = numpy.take(env.grid.compartment_ids, ptsInSphereId, 0) - if self.compNum <= 0: - wrongPt = [cid for cid in compIdsSphere if cid != self.compNum] + if self.compartment_id <= 0: + wrongPt = [cid for cid in compIdsSphere if cid != self.compartment_id] if len(wrongPt): # print wrongPt return True @@ -355,8 +355,10 @@ def checkCylCollisions( ptsInSphereId = numpy.take(pointsInCube, ptsWithinCaps[0], 0) compIdsSphere = numpy.take(env.grid.compartment_ids, ptsInSphereId, 0) # print "compId",compIdsSphere - if self.compNum <= 0: - wrongPt = [cid for cid in compIdsSphere if cid != self.compNum] + if self.compartment_id <= 0: + wrongPt = [ + cid for cid in compIdsSphere if cid != self.compartment_id + ] if len(wrongPt): return True, insidePoints, newDistPoints diff --git a/cellpack/autopack/ingredient/single_sphere.py b/cellpack/autopack/ingredient/single_sphere.py index 23625b926..55bd1f98a 100644 --- a/cellpack/autopack/ingredient/single_sphere.py +++ b/cellpack/autopack/ingredient/single_sphere.py @@ -52,6 +52,7 @@ def __init__( resolution_dictionary=None, rotation_axis=[0.0, 0.0, 0.0], rotation_range=0, + size_options=None, use_orient_bias=False, use_rotation_axis=True, weight=0.2, # use for affinity ie partner.weight @@ -83,6 +84,7 @@ def __init__( representations=representations, rotation_axis=rotation_axis, rotation_range=rotation_range, + size_options=size_options, use_rotation_axis=use_rotation_axis, weight=weight, ) @@ -96,6 +98,44 @@ def __init__( self.encapsulating_radius = radius self.min_radius = radius + @staticmethod + def create_circular_mask( + x_width, y_width, z_width, center=None, radius=None, voxel_size=None + ): + """ + Creates a circular mask of the given shape with the specified center + and radius + """ + if center is None: # use the middle of the image + center = (int(x_width / 2), int(y_width / 2), int(z_width / 2)) + if ( + radius is None + ): # use the smallest distance between the center and image walls + radius = min( + center[0], + center[1], + center[2], + x_width - center[0], + y_width - center[1], + z_width - center[2], + ) + + if voxel_size is None: + voxel_size = numpy.array([1, 1, 1], dtype=int) + + X, Y, Z = numpy.ogrid[:x_width, :y_width, :z_width] + dist_from_center = numpy.sqrt( + ((X - center[0]) * voxel_size[0]) ** 2 + + ((Y - center[1]) * voxel_size[1]) ** 2 + + ((Z - center[2]) * voxel_size[2]) ** 2 + ) + + mask = dist_from_center <= radius + return mask + + def get_radius(self): + return self.radius + def collision_jitter( self, jtrans, @@ -175,8 +215,8 @@ def collides_with_compartment( """ ptsInSphere = env.grid.getPointsInSphere(jtrans, self.radius) # indices compIdsSphere = numpy.take(env.grid.compartment_ids, ptsInSphere, 0) - if self.compNum <= 0: - wrongPt = [cid for cid in compIdsSphere if cid != self.compNum] + if self.compartment_id <= 0: + wrongPt = [cid for cid in compIdsSphere if cid != self.compartment_id] if len(wrongPt): return True return False @@ -253,43 +293,14 @@ def initialize_mesh(self, mesh_store): self.getData() - @staticmethod - def create_circular_mask( - x_width, y_width, z_width, center=None, radius=None, voxel_size=None - ): - """ - Creates a circular mask of the given shape with the specified center - and radius - """ - if center is None: # use the middle of the image - center = (int(x_width / 2), int(y_width / 2), int(z_width / 2)) - if ( - radius is None - ): # use the smallest distance between the center and image walls - radius = min( - center[0], - center[1], - center[2], - x_width - center[0], - y_width - center[1], - z_width - center[2], - ) - - if voxel_size is None: - voxel_size = numpy.array([1, 1, 1], dtype=int) - - X, Y, Z = numpy.ogrid[:x_width, :y_width, :z_width] - dist_from_center = numpy.sqrt( - ((X - center[0]) * voxel_size[0]) ** 2 - + ((Y - center[1]) * voxel_size[1]) ** 2 - + ((Z - center[2]) * voxel_size[2]) ** 2 - ) - - mask = dist_from_center <= radius - return mask - def create_voxelization( - self, image_data, bounding_box, voxel_size, image_size, position, rotation + self, + image_data, + bounding_box, + voxel_size, + image_size, + position, + **kwargs, ): """ Creates a voxelization for the sphere @@ -300,7 +311,7 @@ def create_voxelization( *image_size, center=voxelized_position, radius=self.radius, - voxel_size=voxel_size + voxel_size=voxel_size, ) image_data[mask] = 255 diff --git a/cellpack/autopack/interface_objects/packed_objects.py b/cellpack/autopack/interface_objects/packed_objects.py new file mode 100644 index 000000000..b15079fff --- /dev/null +++ b/cellpack/autopack/interface_objects/packed_objects.py @@ -0,0 +1,78 @@ +import numpy + + +class PackedObject: + def __init__( + self, + position, + rotation, + radius, + pt_index, + ingredient=None, + is_compartment=False, + ) -> None: + self.name = ingredient.name + self.position = position + self.rotation = rotation + self.radius = radius + self.encapsulating_radius = ingredient.encapsulating_radius + self.pt_index = pt_index + self.is_compartment = is_compartment + self.color = ingredient.color + self.ingredient = ingredient + + +class PackedObjects: + def __init__(self): + self._packed_objects = [] + + def add(self, new_object: PackedObject): + self._packed_objects.append(new_object) + + def get_radii(self): + radii = [] + for obj in self._packed_objects: + if not obj.is_compartment: + radii.append(obj.radius) + return radii + + def get_encapsulating_radii(self): + radii = [] + for obj in self._packed_objects: + if not obj.is_compartment: + radii.append(obj.encapsulating_radius) + return radii + + def get_positions(self): + positions = [] + for obj in self._packed_objects: + if not obj.is_compartment: + positions.append(obj.position) + return positions + + def get_positions_for_ingredient(self, ingredient_name): + return numpy.array( + [ + self._packed_objects[i].position + for i in range(len(self._packed_objects)) + if self._packed_objects[i].name == ingredient_name + ] + ) + + def get_rotations_for_ingredient(self, ingredient_name): + return numpy.array( + [ + self._packed_objects[i].rotation + for i in range(len(self._packed_objects)) + if self._packed_objects[i].name == ingredient_name + ] + ) + + def get_ingredients(self): + return [obj for obj in self._packed_objects if obj.is_compartment is False] + + def get_compartment(self): + return [obj for obj in self._packed_objects if obj.is_compartment is True] + + def get_all(self): + return self._packed_objects diff --git a/cellpack/autopack/loaders/config_loader.py b/cellpack/autopack/loaders/config_loader.py index 7b071af21..aaf15ff35 100644 --- a/cellpack/autopack/loaders/config_loader.py +++ b/cellpack/autopack/loaders/config_loader.py @@ -35,7 +35,7 @@ class ConfigLoader(object): "out": "out/", "overwrite_place_method": False, "parallel": False, - "place_method": "jitter", + "place_method": "spheresSST", "randomness_seed": None, "save_analyze_result": False, "save_converted_recipe": False, diff --git a/cellpack/autopack/plotly_result.py b/cellpack/autopack/plotly_result.py index 86268aa4d..c55a901c8 100644 --- a/cellpack/autopack/plotly_result.py +++ b/cellpack/autopack/plotly_result.py @@ -59,12 +59,15 @@ def add_square(self, side_length, pos, rotMat, color, opacity=1): ) def add_ingredient_positions(self, env): - for pos, rot, ingr, ptInd in env.molecules: + for obj in env.packed_objects.get_ingredients(): + ingr = obj.ingredient if hasattr(ingr, "positions") and len(ingr.positions) > 1: + self.add_circle(obj.radius, obj.position, obj.color) + for level in range(len(ingr.positions)): for i in range(len(ingr.positions[level])): position = ingr.apply_rotation( - rot, ingr.positions[level][i], pos + obj.rotation, ingr.positions[level][i], obj.position ) self.add_circle( ingr.radii[level][i], @@ -74,14 +77,16 @@ def add_ingredient_positions(self, env): ) else: if ingr.model_type == "Spheres": - self.add_circle(ingr.encapsulating_radius, pos, ingr.color) - elif ingr.model_type == "Cube": - self.add_square(ingr.radii, pos, rot, ingr.color) + self.add_circle(obj.radius, obj.position, obj.color) + if ingr.model_type == "Cube": + self.add_square(ingr.radii, obj.position, obj.rotation, obj.color) elif ingr.model_type == "Cylinders": length = ingr.length width = 2 * ingr.radii[0][0] side_lengths = [[width, length, 1.0]] - self.add_square(side_lengths, pos, rot, ingr.color) + self.add_square( + side_lengths, obj.position, obj.rotation, ingr.color + ) def make_grid_heatmap(self, env): ids = [] diff --git a/cellpack/autopack/upy/simularium/simularium_helper.py b/cellpack/autopack/upy/simularium/simularium_helper.py index 70a8c03fd..777e75e0f 100644 --- a/cellpack/autopack/upy/simularium/simularium_helper.py +++ b/cellpack/autopack/upy/simularium/simularium_helper.py @@ -484,10 +484,12 @@ def init_scene_with_objects( ): self.time = 0 instance_number = 0 - for position, rotation, ingredient, ptInd in objects: - ingr_name = ingredient.name + for packed_object in objects: + ingr_name = packed_object.name + ingredient = packed_object.ingredient sub_points = None - if self.is_fiber(ingredient.type): + radius = packed_object.radius + if simulariumHelper.is_fiber(ingredient.type): if ingredient.nbCurve == 0: continue # TODO: get sub_points accurately @@ -496,25 +498,25 @@ def init_scene_with_objects( sub_points = ingredient.listePtLinear if ingr_name not in self.display_data: display_type, url = self.get_display_data(ingredient) - self.display_data[ingredient.name] = DisplayData( + self.display_data[packed_object.name] = DisplayData( name=ingr_name, display_type=display_type, url=url, color=simulariumHelper.format_rgb_color(ingredient.color), ) - radius = ingredient.encapsulating_radius if ingredient is not None else 10 + radius = radius if radius is not None else 10 adj_pos = ingredient.representations.get_adjusted_position( - position, rotation + packed_object.position, packed_object.rotation ) self.add_instance( ingr_name, ingredient, - f"{ingr_name}-{ptInd}-{instance_number}", + f"{ingr_name}-{packed_object.pt_index}-{instance_number}", radius, adj_pos, - rotation, + packed_object.rotation, sub_points, ) instance_number += 1 @@ -523,16 +525,18 @@ def init_scene_with_objects( for level in range(len(ingredient.positions)): for i in range(len(ingredient.positions[level])): pos = ingredient.apply_rotation( - rotation, ingredient.positions[level][i], position + packed_object.rotation, + ingredient.positions[level][i], + adj_pos, ) self.add_instance( f"{ingredient.name}-spheres", ingredient, - f"{ingredient.name}-{ptInd}-{i}", + f"{ingredient.name}-{packed_object.pt_index}-{instance_number}-{i}", ingredient.radii[level][i], pos, - rotation, + packed_object.rotation, None, ) diff --git a/cellpack/autopack/utils.py b/cellpack/autopack/utils.py index 745d4e178..1f720d006 100644 --- a/cellpack/autopack/utils.py +++ b/cellpack/autopack/utils.py @@ -181,3 +181,79 @@ def load_object_from_pickle(pickle_file_object): except Exception as e: raise ValueError(f"Error loading saved object: {e}") return output_object + + +def get_min_value_from_distribution(distribution_options, return_int=False): + """ + Returns a low bound on the value from a distribution + """ + value = None + if distribution_options.get("distribution") == "uniform": + value = distribution_options.get("min", 1) + + if distribution_options.get("distribution") == "normal": + value = distribution_options.get("mean", 0) - 2 * distribution_options.get( + "std", 1 + ) + + if distribution_options.get("distribution") == "list": + value = numpy.nanmin(distribution_options.get("list_values", None)) + + if return_int and value is not None: + value = int(numpy.rint(value)) + + return value + + +def get_max_value_from_distribution(distribution_options, return_int=False): + """ + Returns a high bound on the value from a distribution + """ + value = None + if distribution_options.get("distribution") == "uniform": + value = distribution_options.get("max", 1) + + if distribution_options.get("distribution") == "normal": + value = distribution_options.get("mean", 0) + 2 * distribution_options.get( + "std", 1 + ) + + if distribution_options.get("distribution") == "list": + value = numpy.nanmax(distribution_options.get("list_values", None)) + + if return_int and value is not None: + value = int(numpy.rint(value)) + + return value + + +def get_value_from_distribution(distribution_options, return_int=False): + """ + Returns a value from the distribution options + """ + if distribution_options.get("distribution") == "uniform": + if return_int: + return int( + numpy.random.randint( + distribution_options.get("min", 0), + distribution_options.get("max", 1), + ) + ) + else: + return numpy.random.uniform( + distribution_options.get("min", 0), + distribution_options.get("max", 1), + ) + if distribution_options.get("distribution") == "normal": + value = numpy.random.normal( + distribution_options.get("mean", 0), distribution_options.get("std", 1) + ) + elif distribution_options.get("distribution") == "list": + value = numpy.random.choice(distribution_options.get("list_values", None)) + else: + value = None + + if return_int and value is not None: + value = int(numpy.rint(value)) + + return value diff --git a/cellpack/tests/packing-configs/test_variable_size_config.json b/cellpack/tests/packing-configs/test_variable_size_config.json new file mode 100644 index 000000000..d33e9a154 --- /dev/null +++ b/cellpack/tests/packing-configs/test_variable_size_config.json @@ -0,0 +1,17 @@ +{ + "name": "test_variable_size_config", + "format": "simularium", + "inner_grid_method": "trimesh", + "live_packing": false, + "ordered_packing": false, + "out": "cellpack/tests/outputs/", + "overwrite_place_method": true, + "number_of_packings": 10, + "place_method": "spheresSST", + "save_analyze_result": true, + "show_grid_plot": false, + "spacing": null, + "parallel": false, + "use_periodicity": false, + "show_sphere_trees": true +} \ No newline at end of file diff --git a/cellpack/tests/recipes/v2/test_variable_size.json b/cellpack/tests/recipes/v2/test_variable_size.json new file mode 100644 index 000000000..5c3dd30aa --- /dev/null +++ b/cellpack/tests/recipes/v2/test_variable_size.json @@ -0,0 +1,110 @@ +{ + "version": "1.0.0", + "format_version": "2.0", + "name": "test_variable_size", + "bounding_box": [ + [ + 0, + 0, + 0 + ], + [ + 500, + 500, + 1 + ] + ], + "objects": { + "base": { + "jitter_attempts": 10, + "orient_bias_range": [ + -3.1415927, + 3.1415927 + ], + "rotation_range": 6.2831, + "cutoff_boundary": 0, + "max_jitter": [ + 1, + 1, + 0 + ], + "perturb_axis_amplitude": 0.1, + "packing_mode": "random", + "principal_vector": [ + 0, + 0, + 1 + ], + "rejection_threshold": 50, + "place_method": "spheresSST", + "cutoff_surface": 42, + "rotation_axis": [ + 0, + 0, + 1 + ], + "available_regions": { + "interior": {}, + "surface": {}, + "outer_leaflet": {}, + "inner_leaflet": {} + } + }, + "sphere_25": { + "type": "single_sphere", + "inherit": "base", + "color": [ + 0.5, + 0.5, + 0.5 + ], + "radius": 25, + "max_jitter": [ + 1, + 1, + 0 + ] + }, + "sphere_40": { + "type": "single_sphere", + "inherit": "base", + "color": [ + 0.1, + 0.5, + 0.5 + ], + "radius": 40, + "max_jitter": [ + 1, + 1, + 0 + ] + } + }, + "composition": { + "space": { + "regions": { + "interior": [ + "A", "B" + ] + } + }, + "A": { + "object": "sphere_25", + "count": 200, + "size_options": { + "distribution": "uniform", + "min": 20, + "max": 25 + } + }, + "B": { + "object": "sphere_40", + "count": 200, + "size_options": { + "distribution": "list", + "list_values": [35, 40, 41] + } + } + } +} \ No newline at end of file diff --git a/cellpack/tests/test_ingredient.py b/cellpack/tests/test_ingredient.py index 8dd493234..2e5b6344a 100644 --- a/cellpack/tests/test_ingredient.py +++ b/cellpack/tests/test_ingredient.py @@ -52,7 +52,7 @@ "distribution": "invalid_distribution", }, }, - "invalid_distribution is not a valid count distribution", + "invalid_distribution is not a valid distribution", ), ( { diff --git a/examples/recipes/v2/sphere_tree.json b/examples/recipes/v2/sphere_tree.json index 5b45b8991..b54e5e1d5 100644 --- a/examples/recipes/v2/sphere_tree.json +++ b/examples/recipes/v2/sphere_tree.json @@ -11,7 +11,7 @@ [ 1000, 1000, - 1000 + 1 ] ], "objects": { @@ -32,8 +32,8 @@ "place_method": "jitter", "cutoff_surface": 42, "rotation_axis": [ - 1, - 1, + 0, + 0, 1 ], "available_regions": { @@ -62,7 +62,7 @@ "max_jitter": [ 1, 1, - 1 + 0 ] } },