# HG changeset patch # User Matthew Turk # Date 1247673596 14400 # Branch hopwrap # Node ID 60615619173287c363a6f12088c6361c45735deb # Parent cdbd4eedc0b83e628bd95c14aec1df7771135335 # Parent 7098c3bac6551b9e89ebdf10799d328c1624c2fd Merging with the latest hopwrap branch diff -r cdbd4eedc0b83e628bd95c14aec1df7771135335 -r 60615619173287c363a6f12088c6361c45735deb yt/lagos/HaloFinding.py --- a/yt/lagos/HaloFinding.py Wed Jul 15 11:59:15 2009 -0400 +++ b/yt/lagos/HaloFinding.py Wed Jul 15 11:59:56 2009 -0400 @@ -34,7 +34,7 @@ pass from kd import * -import math +import math,sys class Halo(object): """ @@ -428,10 +428,118 @@ else: return False + def _build_task_chain_map(self): + """ + From padded_particles build a mapping that identifies chains on each + task to a different chain on another task. + """ + self.particle_fields["particle_index"] = self._data_source["particle_index"][self._base_indices] + # Figure out our offset + #my_first_id = sum([v for k,v in halo_info.items() if k < mine]) + #after = my_first_id + max(self.tags) + 1 + # build the list of [particle_index, taskID, grpID] + if len(self.padded_particles): + self.particle_index_map = self.particle_fields["particle_index"][self.padded_particles] + tag_map = self.tags[self.padded_particles] + temp_map = [] + for i,particle in enumerate(self.particle_index_map): + temp_map.append([particle,self.mine,tag_map[i]]) # real ID, task ID, chainID + self.particle_index_map = na.array(temp_map) + else: + self.particle_index_map = na.array([]) + # concatenate the map on all tasks, it would be faster to + # intelligently send this to only the processors that need each + # part of the data, but that's harder. Maybe later. + #mylog.info("%s" % str(self.particle_index_map)) + self.particle_index_map = self._mpi_catarray(self.particle_index_map) + if len(self.particle_index_map) == 0: + # none of the chains cross the boundaries, so we're done here. + mylog.info("No chains need to be merged.") + return + # but if there are crossings, move forward. + # On each processor, if one of the particles it owns (in the real + # region!) is in the map, add an entry to the task_chain_map + self.task_chain_map = [] + particle_map_column = self.particle_index_map[:,0] + for index,ownPart in enumerate(self.particle_fields["particle_index"]): + if ownPart in particle_map_column: + remote_index = 0 + while particle_map_column[remote_index] != ownPart: + remote_index += 1 + if self.is_inside[index]: + #mylog.info('point %s bound %s ownPart %s' % \ + #(str(point),str(self.bounds),str(ownPart))) + # find out which local chain this particle lives in + local_chainID = self.particle_fields["tags"][index] + # add a mapping, if it isn't already there + # put lower-indexed task ID first + if self.particle_index_map[remote_index][1] < self.mine: + map = [self.particle_index_map[remote_index][1],\ + self.particle_index_map[remote_index][2], \ + self.mine, local_chainID] + else: + map = [self.mine, local_chainID, \ + self.particle_index_map[remote_index][1],\ + self.particle_index_map[remote_index][2] ] + #mylog.info("map %s" % str(map)) + if map not in self.task_chain_map: + self.task_chain_map.append(map) + # concatenate the task_chain_map on all tasks + self.task_chain_map = self._mpi_catlist(self.task_chain_map) + #mylog.info("end %s" % str(self.task_chain_map)) + + def _build_task_densest_in_chain(self): + """ + 'Upgrade' self.densest_in_chain to include the task each chain is on. + The result will be accessed: self.densest_in_chain[task][chainID]. + """ + temp = {} + temp[self.mine] = self.densest_in_chain + self._mpi_joindict(temp) + self.densest_in_chain = temp + + def _build_task_chain_densest_n(self): + """ + 'Upgrade' self.chain_densest_n to include the task each chain is on. + For now we'll be cheap and make the keys into strings, with the form + 'taskID,chainID' + """ + temp = {} + for chain_high in self.chain_densest_n: + h = '%d,%d' % (self.mine, chain_high) + temp[h] = {} + for chain_low in self.chain_densest_n[chain_high]: + l = '%d,%d' % (self.mine, chain_low) + temp[h][l] = self.chain_densest_n[chain_high][chain_low] + self._mpi_joindict(temp) + self.chain_densest_n = temp + + def _update_chain_densest_n(self): + """ + Update chain_densest_n to reflect the potential neighbors from + chain_potential_n. This involves figuring out the real + particle_index for each of the potential partcles, and communicating + this to the neighboring task. Then, if that task finds that particle + is in a chain on in its domain, adding that information to a list which + is then communicated back to the first. chain_densest_n is updated to + consider the neighboring chains, taking note of identified chains in + task_chain_map. + """ + + def _update_chain_densest_n(self): + """ + Using chain_potential_n as a source of potential mergings between + chains, update chain_densest_n while refering to self.task_chain_map + as to not merge chains that are already identified as being + connected. + """ + + + def _run_finder(self): - xt = self.particle_fields["particle_position_x"] #- 0.961 - yt = self.particle_fields["particle_position_y"] #- 0.369 - zt = self.particle_fields["particle_position_z"] #- 0.710 + xt = self.particle_fields["particle_position_x"] - 0.961 + yt = self.particle_fields["particle_position_y"] - 0.369 + zt = self.particle_fields["particle_position_z"] - 0.710 for i,x in enumerate(xt): if x < 0: xt[i] = 1+x @@ -442,70 +550,22 @@ if z < 0: zt[i] = 1+z obj = RunChainHOP(self.period, self.padding, - self.num_neighbors, self.bounds, - xt, - yt, - zt, + self.num_neighbors, self.bounds,xt,yt,zt, self.particle_fields["ParticleMassMsun"]/self.total_mass, self.threshold) self.densities, self.tags = obj.density, obj.chainID + self.padded_particles = obj.padded_particles + self.is_inside = obj.is_inside + self.densest_in_chain = obj.densest_in_chain + self.chain_densest_n = obj.chain_densest_n self.particle_fields["densities"] = self.densities self.particle_fields["tags"] = self.tags if self._distributed: - self.padded_particles = obj.padded_particles - self.particle_fields["particle_index"] = self._data_source["particle_index"][self._base_indices] - mine, halo_info = self._mpi_info_dict(max(self.tags) + 1) - nhalos = sum(halo_info.values()) - # Figure out our offset - my_first_id = sum([v for k,v in halo_info.items() if k < mine]) - after = my_first_id + max(self.tags) + 1 - # build the list of [particle_index, threadID, grpID] - if len(self.padded_particles): - self.particle_index_map = self.particle_fields["particle_index"][self.padded_particles] - tag_map = self.tags[self.padded_particles] - temp_map = [] - for i,particle in enumerate(self.particle_index_map): - temp_map.append([particle,mine,tag_map[i]]) - self.particle_index_map = na.array(temp_map) - else: - self.particle_index_map = na.array([]) - # concatenate the map on all threads, it would be faster to - # intelligently send this to only the processors that need each - # part of the data, but that's harder. Maybe later. - #mylog.info("%s" % str(self.particle_index_map)) - self.particle_index_map = self._mpi_catarray(self.particle_index_map) - if len(self.particle_index_map) == 0: - # none of the groups cross the boundaries, so we're done here. - mylog.info("No groups need to be merged.") - return - # but if there are crossings, move forward. - # On each processor, if one of the particles it owns (in the real - # region!) is in the map, add an entry to the thread_group_map - self.thread_group_map = [] - particle_map_column = self.particle_index_map[:,0] - for index,ownPart in enumerate(self.particle_fields["particle_index"]): - if ownPart in particle_map_column: - remote_index = 0 - while particle_map_column[remote_index] != ownPart: - remote_index += 1 - point = [xt[index], - yt[index], - zt[index]] - if self._is_inside(point) is True: - #mylog.info('point %s bound %s ownPart %s' % \ - #(str(point),str(self.bounds),str(ownPart))) - # find out which local group this particle lives in - local_groupID = self.particle_fields["tags"][index] - # add a mapping, if it isn't already there - # remote->local, remote has priority! - map = [self.particle_index_map[remote_index][1],\ - self.particle_index_map[remote_index][2], \ - mine, local_groupID] - if map not in self.thread_group_map: - self.thread_group_map.append(map) - # concatenate the thread_group_map on all threads - self.thread_group_map = self._mpi_catlist(self.thread_group_map) - print self.thread_group_map + self.mine, halo_info = self._mpi_info_dict(max(self.tags) + 1) + #nhalos = sum(halo_info.values()) + self._build_task_chain_map() + self._build_task_densest_in_chain() + self._build_task_chain_densest_n() def write_out(self, filename="chainHopAnalysis.out"): HaloList.write_out(self, filename) diff -r cdbd4eedc0b83e628bd95c14aec1df7771135335 -r 60615619173287c363a6f12088c6361c45735deb yt/lagos/chainHOP/chainHOP.py --- a/yt/lagos/chainHOP/chainHOP.py Wed Jul 15 11:59:15 2009 -0400 +++ b/yt/lagos/chainHOP/chainHOP.py Wed Jul 15 11:59:56 2009 -0400 @@ -2,6 +2,7 @@ from bisect import insort from sets import Set +from yt.lagos import * from yt.extensions.kdtree import * from Forthon import * @@ -54,6 +55,8 @@ points = fKD.pos.transpose() self.is_inside = ( (points >= LE).all(axis=1) * \ (points < RE).all(axis=1) ) + # also mark padded particles that are in a right side padding: + self.padded_right = (points >= RE).any(axis=1) def _densestNN(self): """ @@ -142,11 +145,10 @@ by finding the highest boundary density neighbor for each chain. Some chains will have no neighbors! """ - self.chain_densest_n = {} # chainID -> {chainIDs}->boundary dens - self.densest_boundary = {} # chainID -> boundary_density + self.chain_densest_n = {} # chainID -> {chainIDs->boundary dens} + self.chain_potential_n = {} # chainID -> [partIDs] for i in xrange(chain_count): self.chain_densest_n[i] = {} - self.densest_boundary[i] = 0.0 self.reverse_map = {} # chainID -> groupID, one to one for i in xrange(chain_count): self.reverse_map[i] = -1 @@ -170,33 +172,51 @@ # This is probably the same as above, but it's OK. # It can be removed later. if thisNN==i or thisNN==self.densestNN[i]: continue - # If our neighbor is not part of a chain, move on. - if thisNN_chainID < 0: continue - # Find thisNN's chain's max_dens. - thisNN_max_dens = self.densest_in_chain[thisNN_chainID] - # Calculate the two groups boundary density. - boundary_density = (self.density[thisNN] + self.density[i]) / 2. - # Find out who's denser. - if thisNN_max_dens >= part_max_dens: - higher_chain = thisNN_chainID - lower_chain = self.chainID[i] - else: - higher_chain = self.chainID[i] - lower_chain = thisNN_chainID - # See if this boundary density is higher than - # previously recorded for this pair of chains. - # Links only go 'downhill'. - try: - old = self.chain_densest_n[higher_chain][lower_chain] - if old < boundary_density: - # make this the new densest boundary between this pair + # The particles that will be added to chain_potential_n all + # have no chainID, and everything immediately below is for + # particles with a chainID. + if thisNN_chainID >= 0: + # If our neighbor is not part of a chain, move on + #if thisNN_chainID < 0: continue + # Find thisNN's chain's max_dens. + thisNN_max_dens = self.densest_in_chain[thisNN_chainID] + # Calculate the two groups boundary density. + boundary_density = (self.density[thisNN] + self.density[i]) / 2. + # Find out who's denser. + if thisNN_max_dens >= part_max_dens: + higher_chain = thisNN_chainID + lower_chain = self.chainID[i] + else: + higher_chain = self.chainID[i] + lower_chain = thisNN_chainID + # See if this boundary density is higher than + # previously recorded for this pair of chains. + # Links only go 'downhill'. + try: + old = self.chain_densest_n[higher_chain][lower_chain] + if old < boundary_density: + # make this the new densest boundary between this pair + self.chain_densest_n[higher_chain][lower_chain] = \ + boundary_density + except KeyError: + # we haven't seen this pairing before, record this as the + # new densest boundary between chains self.chain_densest_n[higher_chain][lower_chain] = \ boundary_density - except KeyError: - # we haven't seen this pairing before, record this as the - # new densest boundary between chains - self.chain_densest_n[higher_chain][lower_chain] = \ - boundary_density + # now consider right padded, non-assigned particles + elif self.padded_right[thisNN] and thisNN_chainID < 0: + # see if we need to create an entry + try: + potential_n = self.chain_potential_n[self.chainID[i]] + except KeyError: + self.chain_potential_n[self.chainID[i]] = [] + potential_n = [] + # don't add this particle if it's already there + if thisNN not in potential_n: + self.chain_potential_n[self.chainID[i]].append(thisNN) + # Skip non-assigned particles that aren't right padded + else: + continue del self.reverse_map[-1] def _build_groups(self, chain_count): @@ -366,11 +386,11 @@ # Connect the chains into groups. mylog.info('Connecting %d chains into groups...' % chain_count) self._connect_chains(chain_count) - group_count = self._build_groups(chain_count) + #group_count = self._build_groups(chain_count) #mylog.info('Paring %d groups...' % group_count) # Don't pare groups here, it has to happen later. #self._pare_groups_by_max_dens() - self._translate_groupIDs(group_count) - mylog.info('Found %d groups...' % group_count) + #self._translate_groupIDs(group_count) + #mylog.info('Found %d groups...' % group_count)