From cf43a0bec8703b4fad7d44b0a0577cb5216438f0 Mon Sep 17 00:00:00 2001 From: Kiyoshi Masui Date: Wed, 22 Jan 2014 16:54:06 -0800 Subject: [PATCH 01/28] Changes to mapmaker cython to allow subdivision of work. --- map/_mapmaker.pyx | 54 +++++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 52 insertions(+), 2 deletions(-) diff --git a/map/_mapmaker.pyx b/map/_mapmaker.pyx index 6ff00540..5c696393 100644 --- a/map/_mapmaker.pyx +++ b/map/_mapmaker.pyx @@ -194,18 +194,49 @@ def update_map_noise_independant_chan( np.ndarray[np.int_t, ndim=3, mode='c'] pointing_inds not None, np.ndarray[DTYPE_t, ndim=2, mode='c'] pointing_weights not None, f_ind_in, - np.ndarray[DTYPE_t, ndim=4, mode='c'] map_noise_inv not None): + np.ndarray[DTYPE_t, ndim=4, mode='c'] map_noise_inv not None, + ra0_ind_range=None, + dec0_ind_range=None, + ): """Convert noise to map space. This function is only used when ignoring frequency correlations. The noise matrix is update one frequency slice at a time for performance. - """ + The total noise covariance matrix is 5D (freq, ra0, dec0, ra1, dec1). This + function caculates at a since frequency, a range of r0, a range of dec1, + all ra1 and all dec1. + + Note that ``map_noise_inv.shape[0]`` should equal + ``ra0_ind_range[1] - ra0_ind_range[0]`` and likewise for + ``map_noise_inv.shape[1]`` and the dec range. + + """ + + # TODO: Change the order of the arguments to a more sane scheme, and make + # corresponding change in dirty_map.py. + # XXX: Currently the arguments/function is backward compatible with before + # the ra0, dec0 range subdivision. + # Shapes. cdef int n_chan = diagonal_inv.shape[0] cdef int n_time = diagonal_inv.shape[1] cdef int n_pix_per_pointing = pointing_inds.shape[2] cdef int q = time_modes.shape[0] + cdef int ra_ind_start, ra_ind_end + cdef int dec_ind_start, dec_ind_end + if not ra0_ind_range: + ra_ind_start = 0 + ra_ind_end = map_noise_inv.shape[0] + else: + ra_ind_start = ra0_ind_range[0] + ra_ind_end = ra0_ind_range[1] + if not dec0_ind_range: + dec_ind_start = 0 + dec_ind_end = map_noise_inv.shape[1] + else: + dec_ind_start = dec0_ind_range[0] + dec_ind_end = dec0_ind_range[1] # Indecies. cdef int time_ind cdef int f_ind = f_ind_in @@ -227,6 +258,17 @@ def update_map_noise_independant_chan( # Now loop over the pointings. with nogil: for time_ind in xrange(n_time): + # Check if this time touches the ra and dec range we are + # updating. + for jj in xrange(n_pix_per_pointing): + if (pointing_inds_transpose[time_ind,jj,0] >= ra_ind_start + and pointing_inds_transpose[time_ind,jj,0] < ra_ind_end + and pointing_inds_transpose[time_ind,jj,1] >= dec_ind_start + and pointing_inds_transpose[time_ind,jj,1] < dec_ind_end + ): + break + else: + continue # Reset the time noise row to zero. for jj in xrange(n_time): update_term[jj] = 0 @@ -246,6 +288,14 @@ def update_map_noise_independant_chan( for ii in xrange(n_pix_per_pointing): ra_ind = pointing_inds[time_ind,0,ii] dec_ind = pointing_inds[time_ind,1,ii] + if not (ra_ind >= ra_ind_start + and ra_ind < ra_ind_end + and dec_ind >= dec_ind_start + and dec_ind < dec_ind_end + ): + continue + ra_ind = ra_ind - ra_ind_start + dec_ind = dec_ind - dec_ind_start weight = pointing_weights[time_ind,ii] # Loop over the time axes to convert to pixel and accumulate in # the output matrix. From 7b1d73882e9b701953c5c78caa3e56c3b0e92795 Mon Sep 17 00:00:00 2001 From: Chris Anderson Date: Thu, 13 Feb 2014 17:31:32 -0500 Subject: [PATCH 02/28] parallel_dirty_map_tcs.py lets earch thread work on ra subchunk --- input/km/mm_test.ini | 2 +- map/dirty_map.py | 4 +- map/parallel_dirty_map_tcs.py | 1002 +++++++++++++++++++++++++++++++++ 3 files changed, 1005 insertions(+), 3 deletions(-) create mode 100644 map/parallel_dirty_map_tcs.py diff --git a/input/km/mm_test.ini b/input/km/mm_test.ini index 8f7c695c..39c227e8 100644 --- a/input/km/mm_test.ini +++ b/input/km/mm_test.ini @@ -31,7 +31,7 @@ map_spacing = .0627 raid_pro = os.getenv("RAID_PRO") input_data_dir = raid_pro + "kiyo/gbt_out_new/" base_dir = os.getenv("GBT_OUT") -map_dir = base_dir + 'maps/first_lockandwrite_test/' +map_dir = base_dir + 'maps/first_lockandwrite_test/tcs_threading/' #map_root = map_dir + "tmp_test_parallel_freqcorr_" map_root = map_dir + "tmp_test_parallel_inherit_" #map_root = map_dir + "tmp_test_regular_" diff --git a/map/dirty_map.py b/map/dirty_map.py index 76f08fdc..93cdf94b 100644 --- a/map/dirty_map.py +++ b/map/dirty_map.py @@ -1189,7 +1189,7 @@ def noise_to_map_domain(self, Noise, f_ind, ra_ind, map_noise_inv): " `noise_channel_to_map` instead.") raise RuntimeError(msg) - def noise_channel_to_map(self, Noise, f_ind, map_noise_inv): + def noise_channel_to_map(self, Noise, f_ind, map_noise_inv, ra0_range=None, dec0_range=None): """Convert noise to map space. Use this function over `noise_to_map_domain` if the noise has no @@ -1200,7 +1200,7 @@ def noise_channel_to_map(self, Noise, f_ind, map_noise_inv): if not Noise._frequency_correlations: _mapmaker_c.update_map_noise_independant_chan(Noise.diagonal_inv, Noise.time_modes, Noise.time_mode_update, self._pixel_inds, - self._weights, f_ind, map_noise_inv) + self._weights, f_ind, map_noise_inv, ra0_range, dec0_range) else: msg = ("Noise object has frequency correlations. Use " " `noise_to_map_domain` instead.") diff --git a/map/parallel_dirty_map_tcs.py b/map/parallel_dirty_map_tcs.py new file mode 100644 index 00000000..5ee4cffb --- /dev/null +++ b/map/parallel_dirty_map_tcs.py @@ -0,0 +1,1002 @@ +"""Dirty map making module. + +Module converts data in the time domain into noise weighted data in the map +domain, i.e. it creats the dirty map. Module also contains many utilities like +the pointing operator (`Pointing`) and the time domain noise operator +(`Noise`). +""" + +import sys + +from mpi4py import MPI + +import math +import threading +from Queue import Queue +import shelve +import sys +import time as time_mod +import cPickle +import h5py + +import scipy as sp +import numpy.ma as ma +import scipy.fftpack as fft +from scipy import linalg +from scipy import interpolate +import numpy as np +#import kiyopy.pickle_method +import warnings +#import matplotlib.pyplot as plt + +import core.algebra as al +from core import fitsGBT +import utils.misc as utils +import tools +from noise import noise_power +from foreground import ts_measure +import kiyopy.custom_exceptions as ce +import kiyopy.utils +from kiyopy import parse_ini +import _mapmaker as _mapmaker_c +from constants import T_infinity, T_huge, T_large, T_medium, T_small, T_sys +from constants import f_medium, f_large +from utils import misc +from map import dirty_map as dm + +prefix ='dm_' +params_init = {# IO: + 'input_root' : './', + # The unique part of every fname + 'file_middles' : ("testfile_GBTfits",), + 'input_end' : ".fits", + 'output_root' : "./testoutput_", + # Shelve file that holds measurements of the noise. + # What data to include from each file. + 'scans' : (), + 'IFs' : (0,), + 'polarizations' : ('I',), + # Map parameters (Ra (deg), Dec (deg)). + 'field_centre' : (325.0, 0.0), + # In pixels. + 'map_shape' : (5, 5), + 'pixel_spacing' : 0.5, # degrees + # Interpolation between pixel points. Options are 'nearest', + # 'linear' and 'cubic'. + 'interpolation' : 'linear', + # How to treat the data. + # How much data to include at a time (scan by scan or file by + # file) + 'time_block' : 'file', + # Number of files to process at a time. Increasing this number + # reduces IO but takes more memory. + 'n_files_group' : 0, + # What kind of frequency correlations to use. Options are + # 'None', 'mean' and 'measured'. + 'frequency_correlations' : 'None', + # Ignored unless 'frequency_correlations' is 'measured'. + 'number_frequency_modes' : 1, + # Of the above modes, how many to completely discard instead of + # noise weighting. + 'number_frequency_modes_discard' : 0, + # Where to get noise parameters from. + 'noise_parameter_file' : '', + 'deweight_time_mean' : True, + 'deweight_time_slope' : False, + # If there where any foregrounds subtracted in the time stream, + # let the noise know about it. + 'n_ts_foreground_modes' : 0, + 'ts_foreground_mode_file' : '' + } + +comm = MPI.COMM_WORLD + +class DirtyMapMaker(object): + """Dirty map maker. + """ + + def __init__(self, parameter_file_or_dict=None, feedback=2): + # Read in the parameters. + self.params = parse_ini.parse(parameter_file_or_dict, params_init, + prefix=prefix, feedback=feedback) + print self.params['file_middles'] + self.feedback = feedback + self.rank = comm.Get_rank() + self.nproc = comm.Get_size() + + def iterate_data(self, file_middles): + """An iterator over the input data. + + This can either iterate over the file names, processing the whole file + at once, or it can iterate over both the file and the scans. + + Returns the DataBlock objects as well as the file middle for that data + (which is used as a key in several parameter data bases). + """ + + params = self.params + for middle in file_middles: + self.file_middle = middle + fname = params["input_root"] + middle + params["input_end"] + Reader = fitsGBT.Reader(fname, feedback=self.feedback) + Blocks = Reader.read(self.params['scans'], self.band_ind, + force_tuple=True) + if params['time_block'] == 'scan': + for Data in Blocks: + yield (Data,), middle + elif params['time_block'] == 'file': + yield Blocks, middle + else: + msg = "time_block parameter must be 'scan' or 'file'." + raise dm.ValueError(msg) + + + def execute(self, n_processes): + """Driver method.""" + + # n_processes matters when using mpirun from MPI.COMM_WORLD and threading + self.n_processes = n_processes + params = self.params + kiyopy.utils.mkparents(params['output_root']) + parse_ini.write_params(params, params['output_root'] + + 'dirty_map_params.ini', + prefix='dm_') + if self.feedback > 0: + print "Input root is: " + params["input_root"] + print "Output root is: " + params["output_root"] + # Set flag if we are doing independant channels. + if params['frequency_correlations'] == 'None': + self.uncorrelated_channels = True + if self.feedback > 0: + print "Treating frequency channels as being independant." + else: + self.uncorrelated_channels = False + # Open the file that stores noise parameters. + if params['noise_parameter_file']: + self.noise_params = shelve.open(params['noise_parameter_file'], 'r') + else: + self.noise_params = None + # If we are subtracting out foreground modes, open that file. + if params['ts_foreground_mode_file']: + self.foreground_modes = shelve.open( + params['ts_foreground_mode_file'], 'r') + else: + self.foreground_modes = None + # Set the map dimensioning parameters. + self.n_ra = params['map_shape'][0] + self.n_dec = params['map_shape'][1] + spacing = params["pixel_spacing"] + # Negative sign because RA increases from right to left. + ra_spacing = -spacing/sp.cos(params['field_centre'][1]*sp.pi/180.) + # To set some parameters, we have to read the first data file. + first_file_name = (params["input_root"] + params["file_middles"][0] + + params["input_end"]) + Reader = fitsGBT.Reader(first_file_name, feedback=self.feedback) + Blocks = Reader.read(0, (), force_tuple=True) + self.n_chan = Blocks[0].dims[3] + # Figure out which polarization indices to process. + n_pols = Blocks[0].dims[1] + pols = list(Blocks[0].field["CRVAL4"]) + self.pols = pols + if not params['polarizations']: + pol_inds = range(n_pols) + else: + pol_inds = [] + for stokes_par in params['polarizations']: + for ii in range(n_pols): + if stokes_par == utils.polint2str(pols[ii]): + pol_inds.append(ii) + break + else : + msg = "Polarization " + stokes_par + " not in Data files." + raise ce.ValueError(msg) + # Get the centre frequncies of all the bands and figure out which bands + # to process. + n_bands = len(Reader.IF_set) + band_inds = params["IFs"] + if not band_inds: + band_inds = range(n_bands) + delta_freq = Blocks[0].field['CDELT1'] + self.delta_freq = delta_freq + band_centres = [] + for Data in Blocks: + Data.calc_freq() + band_centres.append(Data.freq[self.n_chan//2]) + self.band_centres = band_centres + del Blocks + del Reader + map_shape = (self.n_chan, self.n_ra, self.n_dec) + for ii in pol_inds: + for jj in band_inds: + self.band_ind = jj + band_centre = band_centres[jj] + self.pol_ind = ii + pol = pols[ii] + if self.feedback > 1: + print ("Making map for polarization " + + utils.polint2str(pol) + " and band centred at " + +repr(round(band_centre/1e6)) + "MHz.") + # Work out the file names. + map_filename = (params["output_root"] + "dirty_map_" + + utils.polint2str(pol) + + "_" + str(int(round(band_centre/1e6))) + + '.npy') + cov_filename = (params["output_root"] + "noise_inv_" + + utils.polint2str(pol) + + "_" + str(int(round((band_centre/1e6)))) + + '.hdf5') + self.map_filename = map_filename + self.cov_filename = cov_filename + # Initialization of the outputs. + map = sp.zeros(map_shape, dtype=float) + map = al.make_vect(map, axis_names=('freq', 'ra', 'dec')) + map.set_axis_info('freq', band_centre, delta_freq) + map.set_axis_info('ra', params['field_centre'][0], ra_spacing) + map.set_axis_info('dec', params['field_centre'][1], + params['pixel_spacing']) + '''if self.uncorrelated_channels: + cov_inv = al.open_memmap(cov_filename, mode='w+', + shape=map_shape + map_shape[1:], + dtype=float) + cov_inv = al.make_mat(cov_inv, + axis_names=('freq', 'ra', 'dec', 'ra', 'dec'), + row_axes=(0, 1, 2), col_axes=(0, 3, 4)) + else: + cov_inv = al.open_memmap(cov_filename, mode='w+', + shape=map_shape + map_shape, + dtype=float) + cov_inv = al.make_mat(cov_inv, + axis_names=('freq', 'ra', 'dec', + 'freq', 'ra', 'dec'), + row_axes=(0, 1, 2), col_axes=(3, 4, 5)) + cov_inv.copy_axis_info(map)''' + # The zeroing takes too long. Do it in memory, not hard disk. + #print 'zeroing cov_inv' + #cov_inv[...] = 0 + #print 'zeroed' + self.map = map + #self.cov_inv = cov_inv + # Do work. + try: + print "self.make_map() now" + self.make_map() + except: + # I think yo need to do this to get a sensible error + # message, otherwise it will print cov_inv and fill the + # screen. + #del self.cov_inv, cov_inv + raise + # IO. + # To write the noise_inverse, all we have to do is delete the + # memeory map object. + #del self.cov_inv, cov_inv + al.save(map_filename, map) + if not self.noise_params is None: + self.noise_params.close() + + def make_map(self): + """Makes map for current polarization and band. + + This worker function has been split off for testing reasons. + """ + + params = self.params + map = self.map + #cov_inv = self.cov_inv + # The current polarization index and value. + pol_ind = self.pol_ind + pol_val = self.pols[pol_ind] + # The current band, index and value. + band_ind = self.band_ind + band_centre = self.band_centres[band_ind] + # This outer loop is over groups of files. This is so we don't need to + # hold the noise matrices for all the data in memory at once. + file_middles = params['file_middles'] + n_files = len(file_middles) + n_files_group = params['n_files_group'] + if n_files_group == 0: + n_files_group = n_files + for start_file_ind in range(0, n_files, n_files_group): + all_this_file_middles = file_middles[start_file_ind + :start_file_ind + n_files_group] + # Each process will handle an equal number of input files. + # (Numbers that do not divide are taken care of as evenly as + # possible). This way is most efficient in memory and efficiency. + # This also allows a smaller number of mpi calls when passing + # DataSets around later to fill cov_inv. + + #this_file_middles,junk = split_elems(all_this_file_middles, + #self.nproc) + #this_file_middles = this_file_middles[self.rank] + + this_file_middles = all_this_file_middles + + # Don't need the 'junk' from above right now. That was added + # at a later point for something else. + if self.feedback > 1: + # Will have to change for each process saying something + # unique. + print ("Processing data files %d to %d of %d." + % (start_file_ind, start_file_ind + n_files_group, + n_files)) + # Initialize lists for the data and the noise. + time_stream_list = [] + data_set_list = [] + # Loop an iterator that reads and preprocesses the data. + # This loop can't be threaded without a lot of restructuring, + # because iterate_data, preprocess_data and get_noise_parameter + # all talk to each other through class attributes. + for this_Data, middle in self.iterate_data(this_file_middles): + try: + this_DataSet = dm.DataSet(params, this_Data, map, + self.n_chan, self.delta_freq, self.pols, + self.pol_ind, self.band_centres, + self.band_ind, middle) + except dm.DataSetError: + continue + # Check that the polarization for this data is the correct + # one. + if this_DataSet.pol_val != pol_val: + raise RuntimeError("Data polarizations not consistant.") + # Making the Noise and adding the mask is done in the + # constructor for DataSet. + # + # The thermal part. + thermal_noise = this_DataSet.get_noise_parameter("thermal") + try: + this_DataSet.Noise.add_thermal(thermal_noise) + except dm.NoiseError: + continue + # Get the noise parameter database for this data. + if params['noise_parameter_file']: +# noise_entry = (self.noise_params[middle] + noise_entry = (this_DataSet.noise_params[middle] + [int(round(band_centre/1e6))][pol_val]) + # Should save this into the Data_Set so it doesn't have + # to be kept again there. + else: + noise_entry = None + + # Frequency correlations. + if params['frequency_correlations'] == 'mean': + mean_overf = this_DataSet.get_noise_parameter("mean_over_f") + this_DataSet.Noise.add_correlated_over_f(*mean_overf) + elif params['frequency_correlations'] == 'measured': + # Correlated channel mode noise. + # The first frequency modes are completly deweighted. This + # was seen to be nessisary from looking at the noise + # spectra. + n_modes_discard = params['number_frequency_modes_discard'] + try: + for ii in range(n_modes_discard): + mode_noise_params = this_DataSet.get_noise_parameter( + "over_f_mode_" + repr(ii)) + # Instead of completly discarding, multiply + # amplitude by 2. + #N.deweight_freq_mode(mode_noise_params[4]) + mode_noise_params = ((mode_noise_params[0] * 2,) + + mode_noise_params[1:]) + this_DataSet.Noise.add_over_f_freq_mode(*mode_noise_params) + # The remaining frequency modes are noise weighed + # normally. + n_modes = params['number_frequency_modes'] + for ii in range(n_modes_discard, n_modes): + mode_noise_params = this_DataSet.get_noise_parameter( + "over_f_mode_" + repr(ii)) + this_DataSet.Noise.add_over_f_freq_mode(*mode_noise_params) + # All channel low frequency noise. + all_chan_noise_params = this_DataSet.get_noise_parameter( + 'all_channel') + # In some cases the corner frequncy will be so high that + # this scan will contain no information. + this_DataSet.Noise.add_all_chan_low(*all_chan_noise_params) + except dm.NoiseError: + if self.feedback > 1: + print ("Noise parameters risk numerical" + " instability. Data skipped.") + continue + elif params['frequency_correlations'] == 'None': + pass + else: + raise dm.ValueError("Invalid frequency correlations.") + + +# # From input parameters decide what the noise model should be. +# if params['frequency_correlations'] == 'None': +# model = 'thermal_only' +# else: +# # This must be set to something so code doesn't say +# # it doesn't exist. +# model = None +# # TODO Other options. +# this_DataSet.setup_noise_from_model(model, +# params['deweight_time_mean'], +# params['deweight_time_slope'], noise_entry) + if params['deweight_time_mean']: + this_DataSet.Noise.deweight_time_mean(T_huge**2) + if params['deweight_time_slope']: + this_DataSet.Noise.deweight_time_slope(T_huge**2) + # Set the time stream foregrounds for this data. + if (params['ts_foreground_mode_file'] + and params['n_ts_foreground_modes']) : + foreground_modes_db = self.foreground_modes + key = ts_measure.get_key(file_middle) + # From the database, read the data for this file. + db_freq = foreground_modes_db[key + '.freq'] + db_vects = foreground_modes_db[key + '.vects'] + ts_foregrounds_entry = (db_freq, db_vects) + this_DataSet.setup_ts_foregrounds( + params['n_ts_foreground_modes'], + ts_foregrounds_entry) + # No need for time_stream_list since data in DataSet. + #time_stream_list.append(this_DataSet.time_stream) + data_set_list.append(this_DataSet) + if self.feedback > 1: + print "Finalizing %d noise blocks: " % len(data_set_list) + noise_queue = Queue() + def thread_work(): + while True: + thread_N = noise_queue.get() + #None will be the flag that there is no more work to do. + if thread_N is None: + return + else: + thread_N.Noise.finalize(frequency_correlations=( + not self.uncorrelated_channels), + preserve_matrices=False) + if self.feedback > 1: + sys.stdout.write('.') + sys.stdout.flush() + #Start the worker threads. + thread_list = [] + for ii in range(self.n_processes): + T = threading.Thread(target=thread_work) + T.start() + thread_list.append(T) + #Now put work on the queue for the threads to do. + for a_DataSet in data_set_list: + noise_queue.put(a_DataSet) + #At the end of the queue, tell the threads that they are done. + for ii in range(self.n_processes): + noise_queue.put(None) + #Wait for the threads. + for T in thread_list: + T.join() + if not noise_queue.empty(): + msg = "A thread had an error in Noise finalization." + raise RuntimeError(msg) + + '''# Finalize the noises. Each process has a subset of the total + # noises, so no threading, just a loop. + for a_DataSet in data_set_list: + a_DataSet.Noise.finalize(frequency_correlations= + (not self.uncorrelated_channels), + preserve_matrices=False) + if self.feedback > 1: + sys.stdout.write('.') + sys.stdout.flush()''' + + if self.feedback > 1: + print + print "Noise finalized." + + #temporary code to be sure all noise has time_mode_update attribute + for data in data_set_list: + time_mode=data.Noise.time_mode_update + + # Each process holds a subset of the DataSets. They will be + # passed around "in a circle" with mpi so that each process + # sees all the sets. Right now, find the "previous" and "next" + # process that will make the above statement true. + # Note: This ordering will make others see the package + # in "numerical" order. + if ((self.rank+1) == self.nproc): + self.prev_guy = 0 + else: + self.prev_guy = self.rank + 1 + + if (self.rank == 0): + self.next_guy = self.nproc - 1 + else: + self.next_guy = self.rank - 1 + + # Each process will handle an equal subset of indices of + # the cov_inv to fill. + # This does not have to be rediscovered for every DataSet list + # passed in, so get the indices first. + chan_index_list = range(self.n_chan) + dtype=np.float64 + dsize=MPI.DOUBLE.Get_size() + if self.uncorrelated_channels: + # Using self.rank after gets the info for + # the approproate process. + # 'junk' for now. + index_list,junk = split_elems(chan_index_list,self.nproc) + index_list = index_list[self.rank] + # Allocate hdf5 file + if self.rank == 0: + print self.n_chan + print self.n_ra + print self.n_dec + data_offset, file_size = allocate_hdf5_dataset(self.cov_filename, 'inv_cov', (self.n_chan, self.n_ra, self.n_dec, self.n_ra, self.n_dec), dtype) + else: + ra_index_list = range(self.n_ra) + # cross product = A x B = (a,b) for all a in A, for all b in B + chan_ra_index_list = clacross([chan_index_list,ra_index_list]) + index_list,start_list = split_elems(chan_ra_index_list, + self.nproc) + index_list = index_list[self.rank] + # This is the point in n_chan x n_ra that the first index + # is at. This is needed for knowing which part of the + # total matrix a process has. + f_ra_start_ind = start_list[self.rank] + # Allocate hdf5 file + data_offset, file_size = allocate_hdf5_dataset(self.cov_filename, 'inv_cov', (self.n_chan, self.n_ra, self.n_dec, self.n_chan, self.n_ra, self.n_dec), dtype) + # Wait for file to be written before continuing. + comm.Barrier() + print '\n' + 'Process ' + str(self.rank) + ' Passed first barrier.' + '\n' + # Since the DataSets are split evenly over the processes, + # have to put in the pointing/noise at each index for + # each DataSet list that the processes hold. + + #for run in range(self.nproc): + # The commented loop above is for passing subsections of the data + # between nodes. It currently doesn't work due to pickling issues. + # I replaced it with a loop of just one pass to avoid re-indenting. + for run in range(1): + print '\n' + 'Process' + ' ' + str(self.rank) + ' ' + ',run' + ' ' + str(run) + '\n' + # Each DataSet has to be applied to the dirty map, too. + # Since the '_mpi'dirty map is much smaller than the noise, + # We will just let processor 0 do all the work + # for the dirty map. + # Must be done here to not pass the DataSets around twice. + if self.rank == 0: + for ii in xrange(len(data_set_list)): + print '\n' + 'Process ' + str(self.rank) + ' adding to map data piece' + ' ' + str(ii) + ' of ' + str(len(data_set_list)) +' during run ' + str(run) + '\n' + #time_stream = time_stream_list[ii] + time_stream = data_set_list[ii].time_stream + N = data_set_list[ii].Noise + P = data_set_list[ii].Pointing + w_time_stream = N.weight_time_stream(time_stream) + map += P.apply_to_time_axis(w_time_stream) + + #If first pass, prepare thread_cov_inv_chunk to be worked on by threads + if run == 0 and start_file_ind == 0: + if self.uncorrelated_channels: + thread_cov_inv_chunk = sp.zeros((len(index_list), self.n_ra, self.n_dec, self.n_ra, self.n_dec), dtype=float) + else: + thread_cov_inv_chunk = sp.zeros((len(index_list), self.n_dec, self.n_chan, self.n_ra, self.n_dec), dtype=float) + index_queue = Queue() + def thread_work(): + while True: + thread_inds = index_queue.get() + # None will be the flag that there is no more work to do. + if thread_inds is None: + return + if self.uncorrelated_channels: + thread_f_ind = thread_inds[0] + ra_vals = thread_inds[1] + ra_min = ra_vals[0] + ra_cap = ra_vals[-1]+1 + ra_size = ra_cap - ra_min + thread_cov_inv_block = sp.zeros((ra_size, self.n_dec, self.n_ra, self.n_dec), dtype=float) + for thread_D in data_set_list: + thread_D.Pointing.noise_channel_to_map(thread_D.Noise, thread_f_ind + index_list[0], thread_cov_inv_block, ra0_range = [ra_min, ra_cap]) + if run == 0 and start_file_ind == 0 and ra_min == 0: + #Adding prior to the diagonal. + thread_cov_inv_chunk[thread_f_ind,...].flat[::self.n_ra * self.n_dec + 1] += \ + 1.0 / T_large**2 + #No lock below currently. I don't think it is necessary without the memmap. + thread_cov_inv_chunk[thread_f_ind,ra_min:ra_cap,...] += thread_cov_inv_block + else: + ii = thread_inds + thread_cov_inv_row = sp.zeros((self.n_dec, self.n_chan, self.n_ra, self.n_dec),dtype=float) + thread_f_ind = index_list[ii][0] + thread_ra_ind = index_list[ii][1] + for thread_D in data_set_list: + thread_D.Pointing.noise_to_map_domain(thread_D.Noise, thread_f_ind, thread_ra_ind, thread_cov_inv_row) + if run == 0 and start_file_ind == 0: + #Adding prior to the diagonal. + thread_cov_inv_row_adder = thread_cov_inv_row[:, thread_f_ind, thread_ra_ind, :] + thread_cov_inv_row_adder.flat[:: self.n_dec + 1] += \ + 1.0 / T_large**2 + if (self.feedback > 1 and thread_ra_ind == self.n_ra - 1): + print thread_f_ind, + sys.stdout.flush() + thread_cov_inv_chunk[ii,...] += \ + thread_cov_inv_row + #Now start the worker threads. + thread_list = [] + for ii in range(self.n_processes): + T = threading.Thread(target=thread_work) + T.start() + thread_list.append(T) + print '\n' + 'Process ' + str(self.rank) + ' is using ' + str(self.n_processes) + 'threads.' + '\n' + #Now put work on the queue for the threads to do. + for ii in xrange(len(index_list)): + if self.uncorrelated_channels: + split_ra = split_elems(xrange(self.n_ra),self.n_processes)[0] + for split in split_ra: + index_queue.put((ii, split)) + else: + index_queue.put(ii) + for ii in range(self.n_processes): + index_queue.put(None) + for T in thread_list: + T.join() + if not index_queue.empty(): + msg = "A thread had an error while building map covariance." + raise RuntimeError(msg) + + # Commented section below is without threading. + '''if self.uncorrelated_channels: + # No actual threading going on. + if run == 0 and start_file_ind == 0: + thread_cov_inv_chunk = sp.zeros((len(index_list), + self.n_ra, self.n_dec, + self.n_ra, self.n_dec), dtype=float) + #for thread_f_ind in index_list: + print '\n' + 'Process ' + str(self.rank) + ' just created inv_chunk' + ' during run ' + str(run) + '\n' + for thread_f_ind in index_list: + #for ii in xrange(len(index_list)): + thread_cov_inv_block = sp.zeros((self.n_ra, self.n_dec, + self.n_ra, self.n_dec), dtype=float) + #if run == 0 and start_file_ind == 0: + # thread_cov_inv_block.flat[::self.n_ra * self.n_dec + 1] += \ + # 1.0 / T_large**2 + + for thread_D in data_set_list: + thread_D.Pointing.noise_channel_to_map( + thread_D.Noise, thread_f_ind, thread_cov_inv_block) + + if run == 0 and start_file_ind == 0: + thread_cov_inv_block.flat[::self.n_ra * self.n_dec + 1] += \ + 1.0 / T_large**2 + + #if start_file_ind == 0: + # The first time through the matrix. We 'zero' + # everything by just assigning a value instead of + # setting to zero then += value. + # TODO: writing to same place might hiccup. + #cov_inv[thread_f_ind,...] = thread_cov_inv_block + # Add stuff to diagonal now since it's not + # very fun later. + #thread_cov_inv_block.flat \ + # [::self.n_ra * self.n_dec + 1] += \ + # 1.0 / T_large**2 + #else: + # Not the first time through. Just add in vals. + #cov_inv[thread_f_ind,...] += thread_cov_inv_block + #thread_cov_inv_chunk[thread_f_ind,...] += thread_cov_inv_block + thread_cov_inv_chunk[thread_f_ind-index_list[0],...] += thread_cov_inv_block + print '\n' + 'Process ' + str(self.rank) + ' added to cov_inv block at freq ' + str(thread_f_ind) + ' , using data from pass ' + str(run) + ' of ' + str(self.nproc) + '\n' + if self.feedback > 1: + print thread_f_ind, + sys.stdout.flush() + else: + # Keep all of a process' rows in memory, then use + # MPI and file views to write it out fast. + if run == 0 and start_file_ind == 0: + thread_cov_inv_chunk = sp.zeros((len(index_list), + self.n_dec, self.n_chan, + self.n_ra, self.n_dec), + dtype=float) + #for thread_f_ind,thread_ra_ind in index_list: + for ii in xrange(len(index_list)): + thread_cov_inv_row = sp.zeros((self.n_dec, self.n_chan, + self.n_ra, self.n_dec), + dtype=float) + thread_f_ind = index_list[ii][0] + thread_ra_ind = index_list[ii][1] + if run == 0 and start_file_ind == 0: + thread_cov_inv_row_adder = thread_cov_inv_row[:, thread_f_ind, thread_ra_ind, :] + thread_cov_inv_row_adder.flat[:: self.n_dec + 1] += \ + 1.0 / T_large**2 + for thread_D in data_set_list: + thread_D.Pointing.noise_to_map_domain( + thread_D.Noise, thread_f_ind, + thread_ra_ind, thread_cov_inv_row) + + + #writeout_filename = self.cov_filename \ + # + repr(thread_inds).zfill(20) + #f = open(writeout_filename, 'w') + #np.save(f,thread_cov_inv_row) + #f.close() + if (self.feedback > 1 + and thread_ra_ind == self.n_ra - 1): + print thread_f_ind, + sys.stdout.flush() + thread_cov_inv_chunk[ii,...] += \ + thread_cov_inv_row + # Once done with one list of DataSets, do next. + # Note: No need to pass anything the last time since + # that data was the process' original data and + # has already been processed. + ''' + #if (run != (self.nproc - 1)): + + #Replace if statement below with if statement above to pass data between nodes. + if run == 1: + # Send out the DataSet list I have to next guy. + print '\n' + 'Process ' + str(self.rank) + ' is passing data_set_list to process ' + str(self.next_guy) + ', at the end of run ' + str(run) + '\n' + #Pickling DataSet objects in data_set_list + def pickle_list(list): + pickled_list = [] + for item in list: + pickled_list.append(cPickle.dumps(item)) + return pickled_list + data_set_list_pickled = pickle_list(data_set_list) + comm.send(data_set_list_pickled,dest=self.next_guy) + # Receive the DataSet list being sent to me by prev guy. + data_set_list_pickled = comm.recv(source=self.prev_guy) + #Unpickling DataSet objects + def unpickle_list(list): + unpickled_list = [] + for item in list: + unpickled_list.append(cPickle.loads(item)) + return unpickled_list + data_set_list=unpickle_list(data_set_list_pickled) + ## Same for time_stream_list. + #comm.send(time_stream_list,dest=self.next_guy) + #time_stream_list = comm.recv(source=self.prev_guy) + + # Processes will be relatively in sync here since they + # will block on read until the "previous" process + # sends over its data. + + # Now that thread_cov_inv_chunk is all filled with + # data from all DataSets, write it out. + # Only dealing with correlated channels now. + if not self.uncorrelated_channels: + #total_shape = (self.n_chan*self.n_ra, self.n_dec, + # self.n_chan, self.n_ra, self.n_dec) + #start_ind = (f_ra_start_ind,0,0,0,0) + lock_and_write_buffer(thread_cov_inv_chunk, self.cov_filename, data_offset + dsize*f_ra_start_ind*self.n_dec*self.n_chan*self.n_ra*self.n_dec, dsize*thread_cov_inv_chunk.size) + f = h5py.File(self.cov_filename, 'r+') + cov_inv_dset = f['inv_cov'] + cov_inv_shape = sp.zeros(shape = (self.n_chan, self.n_ra, self.n_dec, self,n_chan, self.n_ra, self.n_dec), dtype = dtype) + cov_inv_info = al.make_mat(cov_inv_shape, axis_names=('freq', 'ra', 'dec','freq', 'ra', 'dec'), row_axes=(0, 1, 2), col_axes=(3, 4, 5)) + cov_inv_info.copy_axis_info(map) + for key, value in cov_inv_info.info.iteritems(): + cov_inv_dset.attrs[key] = repr(value) + # NOTE: using 'float' is not supprted in the saving because + # it has to know if it is 32 or 64 bits. + #dtype = thread_cov_inv_chunk.dtype + #dtype = np.float64 # default for now + #np.save(self.cov_filename+'_mpi_corr_proc'+str(self.rank),thread_cov_inv_chunk) + # Save array. + #mpi_writearray(self.cov_filename+'_mpi', thread_cov_inv_chunk, + #comm, total_shape, start_ind, dtype, + #order='C', displacement=0) + + + if self.uncorrelated_channels: + #total_shape = (self.n_chan, self.n_ra, + # self.n_dec, self.n_ra, self.n_dec) + #start_ind = (index_list[0],0,0,0,0) + lock_and_write_buffer(thread_cov_inv_chunk, self.cov_filename, data_offset + dsize*index_list[0]*self.n_dec*self.n_ra*self.n_dec*self.n_ra, dsize*thread_cov_inv_chunk.size) + f = h5py.File(self.cov_filename, 'r+') + cov_inv_dset = f['inv_cov'] + #cov_inv = al.make_mat(cov_inv_dset, axis_names=('freq', 'ra', 'dec', 'ra', 'dec'), row_axes=(0, 1, 2), col_axes=(0, 3, 4)) + cov_inv_shape = sp.zeros(shape = (self.n_chan, self.n_ra, self.n_dec, self.n_ra, self.n_dec), dtype = dtype) + cov_inv_info = al.make_mat(cov_inv_shape, axis_names=('freq', 'ra', 'dec', 'ra', 'dec'), row_axes=(0, 1, 2), col_axes=(0, 3, 4)) + cov_inv_info.copy_axis_info(map) + for key, value in cov_inv_info.info.iteritems(): + cov_inv_dset.attrs[key] = repr(value) + # NOTE: using 'float' is not supprted in the saving because + # it has to know if it is 32 or 64 bits. + #dtype = thread_cov_inv_chunk.dtype + #dtype = np.float64 # default for now + #np.save(self.cov_filename+'_mpi_uncorr_proc'+str(self.rank),thread_cov_inv_chunk) + # Save array. + #mpi_writearray(self.cov_filename+'_mpi', thread_cov_inv_chunk, + #comm, total_shape, start_ind, dtype, #order='C', displacement=0) + + + +# Close self.noise_params in each DataSet. +# Along those lines, might want to open and close self.foreground_modes, too. + #for a_data in data_set_list: + # if not a_data.noise_params is None: + # a_data.noise_params.close() + # if not a_data.foreground_modes is None: + # a_data.foreground_modes.close( + + + + + + + +#### Utilities #### + +def split_elems(points, num_splits): + '''Split a list of arbitrary type elements, points, into num_splits sets. + Also, returns a list of the index that each set starts at. + Only tested with 'points' being: + list of str + list of int + list of (tuple of int)''' + + size = len(points) + start_list = [] + big_list = [] + for i in range(num_splits): + big_list.append([]) + # Number of elements per split/process. + divdiv = size / num_splits + # Number of processes that get one extra element. + modmod = size % num_splits + # Put in the right number of elements for each process. The 'extra' + # elements from mod get put into the first processes, so the difference + # in the number of elements per process is either 1 or 0. + start = 0 + end = divdiv + for i in range(num_splits): + if i < modmod: + end += 1 + for j in range(start,end): + if type(points[j]) == str: + big_list[i].append(points[j]) + else: + try: + big_list[i].append(tuple(points[j])) + except TypeError: + big_list[i].append(points[j]) + start_list.append(start) + start = end + end += divdiv + return big_list, start_list + + +def cross(set_list): + '''Given a list of sets, return the cross product.''' + # By associativity of cross product, cross the first two sets together + # then cross that with the rest. The big conditional in the list + # comprehension is just to make sure that there are no nested lists + # in the final answer. + if len(set_list) == 1: + ans = [] + for elem in set_list[0]: + # In the 1D case, these elements are not lists. + if type(elem) == list: + ans.append(np.array(elem)) + else: + ans.append(np.array([elem])) + return ans + else: + A = set_list[0] + B = set_list[1] + cross_2 = [a+b if ((type(a) == list) and (type(b) == list)) else \ + (a+[b] if type(a) == list else ([a]+b if type(b) == list else \ + [a]+[b])) for a in A for b in B] + remaining = set_list[2:] + remaining.insert(0,cross_2) + return cross(remaining) + +def lock_and_write_buffer(obj, fname, offset, size): + """Write the contents of a buffer to disk at a given offset, and explicitly + lock the region of the file whilst doing so. + + Parameters + ---------- + obj : buffer + Data to write to disk. + fname : string + Filename to write. + offset : integer + Offset into the file to start writing at. + size : integer + Size of the region to write to (and lock). + """ + import os + #import os.fcntl as fcntl + import fcntl + import h5py + + buf = buffer(obj) + + if len(buf) > size: + raise Exception("Size doesn't match array length.") + + fd = os.open(fname, os.O_RDWR | os.O_CREAT) + #fd = open(fname, 'rw') + #fd = h5py.File(fname, 'r+') + + try: + fcntl.lockf(fd, fcntl.LOCK_EX, size, offset, os.SEEK_SET) + except: + print "Could not obtain lock" + + os.lseek(fd, offset, 0) + + nb = os.write(fd, buf) + + if nb != len(buf): + raise Exception("Something funny happened with the reading.") + + try: + fcntl.lockf(fd, fcntl.LOCK_UN) + except: + print "Could not unlock" + + os.close(fd) + + +def allocate_hdf5_dataset(fname, dsetname, shape, dtype, comm=MPI.COMM_WORLD): + """Create a hdf5 dataset and return its offset and size. + + The dataset will be created contiguously and immediately allocated, + however it will not be filled. + + Parameters + ---------- + fname : string + Name of the file to write. + dsetname : string + Name of the dataset to write (must be at root level). + shape : tuple + Shape of the dataset. + dtype : numpy datatype + Type of the dataset. + comm : MPI communicator + Communicator over which to broadcast results. + + Returns + ------- + offset : integer + Offset into the file at which the dataset starts (in bytes). + size : integer + Size of the dataset in bytes. + + """ + + import h5py + + state = None + + if comm.rank == 0: + + # Create/open file + f = h5py.File(fname, 'a') + + # Create dataspace and HDF5 datatype + sp = h5py.h5s.create_simple(shape, shape) + tp = h5py.h5t.py_create(dtype) + + # Create a new plist and tell it to allocate the space for dataset + # immediately, but don't fill the file with zeros. + plist = h5py.h5p.create(h5py.h5p.DATASET_CREATE) + plist.set_alloc_time(h5py.h5d.ALLOC_TIME_EARLY) + plist.set_fill_time(h5py.h5d.FILL_TIME_NEVER) + + # Create the dataset + dset = h5py.h5d.create(f.id, dsetname, tp, sp, plist) + + # Get the offset of the dataset into the file. + state = dset.get_offset(), dset.get_storage_size() + + f.close() + + state = comm.bcast(state, root=0) + + return state + + + + +# If this file is run from the command line, execute the main function. +if __name__ == "__main__": + import sys + if len(sys.argv) == 2: + par_file = sys.argv[1] + nproc = 1 + DirtyMapMaker(par_file).execute(nproc) + elif len(sys.argv) == 3: + par_file = sys.argv[1] + nproc = int(sys.argv[2]) + DirtyMapMaker(par_file).execute(nproc) + else: + print "Usage: `python map/dirty_map.py parameter_file [num_threads]" + From 1543b474f6573ea00ebca6fe7f3912bc85a1a70d Mon Sep 17 00:00:00 2001 From: Chris Anderson Date: Thu, 6 Mar 2014 16:12:34 -0500 Subject: [PATCH 03/28] parallel_dirty_map.py has option for greater threading. Each thread acts on views of node's covariance chunk. --- input/km/mm_test.ini | 6 ++-- map/parallel_dirty_map.py | 68 ++++++++++++++++++++++++++------------- 2 files changed, 50 insertions(+), 24 deletions(-) diff --git a/input/km/mm_test.ini b/input/km/mm_test.ini index 39c227e8..0fa7ea72 100644 --- a/input/km/mm_test.ini +++ b/input/km/mm_test.ini @@ -31,12 +31,14 @@ map_spacing = .0627 raid_pro = os.getenv("RAID_PRO") input_data_dir = raid_pro + "kiyo/gbt_out_new/" base_dir = os.getenv("GBT_OUT") -map_dir = base_dir + 'maps/first_lockandwrite_test/tcs_threading/' +map_dir = base_dir + 'maps/first_lockandwrite_test/' +#map_dir = base_dir + 'maps/first_lockandwrite_test/tcs_threading/' #map_root = map_dir + "tmp_test_parallel_freqcorr_" map_root = map_dir + "tmp_test_parallel_inherit_" #map_root = map_dir + "tmp_test_regular_" # Map maker parameters. +dm_thread_divide = False dm_input_root = input_data_dir + 'reflagged_sec/' dm_file_middles = file_middles dm_input_end = '.fits' @@ -48,7 +50,7 @@ dm_field_centre = map_centre dm_pixel_spacing = map_spacing dm_map_shape = map_shape dm_time_block = 'scan' -dm_n_files_group = 280 # tpb nodes. +dm_n_files_group = 10 #dm_frequency_correlations = 'measured' dm_frequency_correlations = 'None' dm_number_frequency_modes = 3 diff --git a/map/parallel_dirty_map.py b/map/parallel_dirty_map.py index cb0ef811..ca8d5f09 100644 --- a/map/parallel_dirty_map.py +++ b/map/parallel_dirty_map.py @@ -86,7 +86,8 @@ # If there where any foregrounds subtracted in the time stream, # let the noise know about it. 'n_ts_foreground_modes' : 0, - 'ts_foreground_mode_file' : '' + 'ts_foreground_mode_file' : '', + 'thread_divide' : False } comm = MPI.COMM_WORLD @@ -515,10 +516,12 @@ def thread_work(): index_list = index_list[self.rank] # Allocate hdf5 file if self.rank == 0: + print 'rank zero' + '\n' print self.n_chan print self.n_ra print self.n_dec - data_offset, file_size = allocate_hdf5_dataset(self.cov_filename, 'inv_cov', (self.n_chan, self.n_ra, self.n_dec, self.n_ra, self.n_dec), dtype) + if start_file_ind == 0: + data_offset, file_size = allocate_hdf5_dataset(self.cov_filename, 'inv_cov', (self.n_chan, self.n_ra, self.n_dec, self.n_ra, self.n_dec), dtype) else: ra_index_list = range(self.n_ra) # cross product = A x B = (a,b) for all a in A, for all b in B @@ -531,7 +534,8 @@ def thread_work(): # total matrix a process has. f_ra_start_ind = start_list[self.rank] # Allocate hdf5 file - data_offset, file_size = allocate_hdf5_dataset(self.cov_filename, 'inv_cov', (self.n_chan, self.n_ra, self.n_dec, self.n_chan, self.n_ra, self.n_dec), dtype) + if start_file_ind == 0: + data_offset, file_size = allocate_hdf5_dataset(self.cov_filename, 'inv_cov', (self.n_chan, self.n_ra, self.n_dec, self.n_chan, self.n_ra, self.n_dec), dtype) # Wait for file to be written before continuing. comm.Barrier() print '\n' + 'Process ' + str(self.rank) + ' Passed first barrier.' + '\n' @@ -563,9 +567,13 @@ def thread_work(): #If first pass, prepare thread_cov_inv_chunk to be worked on by threads if run == 0 and start_file_ind == 0: if self.uncorrelated_channels: - thread_cov_inv_chunk = sp.zeros((len(index_list), self.n_ra, self.n_dec, self.n_ra, self.n_dec), dtype=float) + thread_cov_inv_chunk = sp.zeros((len(index_list), self.n_ra, self.n_dec, self.n_ra, self.n_dec), dtype=dtype) + #Adding prior to the diagonal. + for ii in xrange(len(index_list)): + thread_cov_inv_chunk[ii,...].flat[::self.n_ra * self.n_dec + 1] += \ + 1.0 / T_large**2 else: - thread_cov_inv_chunk = sp.zeros((len(index_list), self.n_dec, self.n_chan, self.n_ra, self.n_dec), dtype=float) + thread_cov_inv_chunk = sp.zeros((len(index_list), self.n_dec, self.n_chan, self.n_ra, self.n_dec), dtype=dtype) index_queue = Queue() def thread_work(): while True: @@ -574,33 +582,41 @@ def thread_work(): if thread_inds is None: return if self.uncorrelated_channels: - thread_f_ind = thread_inds - thread_cov_inv_block = sp.zeros((self.n_ra, self.n_dec, self.n_ra, self.n_dec), dtype=float) - for thread_D in data_set_list: - thread_D.Pointing.noise_channel_to_map(thread_D.Noise, thread_f_ind + index_list[0], thread_cov_inv_block) - if run == 0 and start_file_ind == 0: - #Adding prior to the diagonal. - thread_cov_inv_block.flat[::self.n_ra * self.n_dec + 1] += \ - 1.0 / T_large**2 - #No lock below currently. I don't think it is necessary without the memmap. - thread_cov_inv_chunk[thread_f_ind,...] += thread_cov_inv_block + if params['thread_divide'] == False: + thread_f_ind = thread_inds + #thread_cov_inv_block = sp.zeros((self.n_ra, self.n_dec, self.n_ra, self.n_dec), dtype=dtype) + for thread_D in data_set_list: + thread_D.Pointing.noise_channel_to_map(thread_D.Noise, thread_f_ind + index_list[0], thread_cov_inv_chunk[thread_f_ind,...]) + #No lock below currently. I don't think it is necessary without the memmap. + #thread_cov_inv_chunk[thread_f_ind,...] += thread_cov_inv_block + else: + thread_f_ind = thread_inds[0] + ra_vals = thread_inds[1] + ra_min = ra_vals[0] + ra_cap = ra_vals[-1]+1 + ra_size = ra_cap - ra_min + #thread_cov_inv_block = sp.zeros((ra_size, self.n_dec, self.n_ra, self.n_dec), dtype=float) + for thread_D in data_set_list: + thread_D.Pointing.noise_channel_to_map(thread_D.Noise, thread_f_ind + index_list[0], thread_cov_inv_chunk[thread_f_ind,ra_min:ra_cap,...], ra0_range = [ra_min, ra_cap]) + #No lock below currently. I don't think it is necessary without the memmap. + #thread_cov_inv_chunk[thread_f_ind,ra_min:ra_cap,...] += thread_cov_inv_block else: ii = thread_inds - thread_cov_inv_row = sp.zeros((self.n_dec, self.n_chan, self.n_ra, self.n_dec),dtype=float) + #thread_cov_inv_row = sp.zeros((self.n_dec, self.n_chan, self.n_ra, self.n_dec),dtype=dtype) thread_f_ind = index_list[ii][0] thread_ra_ind = index_list[ii][1] for thread_D in data_set_list: - thread_D.Pointing.noise_to_map_domain(thread_D.Noise, thread_f_ind, thread_ra_ind, thread_cov_inv_row) + thread_D.Pointing.noise_to_map_domain(thread_D.Noise, thread_f_ind, thread_ra_ind, thread_cov_inv_chunk[ii,...]) if run == 0 and start_file_ind == 0: #Adding prior to the diagonal. - thread_cov_inv_row_adder = thread_cov_inv_row[:, thread_f_ind, thread_ra_ind, :] - thread_cov_inv_row_adder.flat[:: self.n_dec + 1] += \ + #thread_cov_inv_row_adder = thread_cov_inv_row[:, thread_f_ind, thread_ra_ind, :] + thread_cov_inv_row_chunk[ii,:, thread_f_ind, thread_ra_ind,:].flat[:: self.n_dec + 1] += \ 1.0 / T_large**2 if (self.feedback > 1 and thread_ra_ind == self.n_ra - 1): print thread_f_ind, sys.stdout.flush() - thread_cov_inv_chunk[ii,...] += \ - thread_cov_inv_row + #thread_cov_inv_chunk[ii,...] += \ + # thread_cov_inv_row #Now start the worker threads. thread_list = [] for ii in range(self.n_processes): @@ -610,7 +626,15 @@ def thread_work(): print '\n' + 'Process ' + str(self.rank) + ' is using ' + str(self.n_processes) + 'threads.' + '\n' #Now put work on the queue for the threads to do. for ii in xrange(len(index_list)): - index_queue.put(ii) + if params['thread_divide'] == False: + index_queue.put(ii) + else: + if self.uncorrelated_channels: + split_ra = split_elems(xrange(self.n_ra),self.n_processes)[0] + for split in split_ra: + index_queue.put((ii, split)) + else: + index_queue.put(ii) for ii in range(self.n_processes): index_queue.put(None) for T in thread_list: From 39cde51170159a8be914580f5c91f182b9f25f75 Mon Sep 17 00:00:00 2001 From: Chris Anderson Date: Mon, 10 Mar 2014 14:53:16 -0400 Subject: [PATCH 04/28] Removed parallel_dirty_map_tcs.py. It's functionality is in parallel_dirty_map.py --- map/parallel_dirty_map_tcs.py | 1002 --------------------------------- 1 file changed, 1002 deletions(-) delete mode 100644 map/parallel_dirty_map_tcs.py diff --git a/map/parallel_dirty_map_tcs.py b/map/parallel_dirty_map_tcs.py deleted file mode 100644 index 5ee4cffb..00000000 --- a/map/parallel_dirty_map_tcs.py +++ /dev/null @@ -1,1002 +0,0 @@ -"""Dirty map making module. - -Module converts data in the time domain into noise weighted data in the map -domain, i.e. it creats the dirty map. Module also contains many utilities like -the pointing operator (`Pointing`) and the time domain noise operator -(`Noise`). -""" - -import sys - -from mpi4py import MPI - -import math -import threading -from Queue import Queue -import shelve -import sys -import time as time_mod -import cPickle -import h5py - -import scipy as sp -import numpy.ma as ma -import scipy.fftpack as fft -from scipy import linalg -from scipy import interpolate -import numpy as np -#import kiyopy.pickle_method -import warnings -#import matplotlib.pyplot as plt - -import core.algebra as al -from core import fitsGBT -import utils.misc as utils -import tools -from noise import noise_power -from foreground import ts_measure -import kiyopy.custom_exceptions as ce -import kiyopy.utils -from kiyopy import parse_ini -import _mapmaker as _mapmaker_c -from constants import T_infinity, T_huge, T_large, T_medium, T_small, T_sys -from constants import f_medium, f_large -from utils import misc -from map import dirty_map as dm - -prefix ='dm_' -params_init = {# IO: - 'input_root' : './', - # The unique part of every fname - 'file_middles' : ("testfile_GBTfits",), - 'input_end' : ".fits", - 'output_root' : "./testoutput_", - # Shelve file that holds measurements of the noise. - # What data to include from each file. - 'scans' : (), - 'IFs' : (0,), - 'polarizations' : ('I',), - # Map parameters (Ra (deg), Dec (deg)). - 'field_centre' : (325.0, 0.0), - # In pixels. - 'map_shape' : (5, 5), - 'pixel_spacing' : 0.5, # degrees - # Interpolation between pixel points. Options are 'nearest', - # 'linear' and 'cubic'. - 'interpolation' : 'linear', - # How to treat the data. - # How much data to include at a time (scan by scan or file by - # file) - 'time_block' : 'file', - # Number of files to process at a time. Increasing this number - # reduces IO but takes more memory. - 'n_files_group' : 0, - # What kind of frequency correlations to use. Options are - # 'None', 'mean' and 'measured'. - 'frequency_correlations' : 'None', - # Ignored unless 'frequency_correlations' is 'measured'. - 'number_frequency_modes' : 1, - # Of the above modes, how many to completely discard instead of - # noise weighting. - 'number_frequency_modes_discard' : 0, - # Where to get noise parameters from. - 'noise_parameter_file' : '', - 'deweight_time_mean' : True, - 'deweight_time_slope' : False, - # If there where any foregrounds subtracted in the time stream, - # let the noise know about it. - 'n_ts_foreground_modes' : 0, - 'ts_foreground_mode_file' : '' - } - -comm = MPI.COMM_WORLD - -class DirtyMapMaker(object): - """Dirty map maker. - """ - - def __init__(self, parameter_file_or_dict=None, feedback=2): - # Read in the parameters. - self.params = parse_ini.parse(parameter_file_or_dict, params_init, - prefix=prefix, feedback=feedback) - print self.params['file_middles'] - self.feedback = feedback - self.rank = comm.Get_rank() - self.nproc = comm.Get_size() - - def iterate_data(self, file_middles): - """An iterator over the input data. - - This can either iterate over the file names, processing the whole file - at once, or it can iterate over both the file and the scans. - - Returns the DataBlock objects as well as the file middle for that data - (which is used as a key in several parameter data bases). - """ - - params = self.params - for middle in file_middles: - self.file_middle = middle - fname = params["input_root"] + middle + params["input_end"] - Reader = fitsGBT.Reader(fname, feedback=self.feedback) - Blocks = Reader.read(self.params['scans'], self.band_ind, - force_tuple=True) - if params['time_block'] == 'scan': - for Data in Blocks: - yield (Data,), middle - elif params['time_block'] == 'file': - yield Blocks, middle - else: - msg = "time_block parameter must be 'scan' or 'file'." - raise dm.ValueError(msg) - - - def execute(self, n_processes): - """Driver method.""" - - # n_processes matters when using mpirun from MPI.COMM_WORLD and threading - self.n_processes = n_processes - params = self.params - kiyopy.utils.mkparents(params['output_root']) - parse_ini.write_params(params, params['output_root'] - + 'dirty_map_params.ini', - prefix='dm_') - if self.feedback > 0: - print "Input root is: " + params["input_root"] - print "Output root is: " + params["output_root"] - # Set flag if we are doing independant channels. - if params['frequency_correlations'] == 'None': - self.uncorrelated_channels = True - if self.feedback > 0: - print "Treating frequency channels as being independant." - else: - self.uncorrelated_channels = False - # Open the file that stores noise parameters. - if params['noise_parameter_file']: - self.noise_params = shelve.open(params['noise_parameter_file'], 'r') - else: - self.noise_params = None - # If we are subtracting out foreground modes, open that file. - if params['ts_foreground_mode_file']: - self.foreground_modes = shelve.open( - params['ts_foreground_mode_file'], 'r') - else: - self.foreground_modes = None - # Set the map dimensioning parameters. - self.n_ra = params['map_shape'][0] - self.n_dec = params['map_shape'][1] - spacing = params["pixel_spacing"] - # Negative sign because RA increases from right to left. - ra_spacing = -spacing/sp.cos(params['field_centre'][1]*sp.pi/180.) - # To set some parameters, we have to read the first data file. - first_file_name = (params["input_root"] + params["file_middles"][0] - + params["input_end"]) - Reader = fitsGBT.Reader(first_file_name, feedback=self.feedback) - Blocks = Reader.read(0, (), force_tuple=True) - self.n_chan = Blocks[0].dims[3] - # Figure out which polarization indices to process. - n_pols = Blocks[0].dims[1] - pols = list(Blocks[0].field["CRVAL4"]) - self.pols = pols - if not params['polarizations']: - pol_inds = range(n_pols) - else: - pol_inds = [] - for stokes_par in params['polarizations']: - for ii in range(n_pols): - if stokes_par == utils.polint2str(pols[ii]): - pol_inds.append(ii) - break - else : - msg = "Polarization " + stokes_par + " not in Data files." - raise ce.ValueError(msg) - # Get the centre frequncies of all the bands and figure out which bands - # to process. - n_bands = len(Reader.IF_set) - band_inds = params["IFs"] - if not band_inds: - band_inds = range(n_bands) - delta_freq = Blocks[0].field['CDELT1'] - self.delta_freq = delta_freq - band_centres = [] - for Data in Blocks: - Data.calc_freq() - band_centres.append(Data.freq[self.n_chan//2]) - self.band_centres = band_centres - del Blocks - del Reader - map_shape = (self.n_chan, self.n_ra, self.n_dec) - for ii in pol_inds: - for jj in band_inds: - self.band_ind = jj - band_centre = band_centres[jj] - self.pol_ind = ii - pol = pols[ii] - if self.feedback > 1: - print ("Making map for polarization " - + utils.polint2str(pol) + " and band centred at " - +repr(round(band_centre/1e6)) + "MHz.") - # Work out the file names. - map_filename = (params["output_root"] + "dirty_map_" - + utils.polint2str(pol) - + "_" + str(int(round(band_centre/1e6))) - + '.npy') - cov_filename = (params["output_root"] + "noise_inv_" - + utils.polint2str(pol) - + "_" + str(int(round((band_centre/1e6)))) - + '.hdf5') - self.map_filename = map_filename - self.cov_filename = cov_filename - # Initialization of the outputs. - map = sp.zeros(map_shape, dtype=float) - map = al.make_vect(map, axis_names=('freq', 'ra', 'dec')) - map.set_axis_info('freq', band_centre, delta_freq) - map.set_axis_info('ra', params['field_centre'][0], ra_spacing) - map.set_axis_info('dec', params['field_centre'][1], - params['pixel_spacing']) - '''if self.uncorrelated_channels: - cov_inv = al.open_memmap(cov_filename, mode='w+', - shape=map_shape + map_shape[1:], - dtype=float) - cov_inv = al.make_mat(cov_inv, - axis_names=('freq', 'ra', 'dec', 'ra', 'dec'), - row_axes=(0, 1, 2), col_axes=(0, 3, 4)) - else: - cov_inv = al.open_memmap(cov_filename, mode='w+', - shape=map_shape + map_shape, - dtype=float) - cov_inv = al.make_mat(cov_inv, - axis_names=('freq', 'ra', 'dec', - 'freq', 'ra', 'dec'), - row_axes=(0, 1, 2), col_axes=(3, 4, 5)) - cov_inv.copy_axis_info(map)''' - # The zeroing takes too long. Do it in memory, not hard disk. - #print 'zeroing cov_inv' - #cov_inv[...] = 0 - #print 'zeroed' - self.map = map - #self.cov_inv = cov_inv - # Do work. - try: - print "self.make_map() now" - self.make_map() - except: - # I think yo need to do this to get a sensible error - # message, otherwise it will print cov_inv and fill the - # screen. - #del self.cov_inv, cov_inv - raise - # IO. - # To write the noise_inverse, all we have to do is delete the - # memeory map object. - #del self.cov_inv, cov_inv - al.save(map_filename, map) - if not self.noise_params is None: - self.noise_params.close() - - def make_map(self): - """Makes map for current polarization and band. - - This worker function has been split off for testing reasons. - """ - - params = self.params - map = self.map - #cov_inv = self.cov_inv - # The current polarization index and value. - pol_ind = self.pol_ind - pol_val = self.pols[pol_ind] - # The current band, index and value. - band_ind = self.band_ind - band_centre = self.band_centres[band_ind] - # This outer loop is over groups of files. This is so we don't need to - # hold the noise matrices for all the data in memory at once. - file_middles = params['file_middles'] - n_files = len(file_middles) - n_files_group = params['n_files_group'] - if n_files_group == 0: - n_files_group = n_files - for start_file_ind in range(0, n_files, n_files_group): - all_this_file_middles = file_middles[start_file_ind - :start_file_ind + n_files_group] - # Each process will handle an equal number of input files. - # (Numbers that do not divide are taken care of as evenly as - # possible). This way is most efficient in memory and efficiency. - # This also allows a smaller number of mpi calls when passing - # DataSets around later to fill cov_inv. - - #this_file_middles,junk = split_elems(all_this_file_middles, - #self.nproc) - #this_file_middles = this_file_middles[self.rank] - - this_file_middles = all_this_file_middles - - # Don't need the 'junk' from above right now. That was added - # at a later point for something else. - if self.feedback > 1: - # Will have to change for each process saying something - # unique. - print ("Processing data files %d to %d of %d." - % (start_file_ind, start_file_ind + n_files_group, - n_files)) - # Initialize lists for the data and the noise. - time_stream_list = [] - data_set_list = [] - # Loop an iterator that reads and preprocesses the data. - # This loop can't be threaded without a lot of restructuring, - # because iterate_data, preprocess_data and get_noise_parameter - # all talk to each other through class attributes. - for this_Data, middle in self.iterate_data(this_file_middles): - try: - this_DataSet = dm.DataSet(params, this_Data, map, - self.n_chan, self.delta_freq, self.pols, - self.pol_ind, self.band_centres, - self.band_ind, middle) - except dm.DataSetError: - continue - # Check that the polarization for this data is the correct - # one. - if this_DataSet.pol_val != pol_val: - raise RuntimeError("Data polarizations not consistant.") - # Making the Noise and adding the mask is done in the - # constructor for DataSet. - # - # The thermal part. - thermal_noise = this_DataSet.get_noise_parameter("thermal") - try: - this_DataSet.Noise.add_thermal(thermal_noise) - except dm.NoiseError: - continue - # Get the noise parameter database for this data. - if params['noise_parameter_file']: -# noise_entry = (self.noise_params[middle] - noise_entry = (this_DataSet.noise_params[middle] - [int(round(band_centre/1e6))][pol_val]) - # Should save this into the Data_Set so it doesn't have - # to be kept again there. - else: - noise_entry = None - - # Frequency correlations. - if params['frequency_correlations'] == 'mean': - mean_overf = this_DataSet.get_noise_parameter("mean_over_f") - this_DataSet.Noise.add_correlated_over_f(*mean_overf) - elif params['frequency_correlations'] == 'measured': - # Correlated channel mode noise. - # The first frequency modes are completly deweighted. This - # was seen to be nessisary from looking at the noise - # spectra. - n_modes_discard = params['number_frequency_modes_discard'] - try: - for ii in range(n_modes_discard): - mode_noise_params = this_DataSet.get_noise_parameter( - "over_f_mode_" + repr(ii)) - # Instead of completly discarding, multiply - # amplitude by 2. - #N.deweight_freq_mode(mode_noise_params[4]) - mode_noise_params = ((mode_noise_params[0] * 2,) - + mode_noise_params[1:]) - this_DataSet.Noise.add_over_f_freq_mode(*mode_noise_params) - # The remaining frequency modes are noise weighed - # normally. - n_modes = params['number_frequency_modes'] - for ii in range(n_modes_discard, n_modes): - mode_noise_params = this_DataSet.get_noise_parameter( - "over_f_mode_" + repr(ii)) - this_DataSet.Noise.add_over_f_freq_mode(*mode_noise_params) - # All channel low frequency noise. - all_chan_noise_params = this_DataSet.get_noise_parameter( - 'all_channel') - # In some cases the corner frequncy will be so high that - # this scan will contain no information. - this_DataSet.Noise.add_all_chan_low(*all_chan_noise_params) - except dm.NoiseError: - if self.feedback > 1: - print ("Noise parameters risk numerical" - " instability. Data skipped.") - continue - elif params['frequency_correlations'] == 'None': - pass - else: - raise dm.ValueError("Invalid frequency correlations.") - - -# # From input parameters decide what the noise model should be. -# if params['frequency_correlations'] == 'None': -# model = 'thermal_only' -# else: -# # This must be set to something so code doesn't say -# # it doesn't exist. -# model = None -# # TODO Other options. -# this_DataSet.setup_noise_from_model(model, -# params['deweight_time_mean'], -# params['deweight_time_slope'], noise_entry) - if params['deweight_time_mean']: - this_DataSet.Noise.deweight_time_mean(T_huge**2) - if params['deweight_time_slope']: - this_DataSet.Noise.deweight_time_slope(T_huge**2) - # Set the time stream foregrounds for this data. - if (params['ts_foreground_mode_file'] - and params['n_ts_foreground_modes']) : - foreground_modes_db = self.foreground_modes - key = ts_measure.get_key(file_middle) - # From the database, read the data for this file. - db_freq = foreground_modes_db[key + '.freq'] - db_vects = foreground_modes_db[key + '.vects'] - ts_foregrounds_entry = (db_freq, db_vects) - this_DataSet.setup_ts_foregrounds( - params['n_ts_foreground_modes'], - ts_foregrounds_entry) - # No need for time_stream_list since data in DataSet. - #time_stream_list.append(this_DataSet.time_stream) - data_set_list.append(this_DataSet) - if self.feedback > 1: - print "Finalizing %d noise blocks: " % len(data_set_list) - noise_queue = Queue() - def thread_work(): - while True: - thread_N = noise_queue.get() - #None will be the flag that there is no more work to do. - if thread_N is None: - return - else: - thread_N.Noise.finalize(frequency_correlations=( - not self.uncorrelated_channels), - preserve_matrices=False) - if self.feedback > 1: - sys.stdout.write('.') - sys.stdout.flush() - #Start the worker threads. - thread_list = [] - for ii in range(self.n_processes): - T = threading.Thread(target=thread_work) - T.start() - thread_list.append(T) - #Now put work on the queue for the threads to do. - for a_DataSet in data_set_list: - noise_queue.put(a_DataSet) - #At the end of the queue, tell the threads that they are done. - for ii in range(self.n_processes): - noise_queue.put(None) - #Wait for the threads. - for T in thread_list: - T.join() - if not noise_queue.empty(): - msg = "A thread had an error in Noise finalization." - raise RuntimeError(msg) - - '''# Finalize the noises. Each process has a subset of the total - # noises, so no threading, just a loop. - for a_DataSet in data_set_list: - a_DataSet.Noise.finalize(frequency_correlations= - (not self.uncorrelated_channels), - preserve_matrices=False) - if self.feedback > 1: - sys.stdout.write('.') - sys.stdout.flush()''' - - if self.feedback > 1: - print - print "Noise finalized." - - #temporary code to be sure all noise has time_mode_update attribute - for data in data_set_list: - time_mode=data.Noise.time_mode_update - - # Each process holds a subset of the DataSets. They will be - # passed around "in a circle" with mpi so that each process - # sees all the sets. Right now, find the "previous" and "next" - # process that will make the above statement true. - # Note: This ordering will make others see the package - # in "numerical" order. - if ((self.rank+1) == self.nproc): - self.prev_guy = 0 - else: - self.prev_guy = self.rank + 1 - - if (self.rank == 0): - self.next_guy = self.nproc - 1 - else: - self.next_guy = self.rank - 1 - - # Each process will handle an equal subset of indices of - # the cov_inv to fill. - # This does not have to be rediscovered for every DataSet list - # passed in, so get the indices first. - chan_index_list = range(self.n_chan) - dtype=np.float64 - dsize=MPI.DOUBLE.Get_size() - if self.uncorrelated_channels: - # Using self.rank after gets the info for - # the approproate process. - # 'junk' for now. - index_list,junk = split_elems(chan_index_list,self.nproc) - index_list = index_list[self.rank] - # Allocate hdf5 file - if self.rank == 0: - print self.n_chan - print self.n_ra - print self.n_dec - data_offset, file_size = allocate_hdf5_dataset(self.cov_filename, 'inv_cov', (self.n_chan, self.n_ra, self.n_dec, self.n_ra, self.n_dec), dtype) - else: - ra_index_list = range(self.n_ra) - # cross product = A x B = (a,b) for all a in A, for all b in B - chan_ra_index_list = clacross([chan_index_list,ra_index_list]) - index_list,start_list = split_elems(chan_ra_index_list, - self.nproc) - index_list = index_list[self.rank] - # This is the point in n_chan x n_ra that the first index - # is at. This is needed for knowing which part of the - # total matrix a process has. - f_ra_start_ind = start_list[self.rank] - # Allocate hdf5 file - data_offset, file_size = allocate_hdf5_dataset(self.cov_filename, 'inv_cov', (self.n_chan, self.n_ra, self.n_dec, self.n_chan, self.n_ra, self.n_dec), dtype) - # Wait for file to be written before continuing. - comm.Barrier() - print '\n' + 'Process ' + str(self.rank) + ' Passed first barrier.' + '\n' - # Since the DataSets are split evenly over the processes, - # have to put in the pointing/noise at each index for - # each DataSet list that the processes hold. - - #for run in range(self.nproc): - # The commented loop above is for passing subsections of the data - # between nodes. It currently doesn't work due to pickling issues. - # I replaced it with a loop of just one pass to avoid re-indenting. - for run in range(1): - print '\n' + 'Process' + ' ' + str(self.rank) + ' ' + ',run' + ' ' + str(run) + '\n' - # Each DataSet has to be applied to the dirty map, too. - # Since the '_mpi'dirty map is much smaller than the noise, - # We will just let processor 0 do all the work - # for the dirty map. - # Must be done here to not pass the DataSets around twice. - if self.rank == 0: - for ii in xrange(len(data_set_list)): - print '\n' + 'Process ' + str(self.rank) + ' adding to map data piece' + ' ' + str(ii) + ' of ' + str(len(data_set_list)) +' during run ' + str(run) + '\n' - #time_stream = time_stream_list[ii] - time_stream = data_set_list[ii].time_stream - N = data_set_list[ii].Noise - P = data_set_list[ii].Pointing - w_time_stream = N.weight_time_stream(time_stream) - map += P.apply_to_time_axis(w_time_stream) - - #If first pass, prepare thread_cov_inv_chunk to be worked on by threads - if run == 0 and start_file_ind == 0: - if self.uncorrelated_channels: - thread_cov_inv_chunk = sp.zeros((len(index_list), self.n_ra, self.n_dec, self.n_ra, self.n_dec), dtype=float) - else: - thread_cov_inv_chunk = sp.zeros((len(index_list), self.n_dec, self.n_chan, self.n_ra, self.n_dec), dtype=float) - index_queue = Queue() - def thread_work(): - while True: - thread_inds = index_queue.get() - # None will be the flag that there is no more work to do. - if thread_inds is None: - return - if self.uncorrelated_channels: - thread_f_ind = thread_inds[0] - ra_vals = thread_inds[1] - ra_min = ra_vals[0] - ra_cap = ra_vals[-1]+1 - ra_size = ra_cap - ra_min - thread_cov_inv_block = sp.zeros((ra_size, self.n_dec, self.n_ra, self.n_dec), dtype=float) - for thread_D in data_set_list: - thread_D.Pointing.noise_channel_to_map(thread_D.Noise, thread_f_ind + index_list[0], thread_cov_inv_block, ra0_range = [ra_min, ra_cap]) - if run == 0 and start_file_ind == 0 and ra_min == 0: - #Adding prior to the diagonal. - thread_cov_inv_chunk[thread_f_ind,...].flat[::self.n_ra * self.n_dec + 1] += \ - 1.0 / T_large**2 - #No lock below currently. I don't think it is necessary without the memmap. - thread_cov_inv_chunk[thread_f_ind,ra_min:ra_cap,...] += thread_cov_inv_block - else: - ii = thread_inds - thread_cov_inv_row = sp.zeros((self.n_dec, self.n_chan, self.n_ra, self.n_dec),dtype=float) - thread_f_ind = index_list[ii][0] - thread_ra_ind = index_list[ii][1] - for thread_D in data_set_list: - thread_D.Pointing.noise_to_map_domain(thread_D.Noise, thread_f_ind, thread_ra_ind, thread_cov_inv_row) - if run == 0 and start_file_ind == 0: - #Adding prior to the diagonal. - thread_cov_inv_row_adder = thread_cov_inv_row[:, thread_f_ind, thread_ra_ind, :] - thread_cov_inv_row_adder.flat[:: self.n_dec + 1] += \ - 1.0 / T_large**2 - if (self.feedback > 1 and thread_ra_ind == self.n_ra - 1): - print thread_f_ind, - sys.stdout.flush() - thread_cov_inv_chunk[ii,...] += \ - thread_cov_inv_row - #Now start the worker threads. - thread_list = [] - for ii in range(self.n_processes): - T = threading.Thread(target=thread_work) - T.start() - thread_list.append(T) - print '\n' + 'Process ' + str(self.rank) + ' is using ' + str(self.n_processes) + 'threads.' + '\n' - #Now put work on the queue for the threads to do. - for ii in xrange(len(index_list)): - if self.uncorrelated_channels: - split_ra = split_elems(xrange(self.n_ra),self.n_processes)[0] - for split in split_ra: - index_queue.put((ii, split)) - else: - index_queue.put(ii) - for ii in range(self.n_processes): - index_queue.put(None) - for T in thread_list: - T.join() - if not index_queue.empty(): - msg = "A thread had an error while building map covariance." - raise RuntimeError(msg) - - # Commented section below is without threading. - '''if self.uncorrelated_channels: - # No actual threading going on. - if run == 0 and start_file_ind == 0: - thread_cov_inv_chunk = sp.zeros((len(index_list), - self.n_ra, self.n_dec, - self.n_ra, self.n_dec), dtype=float) - #for thread_f_ind in index_list: - print '\n' + 'Process ' + str(self.rank) + ' just created inv_chunk' + ' during run ' + str(run) + '\n' - for thread_f_ind in index_list: - #for ii in xrange(len(index_list)): - thread_cov_inv_block = sp.zeros((self.n_ra, self.n_dec, - self.n_ra, self.n_dec), dtype=float) - #if run == 0 and start_file_ind == 0: - # thread_cov_inv_block.flat[::self.n_ra * self.n_dec + 1] += \ - # 1.0 / T_large**2 - - for thread_D in data_set_list: - thread_D.Pointing.noise_channel_to_map( - thread_D.Noise, thread_f_ind, thread_cov_inv_block) - - if run == 0 and start_file_ind == 0: - thread_cov_inv_block.flat[::self.n_ra * self.n_dec + 1] += \ - 1.0 / T_large**2 - - #if start_file_ind == 0: - # The first time through the matrix. We 'zero' - # everything by just assigning a value instead of - # setting to zero then += value. - # TODO: writing to same place might hiccup. - #cov_inv[thread_f_ind,...] = thread_cov_inv_block - # Add stuff to diagonal now since it's not - # very fun later. - #thread_cov_inv_block.flat \ - # [::self.n_ra * self.n_dec + 1] += \ - # 1.0 / T_large**2 - #else: - # Not the first time through. Just add in vals. - #cov_inv[thread_f_ind,...] += thread_cov_inv_block - #thread_cov_inv_chunk[thread_f_ind,...] += thread_cov_inv_block - thread_cov_inv_chunk[thread_f_ind-index_list[0],...] += thread_cov_inv_block - print '\n' + 'Process ' + str(self.rank) + ' added to cov_inv block at freq ' + str(thread_f_ind) + ' , using data from pass ' + str(run) + ' of ' + str(self.nproc) + '\n' - if self.feedback > 1: - print thread_f_ind, - sys.stdout.flush() - else: - # Keep all of a process' rows in memory, then use - # MPI and file views to write it out fast. - if run == 0 and start_file_ind == 0: - thread_cov_inv_chunk = sp.zeros((len(index_list), - self.n_dec, self.n_chan, - self.n_ra, self.n_dec), - dtype=float) - #for thread_f_ind,thread_ra_ind in index_list: - for ii in xrange(len(index_list)): - thread_cov_inv_row = sp.zeros((self.n_dec, self.n_chan, - self.n_ra, self.n_dec), - dtype=float) - thread_f_ind = index_list[ii][0] - thread_ra_ind = index_list[ii][1] - if run == 0 and start_file_ind == 0: - thread_cov_inv_row_adder = thread_cov_inv_row[:, thread_f_ind, thread_ra_ind, :] - thread_cov_inv_row_adder.flat[:: self.n_dec + 1] += \ - 1.0 / T_large**2 - for thread_D in data_set_list: - thread_D.Pointing.noise_to_map_domain( - thread_D.Noise, thread_f_ind, - thread_ra_ind, thread_cov_inv_row) - - - #writeout_filename = self.cov_filename \ - # + repr(thread_inds).zfill(20) - #f = open(writeout_filename, 'w') - #np.save(f,thread_cov_inv_row) - #f.close() - if (self.feedback > 1 - and thread_ra_ind == self.n_ra - 1): - print thread_f_ind, - sys.stdout.flush() - thread_cov_inv_chunk[ii,...] += \ - thread_cov_inv_row - # Once done with one list of DataSets, do next. - # Note: No need to pass anything the last time since - # that data was the process' original data and - # has already been processed. - ''' - #if (run != (self.nproc - 1)): - - #Replace if statement below with if statement above to pass data between nodes. - if run == 1: - # Send out the DataSet list I have to next guy. - print '\n' + 'Process ' + str(self.rank) + ' is passing data_set_list to process ' + str(self.next_guy) + ', at the end of run ' + str(run) + '\n' - #Pickling DataSet objects in data_set_list - def pickle_list(list): - pickled_list = [] - for item in list: - pickled_list.append(cPickle.dumps(item)) - return pickled_list - data_set_list_pickled = pickle_list(data_set_list) - comm.send(data_set_list_pickled,dest=self.next_guy) - # Receive the DataSet list being sent to me by prev guy. - data_set_list_pickled = comm.recv(source=self.prev_guy) - #Unpickling DataSet objects - def unpickle_list(list): - unpickled_list = [] - for item in list: - unpickled_list.append(cPickle.loads(item)) - return unpickled_list - data_set_list=unpickle_list(data_set_list_pickled) - ## Same for time_stream_list. - #comm.send(time_stream_list,dest=self.next_guy) - #time_stream_list = comm.recv(source=self.prev_guy) - - # Processes will be relatively in sync here since they - # will block on read until the "previous" process - # sends over its data. - - # Now that thread_cov_inv_chunk is all filled with - # data from all DataSets, write it out. - # Only dealing with correlated channels now. - if not self.uncorrelated_channels: - #total_shape = (self.n_chan*self.n_ra, self.n_dec, - # self.n_chan, self.n_ra, self.n_dec) - #start_ind = (f_ra_start_ind,0,0,0,0) - lock_and_write_buffer(thread_cov_inv_chunk, self.cov_filename, data_offset + dsize*f_ra_start_ind*self.n_dec*self.n_chan*self.n_ra*self.n_dec, dsize*thread_cov_inv_chunk.size) - f = h5py.File(self.cov_filename, 'r+') - cov_inv_dset = f['inv_cov'] - cov_inv_shape = sp.zeros(shape = (self.n_chan, self.n_ra, self.n_dec, self,n_chan, self.n_ra, self.n_dec), dtype = dtype) - cov_inv_info = al.make_mat(cov_inv_shape, axis_names=('freq', 'ra', 'dec','freq', 'ra', 'dec'), row_axes=(0, 1, 2), col_axes=(3, 4, 5)) - cov_inv_info.copy_axis_info(map) - for key, value in cov_inv_info.info.iteritems(): - cov_inv_dset.attrs[key] = repr(value) - # NOTE: using 'float' is not supprted in the saving because - # it has to know if it is 32 or 64 bits. - #dtype = thread_cov_inv_chunk.dtype - #dtype = np.float64 # default for now - #np.save(self.cov_filename+'_mpi_corr_proc'+str(self.rank),thread_cov_inv_chunk) - # Save array. - #mpi_writearray(self.cov_filename+'_mpi', thread_cov_inv_chunk, - #comm, total_shape, start_ind, dtype, - #order='C', displacement=0) - - - if self.uncorrelated_channels: - #total_shape = (self.n_chan, self.n_ra, - # self.n_dec, self.n_ra, self.n_dec) - #start_ind = (index_list[0],0,0,0,0) - lock_and_write_buffer(thread_cov_inv_chunk, self.cov_filename, data_offset + dsize*index_list[0]*self.n_dec*self.n_ra*self.n_dec*self.n_ra, dsize*thread_cov_inv_chunk.size) - f = h5py.File(self.cov_filename, 'r+') - cov_inv_dset = f['inv_cov'] - #cov_inv = al.make_mat(cov_inv_dset, axis_names=('freq', 'ra', 'dec', 'ra', 'dec'), row_axes=(0, 1, 2), col_axes=(0, 3, 4)) - cov_inv_shape = sp.zeros(shape = (self.n_chan, self.n_ra, self.n_dec, self.n_ra, self.n_dec), dtype = dtype) - cov_inv_info = al.make_mat(cov_inv_shape, axis_names=('freq', 'ra', 'dec', 'ra', 'dec'), row_axes=(0, 1, 2), col_axes=(0, 3, 4)) - cov_inv_info.copy_axis_info(map) - for key, value in cov_inv_info.info.iteritems(): - cov_inv_dset.attrs[key] = repr(value) - # NOTE: using 'float' is not supprted in the saving because - # it has to know if it is 32 or 64 bits. - #dtype = thread_cov_inv_chunk.dtype - #dtype = np.float64 # default for now - #np.save(self.cov_filename+'_mpi_uncorr_proc'+str(self.rank),thread_cov_inv_chunk) - # Save array. - #mpi_writearray(self.cov_filename+'_mpi', thread_cov_inv_chunk, - #comm, total_shape, start_ind, dtype, #order='C', displacement=0) - - - -# Close self.noise_params in each DataSet. -# Along those lines, might want to open and close self.foreground_modes, too. - #for a_data in data_set_list: - # if not a_data.noise_params is None: - # a_data.noise_params.close() - # if not a_data.foreground_modes is None: - # a_data.foreground_modes.close( - - - - - - - -#### Utilities #### - -def split_elems(points, num_splits): - '''Split a list of arbitrary type elements, points, into num_splits sets. - Also, returns a list of the index that each set starts at. - Only tested with 'points' being: - list of str - list of int - list of (tuple of int)''' - - size = len(points) - start_list = [] - big_list = [] - for i in range(num_splits): - big_list.append([]) - # Number of elements per split/process. - divdiv = size / num_splits - # Number of processes that get one extra element. - modmod = size % num_splits - # Put in the right number of elements for each process. The 'extra' - # elements from mod get put into the first processes, so the difference - # in the number of elements per process is either 1 or 0. - start = 0 - end = divdiv - for i in range(num_splits): - if i < modmod: - end += 1 - for j in range(start,end): - if type(points[j]) == str: - big_list[i].append(points[j]) - else: - try: - big_list[i].append(tuple(points[j])) - except TypeError: - big_list[i].append(points[j]) - start_list.append(start) - start = end - end += divdiv - return big_list, start_list - - -def cross(set_list): - '''Given a list of sets, return the cross product.''' - # By associativity of cross product, cross the first two sets together - # then cross that with the rest. The big conditional in the list - # comprehension is just to make sure that there are no nested lists - # in the final answer. - if len(set_list) == 1: - ans = [] - for elem in set_list[0]: - # In the 1D case, these elements are not lists. - if type(elem) == list: - ans.append(np.array(elem)) - else: - ans.append(np.array([elem])) - return ans - else: - A = set_list[0] - B = set_list[1] - cross_2 = [a+b if ((type(a) == list) and (type(b) == list)) else \ - (a+[b] if type(a) == list else ([a]+b if type(b) == list else \ - [a]+[b])) for a in A for b in B] - remaining = set_list[2:] - remaining.insert(0,cross_2) - return cross(remaining) - -def lock_and_write_buffer(obj, fname, offset, size): - """Write the contents of a buffer to disk at a given offset, and explicitly - lock the region of the file whilst doing so. - - Parameters - ---------- - obj : buffer - Data to write to disk. - fname : string - Filename to write. - offset : integer - Offset into the file to start writing at. - size : integer - Size of the region to write to (and lock). - """ - import os - #import os.fcntl as fcntl - import fcntl - import h5py - - buf = buffer(obj) - - if len(buf) > size: - raise Exception("Size doesn't match array length.") - - fd = os.open(fname, os.O_RDWR | os.O_CREAT) - #fd = open(fname, 'rw') - #fd = h5py.File(fname, 'r+') - - try: - fcntl.lockf(fd, fcntl.LOCK_EX, size, offset, os.SEEK_SET) - except: - print "Could not obtain lock" - - os.lseek(fd, offset, 0) - - nb = os.write(fd, buf) - - if nb != len(buf): - raise Exception("Something funny happened with the reading.") - - try: - fcntl.lockf(fd, fcntl.LOCK_UN) - except: - print "Could not unlock" - - os.close(fd) - - -def allocate_hdf5_dataset(fname, dsetname, shape, dtype, comm=MPI.COMM_WORLD): - """Create a hdf5 dataset and return its offset and size. - - The dataset will be created contiguously and immediately allocated, - however it will not be filled. - - Parameters - ---------- - fname : string - Name of the file to write. - dsetname : string - Name of the dataset to write (must be at root level). - shape : tuple - Shape of the dataset. - dtype : numpy datatype - Type of the dataset. - comm : MPI communicator - Communicator over which to broadcast results. - - Returns - ------- - offset : integer - Offset into the file at which the dataset starts (in bytes). - size : integer - Size of the dataset in bytes. - - """ - - import h5py - - state = None - - if comm.rank == 0: - - # Create/open file - f = h5py.File(fname, 'a') - - # Create dataspace and HDF5 datatype - sp = h5py.h5s.create_simple(shape, shape) - tp = h5py.h5t.py_create(dtype) - - # Create a new plist and tell it to allocate the space for dataset - # immediately, but don't fill the file with zeros. - plist = h5py.h5p.create(h5py.h5p.DATASET_CREATE) - plist.set_alloc_time(h5py.h5d.ALLOC_TIME_EARLY) - plist.set_fill_time(h5py.h5d.FILL_TIME_NEVER) - - # Create the dataset - dset = h5py.h5d.create(f.id, dsetname, tp, sp, plist) - - # Get the offset of the dataset into the file. - state = dset.get_offset(), dset.get_storage_size() - - f.close() - - state = comm.bcast(state, root=0) - - return state - - - - -# If this file is run from the command line, execute the main function. -if __name__ == "__main__": - import sys - if len(sys.argv) == 2: - par_file = sys.argv[1] - nproc = 1 - DirtyMapMaker(par_file).execute(nproc) - elif len(sys.argv) == 3: - par_file = sys.argv[1] - nproc = int(sys.argv[2]) - DirtyMapMaker(par_file).execute(nproc) - else: - print "Usage: `python map/dirty_map.py parameter_file [num_threads]" - From 15ca1a39653996144e5c708fc0c1740e8cb3ad17 Mon Sep 17 00:00:00 2001 From: Tabitha Voytek Date: Mon, 24 Mar 2014 16:24:58 -0400 Subject: [PATCH 05/28] usage edits --- input/tcv/pointcorr_mapmaker.pipe | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/input/tcv/pointcorr_mapmaker.pipe b/input/tcv/pointcorr_mapmaker.pipe index 90333883..967f4f05 100644 --- a/input/tcv/pointcorr_mapmaker.pipe +++ b/input/tcv/pointcorr_mapmaker.pipe @@ -152,17 +152,17 @@ pipe_modules = [] #from time_stream import reflag #pipe_modules.append((reflag.ReFlag, ('sf2_', 'sf_'))) -from noise import measure_noise -pipe_modules.append((measure_noise.Measure, ('mn2_', 'mn_'))) +#from noise import measure_noise +#pipe_modules.append((measure_noise.Measure, ('mn2_', 'mn_'))) -from map import dirty_map -pipe_modules.append((dirty_map.DirtyMapMaker, ('dmA_', 'dm_'))) +#from map import dirty_map +#pipe_modules.append((dirty_map.DirtyMapMaker, ('dmA_', 'dm_'))) #pipe_modules.append((dirty_map.DirtyMapMaker, ('dmB_', 'dm_'))) #pipe_modules.append((dirty_map.DirtyMapMaker, ('dmC_', 'dm_'))) #pipe_modules.append((dirty_map.DirtyMapMaker, ('dmD_', 'dm_'))) -#from map import clean_map -#pipe_modules.append((clean_map.CleanMapMaker, ('cmA_', 'cm_'))) +from map import clean_map +pipe_modules.append((clean_map.CleanMapMaker, ('cmA_', 'cm_'))) #pipe_modules.append((clean_map.CleanMapMaker, ('cmB_', 'cm_'))) #pipe_modules.append((clean_map.CleanMapMaker, ('cmC_', 'cm_'))) #pipe_modules.append((clean_map.CleanMapMaker, ('cmD_', 'cm_'))) @@ -413,7 +413,7 @@ mn2_input_root = sf2_subtracted_output_root mn2_file_middles = new_file_middles mn2_input_end = '.fits' mn2_output_root = base_dir + 'noise_measurments_sec/' -mn2_output_filename = "noise_parameters.shelve" +mn2_output_filename = "noise_parameters_take2.shelve" #mn2_save_spectra_plots = True mn2_save_spectra_plots = False mn2_time_block = 'scan' From 6a2cc95aeb7d9244a38f2613e899e73f87c3ba65 Mon Sep 17 00:00:00 2001 From: Christopher Anderson Date: Tue, 15 Apr 2014 22:53:59 -0400 Subject: [PATCH 06/28] Trying to fix write to hdf5 file --- map/parallel_dirty_map.py | 31 +++++++++++++++++++++++++------ 1 file changed, 25 insertions(+), 6 deletions(-) diff --git a/map/parallel_dirty_map.py b/map/parallel_dirty_map.py index ca8d5f09..20bc77a5 100644 --- a/map/parallel_dirty_map.py +++ b/map/parallel_dirty_map.py @@ -767,7 +767,8 @@ def unpickle_list(list): #total_shape = (self.n_chan*self.n_ra, self.n_dec, # self.n_chan, self.n_ra, self.n_dec) #start_ind = (f_ra_start_ind,0,0,0,0) - lock_and_write_buffer(thread_cov_inv_chunk, self.cov_filename, data_offset + dsize*f_ra_start_ind*self.n_dec*self.n_chan*self.n_ra*self.n_dec, dsize*thread_cov_inv_chunk.size) + lock_and_write_buffer(thread_cov_inv_chunk, self.cov_filename, data_offset + dsize*f_ra_start_ind*self.n_dec*self.n_chan*self.n_ra*self.n_dec, dsize*thread_cov_inv_chunk.size, self.rank) + del thread_cov_inv_chunk f = h5py.File(self.cov_filename, 'r+') cov_inv_dset = f['inv_cov'] cov_inv_shape = sp.zeros(shape = (self.n_chan, self.n_ra, self.n_dec, self,n_chan, self.n_ra, self.n_dec), dtype = dtype) @@ -790,11 +791,13 @@ def unpickle_list(list): #total_shape = (self.n_chan, self.n_ra, # self.n_dec, self.n_ra, self.n_dec) #start_ind = (index_list[0],0,0,0,0) - lock_and_write_buffer(thread_cov_inv_chunk, self.cov_filename, data_offset + dsize*index_list[0]*self.n_dec*self.n_ra*self.n_dec*self.n_ra, dsize*thread_cov_inv_chunk.size) + lock_and_write_buffer(thread_cov_inv_chunk, self.cov_filename, data_offset + dsize*index_list[0]*self.n_dec*self.n_ra*self.n_dec*self.n_ra, dsize*thread_cov_inv_chunk.size, self.rank) + #del thread_cov_inv_chunk f = h5py.File(self.cov_filename, 'r+') cov_inv_dset = f['inv_cov'] #cov_inv = al.make_mat(cov_inv_dset, axis_names=('freq', 'ra', 'dec', 'ra', 'dec'), row_axes=(0, 1, 2), col_axes=(0, 3, 4)) - cov_inv_shape = sp.zeros(shape = (self.n_chan, self.n_ra, self.n_dec, self.n_ra, self.n_dec), dtype = dtype) + #cov_inv_shape = sp.zeros(shape = (self.n_chan, self.n_ra, self.n_dec, self.n_ra, self.n_dec), dtype = dtype) + #cov_inv_info = al.make_mat(cov_inv_shape, axis_names=('freq', 'ra', 'dec', 'ra', 'dec'), row_axes=(0, 1, 2), col_axes=(0, 3, 4)) cov_inv_info = al.make_mat(cov_inv_shape, axis_names=('freq', 'ra', 'dec', 'ra', 'dec'), row_axes=(0, 1, 2), col_axes=(0, 3, 4)) cov_inv_info.copy_axis_info(map) for key, value in cov_inv_info.info.iteritems(): @@ -890,7 +893,7 @@ def cross(set_list): remaining.insert(0,cross_2) return cross(remaining) -def lock_and_write_buffer(obj, fname, offset, size): +def lock_and_write_buffer(obj, fname, offset, size, proc): """Write the contents of a buffer to disk at a given offset, and explicitly lock the region of the file whilst doing so. @@ -925,12 +928,28 @@ def lock_and_write_buffer(obj, fname, offset, size): print "Could not obtain lock" os.lseek(fd, offset, 0) + + a=0 + while(True): + #while(len(buf)>0): + #nb = os.write(fd,buf) + nb = os.write(fd, buf[a:]) + print "The buffer save started at byte " + str(a) + " for process " + str(proc) + if nb < 0: + raise Exception("Failed write") + + if nb == len(buf[a:]): + break + else: + a += nb + #buf = buf[nb:] + ''' nb = os.write(fd, buf) if nb != len(buf): - raise Exception("Something funny happened with the reading.") - + #raise Exception("Something funny happened with the reading.") + raise Exception("Something funny happened with the reading." + " The buffer is length " + str(len(buf)) + " but the number of bytes written was " + str(nb) + " for process " + str(proc) + "." ) ''' try: fcntl.lockf(fd, fcntl.LOCK_UN) except: From d3d94ffac9310cb6e759ac2fdcff0187dd00fa4c Mon Sep 17 00:00:00 2001 From: Christopher Anderson Date: Thu, 17 Apr 2014 21:49:08 -0400 Subject: [PATCH 07/28] Trying to fix MemoryError. --- map/_mapmaker.pyx | 51 +++++++++++++++++++++++++++++++++++++++ map/parallel_dirty_map.py | 36 +++++++++++++-------------- 2 files changed, 69 insertions(+), 18 deletions(-) diff --git a/map/_mapmaker.pyx b/map/_mapmaker.pyx index 5c696393..4ca0f476 100644 --- a/map/_mapmaker.pyx +++ b/map/_mapmaker.pyx @@ -9,8 +9,13 @@ cimport cython # We will do everything in double precision. DTYPE = np.float +DTYPE_num = np.NPY_FLOAT64 ctypedef np.float_t DTYPE_t +cdef extern from "stdlib.h": + void free(void* ptr) + void* malloc(size_t size) + @cython.boundscheck(False) @cython.wraparound(False) @@ -131,6 +136,52 @@ def update_map_noise_chan_ra_row( tmp2 = tmp1 * pointing_weights[kk,pp] map_noise_inv[dec_ind,jj,this_ra,this_dec] += tmp2 +cdef class memory_container: + """Creates and holds a memory buffer. +Deallocates the memory when instance goes out of scope, even if other +objects are using the buffer. +""" + + cdef void * buffer + + def __cinit__(self, size): + self.buffer = malloc(size) + if not buffer: + raise MemoryError() + + cdef void * get_buffer(self): + return self.buffer + + def __dealloc__(self): + free(self.buffer) + +class buffer_array(np.ndarray): + _memory_handler = None + +def large_empty(shape): + + # Get all the dimensions in the right format. + cdef int nd = len(shape) + # Very important to use this type. + cdef np.npy_intp * dims = malloc(nd * sizeof(np.npy_intp)) + cdef long size = 1 + for ii, dim in enumerate(shape): + dims[ii] = dim + size *= dim + # Allocate all the memory we need in a buffer. + cdef memory_container mem + mem = memory_container(size * sizeof(DTYPE_t)) + # Get the buffer and convert it to a numpy array. + cdef void * buf = mem.get_buffer() + arr = np.PyArray_SimpleNewFromData(nd, dims, DTYPE_num, buf) + # Store the only living reference (once this function returns) to the + # memory container on the array. + arr = arr.view(buffer_array) + arr._memory_handler = mem + + free(dims) + return arr + def get_noise_inv_diag( np.ndarray[DTYPE_t, ndim=2, mode='c'] diagonal_inv not None, np.ndarray[DTYPE_t, ndim=2, mode='c'] freq_modes not None, diff --git a/map/parallel_dirty_map.py b/map/parallel_dirty_map.py index 20bc77a5..0af7aec3 100644 --- a/map/parallel_dirty_map.py +++ b/map/parallel_dirty_map.py @@ -768,14 +768,14 @@ def unpickle_list(list): # self.n_chan, self.n_ra, self.n_dec) #start_ind = (f_ra_start_ind,0,0,0,0) lock_and_write_buffer(thread_cov_inv_chunk, self.cov_filename, data_offset + dsize*f_ra_start_ind*self.n_dec*self.n_chan*self.n_ra*self.n_dec, dsize*thread_cov_inv_chunk.size, self.rank) - del thread_cov_inv_chunk - f = h5py.File(self.cov_filename, 'r+') - cov_inv_dset = f['inv_cov'] - cov_inv_shape = sp.zeros(shape = (self.n_chan, self.n_ra, self.n_dec, self,n_chan, self.n_ra, self.n_dec), dtype = dtype) - cov_inv_info = al.make_mat(cov_inv_shape, axis_names=('freq', 'ra', 'dec','freq', 'ra', 'dec'), row_axes=(0, 1, 2), col_axes=(3, 4, 5)) - cov_inv_info.copy_axis_info(map) - for key, value in cov_inv_info.info.iteritems(): - cov_inv_dset.attrs[key] = repr(value) + if self.rank == 0: + f = h5py.File(self.cov_filename, 'r+') + cov_inv_dset = f['inv_cov'] + cov_inv_shape = sp.zeros(shape = (self.n_chan, self.n_ra, self.n_dec, self,n_chan, self.n_ra, self.n_dec), dtype = dtype) + cov_inv_info = al.make_mat(cov_inv_shape, axis_names=('freq', 'ra', 'dec','freq', 'ra', 'dec'), row_axes=(0, 1, 2), col_axes=(3, 4, 5)) + cov_inv_info.copy_axis_info(map) + for key, value in cov_inv_info.info.iteritems(): + cov_inv_dset.attrs[key] = repr(value) # NOTE: using 'float' is not supprted in the saving because # it has to know if it is 32 or 64 bits. #dtype = thread_cov_inv_chunk.dtype @@ -792,16 +792,16 @@ def unpickle_list(list): # self.n_dec, self.n_ra, self.n_dec) #start_ind = (index_list[0],0,0,0,0) lock_and_write_buffer(thread_cov_inv_chunk, self.cov_filename, data_offset + dsize*index_list[0]*self.n_dec*self.n_ra*self.n_dec*self.n_ra, dsize*thread_cov_inv_chunk.size, self.rank) - #del thread_cov_inv_chunk - f = h5py.File(self.cov_filename, 'r+') - cov_inv_dset = f['inv_cov'] - #cov_inv = al.make_mat(cov_inv_dset, axis_names=('freq', 'ra', 'dec', 'ra', 'dec'), row_axes=(0, 1, 2), col_axes=(0, 3, 4)) - #cov_inv_shape = sp.zeros(shape = (self.n_chan, self.n_ra, self.n_dec, self.n_ra, self.n_dec), dtype = dtype) - #cov_inv_info = al.make_mat(cov_inv_shape, axis_names=('freq', 'ra', 'dec', 'ra', 'dec'), row_axes=(0, 1, 2), col_axes=(0, 3, 4)) - cov_inv_info = al.make_mat(cov_inv_shape, axis_names=('freq', 'ra', 'dec', 'ra', 'dec'), row_axes=(0, 1, 2), col_axes=(0, 3, 4)) - cov_inv_info.copy_axis_info(map) - for key, value in cov_inv_info.info.iteritems(): - cov_inv_dset.attrs[key] = repr(value) + if self.rank == 0: + f = h5py.File(self.cov_filename, 'r+') + cov_inv_dset = f['inv_cov'] + #cov_inv = al.make_mat(cov_inv_dset, axis_names=('freq', 'ra', 'dec', 'ra', 'dec'), row_axes=(0, 1, 2), col_axes=(0, 3, 4)) + #cov_inv_shape = sp.zeros(shape = (self.n_chan, self.n_ra, self.n_dec, self.n_ra, self.n_dec), dtype = dtype) + cov_inv_shape = _mapmaker_c.large_empty((self.n_chan, self.n_ra, self.n_dec, self.n_ra, self.n_dec)) + cov_inv_info = al.make_mat(cov_inv_shape, axis_names=('freq', 'ra', 'dec', 'ra', 'dec'), row_axes=(0, 1, 2), col_axes=(0, 3, 4)) + cov_inv_info.copy_axis_info(map) + for key, value in cov_inv_info.info.iteritems(): + cov_inv_dset.attrs[key] = repr(value) # NOTE: using 'float' is not supprted in the saving because # it has to know if it is 32 or 64 bits. #dtype = thread_cov_inv_chunk.dtype From 59930fa1d293d2b30ae5e8b226786a901f143345 Mon Sep 17 00:00:00 2001 From: Christopher Anderson Date: Fri, 18 Apr 2014 18:42:46 -0400 Subject: [PATCH 08/28] Putting in dataste attributes by hand, for no freq. corr. case.. --- map/parallel_dirty_map.py | 20 ++++++++++++++++---- 1 file changed, 16 insertions(+), 4 deletions(-) diff --git a/map/parallel_dirty_map.py b/map/parallel_dirty_map.py index 0af7aec3..19c8df4b 100644 --- a/map/parallel_dirty_map.py +++ b/map/parallel_dirty_map.py @@ -169,6 +169,7 @@ def execute(self, n_processes): spacing = params["pixel_spacing"] # Negative sign because RA increases from right to left. ra_spacing = -spacing/sp.cos(params['field_centre'][1]*sp.pi/180.) + self.ra_spacing = ra_spacing # To set some parameters, we have to read the first data file. first_file_name = (params["input_root"] + params["file_middles"][0] + params["input_end"]) @@ -795,12 +796,23 @@ def unpickle_list(list): if self.rank == 0: f = h5py.File(self.cov_filename, 'r+') cov_inv_dset = f['inv_cov'] + attrs = con_inv_dset.attrs + attrs.__setitem__('rows', (0, 1, 2)) + attrs.__setitem__('cols', (0, 3, 4)) + attrs.__setitem__('dec_centre', self.params['field_centre'][1]) + attrs.__setitem__('ra_centre', self.params['field_centre'][0]) + attrs.__setitem__('type', "'mat'") + attrs.__setitem__('axes', "('freq', 'ra', 'dec', 'ra', 'dec')") + attrs.__setitem__('freq_centre', band_centre) + attrs.__setitem__('freq_delta', self.delta_freq) + attrs.__setitem__('ra_delta', self.ra_spacing) + attrs.__setitem__('dec_delta', self.params['pixel_spacing']) #cov_inv = al.make_mat(cov_inv_dset, axis_names=('freq', 'ra', 'dec', 'ra', 'dec'), row_axes=(0, 1, 2), col_axes=(0, 3, 4)) #cov_inv_shape = sp.zeros(shape = (self.n_chan, self.n_ra, self.n_dec, self.n_ra, self.n_dec), dtype = dtype) - cov_inv_shape = _mapmaker_c.large_empty((self.n_chan, self.n_ra, self.n_dec, self.n_ra, self.n_dec)) - cov_inv_info = al.make_mat(cov_inv_shape, axis_names=('freq', 'ra', 'dec', 'ra', 'dec'), row_axes=(0, 1, 2), col_axes=(0, 3, 4)) - cov_inv_info.copy_axis_info(map) - for key, value in cov_inv_info.info.iteritems(): + #cov_inv_shape = _mapmaker_c.large_empty((self.n_chan, self.n_ra, self.n_dec, self.n_ra, self.n_dec)) + #cov_inv_info = al.make_mat(cov_inv_shape, axis_names=('freq', 'ra', 'dec', 'ra', 'dec'), row_axes=(0, 1, 2), col_axes=(0, 3, 4)) + #cov_inv_info.copy_axis_info(map) + #for key, value in cov_inv_info.info.iteritems(): cov_inv_dset.attrs[key] = repr(value) # NOTE: using 'float' is not supprted in the saving because # it has to know if it is 32 or 64 bits. From 07526a7b3a3a63334e64b3dd92ed79c8e8cde1db Mon Sep 17 00:00:00 2001 From: Christopher Anderson Date: Sat, 19 Apr 2014 02:32:09 -0400 Subject: [PATCH 09/28] Appears to be working with larger maps. --- map/parallel_dirty_map.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/map/parallel_dirty_map.py b/map/parallel_dirty_map.py index 19c8df4b..0f90a191 100644 --- a/map/parallel_dirty_map.py +++ b/map/parallel_dirty_map.py @@ -272,7 +272,8 @@ def execute(self, n_processes): # To write the noise_inverse, all we have to do is delete the # memeory map object. #del self.cov_inv, cov_inv - al.save(map_filename, map) + if self.rank == 0: + al.save(map_filename, map) if not self.noise_params is None: self.noise_params.close() @@ -796,7 +797,7 @@ def unpickle_list(list): if self.rank == 0: f = h5py.File(self.cov_filename, 'r+') cov_inv_dset = f['inv_cov'] - attrs = con_inv_dset.attrs + attrs = cov_inv_dset.attrs attrs.__setitem__('rows', (0, 1, 2)) attrs.__setitem__('cols', (0, 3, 4)) attrs.__setitem__('dec_centre', self.params['field_centre'][1]) @@ -813,7 +814,7 @@ def unpickle_list(list): #cov_inv_info = al.make_mat(cov_inv_shape, axis_names=('freq', 'ra', 'dec', 'ra', 'dec'), row_axes=(0, 1, 2), col_axes=(0, 3, 4)) #cov_inv_info.copy_axis_info(map) #for key, value in cov_inv_info.info.iteritems(): - cov_inv_dset.attrs[key] = repr(value) + # cov_inv_dset.attrs[key] = repr(value) # NOTE: using 'float' is not supprted in the saving because # it has to know if it is 32 or 64 bits. #dtype = thread_cov_inv_chunk.dtype From 714e8fd315d9dbed608568bddecce92b21374589 Mon Sep 17 00:00:00 2001 From: Tabitha Voytek Date: Mon, 21 Apr 2014 19:08:39 -0400 Subject: [PATCH 10/28] usage edits --- input/tcv/mm_test.ini | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/input/tcv/mm_test.ini b/input/tcv/mm_test.ini index cdd8a925..5a075955 100644 --- a/input/tcv/mm_test.ini +++ b/input/tcv/mm_test.ini @@ -10,7 +10,7 @@ from core import dir_data file_middles = tuple(dir_data.get_data_files(range(80,91), field='1hr', project='GBT10B_036', type='ralongmap')+dir_data.get_data_files(range(0,14), field='1hr', project='GBT11B_055',type='ralongmap')+dir_data.get_data_files(range(15,19), field='1hr', project='GBT11B_055', type='ralongmap')) # Use this line to shorten amount of data to use. -file_middles = file_middles[:10] +#file_middles = file_middles[:10] # file_middles = file_middles[:1] map_centre = (13.0, 1.85) @@ -49,7 +49,7 @@ dm_field_centre = map_centre dm_pixel_spacing = map_spacing dm_map_shape = map_shape dm_time_block = 'scan' -dm_n_files_group = 280 # tpb nodes. +dm_n_files_group = 10 # tpb nodes. dm_frequency_correlations = 'None' #dm_frequency_correlations = 'None' dm_number_frequency_modes = 0 From 98df3bedf6f8ec25b5fbbeab5a8f93d07d4de83c Mon Sep 17 00:00:00 2001 From: Christopher Anderson Date: Tue, 29 Apr 2014 00:07:57 -0400 Subject: [PATCH 11/28] clean_map makes .npy memmap for .hdf5 inv_cov. --- core/algebra.py | 78 +++++++++++++++++++++++++++++++++++++++ map/clean_map.py | 3 +- map/parallel_dirty_map.py | 3 ++ 3 files changed, 83 insertions(+), 1 deletion(-) diff --git a/core/algebra.py b/core/algebra.py index 31968287..617aa514 100644 --- a/core/algebra.py +++ b/core/algebra.py @@ -563,6 +563,84 @@ def load_h5(h5obj, path): iarray = info_array(iarray, info) return iarray +def load_h5_memmap(h5obj, path, fname, gb_lim): + """Load a info memmap from an hdf5 file. + + Parameters + ---------- + h5obj : h5py File or Group object + File from which the info array will be read from. + path : string + Path within `h5obj` to read the array. + + Returns + ------- + memarray : info_memmap + Array loaded from file. + """ + + # TODO: Allow `h5obj` to be a string with a path to a file to be opened + # and then closed. + #data = h5obj[path] + dset = h5obj[path] + memarray = sp.memmap(fname, shape=dset.shape, dtype=dset.dtype, mode='w+') + memarray = fill_memmap( memarray, 1000000000*gb_lim, dset) + #memarray = data.value + info = {} + for key, value in dset.attrs.iteritems(): + info[key] = safe_eval(value) + memarray = info_memmap(memarray, info) + return memarray + +def fill_memmap(memmap, mem_lim, dset): + chunk_size = dset.size*dset.dtype.itemsize + for ind in range(len(dset.shape)): + if chunk_size <= mem_lim: + i = ind + print str(i) + break + else: + chunk_size = chunk_size/dset.shape[ind] + if i == 0: + memmap = dset.value + else: + index_list = [] + for num in range(i): + index_list.append(range(dset.shape[num])) + indexes = cross(index_list) + for index in indexes: + print index + memmap[tuple(index)] = dset[tuple(index)] + return memmap + + +def cross(set_list): + '''Given a list of sets, return the cross product.''' + # By associativity of cross product, cross the first two sets together + # then cross that with the rest. The big conditional in the list + # comprehension is just to make sure that there are no nested lists + # in the final answer. + if len(set_list) == 1: + ans = [] + for elem in set_list[0]: + # In the 1D case, these elements are not lists. + if type(elem) == list: + ans.append(np.array(elem)) + else: + ans.append(np.array([elem])) + return ans + else: + A = set_list[0] + B = set_list[1] + cross_2 = [a+b if ((type(a) == list) and (type(b) == list)) else \ + (a+[b] if type(a) == list else ([a]+b if type(b) == list else \ + [a]+[b])) for a in A for b in B] + remaining = set_list[2:] + remaining.insert(0,cross_2) + return cross(remaining) + + + # ---- Functions for manipulating above arrays as matrices and vectors. ------- diff --git a/map/clean_map.py b/map/clean_map.py index 509eba23..5993c9d6 100644 --- a/map/clean_map.py +++ b/map/clean_map.py @@ -25,6 +25,7 @@ 'save_cholesky' : False, 'from_eig' : False, 'bands' : () + 'mem_lim' : () } prefix = 'cm_' @@ -130,7 +131,7 @@ def execute(self, nprocesses=1) : noise_inv = algebra.make_mat(noise_inv) elif noise_fname.split('.')[-1] == 'hdf5': noise_h5 = h5py.File(noise_fname, 'r') - noise_inv = algebra.load_h5(noise_h5, 'inv_cov') + noise_inv = algebra.load_h5_memmap(noise_h5, 'inv_cov', noisefname + ".npy" , params["mem_lim"]) noise_inv = algebra.make_mat(noise_inv) else: raise ValueError("Noise file is of unsupported type, neither .npy or .hdf5.") diff --git a/map/parallel_dirty_map.py b/map/parallel_dirty_map.py index 0f90a191..5c0309e0 100644 --- a/map/parallel_dirty_map.py +++ b/map/parallel_dirty_map.py @@ -10,6 +10,8 @@ from mpi4py import MPI +#from memory_profiler import profile + import math import threading from Queue import Queue @@ -906,6 +908,7 @@ def cross(set_list): remaining.insert(0,cross_2) return cross(remaining) +#@profile def lock_and_write_buffer(obj, fname, offset, size, proc): """Write the contents of a buffer to disk at a given offset, and explicitly lock the region of the file whilst doing so. From b08f702065f65e500ff300cb7be650e9c85a4b0d Mon Sep 17 00:00:00 2001 From: Tabitha Voytek Date: Tue, 29 Apr 2014 20:31:35 -0400 Subject: [PATCH 12/28] fixed small typos --- input/tcv/pointcorr_mapmaker.pipe | 27 +++++++++++++++------------ map/clean_map.py | 4 ++-- 2 files changed, 17 insertions(+), 14 deletions(-) diff --git a/input/tcv/pointcorr_mapmaker.pipe b/input/tcv/pointcorr_mapmaker.pipe index 967f4f05..9d549a70 100644 --- a/input/tcv/pointcorr_mapmaker.pipe +++ b/input/tcv/pointcorr_mapmaker.pipe @@ -7,7 +7,8 @@ import scipy as sp # What data to process and how to split it up. #field = '11hr' -field = '15hr' +#field = '15hr' +field = '1hr' #sessions = range(41, 90) # According to Tabitha, cal unstable in sessions 52 to 58. @@ -67,11 +68,11 @@ middles_temp = tuple(dir_data.get_data_files([47], field='15hr', # adjusted to (140,81) due to memory limits. # 11hr: RA center = 162.5, DEC center = 3.25 (Total is 18.5,9.5). Try pixel spacing = 0.25, map shape =(74,38) -map_centre = (217.8688,2.0) +map_centre = (13.0,1.85) # For maps where performace isn't important. #map0_shape = (100, 60) #map0_spacing = .040 -map0_shape = (78, 43) +map0_shape = (161, 83) map0_spacing = .0627 #map0_shape = (78, 43) #map0_spacing = .0627 @@ -88,11 +89,11 @@ base_dir = base_dir_main #base_tcv = os.getenv('GBT_TCV') # Which set of maps we are working on. #map_base = base_dir + 'maps/' -map_base = base_dir_main + 'maps/'+'15hr_41-80_pointcorr/' +map_base = base_dir_main + 'maps/'+'1hr_parallel_test/' # IO directory and file prefixes. prefix = '' -map_prefix = '' + field + '_' + '41-80_pointcorr' + '_' +map_prefix = '' + field + '_' + '80-18_ptcorr' + '_' # Maximum number of processes to use. pipe_processes = 1 @@ -137,8 +138,8 @@ pipe_modules = [] #from map import dirty_map #pipe_modules.append(dirty_map.DirtyMapMaker) -#from map import clean_map -#pipe_modules.append(clean_map.CleanMapMaker) +from map import clean_map +pipe_modules.append(clean_map.CleanMapMaker) #### Second Map making iteration #### @@ -161,8 +162,8 @@ pipe_modules = [] #pipe_modules.append((dirty_map.DirtyMapMaker, ('dmC_', 'dm_'))) #pipe_modules.append((dirty_map.DirtyMapMaker, ('dmD_', 'dm_'))) -from map import clean_map -pipe_modules.append((clean_map.CleanMapMaker, ('cmA_', 'cm_'))) +#from map import clean_map +#pipe_modules.append((clean_map.CleanMapMaker, ('cmA_', 'cm_'))) #pipe_modules.append((clean_map.CleanMapMaker, ('cmB_', 'cm_'))) #pipe_modules.append((clean_map.CleanMapMaker, ('cmC_', 'cm_'))) #pipe_modules.append((clean_map.CleanMapMaker, ('cmD_', 'cm_'))) @@ -321,7 +322,8 @@ dm_input_root = rp_output_root dm_file_middles = file_middles #dm_file_middles = middles_firs dm_input_end = '.fits' -dm_output_root = map_base + 'fir_' + map_prefix +dm_output_root = map_base + map_prefix +#dm_output_root = map_base + 'fir_' + map_prefix dm_scans = () dm_IFs = () @@ -342,10 +344,11 @@ dm_interpolation = 'cubic' # clean_map cm_input_root = dm_output_root cm_output_root = cm_input_root -#cm_polarizations = ('V',) -cm_polarizations = ('I','Q','U','V') +cm_polarizations = ('I',) +#cm_polarizations = ('I','Q','U','V') cm_bands = (800,) cm_save_noise_diag = True +cm_mem_lim = (2) #### Second round of map making ##### diff --git a/map/clean_map.py b/map/clean_map.py index 5993c9d6..53d9b1cc 100644 --- a/map/clean_map.py +++ b/map/clean_map.py @@ -24,7 +24,7 @@ 'save_noise_inv_diag' : False, 'save_cholesky' : False, 'from_eig' : False, - 'bands' : () + 'bands' : (), 'mem_lim' : () } prefix = 'cm_' @@ -131,7 +131,7 @@ def execute(self, nprocesses=1) : noise_inv = algebra.make_mat(noise_inv) elif noise_fname.split('.')[-1] == 'hdf5': noise_h5 = h5py.File(noise_fname, 'r') - noise_inv = algebra.load_h5_memmap(noise_h5, 'inv_cov', noisefname + ".npy" , params["mem_lim"]) + noise_inv = algebra.load_h5_memmap(noise_h5, 'inv_cov', noise_fname + ".npy" , params["mem_lim"]) noise_inv = algebra.make_mat(noise_inv) else: raise ValueError("Noise file is of unsupported type, neither .npy or .hdf5.") From 99f86f9b0ef594eda8a3253962af7ce3acf34ec5 Mon Sep 17 00:00:00 2001 From: Christopher Anderson Date: Wed, 30 Apr 2014 23:44:39 -0400 Subject: [PATCH 13/28] Dirty map should save inv_cov metadata as strings now. --- map/parallel_dirty_map.py | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/map/parallel_dirty_map.py b/map/parallel_dirty_map.py index 5c0309e0..66510ee9 100644 --- a/map/parallel_dirty_map.py +++ b/map/parallel_dirty_map.py @@ -10,7 +10,7 @@ from mpi4py import MPI -#from memory_profiler import profile +from memory_profiler import profile import math import threading @@ -800,16 +800,16 @@ def unpickle_list(list): f = h5py.File(self.cov_filename, 'r+') cov_inv_dset = f['inv_cov'] attrs = cov_inv_dset.attrs - attrs.__setitem__('rows', (0, 1, 2)) - attrs.__setitem__('cols', (0, 3, 4)) - attrs.__setitem__('dec_centre', self.params['field_centre'][1]) - attrs.__setitem__('ra_centre', self.params['field_centre'][0]) + attrs.__setitem__('rows', str((0, 1, 2))) + attrs.__setitem__('cols', str((0, 3, 4))) + attrs.__setitem__('dec_centre', str(self.params['field_centre'][1])) + attrs.__setitem__('ra_centre', str(self.params['field_centre'][0])) attrs.__setitem__('type', "'mat'") attrs.__setitem__('axes', "('freq', 'ra', 'dec', 'ra', 'dec')") - attrs.__setitem__('freq_centre', band_centre) - attrs.__setitem__('freq_delta', self.delta_freq) - attrs.__setitem__('ra_delta', self.ra_spacing) - attrs.__setitem__('dec_delta', self.params['pixel_spacing']) + attrs.__setitem__('freq_centre', str(band_centre)) + attrs.__setitem__('freq_delta', str(self.delta_freq)) + attrs.__setitem__('ra_delta', str(self.ra_spacing)) + attrs.__setitem__('dec_delta', str(self.params['pixel_spacing'])) #cov_inv = al.make_mat(cov_inv_dset, axis_names=('freq', 'ra', 'dec', 'ra', 'dec'), row_axes=(0, 1, 2), col_axes=(0, 3, 4)) #cov_inv_shape = sp.zeros(shape = (self.n_chan, self.n_ra, self.n_dec, self.n_ra, self.n_dec), dtype = dtype) #cov_inv_shape = _mapmaker_c.large_empty((self.n_chan, self.n_ra, self.n_dec, self.n_ra, self.n_dec)) @@ -908,7 +908,7 @@ def cross(set_list): remaining.insert(0,cross_2) return cross(remaining) -#@profile +@profile def lock_and_write_buffer(obj, fname, offset, size, proc): """Write the contents of a buffer to disk at a given offset, and explicitly lock the region of the file whilst doing so. From 8227f3a500361f26e08f64cc1ac2a65daf626074 Mon Sep 17 00:00:00 2001 From: Christopher Anderson Date: Mon, 5 May 2014 21:48:11 -0400 Subject: [PATCH 14/28] Much simpler fix for clean map with hdf5 noise inverse. --- input/km/mm_test_clean.ini | 67 +++++--------------------------------- map/clean_map.py | 19 +++++++---- map/parallel_dirty_map.py | 4 +-- 3 files changed, 23 insertions(+), 67 deletions(-) diff --git a/input/km/mm_test_clean.ini b/input/km/mm_test_clean.ini index f6b30481..34756646 100644 --- a/input/km/mm_test_clean.ini +++ b/input/km/mm_test_clean.ini @@ -4,69 +4,20 @@ import os from core import dir_data - -sessions = range(41, 81) -file_middles = tuple(dir_data.get_data_files(sessions, field='15hr', - project="GBT10B_036", - type='ralongmap')) -# Use this line to shorten amount of data to use. -file_middles = file_middles[:10] -# file_middles = file_middles[:1] - -map_centre = (217.87, 2.0) -# Large problem size, full map. -#map_shape = (90, 48) -#map_spacing = .0627 -# Standard problem size. -#map_shape = (72, 38) -#map_spacing = .0627 -# Small problem size. -#map_shape = (44, 24) -#map_spacing = .0627 -# Tiny problem size. -map_shape = (24, 14) -map_spacing = .0627 +map_centre = (30, -30) +map_shape = (10, 10) +map_spacing = 1 # Data paths -raid_pro = os.getenv("RAID_PRO") -input_data_dir = raid_pro + "kiyo/gbt_out_new/" -base_dir = os.getenv("GBT_OUT") -map_dir = base_dir + 'maps/first_lockandwrite_test/' -#map_dir = base_dir + 'maps/' -#map_root = map_dir + "tmp_test_parallel_freqcorr_" -map_root = map_dir + "tmp_test_parallel_" -#map_root = map_dir + "tmp_test_mustbesure_" -#map_root = map_dir + "tmp_test_regular_" -''' -# Map maker parameters. -dm_input_root = input_data_dir + 'reflagged_sec/' -dm_file_middles = file_middles -dm_input_end = '.fits' -dm_output_root = map_root -dm_scans = () -dm_IFs = (0,) -dm_polarizations = ('I',) -dm_field_centre = map_centre -dm_pixel_spacing = map_spacing -dm_map_shape = map_shape -dm_time_block = 'scan' -dm_n_files_group = 280 # tpb nodes. -#dm_frequency_correlations = 'measured' -dm_frequency_correlations = 'None' -dm_number_frequency_modes = 3 -dm_number_frequency_modes_discard = 1 -dm_noise_parameter_file = (input_data_dir + 'noise_measurments_sec/' - 'noise_parameters.shelve') -dm_deweight_time_mean = True -dm_deweight_time_slope = True -dm_interpolation = 'cubic' -dm_ts_foreground_mode_file = '' -dm_n_ts_foreground_modes = 0 -''' +input_data_dir = '/scratch/p/pen/andersoc/first_parkes_pipe/rotated_to_I/' +base_dir = '/scratch/p/pen/andersoc/first_parkes_pipe/' +map_dir = base_dir + 'parallel_maps/' +map_root = map_dir + "parkes_test_parallel_thread_ra30decn30_small_" + cm_input_root = map_root cm_output_root = cm_input_root cm_polarizations = ('I',) -cm_bands = (762,) +cm_bands = (1315,) cm_save_noise_diag = True diff --git a/map/clean_map.py b/map/clean_map.py index 53d9b1cc..405119a1 100644 --- a/map/clean_map.py +++ b/map/clean_map.py @@ -131,14 +131,15 @@ def execute(self, nprocesses=1) : noise_inv = algebra.make_mat(noise_inv) elif noise_fname.split('.')[-1] == 'hdf5': noise_h5 = h5py.File(noise_fname, 'r') - noise_inv = algebra.load_h5_memmap(noise_h5, 'inv_cov', noise_fname + ".npy" , params["mem_lim"]) - noise_inv = algebra.make_mat(noise_inv) + #noise_inv = algebra.load_h5_memmap(noise_h5, 'inv_cov', noise_fname + ".npy" , params["mem_lim"]) + #noise_inv = algebra.make_mat(noise_inv) + noise_inv = noise_h5['inv_cov'] else: raise ValueError("Noise file is of unsupported type, neither .npy or .hdf5.") # Two cases for the noise. If its the same shape as the map # then the noise is diagonal. Otherwise, it should be # block diagonal in frequency. - if noise_inv.ndim == 3 : + if len(noise_inv.shape) == 3 : if noise_inv.axes != ('freq', 'ra', 'dec') : msg = ("Expeced noise matrix to have axes " "('freq', 'ra', 'dec'), but it has: " @@ -155,12 +156,16 @@ def execute(self, nprocesses=1) : if save_noise_diag : noise_diag[good_data] = \ 1/noise_inv_memory[good_data] - elif noise_inv.ndim == 5 : - if noise_inv.axes != ('freq', 'ra', 'dec', 'ra', + elif len(noise_inv.shape) == 5 : + try: + axis_labels = noise_inv.axes + except: + axis_labels = eval(noise_inv.attrs['axes']) + if axis_labels != ('freq', 'ra', 'dec', 'ra', 'dec'): msg = ("Expeced noise matrix to have axes " "('freq', 'ra', 'dec', 'ra', 'dec'), " - "but it has: " + str(noise_inv.axes)) + "but it has: " + str(axis_labels)) raise ce.DataError(msg) # Arrange the dirty map as a vector. dirty_map_vect = sp.array(dirty_map) # A view. @@ -234,7 +239,7 @@ def execute(self, nprocesses=1) : if self.feedback > 1: print "" sys.stdout.flush() - elif noise_inv.ndim == 6 : + elif len(noise_inv.shape) == 6 : if save_noise_diag: # OLD WAY. #clean_map, noise_diag, chol = solve(noise_inv, diff --git a/map/parallel_dirty_map.py b/map/parallel_dirty_map.py index 66510ee9..8b910d59 100644 --- a/map/parallel_dirty_map.py +++ b/map/parallel_dirty_map.py @@ -10,7 +10,7 @@ from mpi4py import MPI -from memory_profiler import profile +#from memory_profiler import profile import math import threading @@ -908,7 +908,7 @@ def cross(set_list): remaining.insert(0,cross_2) return cross(remaining) -@profile +#@profile def lock_and_write_buffer(obj, fname, offset, size, proc): """Write the contents of a buffer to disk at a given offset, and explicitly lock the region of the file whilst doing so. From 04b231a8dd23a5d5e536c01845c72e7c864b2ba1 Mon Sep 17 00:00:00 2001 From: Christopher Anderson Date: Tue, 3 Jun 2014 22:35:47 -0400 Subject: [PATCH 15/28] Added MPI parallel clean mapper for no freq. correlation case --- map/parallel_clean_map.py | 574 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 574 insertions(+) create mode 100644 map/parallel_clean_map.py diff --git a/map/parallel_clean_map.py b/map/parallel_clean_map.py new file mode 100644 index 00000000..ebde75e2 --- /dev/null +++ b/map/parallel_clean_map.py @@ -0,0 +1,574 @@ +"""Converter from a dirty map to a clean map.""" + +from mpi4py import MPI +from map import parallel_dirty_map as pdm + +import sys +import glob +import time + +import numpy as np +import scipy as sp +import scipy.linalg as linalg +from ast import literal_eval as ev + +from core import algebra, hist +from kiyopy import parse_ini +import kiyopy.utils +import kiyopy.custom_exceptions as ce +import constants +import h5py + + +params_init = {'input_root' : './', + 'polarizations' : ('I',), + 'output_root' : './', + 'save_noise_diag' : False, + 'save_noise_inv_diag' : False, + 'save_cholesky' : False, + 'from_eig' : False, + 'bands' : (), + 'mem_lim' : () + } +prefix = 'cm_' + +comm = MPI.COMM_WORLD + +class CleanMapMaker(object) : + """Converts a Dirty map to a clean map.""" + + def __init__(self, parameter_file_or_dict=None, feedback=2) : + # Read in the parameters. + self.params = parse_ini.parse(parameter_file_or_dict, params_init, + prefix=prefix, feedback=feedback) + self.feedback = feedback + self.rank = comm.Get_rank() + self.nproc = comm.Get_size() + + def execute(self, nprocesses=1) : + """Worker funciton.""" + params = self.params + # Make parent directory and write parameter file. + kiyopy.utils.mkparents(params['output_root']) + parse_ini.write_params(params, params['output_root'] + 'params.ini', + prefix=prefix) + save_noise_diag = params['save_noise_diag'] + in_root = params['input_root'] + all_out_fname_list = [] + all_in_fname_list = [] + # Figure out what the band names are. + bands = params['bands'] + if not bands: + map_files = glob.glob(in_root + 'dirty_map_' + pol_str + "_*.npy") + bands = [] + root_len = len(in_root + 'dirty_map_') + for file_name in map_files: + bands.append(file_name[root_len:-4]) + # Loop over files to process. + for pol_str in params['polarizations']: + for band in bands: + if band == -1: + band_str = '' + else: + band_str = "_" + repr(band) + dmap_fname = (in_root + 'dirty_map_' + pol_str + + band_str + '.npy') + all_in_fname_list.append( + kiyopy.utils.abbreviate_file_path(dmap_fname)) + # Load the dirty map and the noise matrix. + dirty_map = algebra.load(dmap_fname) + dirty_map = algebra.make_vect(dirty_map) + if dirty_map.axes != ('freq', 'ra', 'dec') : + msg = ("Expeced dirty map to have axes ('freq'," + "'ra', 'dec'), but it has axes: " + + str(dirty_map.axes)) + raise ce.DataError(msg) + shape = dirty_map.shape + # Initialize the clean map. + clean_map = algebra.info_array(sp.zeros(dirty_map.shape)) + clean_map.info = dict(dirty_map.info) + clean_map = algebra.make_vect(clean_map) + # If needed, initialize a map for the noise diagonal. + if save_noise_diag : + noise_diag = algebra.zeros_like(clean_map) + if params["from_eig"]: + # Solving from eigen decomposition of the noise instead of + # the noise itself. + # Load in the decomposition. + evects_fname = (in_root + 'noise_evects_' + pol_str + + + band_str + '.npy') + if self.feedback > 1: + print "Using dirty map: " + dmap_fname + print "Using eigenvectors: " + evects_fname + evects = algebra.open_memmap(evects_fname, 'r') + evects = algebra.make_mat(evects) + evals_inv_fname = (in_root + 'noise_evalsinv_' + pol_str + + "_" + repr(band) + '.npy') + evals_inv = algebra.load(evals_inv_fname) + evals_inv = algebra.make_mat(evals_inv) + # Solve for the map. + if params["save_noise_diag"]: + clean_map, noise_diag = solve_from_eig(evals_inv, + evects, dirty_map, True, self.feedback) + else: + clean_map = solve_from_eig(evals_inv, + evects, dirty_map, False, self.feedback) + # Delete the eigen vectors to recover memory. + del evects + else: + # Solving from the noise. + noise_fname1 = (in_root + 'noise_inv_' + pol_str + + band_str + '.npy') + noise_fname2 = (in_root + 'noise_inv_' + pol_str + + band_str + '.hdf5') + noise_fnames = glob.glob(noise_fname1)+glob.glob(noise_fname2) + if len(noise_fnames) == 0: + raise ValueError("Couldn't find noise_inv file.") + if len(noise_fnames) > 1: + raise ValueError("Multiple files have naming pattern of noise_inv") + noise_fname = noise_fnames[0] + if self.feedback > 1: + print "Using dirty map: " + dmap_fname + print "Using noise inverse: " + noise_fname + all_in_fname_list.append( + kiyopy.utils.abbreviate_file_path(noise_fname)) + + if noise_fname.split('.')[-1] == 'npy': + noise_inv = algebra.open_memmap(noise_fname, 'r') + noise_inv = algebra.make_mat(noise_inv) + elif noise_fname.split('.')[-1] == 'hdf5': + noise_h5 = h5py.File(noise_fname, 'r') + #noise_inv = algebra.load_h5_memmap(noise_h5, 'inv_cov', noise_fname + ".npy" , params["mem_lim"]) + #noise_inv = algebra.make_mat(noise_inv) + noise_inv = noise_h5['inv_cov'] + else: + raise ValueError("Noise file is of unsupported type, neither .npy or .hdf5.") + # Two cases for the noise. If its the same shape as the map + # then the noise is diagonal. Otherwise, it should be + # block diagonal in frequency. + if len(noise_inv.shape) == 3 : + if noise_inv.axes != ('freq', 'ra', 'dec') : + msg = ("Expeced noise matrix to have axes " + "('freq', 'ra', 'dec'), but it has: " + + str(noise_inv.axes)) + raise ce.DataError(msg) + # Noise inverse can fit in memory, so copy it. + noise_inv_memory = sp.array(noise_inv, copy=True) + # Find the non-singular (covered) pixels. + max_information = noise_inv_memory.max() + good_data = noise_inv_memory < 1.0e-10*max_information + # Make the clean map. + clean_map[good_data] = (dirty_map[good_data] + / noise_inv_memory[good_data]) + if save_noise_diag : + noise_diag[good_data] = \ + 1/noise_inv_memory[good_data] + elif len(noise_inv.shape) == 5 : + try: + axis_labels = noise_inv.axes + except: + axis_labels = eval(noise_inv.attrs['axes']) + if axis_labels != ('freq', 'ra', 'dec', 'ra', + 'dec'): + msg = ("Expeced noise matrix to have axes " + "('freq', 'ra', 'dec', 'ra', 'dec'), " + "but it has: " + str(axis_labels)) + raise ce.DataError(msg) + # Arrange the dirty map as a vector. + dirty_map_vect = sp.array(dirty_map) # A view. + dirty_map_vect.shape = (shape[0], shape[1]*shape[2]) + frequencies = dirty_map.get_axis('freq')/1.0e6 + # Allowcate memory only once. + noise_inv_freq = sp.empty((shape[1], shape[2], + shape[1], shape[2]), dtype=float) + if self.feedback > 1 : + print "Inverting noise matrix." + chan_index_list = range(noise_inv.shape[0]) + index_list, junk = pdm.split_elems(chan_index_list, self.nproc) + index_list = index_list(self.rank) + # Block diagonal in frequency so loop over frequencies. + for ii in index_list: + if self.feedback > 1: + print "Frequency: ", "%5.1f"%(frequencies[ii]), + if self.feedback > 2: + print ", start mmap read:", + sys.stdout.flush() + noise_inv_freq[...] = noise_inv[ii, ...] + if self.feedback > 2: + print "done, start eig:", + sys.stdout.flush() + noise_inv_freq.shape = (shape[1]*shape[2], + shape[1]*shape[2]) + # Solve the map making equation by diagonalization. + noise_inv_diag, Rot = sp.linalg.eigh( + noise_inv_freq, overwrite_a=True) + if self.feedback > 2: + print "done", + map_rotated = sp.dot(Rot.T, dirty_map_vect[ii]) + # Zero out infinite noise modes. + bad_modes = (noise_inv_diag + < 1.0e-5 * noise_inv_diag.max()) + if self.feedback > 1: + print ", discarded: ", + print "%4.1f" % (100.0 * sp.sum(bad_modes) + / bad_modes.size), + print "% of modes", + if self.feedback > 2: + print ", start rotations:", + sys.stdout.flush() + map_rotated[bad_modes] = 0. + noise_inv_diag[bad_modes] = 1.0 + # Solve for the clean map and rotate back. + map_rotated /= noise_inv_diag + map = sp.dot(Rot, map_rotated) + if self.feedback > 2: + print "done", + sys.stdout.flush() + # Fill the clean array. + map.shape = (shape[1], shape[2]) + clean_map[ii, ...] = map + if save_noise_diag : + # Using C = R Lambda R^T + # where Lambda = diag(1/noise_inv_diag). + temp_noise_diag = 1/noise_inv_diag + temp_noise_diag[bad_modes] = 0 + # Multiply R by the diagonal eigenvalue matrix. + # Broadcasting does equivalent of mult by diag + # matrix. + temp_mat = Rot*temp_noise_diag + # Multiply by R^T, but only calculate the + # diagonal elements. + for jj in range(shape[1]*shape[2]) : + temp_noise_diag[jj] = sp.dot( + temp_mat[jj,:], Rot[jj,:]) + temp_noise_diag.shape = (shape[1], shape[2]) + noise_diag[ii, ...] = temp_noise_diag + # Return workspace memory to origional shape. + noise_inv_freq.shape = (shape[1], shape[2], + shape[1], shape[2]) + if self.feedback > 1: + print "" + sys.stdout.flush() + elif len(noise_inv.shape) == 6 : + if save_noise_diag: + # OLD WAY. + #clean_map, noise_diag, chol = solve(noise_inv, + # dirty_map, True, feedback=self.feedback) + # NEW WAY. + clean_map, noise_diag, noise_inv_diag, chol = \ + solve(noise_fname, noise_inv, dirty_map, + True, feedback=self.feedback) + else: + # OLD WAY. + #clean_map, chol = solve(noise_inv, dirty_map, + # False, feedback=self.feedback) + # NEW WAY. + clean_map, noise_inv_diag, chol = \ + solve(noise_fname, noise_inv, dirty_map, + False, feedback=self.feedback) + if params['save_cholesky']: + chol_fname = (params['output_root'] + 'chol_' + + pol_str + band_str + '.npy') + sp.save(chol_fname, chol) + if params['save_noise_inv_diag']: + noise_inv_diag_fname = (params['output_root'] + + 'noise_inv_diag_' + pol_str + band_str + + '.npy') + algebra.save(noise_inv_diag_fname, noise_inv_diag) + # Delete the cholesky to recover memory. + del chol + else : + raise ce.DataError("Noise matrix has bad shape.") + # In all cases delete the noise object to recover memeory. + del noise_inv + # Write the clean map to file. + out_fname = (params['output_root'] + 'clean_map_' + + pol_str + band_str + '.npy') + if self.feedback > 1: + print "Writing clean map to: " + out_fname + algebra.save(out_fname, clean_map) + all_out_fname_list.append( + kiyopy.utils.abbreviate_file_path(out_fname)) + if save_noise_diag : + noise_diag_fname = (params['output_root'] + 'noise_diag_' + + pol_str + band_str + '.npy') + algebra.save(noise_diag_fname, noise_diag) + all_out_fname_list.append( + kiyopy.utils.abbreviate_file_path(noise_diag_fname)) + # Check the clean map for faileur. + if not sp.alltrue(sp.isfinite(clean_map)): + n_bad = sp.sum(sp.logical_not(sp.isfinite(clean_map))) + msg = ("Non finite entries found in clean map. Solve" + " failed. %d out of %d entries bad" + % (n_bad, clean_map.size)) + raise RuntimeError(msg) + # This needs to be added to the new dirty map maker before I can + # add it here. + # Finally update the history object. + #history = hist.read(in_root + 'history.hist') + #history.add('Read map and noise files:', all_in_fname_list) + #history.add('Converted dirty map to clean map.', + # all_out_fname_list) + #h_fname = params['output_root'] + "history.hist" + #history.write(h_fname) + +def solve(noise_inv_filename, noise_inv, dirty_map, + return_noise_diag=False, feedback=0): + """Solve for the clean map. + + Matrix and vector passed in as algebra objects with no block diagonality. + The system is solved using a GPU accelerated cholesky decomposition. + Optionally, the diagonal on the noise matrix is returned. + """ + + # Import the cython stuff locally so some clean maps can be made on any + # machine. + import _cholesky as _c + # Put into the 2D matrix shape. + expanded = noise_inv.view() + side_size = noise_inv.shape[0] * noise_inv.shape[1] * noise_inv.shape[2] + expanded.shape = (side_size,) * 2 + # Allowcate memory for the cholesky and copy the upper triangular data. + # Instead of copying, open the matrix in copy on write mode. + if feedback > 1: + print "Copying matrix." + time_before = time.time() / 60. + # OLD WAY. + # tri_copy = _c.up_tri_copy(expanded) + # NEW WAY. + tri_copy = up_tri_copy_from_file(noise_inv_filename) + time_after = time.time() / 60. + if feedback > 1: + print "\nMatrix copying time: %.2f minutes." % (time_after-time_before) + #tri_copy = expanded + # Get the diagonal of the giant noise inverse. + if feedback > 1: + print "Getting noise inverse diagonal." + time_before = time.time() / 60. + tri_copy.shape = side_size*side_size + noise_inv_diag = tri_copy[::side_size+1] + tri_copy.shape = (side_size,side_size) + noise_inv_diag.shape = dirty_map.shape + noise_inv_diag = algebra.as_alg_like(noise_inv_diag, dirty_map) + time_after = time.time() / 60. + if feedback > 1: + print "\nNoise inverse diagonal gotten in: %.2f minutes." \ + % (time_after-time_before) + # Cholesky decompose it. + if feedback > 1: + print "Cholesky decomposition." + time_before = time.time() / 60. + _c.call_cholesky(tri_copy) + time_after = time.time() / 60. + if feedback > 1: + print "\nCholesky decomposition time: %.2f minutes." \ + % (time_after-time_before) + # Solve for the clean map. + flat_map = dirty_map.view() + flat_map.shape = (flat_map.size,) + if feedback > 1: + print "Solving for clean map." + time_before = time.time() / 60. + clean_map = _c.cho_solve(tri_copy, flat_map) + time_after = time.time() / 60. + if feedback > 1: + print "\nClean map solving time: %.2f minutes." % (time_after-time_before) + # Reshape and cast as a map. + clean_map.shape = dirty_map.shape + clean_map = algebra.as_alg_like(clean_map, dirty_map) + if not return_noise_diag: + return clean_map, noise_inv_diag, tri_copy + else: + if feedback > 1: + print "Getting noise diagonal." + noise_diag = sp.empty(side_size, dtype=float) + time_before = time.time() / 60. + _c.inv_diag_from_chol(tri_copy, noise_diag) + time_after = time.time() / 60. + if feedback > 1: + print "\nNoise diagonal gotten in: %.2f minutes." \ + % (time_after-time_before) + noise_diag.shape = dirty_map.shape + noise_diag = algebra.as_alg_like(noise_diag, dirty_map) + return clean_map, noise_diag, noise_inv_diag, tri_copy + +def solve_from_eig(noise_evalsinv, noise_evects, dirty_map, + return_noise_diag=False, feedback=0): + """Converts a dirty map to a clean map using the eigen decomposition of the + noise inverse. + """ + + # Check the shapes. + if noise_evects.ndim != 4: + raise ValueError("Expected 4D array for 'noise_evects`.") + if noise_evalsinv.shape != (noise_evects.shape[-1],): + raise ValueError("Wrong number of eigenvalues.") + if dirty_map.shape != noise_evects.shape[:-1]: + raise ValueError("Dirty map and noise don't have matching dimensions.") + if dirty_map.size != noise_evects.shape[-1]: + raise ValueError("Eigen space not the same total size as map space.") + n = noise_evalsinv.shape[0] + nf = dirty_map.shape[0] + nr = dirty_map.shape[1] + nd = dirty_map.shape[2] + # Copy the eigenvalues. + noise_evalsinv = noise_evalsinv.copy() + # Find poorly constrained modes and zero them out. + bad_inds = noise_evalsinv < 1./constants.T_huge**2 + n_bad = sp.sum(bad_inds) + if feedback > 1: + print ("Discarding %d modes of %d. %f percent." + % (n_bad, n, 100. * n_bad / n)) + noise_evalsinv[bad_inds] = 1. + # Rotate the dirty map into the diagonal noise space. + if feedback > 1: + print "Rotating map to eigenspace." + map_rot = sp.zeros(n, dtype=sp.float64) + for ii in xrange(nf): + for jj in xrange(nr): + for kk in xrange(nd): + tmp = noise_evects[ii,jj,kk,:].copy() + map_rot += dirty_map[ii,jj,kk] * tmp + # Multiply by the (diagonal) noise (inverse, inverse). Zero out any poorly + # constrained modes. + map_rot[bad_inds] = 0 + # Take inverse and multiply. + map_rot = map_rot / noise_evalsinv + # Now rotate back to the origional space. + if feedback > 1: + print "Rotating back to map space." + clean_map = algebra.zeros_like(dirty_map) + for ii in xrange(nf): + for jj in xrange(nr): + for kk in xrange(nd): + tmp = noise_evects[ii,jj,kk,:].copy() + clean_map[ii,jj,kk] = sp.sum(map_rot * tmp) + if return_noise_diag: + if feedback > 1: + print "Getting noise diagonal." + noise_diag = algebra.zeros_like(dirty_map) + noise_evals = 1. / noise_evalsinv + noise_evals[bad_inds] = constants.T_huge**2 + for ii in xrange(nf): + for jj in xrange(nr): + for kk in xrange(nd): + tmp = noise_evects[ii,jj,kk,:].copy() + noise_diag[ii,jj,kk] = sp.sum(tmp**2 * noise_evals) + return clean_map, noise_diag + else: + return clean_map + + +def read_n_numbers(f,arr,dt,n,row_len,row_pos): + '''Read n numbers from an open file stream, f, and add them to the next + available row, row_pos, in the 2D array, arr. The data type, dt, is needed + to convert from bytes to usable numbers. row_len and row_pos are needed + so that only the upper triangular part of a row is added to arr.''' + bytes_per_num = 0 + # Assign how long a number is in memory. + if ((dt.name == 'int64') or (dt.name == 'float64')): + bytes_per_num = 8 + if ((dt.name == 'int32') or (dt.name == 'float32')): + bytes_per_num = 4 + # Check if possible to run. + if (bytes_per_num == 0): + msg = "%s is not a supported number type right now." % (dt.name) + raise ce.DataError(msg) + # Read in n numbers. + raw_nums = f.read(bytes_per_num*n) + # Convert numbers in bytes to useable numbers all at once. + new_slice = np.fromstring(raw_nums,dtype=dt,count=n) + # n is a multiple of rows. + rows_to_read = int(n/row_len) + # Copy each row from the buffer (in memory) to the array (in memory). + # We are copying in only the upper triangular part of the array so + # we must only write in the appropriate values. + for i in range(rows_to_read): + # Every row down the matrix, we want to ignore one more number + # at the beginning of a row. + arr[row_pos+i,row_pos+i:row_len] = \ + new_slice[i*row_len+row_pos+i:i*row_len+row_len] + + +def up_tri_copy_from_file(filename): + '''Return the upper triagular part of the array in filename. The file + must be a binary .npy file.''' + f = open(filename, 'rb') + # The first 6 bytes are a magic string: exactly "\x93NUMPY". + numpy_string = f.read(6) + # The next 1 byte is an unsigned byte: the major version number + # of the file format, e.g. \x01. + major_ver = ord(f.read(1)) + # The next 1 byte is an unsigned byte: the minor version number + # of the file format, e.g. \x00. + minor_ver = ord(f.read(1)) + # Check that the version of the file used is what this code can handle. + if ((numpy_string != "\x93NUMPY") or (major_ver != 1) or (minor_ver != 0)): + msg = "Array can only be read from a '\93NUMPY' version 1.0 .npy file " + msg += "not %s %d.%d" % (numpy_string, major_ver, minor_ver) + raise ce.DataError(msg) + # The next 2 bytes form a little-endian unsigned short int: the + # length of the header data HEADER_LEN. + byte2,byte1 = f.read(2) + # Get the value from a short int (2 bytes) to decimal. + header_len = (16**2)*ord(byte1) + ord(byte2) + # The next HEADER_LEN bytes form the header data describing the + # array's format. It is an ASCII string which contains a Python + # literal expression of a dictionary. It is terminated by a newline + # ('\n') and padded with spaces ('\x20') to make the total length of + # the magic string + 4 + HEADER_LEN be evenly divisible by 16 for + # alignment purposes. + raw_dict = f.read(header_len) + # Take off trailing uselessness. + raw_dict = raw_dict.rstrip('\n') + raw_dict = raw_dict.rstrip() + header_dict = ev(raw_dict) + # Check that it is possible to read array. + # "fortran_order" : bool + # Whether the array data is Fortran-contiguous or not. + if (header_dict['fortran_order']): + msg = "Array must be in C order, not fortran order." + raise ce.DataError(msg) + # "shape" : tuple of int + # The shape of the array. + arr_shape = header_dict['shape'] + # "descr" : dtype.descr + # An object that can be passed as an argument to the + # numpy.dtype() constructor to create the array's dtype. + dt = np.dtype(header_dict['descr']) + # Where to save all the data. Note this does not make a new array in memory. + num_nums = np.product(arr_shape) + # Assume array is square + row_len = int(np.sqrt(num_nums)) + num_rows = row_len + loaded_arr = np.empty((row_len,row_len)) + print "Rows copied:" + # Load the array reading n numbers at a time. Reading one number at + # a time is best for memory but requires too many disk seeks for a + # memory map. Loading too many numbers at once requires too much memory, + # so choose a happy medium. + # Assume array is square. Reads n rows of numbers at a time. + rows_to_read = 1000 + n = rows_to_read*row_len + row_pos = 0 + while (row_pos < num_rows): + # n may not divide evenly into the total number of numbers. + # This statement can only be True at the end of the array. + if ((num_rows-row_pos) < rows_to_read): + rows_to_read = num_rows-row_pos + n = rows_to_read*row_len + read_n_numbers(f,loaded_arr,dt,n,row_len,row_pos) + # Changing row_pos in the function called does not change it here. + row_pos += rows_to_read + sys.stderr.write(str(row_pos)+' ') + print + f.close() + return loaded_arr + +# If this file is run from the command line, execute the main function. +if __name__ == "__main__": + import sys + if len(sys.argv) == 2: + par_file = sys.argv[1] + nproc = 1 + CleanMapMaker(par_file).execute(nproc) From 811e0be6209f02a1715909e12fc3ba366d108c6e Mon Sep 17 00:00:00 2001 From: Christopher Anderson Date: Wed, 4 Jun 2014 11:48:09 -0400 Subject: [PATCH 16/28] Fixed I/O for parallel clean map. Hopefully works now. --- map/parallel_clean_map.py | 39 ++++++++++++++++++++++----------------- 1 file changed, 22 insertions(+), 17 deletions(-) diff --git a/map/parallel_clean_map.py b/map/parallel_clean_map.py index ebde75e2..944ce93f 100644 --- a/map/parallel_clean_map.py +++ b/map/parallel_clean_map.py @@ -184,8 +184,8 @@ def execute(self, nprocesses=1) : if self.feedback > 1 : print "Inverting noise matrix." chan_index_list = range(noise_inv.shape[0]) - index_list, junk = pdm.split_elems(chan_index_list, self.nproc) - index_list = index_list(self.rank) + index_list_full, junk = pdm.split_elems(chan_index_list, self.nproc) + index_list = index_list_full(self.rank) # Block diagonal in frequency so loop over frequencies. for ii in index_list: if self.feedback > 1: @@ -249,6 +249,10 @@ def execute(self, nprocesses=1) : if self.feedback > 1: print "" sys.stdout.flush() + for process_n in self.nproc: + for ii in index_list_full[process_n]: + clean_map[ii, ...] = comm.bcast(clean_map[ii, ...], root = process_n) + noise_diag[ii, ...] = comm.bcast(noise_diag[ii, ...], root = process_n) elif len(noise_inv.shape) == 6 : if save_noise_diag: # OLD WAY. @@ -284,24 +288,25 @@ def execute(self, nprocesses=1) : # Write the clean map to file. out_fname = (params['output_root'] + 'clean_map_' + pol_str + band_str + '.npy') - if self.feedback > 1: - print "Writing clean map to: " + out_fname - algebra.save(out_fname, clean_map) - all_out_fname_list.append( - kiyopy.utils.abbreviate_file_path(out_fname)) - if save_noise_diag : - noise_diag_fname = (params['output_root'] + 'noise_diag_' - + pol_str + band_str + '.npy') - algebra.save(noise_diag_fname, noise_diag) + if self.rank == 0: + if self.feedback > 1: + print "Writing clean map to: " + out_fname + algebra.save(out_fname, clean_map) all_out_fname_list.append( - kiyopy.utils.abbreviate_file_path(noise_diag_fname)) - # Check the clean map for faileur. - if not sp.alltrue(sp.isfinite(clean_map)): - n_bad = sp.sum(sp.logical_not(sp.isfinite(clean_map))) - msg = ("Non finite entries found in clean map. Solve" + kiyopy.utils.abbreviate_file_path(out_fname)) + if save_noise_diag : + noise_diag_fname = (params['output_root'] + 'noise_diag_' + + pol_str + band_str + '.npy') + algebra.save(noise_diag_fname, noise_diag) + all_out_fname_list.append( + kiyopy.utils.abbreviate_file_path(noise_diag_fname)) + # Check the clean map for faileur. + if not sp.alltrue(sp.isfinite(clean_map)): + n_bad = sp.sum(sp.logical_not(sp.isfinite(clean_map))) + msg = ("Non finite entries found in clean map. Solve" " failed. %d out of %d entries bad" % (n_bad, clean_map.size)) - raise RuntimeError(msg) + raise RuntimeError(msg) # This needs to be added to the new dirty map maker before I can # add it here. # Finally update the history object. From 164a46d7594426328c94e25a6aff2f5d48a56fdf Mon Sep 17 00:00:00 2001 From: Christopher Anderson Date: Tue, 24 Jun 2014 14:57:37 -0400 Subject: [PATCH 17/28] Fixed parallel clean mapmaker --- input/km/dm_parkes_thread_test.ini | 64 ++++++++++++++++++++++++++++++ input/km/mm_test_clean.ini | 2 +- map/parallel_clean_map.py | 19 +++++++-- 3 files changed, 80 insertions(+), 5 deletions(-) create mode 100755 input/km/dm_parkes_thread_test.ini diff --git a/input/km/dm_parkes_thread_test.ini b/input/km/dm_parkes_thread_test.ini new file mode 100755 index 00000000..9c4f7f82 --- /dev/null +++ b/input/km/dm_parkes_thread_test.ini @@ -0,0 +1,64 @@ +# Input file for map maker testing. + +import os + +from core import dir_data + + +sessions = range(41, 81) +file_list = open('scripts/datalist_2008_center.txt', 'r').read().splitlines() +#file_list = file_list[:4] +file_middles = [] +for file in file_list: + file_middles.append(file[72:-7]) + + +map_centre = (30, -30) +# Large problem size, full map. +#map_shape = (90, 48) +#map_spacing = .0627 +# Standard problem size. +#map_shape = (72, 38) +#map_spacing = .0627 +# Small problem size. +#map_shape = (44, 24) +#map_spacing = .0627 +# Tiny problem size. +map_shape = (10, 10) +map_spacing = 1 + +# Data paths +#raid_pro = os.getenv("RAID_PRO") +input_data_dir = '/scratch/p/pen/andersoc/first_parkes_pipe/rotated_to_I/' +base_dir = '/scratch/p/pen/andersoc/first_parkes_pipe/' +map_dir = base_dir + 'parallel_maps/' +map_root = map_dir + "parkes_test_parallel_thread_ra30decn30_small_" + + +# Map maker parameters. +dm_input_root = '/scratch/p/pen/andersoc/first_parkes_pipe/rotated_to_I/' +dm_file_middles = file_middles +dm_input_end = '.fits' +dm_output_root = map_root +dm_scans = () +dm_IFs = (0,) +dm_beams = (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12) + +dm_thread_divide = True +dm_polarizations = ('I',) +dm_field_centre = map_centre +dm_pixel_spacing = map_spacing +dm_map_shape = map_shape +dm_time_block = 'scan' +dm_n_files_group = 10 +#dm_frequency_correlations = 'measured' +dm_frequency_correlations = 'None' +dm_number_frequency_modes = 0 +#dm_number_frequency_modes_discard = 1 +#dm_noise_parmeters_input_root = 'None' +dm_noise_parameter_file = '' +dm_deweight_time_mean = True +dm_deweight_time_slope = True +dm_interpolation = 'linear' +#dm_ts_foreground_mode_file = '' +#dm_n_ts_foreground_modes = 0 diff --git a/input/km/mm_test_clean.ini b/input/km/mm_test_clean.ini index 34756646..270cd5af 100644 --- a/input/km/mm_test_clean.ini +++ b/input/km/mm_test_clean.ini @@ -17,7 +17,7 @@ map_root = map_dir + "parkes_test_parallel_thread_ra30decn30_small_" cm_input_root = map_root -cm_output_root = cm_input_root +cm_output_root = cm_input_root + '_not_parallel' cm_polarizations = ('I',) cm_bands = (1315,) cm_save_noise_diag = True diff --git a/map/parallel_clean_map.py b/map/parallel_clean_map.py index 944ce93f..48c8617f 100644 --- a/map/parallel_clean_map.py +++ b/map/parallel_clean_map.py @@ -185,7 +185,7 @@ def execute(self, nprocesses=1) : print "Inverting noise matrix." chan_index_list = range(noise_inv.shape[0]) index_list_full, junk = pdm.split_elems(chan_index_list, self.nproc) - index_list = index_list_full(self.rank) + index_list = index_list_full[self.rank] # Block diagonal in frequency so loop over frequencies. for ii in index_list: if self.feedback > 1: @@ -249,10 +249,21 @@ def execute(self, nprocesses=1) : if self.feedback > 1: print "" sys.stdout.flush() - for process_n in self.nproc: + for process_n in range(self.nproc): for ii in index_list_full[process_n]: - clean_map[ii, ...] = comm.bcast(clean_map[ii, ...], root = process_n) - noise_diag[ii, ...] = comm.bcast(noise_diag[ii, ...], root = process_n) + print clean_map[ii, ...].shape + if self.rank == process_n: + comm.Send(clean_map[ii, ...], dest=0, tag=13) + if self.rank == 0: + comm.Recv(clean_map[ii, ...], source=process_n, tag=13) + comm.Barrier() + if self.rank == process_n: + comm.Send(noise_diag[ii, ...], dest=0, tag=13) + if self.rank == 0: + comm.Recv(noise_diag[ii, ...], source=process_n, tag=13) + comm.Barrier() + #clean_map[ii, ...] = comm.bcast(clean_map[ii, ...], root = process_n) + #noise_diag[ii, ...] = comm.bcast(noise_diag[ii, ...], root = process_n) elif len(noise_inv.shape) == 6 : if save_noise_diag: # OLD WAY. From f144100ab8dc3b299f15cde97b1b764b204e7461 Mon Sep 17 00:00:00 2001 From: Tabitha Voytek Date: Mon, 30 Jun 2014 19:09:42 -0400 Subject: [PATCH 18/28] usage change --- input/tcv/mm_test.ini | 17 +++++++++++++---- 1 file changed, 13 insertions(+), 4 deletions(-) diff --git a/input/tcv/mm_test.ini b/input/tcv/mm_test.ini index cdd8a925..1fc6424f 100644 --- a/input/tcv/mm_test.ini +++ b/input/tcv/mm_test.ini @@ -7,10 +7,19 @@ from core import dir_data #sessions = range(41, 81) #file_middles = tuple(dir_data.get_data_files(sessions, field='15hr', project="GBT10B_036", type = 'ralongmap')) -file_middles = tuple(dir_data.get_data_files(range(80,91), field='1hr', project='GBT10B_036', type='ralongmap')+dir_data.get_data_files(range(0,14), field='1hr', project='GBT11B_055',type='ralongmap')+dir_data.get_data_files(range(15,19), field='1hr', project='GBT11B_055', type='ralongmap')) + +first_set = dir_data.get_data_files(range(80,91),field='1hr',project='GBT10B_036', type='ralongmap') +second_set = dir_data.get_data_files(range(0,14),field='1hr',project='GBT11B_055', type='ralongmap') +third_set = dir_data.get_data_files(range(15,19),field='1hr',project='GBT11B_055',type='ralongmap') +fourth_set = dir_data.get_data_files(range(1,28),field='1hr',project='GBT13B_352',type='ralongmap') +fourth_set = fourth_set[1:] +file_middles = tuple(first_set+second_set+third_set+fourth_set) +print len(file_middles) + +#file_middles = tuple(dir_data.get_data_files(range(80,91), field='1hr', project='GBT10B_036', type='ralongmap')+dir_data.get_data_files(range(0,14), field='1hr', project='GBT11B_055',type='ralongmap')+dir_data.get_data_files(range(15,19), field='1hr', project='GBT11B_055', type='ralongmap')) # Use this line to shorten amount of data to use. -file_middles = file_middles[:10] +#file_middles = file_middles[:10] # file_middles = file_middles[:1] map_centre = (13.0, 1.85) @@ -34,7 +43,7 @@ base_dir = os.getenv("GBT_OUT") input_data_dir = base_dir+"rotated_to_I_Q_avg_fdgp_ptcorr/" #map_dir = base_dir + 'maps/first_lockandwrite_test/' #map_root = map_dir + "tmp_test_parallel_freqcorr_" -map_root = base_dir+"maps/1hr_parallel_test/1hr_80-18_ptcorr_" +map_root = base_dir+"maps/1hr_parallel_test/1hr_80-28_ptcorr_" #map_root = map_dir + "tmp_test_parallel_" # Map maker parameters. @@ -49,7 +58,7 @@ dm_field_centre = map_centre dm_pixel_spacing = map_spacing dm_map_shape = map_shape dm_time_block = 'scan' -dm_n_files_group = 280 # tpb nodes. +dm_n_files_group = 10 # tpb nodes. dm_frequency_correlations = 'None' #dm_frequency_correlations = 'None' dm_number_frequency_modes = 0 From b133e785c6db8d926d89365bb765a9d866567ea9 Mon Sep 17 00:00:00 2001 From: Christopher Anderson Date: Fri, 18 Jul 2014 17:53:21 -0400 Subject: [PATCH 19/28] Dirty mapper can make beam maps for parkes. --- map/parallel_clean_map.py | 2 ++ map/parallel_dirty_map.py | 17 +++++++++++------ 2 files changed, 13 insertions(+), 6 deletions(-) diff --git a/map/parallel_clean_map.py b/map/parallel_clean_map.py index 48c8617f..3619d030 100644 --- a/map/parallel_clean_map.py +++ b/map/parallel_clean_map.py @@ -45,6 +45,7 @@ def __init__(self, parameter_file_or_dict=None, feedback=2) : self.rank = comm.Get_rank() self.nproc = comm.Get_size() + #@profile def execute(self, nprocesses=1) : """Worker funciton.""" params = self.params @@ -141,6 +142,7 @@ def execute(self, nprocesses=1) : #noise_inv = algebra.load_h5_memmap(noise_h5, 'inv_cov', noise_fname + ".npy" , params["mem_lim"]) #noise_inv = algebra.make_mat(noise_inv) noise_inv = noise_h5['inv_cov'] + noise_h5.close() else: raise ValueError("Noise file is of unsupported type, neither .npy or .hdf5.") # Two cases for the noise. If its the same shape as the map diff --git a/map/parallel_dirty_map.py b/map/parallel_dirty_map.py index 8b910d59..f3458b2a 100644 --- a/map/parallel_dirty_map.py +++ b/map/parallel_dirty_map.py @@ -5,13 +5,12 @@ the pointing operator (`Pointing`) and the time domain noise operator (`Noise`). """ - +#import os +#os.environ['PYTHON_EGG_CACHE'] = '/scratch/p/pen/andersoc/.python-eggs' import sys from mpi4py import MPI -#from memory_profiler import profile - import math import threading from Queue import Queue @@ -89,7 +88,9 @@ # let the noise know about it. 'n_ts_foreground_modes' : 0, 'ts_foreground_mode_file' : '', - 'thread_divide' : False + 'thread_divide' : False, + # Which beam to use to make map. Choose zero to use all beams. + 'beam' : 0 } comm = MPI.COMM_WORLD @@ -126,9 +127,13 @@ def iterate_data(self, file_middles): force_tuple=True) if params['time_block'] == 'scan': for Data in Blocks: - yield (Data,), middle + if params["beam"] == 0 or params["beam"] == Data.field['BEAM']: + yield (Data,), middle elif params['time_block'] == 'file': - yield Blocks, middle + if params["beam"] == 0: + yield Blocks, middle + else: + yield tuple([data for data in Blocks if data.field['BEAM'] == params["beam"]]) else: msg = "time_block parameter must be 'scan' or 'file'." raise dm.ValueError(msg) From c8123583bbb5ab0288dc781d88d150c6cfee6596 Mon Sep 17 00:00:00 2001 From: Christopher Anderson Date: Tue, 22 Jul 2014 23:16:10 -0400 Subject: [PATCH 20/28] memory_profiler module --- map/parallel_dirty_map.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/map/parallel_dirty_map.py b/map/parallel_dirty_map.py index f3458b2a..38bc0758 100644 --- a/map/parallel_dirty_map.py +++ b/map/parallel_dirty_map.py @@ -5,9 +5,10 @@ the pointing operator (`Pointing`) and the time domain noise operator (`Noise`). """ -#import os -#os.environ['PYTHON_EGG_CACHE'] = '/scratch/p/pen/andersoc/.python-eggs' +import os +os.environ['PYTHON_EGG_CACHE'] = '/scratch/p/pen/andersoc/.python-eggs' import sys +from memory_profiler import profile from mpi4py import MPI @@ -913,7 +914,7 @@ def cross(set_list): remaining.insert(0,cross_2) return cross(remaining) -#@profile +@profile def lock_and_write_buffer(obj, fname, offset, size, proc): """Write the contents of a buffer to disk at a given offset, and explicitly lock the region of the file whilst doing so. From 86cc16a8a620de351416faf8ac999bda6adb46bc Mon Sep 17 00:00:00 2001 From: Tabitha Voytek Date: Wed, 23 Jul 2014 11:38:58 -0400 Subject: [PATCH 21/28] Usage edits --- input/tcv/mm_test.ini | 184 +++++++++++++++++++++++++++++- input/tcv/pointcorr_mapmaker.pipe | 14 ++- 2 files changed, 191 insertions(+), 7 deletions(-) diff --git a/input/tcv/mm_test.ini b/input/tcv/mm_test.ini index 70440e7c..e86addd0 100644 --- a/input/tcv/mm_test.ini +++ b/input/tcv/mm_test.ini @@ -3,7 +3,7 @@ import os from core import dir_data - +import scipy as sp #sessions = range(41, 81) #file_middles = tuple(dir_data.get_data_files(sessions, field='15hr', project="GBT10B_036", type = 'ralongmap')) @@ -14,8 +14,33 @@ third_set = dir_data.get_data_files(range(15,19),field='1hr',project='GBT11B_055 fourth_set = dir_data.get_data_files(range(1,28),field='1hr',project='GBT13B_352',type='ralongmap') fourth_set = fourth_set[1:] file_middles = tuple(first_set+second_set+third_set+fourth_set) +file_middles = file_middles[:-3] print len(file_middles) +n_files = len(file_middles) +new_file_middles = [] +for i in range(0,n_files): +# if file_middles[i]=='GBT10B_036/44_wigglez15hrst_ralongmap_169-176': +# print file_middles[i] +# elif file_middles[i]=='GBT10B_036/47_wigglez15hrst_ralongmap_230-237': +# print file_middles[i] +# elif file_middles[i]=='GBT10B_036/47_wigglez15hrst_ralongmap_58-65': +# print file_middles[i] +# elif file_middles[i]=='GBT10B_036/47_wigglez15hrst_ralongmap_10-17': +# print file_middles[i] +# elif file_middles[i]=='GBT10B_036/47_wigglez15hrst_ralongmap_110-117': +# print file_middles[i] +# else : + new_file_middles.append(file_middles[i]) +new_file_middles = tuple(new_file_middles) +n_files = len(new_file_middles) + +middles_a = file_middles[:n_files//4] +middles_b = file_middles[n_files//4:2*n_files//4] +middles_c = file_middles[2*n_files//4:3*n_files//4] +middles_d = file_middles[3*n_files//4:n_files] + + #file_middles = tuple(dir_data.get_data_files(range(80,91), field='1hr', project='GBT10B_036', type='ralongmap')+dir_data.get_data_files(range(0,14), field='1hr', project='GBT11B_055',type='ralongmap')+dir_data.get_data_files(range(15,19), field='1hr', project='GBT11B_055', type='ralongmap')) # Use this line to shorten amount of data to use. @@ -45,6 +70,8 @@ input_data_dir = base_dir+"rotated_to_I_Q_avg_fdgp_ptcorr/" #map_root = map_dir + "tmp_test_parallel_freqcorr_" map_root = base_dir+"maps/1hr_parallel_test/1hr_80-28_ptcorr_" #map_root = map_dir + "tmp_test_parallel_" +map_base = base_dir+'maps/1hr_parallel_test/' +map_prefix = '_1hr_80-28_ptcorr_' # Map maker parameters. dm_input_root = input_data_dir @@ -58,11 +85,7 @@ dm_field_centre = map_centre dm_pixel_spacing = map_spacing dm_map_shape = map_shape dm_time_block = 'scan' -<<<<<<< HEAD dm_n_files_group = 10 # tpb nodes. -======= -dm_n_files_group = 10 # tpb nodes. ->>>>>>> 164a46d7594426328c94e25a6aff2f5d48a56fdf dm_frequency_correlations = 'None' #dm_frequency_correlations = 'None' dm_number_frequency_modes = 0 @@ -73,3 +96,154 @@ dm_deweight_time_slope = True dm_interpolation = 'cubic' #dm_ts_foreground_mode_file = '' #dm_n_ts_foreground_modes = 0 + +cm_input_root = dm_output_root +cm_output_root = cm_input_root +cm_polarizations = ('I',) +#cm_polarizations = ('I','Q','U','V') +cm_bands = (800,) +cm_save_noise_diag = True +#cm_mem_lim = (2) + +# dirty_map +# Map A +dmA_input_root = input_data_dir +dmA_file_middles = middles_a +dmA_input_end = '.fits' +dmA_output_root = map_base + 'secA_' + map_prefix +dmA_scans = () +dmA_IFs = () + +#dmA_polarizations = ('I','Q','U','V') +dmA_polarizations = ('I',) +dmA_field_centre = map_centre +dmA_pixel_spacing = map_spacing +dmA_map_shape = map_shape +dmA_time_block = 'scan' +#dmA_n_files_group = 420 # prawn +dmA_n_files_group = 10 # tpb nodes. +#dmA_n_files_group = 120 +dmA_frequency_correlations = 'None' +#dmA_number_frequency_modes = 3 +dmA_number_frequency_modes = 3 # Probably most appropriate. +dmA_number_frequency_modes_discard = 1 +#dmA_noise_parameter_file = '' +dmA_noise_parameter_file = base_dir+'noise_measurements_sec/noise_parameters_take2.shelve' +dmA_deweight_time_mean = True +dmA_deweight_time_slope = True +dmA_interpolation = 'cubic' +dmA_ts_foreground_mode_file = '' +dmA_n_ts_foreground_modes = 0 + +# Other maps mostly copy parameters of map A. +dmB_file_middles = middles_b +dmB_output_root = map_base + 'secB_' + map_prefix + +dmC_file_middles = middles_c +dmC_output_root = map_base + 'secC_' + map_prefix + +dmD_file_middles = middles_d +dmD_output_root = map_base + 'secD_' + map_prefix + +# clean_map +# Map A +cmA_input_root = dmA_output_root +cmA_output_root = cmA_input_root + +cmA_polarizations = ('I',) +#cmA_polarizations = ('I','Q','U','V') +cmA_bands = (800,) +cmA_save_noise_diag = True +cmA_save_cholesky = False +cmA_from_eig = False + +# Other maps +cmB_input_root = dmB_output_root +cmB_output_root = cmB_input_root + +cmC_input_root = dmC_output_root +cmC_output_root = cmC_input_root + +cmD_input_root = dmD_output_root +cmD_output_root = cmD_input_root + +# Parameters that are set to be the same as the A case. +dmB_input_root = dmA_input_root +dmB_input_end = dmA_input_end +dmB_scans = dmA_scans +dmB_IFs = dmA_IFs +dmB_polarizations = dmA_polarizations +dmB_field_centre = dmA_field_centre +dmB_map_shape = dmA_map_shape +dmB_pixel_spacing = dmA_pixel_spacing +dmB_time_block = dmA_time_block +dmB_n_files_group = dmA_n_files_group +dmB_frequency_correlations = dmA_frequency_correlations +dmB_number_frequency_modes = dmA_number_frequency_modes +dmB_number_frequency_modes_discard = dmA_number_frequency_modes_discard +dmB_noise_parameter_file = dmA_noise_parameter_file +dmB_deweight_time_mean = dmA_deweight_time_mean +dmB_deweight_time_slope = dmA_deweight_time_slope +dmB_interpolation = dmA_interpolation +dmB_ts_foreground_mode_file = dmA_ts_foreground_mode_file +dmB_n_ts_foreground_modes = dmA_n_ts_foreground_modes + +dmC_input_root = dmA_input_root +dmC_input_end = dmA_input_end +dmC_scans = dmA_scans +dmC_IFs = dmA_IFs +dmC_polarizations = dmA_polarizations +dmC_field_centre = dmA_field_centre +dmC_map_shape = dmA_map_shape +dmC_pixel_spacing = dmA_pixel_spacing +dmC_time_block = dmA_time_block +dmC_n_files_group = dmA_n_files_group +dmC_frequency_correlations = dmA_frequency_correlations +dmC_number_frequency_modes = dmA_number_frequency_modes +dmC_number_frequency_modes_discard = dmA_number_frequency_modes_discard +dmC_noise_parameter_file = dmA_noise_parameter_file +dmC_deweight_time_mean = dmA_deweight_time_mean +dmC_deweight_time_slope = dmA_deweight_time_slope +dmC_interpolation = dmA_interpolation +dmC_ts_foreground_mode_file = dmA_ts_foreground_mode_file +dmC_n_ts_foreground_modes = dmA_n_ts_foreground_modes + +dmD_input_root = dmA_input_root +dmD_input_end = dmA_input_end +dmD_scans = dmA_scans +dmD_IFs = dmA_IFs +dmD_polarizations = dmA_polarizations +dmD_field_centre = dmA_field_centre +dmD_map_shape = dmA_map_shape +dmD_pixel_spacing = dmA_pixel_spacing +dmD_time_block = dmA_time_block +dmD_n_files_group = dmA_n_files_group +dmD_frequency_correlations = dmA_frequency_correlations +dmD_number_frequency_modes = dmA_number_frequency_modes +dmD_number_frequency_modes_discard = dmA_number_frequency_modes_discard +dmD_noise_parameter_file = dmA_noise_parameter_file +dmD_deweight_time_mean = dmA_deweight_time_mean +dmD_deweight_time_slope = dmA_deweight_time_slope +dmD_interpolation = dmA_interpolation +dmD_ts_foreground_mode_file = dmA_ts_foreground_mode_file +dmD_n_ts_foreground_modes = dmA_n_ts_foreground_modes + + +cmB_polarizations = cmA_polarizations +cmB_bands = cmA_bands +cmB_save_noise_diag = cmA_save_noise_diag +cmB_save_cholesky = cmA_save_cholesky +cmB_from_eig = cmA_from_eig + +cmC_polarizations = cmA_polarizations +cmC_bands = cmA_bands +cmC_save_noise_diag = cmA_save_noise_diag +cmC_save_cholesky = cmA_save_cholesky +cmC_from_eig = cmA_from_eig + +cmD_polarizations = cmA_polarizations +cmD_bands = cmA_bands +cmD_save_noise_diag = cmA_save_noise_diag +cmD_save_cholesky = cmA_save_cholesky +cmD_from_eig = cmA_from_eig + diff --git a/input/tcv/pointcorr_mapmaker.pipe b/input/tcv/pointcorr_mapmaker.pipe index 9d549a70..a7f37876 100644 --- a/input/tcv/pointcorr_mapmaker.pipe +++ b/input/tcv/pointcorr_mapmaker.pipe @@ -28,11 +28,21 @@ field = '1hr' #file_middles = tuple(dir_data.get_data_files(range(01,19),field='11hr',project='GBT12A_418',type='ralongmap')) #file_middles = tuple(dir_data.get_data_files(range(80,91), field='1hr', project='GBT10B_036', type='ralongmap')+dir_data.get_data_files(range(0,14), field='1hr', project='GBT11B_055',type='ralongmap')+dir_data.get_data_files(range(15,19), field='1hr', project='GBT11B_055', type='ralongmap')) #file_middles = tuple(dir_data.get_data_files(range(17,19),field='1hr',project='GBT11B_055', type='ralongmap')) -file_middles = tuple(dir_data.get_data_files(range(41,57),field='15hr',project='GBT10B_036',type='ralongmap')+dir_data.get_data_files(range(58,80),field='15hr',project='GBT10B_036',type='ralongmap')) +#file_middles = tuple(dir_data.get_data_files(range(41,57),field='15hr',project='GBT10B_036',type='ralongmap')+dir_data.get_data_files(range(58,80),field='15hr',project='GBT10B_036',type='ralongmap')) #bad_ind = sp.where(file_middles=='GBT10B_036/44_wigglez15hr_ralongmap_169-176')[0] #print bad_ind #file_middles[bad_ind] = '' +first_set = dir_data.get_data_files(range(80,91),field='1hr',project='GBT10B_036', type='ralongmap') +second_set = dir_data.get_data_files(range(0,14),field='1hr',project='GBT11B_055', type='ralongmap') +third_set = dir_data.get_data_files(range(15,19),field='1hr',project='GBT11B_055',type='ralongmap') +fourth_set = dir_data.get_data_files(range(1,28),field='1hr',project='GBT13B_352',type='ralongmap') +fourth_set = fourth_set[1:] +file_middles = tuple(first_set+second_set+third_set+fourth_set) +file_middles = file_middles[:-3] +print len(file_middles) + + n_files = len(file_middles) new_file_middles = [] for i in range(0,n_files): @@ -93,7 +103,7 @@ map_base = base_dir_main + 'maps/'+'1hr_parallel_test/' # IO directory and file prefixes. prefix = '' -map_prefix = '' + field + '_' + '80-18_ptcorr' + '_' +map_prefix = '' + field + '_' + '80-28_ptcorr' + '_' # Maximum number of processes to use. pipe_processes = 1 From 47920b8c10ebe3bcbd71d695f5a2ba82f68033dd Mon Sep 17 00:00:00 2001 From: Chris Anderson Date: Fri, 25 Jul 2014 14:18:55 -0400 Subject: [PATCH 22/28] Little fix for beam dirty map maker. --- map/parallel_dirty_map.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/map/parallel_dirty_map.py b/map/parallel_dirty_map.py index f3458b2a..7de94c03 100644 --- a/map/parallel_dirty_map.py +++ b/map/parallel_dirty_map.py @@ -133,7 +133,8 @@ def iterate_data(self, file_middles): if params["beam"] == 0: yield Blocks, middle else: - yield tuple([data for data in Blocks if data.field['BEAM'] == params["beam"]]) + yield tuple([data for data in Blocks if data.field['BEAM'] +== params["beam"]]), middle else: msg = "time_block parameter must be 'scan' or 'file'." raise dm.ValueError(msg) From 15f074ddcc0392740aff51182d732ef81b80e41a Mon Sep 17 00:00:00 2001 From: Chris Anderson Date: Fri, 25 Jul 2014 17:54:17 -0400 Subject: [PATCH 23/28] Parallel clean mapper closes hdf5 file after using. --- map/parallel_clean_map.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/map/parallel_clean_map.py b/map/parallel_clean_map.py index 3619d030..e17899b9 100644 --- a/map/parallel_clean_map.py +++ b/map/parallel_clean_map.py @@ -142,7 +142,6 @@ def execute(self, nprocesses=1) : #noise_inv = algebra.load_h5_memmap(noise_h5, 'inv_cov', noise_fname + ".npy" , params["mem_lim"]) #noise_inv = algebra.make_mat(noise_inv) noise_inv = noise_h5['inv_cov'] - noise_h5.close() else: raise ValueError("Noise file is of unsupported type, neither .npy or .hdf5.") # Two cases for the noise. If its the same shape as the map @@ -298,6 +297,12 @@ def execute(self, nprocesses=1) : raise ce.DataError("Noise matrix has bad shape.") # In all cases delete the noise object to recover memeory. del noise_inv + try: + noise_h5 + except: + print "Not using hdf5 noise, no need to close." + else: + noise_h5.close() # Write the clean map to file. out_fname = (params['output_root'] + 'clean_map_' + pol_str + band_str + '.npy') From 41dfdb735c63c4a110fd590a9833923461bae3b1 Mon Sep 17 00:00:00 2001 From: Chris Anderson Date: Tue, 29 Jul 2014 01:29:10 -0400 Subject: [PATCH 24/28] Fixes for freq correl parallel dirty map. --- map/parallel_dirty_map.py | 25 ++++++++++++++++++------- 1 file changed, 18 insertions(+), 7 deletions(-) diff --git a/map/parallel_dirty_map.py b/map/parallel_dirty_map.py index 7de94c03..37cbabd7 100644 --- a/map/parallel_dirty_map.py +++ b/map/parallel_dirty_map.py @@ -535,7 +535,7 @@ def thread_work(): else: ra_index_list = range(self.n_ra) # cross product = A x B = (a,b) for all a in A, for all b in B - chan_ra_index_list = clacross([chan_index_list,ra_index_list]) + chan_ra_index_list = cross([chan_index_list,ra_index_list]) index_list,start_list = split_elems(chan_ra_index_list, self.nproc) index_list = index_list[self.rank] @@ -620,7 +620,7 @@ def thread_work(): if run == 0 and start_file_ind == 0: #Adding prior to the diagonal. #thread_cov_inv_row_adder = thread_cov_inv_row[:, thread_f_ind, thread_ra_ind, :] - thread_cov_inv_row_chunk[ii,:, thread_f_ind, thread_ra_ind,:].flat[:: self.n_dec + 1] += \ + thread_cov_inv_chunk[ii,:, thread_f_ind, thread_ra_ind,:].flat[:: self.n_dec + 1] += \ 1.0 / T_large**2 if (self.feedback > 1 and thread_ra_ind == self.n_ra - 1): print thread_f_ind, @@ -781,11 +781,22 @@ def unpickle_list(list): if self.rank == 0: f = h5py.File(self.cov_filename, 'r+') cov_inv_dset = f['inv_cov'] - cov_inv_shape = sp.zeros(shape = (self.n_chan, self.n_ra, self.n_dec, self,n_chan, self.n_ra, self.n_dec), dtype = dtype) - cov_inv_info = al.make_mat(cov_inv_shape, axis_names=('freq', 'ra', 'dec','freq', 'ra', 'dec'), row_axes=(0, 1, 2), col_axes=(3, 4, 5)) - cov_inv_info.copy_axis_info(map) - for key, value in cov_inv_info.info.iteritems(): - cov_inv_dset.attrs[key] = repr(value) + attrs = cov_inv_dset.attrs + attrs.__setitem__('rows', str((0, 1, 2))) + attrs.__setitem__('cols', str((3, 4, 5))) + attrs.__setitem__('dec_centre', str(self.params['field_centre'][1])) + attrs.__setitem__('ra_centre', str(self.params['field_centre'][0])) + attrs.__setitem__('type', "'mat'") + attrs.__setitem__('axes', "('freq', 'ra', 'dec', 'freq', 'ra', 'dec')") + attrs.__setitem__('freq_centre', str(band_centre)) + attrs.__setitem__('freq_delta', str(self.delta_freq)) + attrs.__setitem__('ra_delta', str(self.ra_spacing)) + attrs.__setitem__('dec_delta', str(self.params['pixel_spacing'])) + #cov_inv_shape = sp.zeros(shape = (self.n_chan, self.n_ra, self.n_dec, self.n_chan, self.n_ra, self.n_dec), dtype = dtype) + #cov_inv_info = al.make_mat(cov_inv_shape, axis_names=('freq', 'ra', 'dec','freq', 'ra', 'dec'), row_axes=(0, 1, 2), col_axes=(3, 4, 5)) + #cov_inv_info.copy_axis_info(map) + #for key, value in cov_inv_info.info.iteritems(): + # cov_inv_dset.attrs[key] = repr(value) # NOTE: using 'float' is not supprted in the saving because # it has to know if it is 32 or 64 bits. #dtype = thread_cov_inv_chunk.dtype From 048733d8248a8bd1fb78b6ffaae5b3ffdcda6ccb Mon Sep 17 00:00:00 2001 From: Christopher Anderson Date: Thu, 2 Oct 2014 14:57:04 -0400 Subject: [PATCH 25/28] Minor changes. --- input/km/mm_test_clean.ini | 2 +- map/parallel_clean_map.py | 9 +++++++-- map/parallel_dirty_map.py | 1 + 3 files changed, 9 insertions(+), 3 deletions(-) diff --git a/input/km/mm_test_clean.ini b/input/km/mm_test_clean.ini index 270cd5af..34756646 100644 --- a/input/km/mm_test_clean.ini +++ b/input/km/mm_test_clean.ini @@ -17,7 +17,7 @@ map_root = map_dir + "parkes_test_parallel_thread_ra30decn30_small_" cm_input_root = map_root -cm_output_root = cm_input_root + '_not_parallel' +cm_output_root = cm_input_root cm_polarizations = ('I',) cm_bands = (1315,) cm_save_noise_diag = True diff --git a/map/parallel_clean_map.py b/map/parallel_clean_map.py index e17899b9..a4cf1626 100644 --- a/map/parallel_clean_map.py +++ b/map/parallel_clean_map.py @@ -1,9 +1,14 @@ """Converter from a dirty map to a clean map.""" +import os +os.environ['PYTHON_EGG_CACHE'] = '/scratch/p/pen/andersoc/.python-eggs' +from memory_profiler import profile +import sys + from mpi4py import MPI from map import parallel_dirty_map as pdm -import sys +#import sys import glob import time @@ -45,7 +50,7 @@ def __init__(self, parameter_file_or_dict=None, feedback=2) : self.rank = comm.Get_rank() self.nproc = comm.Get_size() - #@profile + @profile def execute(self, nprocesses=1) : """Worker funciton.""" params = self.params diff --git a/map/parallel_dirty_map.py b/map/parallel_dirty_map.py index bbfec4d1..2d695845 100644 --- a/map/parallel_dirty_map.py +++ b/map/parallel_dirty_map.py @@ -286,6 +286,7 @@ def execute(self, n_processes): if not self.noise_params is None: self.noise_params.close() + @profile def make_map(self): """Makes map for current polarization and band. From 8cf6fb2631d12741405d2c5edac180757bdf0aeb Mon Sep 17 00:00:00 2001 From: Christopher Anderson Date: Wed, 29 Oct 2014 18:23:48 -0400 Subject: [PATCH 26/28] Added .ini files in input/km --- input/km/cm_parkes_parallel2.ini | 103 ++++++++++++++++++++++++++++ input/km/dm_parkes_thread2.ini | 111 +++++++++++++++++++++++++++++++ 2 files changed, 214 insertions(+) create mode 100755 input/km/cm_parkes_parallel2.ini create mode 100755 input/km/dm_parkes_thread2.ini diff --git a/input/km/cm_parkes_parallel2.ini b/input/km/cm_parkes_parallel2.ini new file mode 100755 index 00000000..367d6cdf --- /dev/null +++ b/input/km/cm_parkes_parallel2.ini @@ -0,0 +1,103 @@ +# Input file for map maker testing. + +import os +import glob + +#from core import dir_data + + +#sessions = range(41, 81) +#file_list = open('scripts/datalist_2008_center.txt', 'r').read().splitlines() +#file_list = file_list[:4] +file_middles = [] +#for file in file_list: +# file_middles.append(file[72:-7]) + +data_dir = '/scratch/p/pen/andersoc/second_parkes_pipe/rebinned/' +dir_len = len(data_dir) +for el in os.walk(data_dir): + dir_files = glob.glob(el[0] + '/*' + '2df1' + '*.fits') + for file in dir_files: + file_middles.append(file.split('.')[0][dir_len:]) + + +map_centre = (10, -30) +# Large problem size, full map. +#map_shape = (90, 48) +#map_spacing = .0627 +# Standard problem size. +#map_shape = (72, 38) +#map_spacing = .0627 +# Small problem size. +#map_shape = (44, 24) +#map_spacing = .0627 +# Tiny problem size. +map_shape = (188, 88) +map_spacing = .08 + +# Data paths +#raid_pro = os.getenv("RAID_PRO") +input_data_dir = data_dir +#base_dir = '/scratch/p/pen/andersoc/first_parkes_pipe/' +base_dir = '/scratch/p/pen/andersoc/second_parkes_pipe/' +map_dir = base_dir + 'maps/' +#map_dir = base_dir + 'parallel_maps/' +#map_root = map_dir + "parkes_test_parallel_thread_ra30decn30_small_" +map_root = map_dir + "parkes_parallel_thread_ra20decn30_p08_125by88_beam1_" +#map_root = map_dir + "parkes_parallel_thread_ra10decn30_p08_125by88_beam2_" +#map_root = map_dir + "parkes_parallel_thread_ra10decn30_p08_125by88_beam3_" +#map_root = map_dir + "parkes_parallel_thread_ra10decn30_p08_125by88_beam4_" +#map_root = map_dir + "parkes_parallel_thread_ra10decn30_p08_125by88_beam5_" +#map_root = map_dir + "parkes_parallel_thread_ra10decn30_p08_125by88_beam6_" +#map_root = map_dir + "parkes_parallel_thread_ra10decn30_p08_125by88_beam7_" +#map_root = map_dir + "parkes_parallel_thread_ra10decn30_p08_125by88_beam8_" +#map_root = map_dir + "parkes_parallel_thread_ra10decn30_p08_125by88_beam9_" +#map_root = map_dir + "parkes_parallel_thread_ra10decn30_p08_125by88_beam10_" +#map_root = map_dir + "parkes_parallel_thread_ra10decn30_p08_125by88_beam11_" +#map_root = map_dir + "parkes_parallel_thread_ra10decn30_p08_125by88_beam12_" +#map_root = map_dir + "parkes_parallel_thread_ra10decn30_p08_125by88_beam13_" +#map_root = map_dir + "parkes_parallel_thread_ra10decn30_p08_312by88_beam2_" +#map_root = map_dir + "parkes_parallel_thread_ra10decn30_p08_312by88_beam3_" + +# Map maker parameters. +#dm_input_root = '/scratch/p/pen/andersoc/first_parkes_pipe/rotated_to_I/' +dm_input_root = input_data_dir +dm_file_middles = file_middles +dm_input_end = '.fits' +dm_output_root = map_root +dm_scans = () +dm_IFs = (0,) +#dm_beams = (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12) +#dm_beam = 1 +#dm_beam = 2 +dm_beam = 3 + +dm_thread_divide = True +dm_polarizations = ('XX', 'YY') +#dm_polarizations = ('XX',) +dm_field_centre = map_centre +dm_pixel_spacing = map_spacing +dm_map_shape = map_shape +dm_time_block = 'scan' +dm_n_files_group = 10 +#dm_frequency_correlations = 'measured' +dm_frequency_correlations = 'None' +dm_number_frequency_modes = 0 +#dm_number_frequency_modes_discard = 1 +#dm_noise_parmeters_input_root = 'None' +dm_noise_parameter_file = '' +dm_deweight_time_mean = True +dm_deweight_time_slope = True +dm_interpolation = 'linear' +#dm_ts_foreground_mode_file = '' +#dm_n_ts_foreground_modes = 0 + +cm_input_root = dm_output_root +cm_output_root = cm_input_root +cm_polarizations = ('XX','YY') +#cm_polarizations = ('XX',) +#cm_polarizations = ('YY',) +#cm_polarizations = ('I',) +cm_save_noise_diag = True +cm_bands = (1316,) +#cm_bands = (1315,) diff --git a/input/km/dm_parkes_thread2.ini b/input/km/dm_parkes_thread2.ini new file mode 100755 index 00000000..e846378f --- /dev/null +++ b/input/km/dm_parkes_thread2.ini @@ -0,0 +1,111 @@ +# Input file for map maker testing. + +import os +import glob + +#from core import dir_data + + +#sessions = range(41, 81) +#file_list = open('scripts/datalist_2008_center.txt', 'r').read().splitlines() +#file_list = file_list[:4] +file_middles = [] +#for file in file_list: +# file_middles.append(file[72:-7]) + +data_dir = '/scratch/p/pen/andersoc/second_parkes_pipe/rebinned/' +dir_len = len(data_dir) +for el in os.walk(data_dir): + dir_files = glob.glob(el[0] + '/*' + '2df1' + '*.fits') + for file in dir_files: + file_middles.append(file.split('.')[0][dir_len:]) + + +#map_centre = (10, -30) +map_centre = (20, -30) +# Large problem size, full map. +#map_shape = (90, 48) +#map_spacing = .0627 +# Standard problem size. +#map_shape = (72, 38) +#map_spacing = .0627 +# Small problem size. +#map_shape = (44, 24) +#map_spacing = .0627 +# Tiny problem size. +map_shape = (125, 88) +map_spacing = .08 + +# Data paths +#raid_pro = os.getenv("RAID_PRO") +input_data_dir = data_dir +#base_dir = '/scratch/p/pen/andersoc/first_parkes_pipe/' +base_dir = '/scratch/p/pen/andersoc/second_parkes_pipe/' +map_dir = base_dir + 'maps/' +#map_dir = base_dir + 'parallel_maps/' +#map_root = map_dir + "parkes_test_parallel_thread_ra30decn30_small_" +#map_root = map_dir + "parkes_parallel_thread_ra20decn30_p08_125by88_beam1_" +map_root = map_dir + "parkes_parallel_thread_ra20decn30_p08_125by88_beam2_" +#map_root = map_dir + "parkes_parallel_thread_ra10decn30_p08_125by88_beam3_" +#map_root = map_dir + "parkes_parallel_thread_ra10decn30_p08_125by88_beam4_" +#map_root = map_dir + "parkes_parallel_thread_ra10decn30_p08_125by88_beam5_" +#map_root = map_dir + "parkes_parallel_thread_ra10decn30_p08_125by88_beam6_" +#map_root = map_dir + "parkes_parallel_thread_ra10decn30_p08_125by88_beam7_" +#map_root = map_dir + "parkes_parallel_thread_ra10decn30_p08_125by88_beam8_" +#map_root = map_dir + "parkes_parallel_thread_ra10decn30_p08_125by88_beam9_" +#map_root = map_dir + "parkes_parallel_thread_ra10decn30_p08_125by88_beam10_" +#map_root = map_dir + "parkes_parallel_thread_ra10decn30_p08_125by88_beam11_" +#map_root = map_dir + "parkes_parallel_thread_ra10decn30_p08_125by88_beam12_" +#map_root = map_dir + "parkes_parallel_thread_ra10decn30_p08_125by88_beam13_" + +# Map maker parameters. +#dm_input_root = '/scratch/p/pen/andersoc/first_parkes_pipe/rotated_to_I/' +dm_input_root = input_data_dir +dm_file_middles = file_middles +dm_input_end = '.fits' +dm_output_root = map_root +dm_scans = () +dm_IFs = (0,) +#dm_beams = (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12) +#dm_beam = 1 +dm_beam = 2 +#dm_beam = 3 +#dm_beam = 4 +#dm_beam = 5 +#dm_beam = 6 +#dm_beam = 7 +#dm_beam = 8 +#dm_beam = 9 +#dm_beam = 10 +#dm_beam = 11 +#dm_beam = 12 +#dm_beam = 13 + +dm_thread_divide = True +#dm_polarizations = ('XX', 'YY') +#dm_polarizations = ('XX',) +dm_polarizations = ('YY',) +dm_field_centre = map_centre +dm_pixel_spacing = map_spacing +dm_map_shape = map_shape +dm_time_block = 'scan' +dm_n_files_group = 10 +#dm_frequency_correlations = 'measured' +dm_frequency_correlations = 'None' +dm_number_frequency_modes = 0 +#dm_number_frequency_modes_discard = 1 +#dm_noise_parmeters_input_root = 'None' +dm_noise_parameter_file = '' +dm_deweight_time_mean = True +dm_deweight_time_slope = True +dm_interpolation = 'linear' +#dm_ts_foreground_mode_file = '' +#dm_n_ts_foreground_modes = 0 + +cm_input_root = dm_output_root +cm_output_root = cm_input_root +cm_polarizations = ('XX','YY') +#cm_polarizations = ('I',) +cm_save_noise_diag = True +cm_bands = (1316,) +#cm_bands = (1315,) From a9404bda9e9cc656963274d16cb1f6e54a448223 Mon Sep 17 00:00:00 2001 From: Christopher Anderson Date: Tue, 13 Jan 2015 15:58:37 -0500 Subject: [PATCH 27/28] Includes ra modulo fix that works for dirt_map and parallel_dirty_map. --- map/dirty_map.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/map/dirty_map.py b/map/dirty_map.py index 8db14f24..89864726 100644 --- a/map/dirty_map.py +++ b/map/dirty_map.py @@ -795,7 +795,8 @@ def preprocess_data(self, Blocks, map): mask[:,tmp_time_ind:tmp_time_ind + this_nt] = this_mask # Copy the other fields. Data.calc_pointing() - ra[tmp_time_ind:tmp_time_ind + this_nt] = Data.ra + map_centre_ra = np.mean(map.get_axis('ra')) + ra[tmp_time_ind:tmp_time_ind + this_nt] = ((Data.ra - map_centre_ra + 180) % 360) + map_centre_ra - 180 dec[tmp_time_ind:tmp_time_ind + this_nt] = Data.dec Data.calc_time() time[tmp_time_ind:tmp_time_ind + this_nt] = Data.time From 2a6a61e1f71d38762a3722542a48cec17428f0fb Mon Sep 17 00:00:00 2001 From: Christopher Anderson Date: Wed, 14 Jan 2015 23:46:42 -0500 Subject: [PATCH 28/28] Fixed memory leak in parallel dirty map writes over 2 GB per node. --- input/km/dm_parkes_thread2.ini | 18 ++++++++++-------- map/parallel_dirty_map.py | 22 ++++++++++++---------- 2 files changed, 22 insertions(+), 18 deletions(-) diff --git a/input/km/dm_parkes_thread2.ini b/input/km/dm_parkes_thread2.ini index e846378f..7a7b39a8 100755 --- a/input/km/dm_parkes_thread2.ini +++ b/input/km/dm_parkes_thread2.ini @@ -22,7 +22,8 @@ for el in os.walk(data_dir): #map_centre = (10, -30) -map_centre = (20, -30) +#map_centre = (20, -30) +map_centre = (-10, -30) # Large problem size, full map. #map_shape = (90, 48) #map_spacing = .0627 @@ -33,6 +34,7 @@ map_centre = (20, -30) #map_shape = (44, 24) #map_spacing = .0627 # Tiny problem size. +#map_shape = (375, 100) map_shape = (125, 88) map_spacing = .08 @@ -44,8 +46,8 @@ base_dir = '/scratch/p/pen/andersoc/second_parkes_pipe/' map_dir = base_dir + 'maps/' #map_dir = base_dir + 'parallel_maps/' #map_root = map_dir + "parkes_test_parallel_thread_ra30decn30_small_" -#map_root = map_dir + "parkes_parallel_thread_ra20decn30_p08_125by88_beam1_" -map_root = map_dir + "parkes_parallel_thread_ra20decn30_p08_125by88_beam2_" +#map_root = map_dir + "parkes_parallel_thread_ra20decn30_p08_375by100_beam1_" +#map_root = map_dir + "parkes_parallel_thread_ra20decn30_p08_125by88_beam2_" #map_root = map_dir + "parkes_parallel_thread_ra10decn30_p08_125by88_beam3_" #map_root = map_dir + "parkes_parallel_thread_ra10decn30_p08_125by88_beam4_" #map_root = map_dir + "parkes_parallel_thread_ra10decn30_p08_125by88_beam5_" @@ -53,7 +55,7 @@ map_root = map_dir + "parkes_parallel_thread_ra20decn30_p08_125by88_beam2_" #map_root = map_dir + "parkes_parallel_thread_ra10decn30_p08_125by88_beam7_" #map_root = map_dir + "parkes_parallel_thread_ra10decn30_p08_125by88_beam8_" #map_root = map_dir + "parkes_parallel_thread_ra10decn30_p08_125by88_beam9_" -#map_root = map_dir + "parkes_parallel_thread_ra10decn30_p08_125by88_beam10_" +map_root = map_dir + "parkes_parallel_thread_ra10decn30_p08_125by88_beam10_testsave_" #map_root = map_dir + "parkes_parallel_thread_ra10decn30_p08_125by88_beam11_" #map_root = map_dir + "parkes_parallel_thread_ra10decn30_p08_125by88_beam12_" #map_root = map_dir + "parkes_parallel_thread_ra10decn30_p08_125by88_beam13_" @@ -68,7 +70,7 @@ dm_scans = () dm_IFs = (0,) #dm_beams = (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12) #dm_beam = 1 -dm_beam = 2 +#dm_beam = 2 #dm_beam = 3 #dm_beam = 4 #dm_beam = 5 @@ -76,15 +78,15 @@ dm_beam = 2 #dm_beam = 7 #dm_beam = 8 #dm_beam = 9 -#dm_beam = 10 +dm_beam = 10 #dm_beam = 11 #dm_beam = 12 #dm_beam = 13 dm_thread_divide = True #dm_polarizations = ('XX', 'YY') -#dm_polarizations = ('XX',) -dm_polarizations = ('YY',) +dm_polarizations = ('XX',) +#dm_polarizations = ('YY',) dm_field_centre = map_centre dm_pixel_spacing = map_spacing dm_map_shape = map_shape diff --git a/map/parallel_dirty_map.py b/map/parallel_dirty_map.py index 2d695845..73958a6b 100644 --- a/map/parallel_dirty_map.py +++ b/map/parallel_dirty_map.py @@ -964,20 +964,22 @@ def lock_and_write_buffer(obj, fname, offset, size, proc): os.lseek(fd, offset, 0) - a=0 - while(True): - #while(len(buf)>0): - #nb = os.write(fd,buf) - nb = os.write(fd, buf[a:]) - print "The buffer save started at byte " + str(a) + " for process " + str(proc) + nb = 0 + #a=0 + #while(True): + while(len(buf)>0): + nb += os.write(fd,buf) + #nb = os.write(fd, buf[a:]) + #print "The buffer save started at byte " + str(a) + " for process " + str(proc) if nb < 0: raise Exception("Failed write") - if nb == len(buf[a:]): - break - else: - a += nb + #if nb == len(buf[a:]): + # break + #else: + # a += nb #buf = buf[nb:] + buf = buffer(obj, nb) ''' nb = os.write(fd, buf)