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/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..7a7b39a8
--- /dev/null
+++ b/input/km/dm_parkes_thread2.ini
@@ -0,0 +1,113 @@
+# 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)
+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 = (375, 100)
+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_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_"
+#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_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_"
+
+# 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,)
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.ini b/input/km/mm_test.ini
index 8f7c695c..0fa7ea72 100644
--- a/input/km/mm_test.ini
+++ b/input/km/mm_test.ini
@@ -32,11 +32,13 @@ 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_"
 
 # 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/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/input/tcv/mm_test.ini b/input/tcv/mm_test.ini
index cdd8a925..e86addd0 100644
--- a/input/tcv/mm_test.ini
+++ b/input/tcv/mm_test.ini
@@ -3,14 +3,48 @@
 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'))
-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)
+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.
-file_middles = file_middles[:10]
+#file_middles = file_middles[:10]
 # file_middles = file_middles[:1]
 
 map_centre = (13.0, 1.85)
@@ -34,8 +68,10 @@ 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_base = base_dir+'maps/1hr_parallel_test/'
+map_prefix = '_1hr_80-28_ptcorr_'
 
 # Map maker parameters.
 dm_input_root = input_data_dir 
@@ -49,7 +85,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
@@ -60,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 90333883..a7f37876 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.
@@ -27,11 +28,21 @@ field = '15hr'
 #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):
@@ -67,11 +78,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 +99,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-28_ptcorr' + '_'
 
 # Maximum number of processes to use.
 pipe_processes = 1
@@ -137,8 +148,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 ####
@@ -152,11 +163,11 @@ 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_')))
@@ -321,7 +332,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 +354,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 #####
 
@@ -413,7 +426,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'
diff --git a/map/_mapmaker.pyx b/map/_mapmaker.pyx
index 6ff00540..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 = <DTYPE_t *> 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 = <np.npy_intp *> 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, 
@@ -194,18 +245,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 +309,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 +339,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.
diff --git a/map/clean_map.py b/map/clean_map.py
index 509eba23..405119a1 100644
--- a/map/clean_map.py
+++ b/map/clean_map.py
@@ -24,7 +24,8 @@
                'save_noise_inv_diag' : False,
                'save_cholesky' : False,
                'from_eig' : False,
-               'bands' : ()
+               'bands' : (),
+               'mem_lim' : ()
                }
 prefix = 'cm_'
 
@@ -130,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(noise_h5, 'inv_cov')
-                        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: "
@@ -154,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.
@@ -233,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/dirty_map.py b/map/dirty_map.py
index c6896d0a..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
@@ -1195,7 +1196,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
@@ -1206,7 +1207,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_clean_map.py b/map/parallel_clean_map.py
new file mode 100644
index 00000000..a4cf1626
--- /dev/null
+++ b/map/parallel_clean_map.py
@@ -0,0 +1,602 @@
+"""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 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()
+
+    @profile
+    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_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:
+                                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()
+                        for process_n in range(self.nproc):
+                            for ii in index_list_full[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.
+                            #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
+                    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')
+                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(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)
diff --git a/map/parallel_dirty_map.py b/map/parallel_dirty_map.py
index eecae4c8..73958a6b 100644
--- a/map/parallel_dirty_map.py
+++ b/map/parallel_dirty_map.py
@@ -5,8 +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 sys
+from memory_profiler import profile
 
 from mpi4py import MPI
 
@@ -86,7 +88,10 @@
                # 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,
+               # Which beam to use to make map.  Choose zero to use all beams.
+               'beam' : 0
                }
 
 comm = MPI.COMM_WORLD
@@ -123,9 +128,14 @@ 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"]]), middle
             else:
                 msg = "time_block parameter must be 'scan' or 'file'."
                 raise dm.ValueError(msg)
@@ -168,6 +178,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"])
@@ -270,10 +281,12 @@ 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()
     
+    @profile    
     def make_map(self):
         """Makes map for current polarization and band.
         
@@ -524,7 +537,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]
@@ -567,6 +580,10 @@ def thread_work():
                 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=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=dtype)
                 index_queue = Queue()
@@ -577,35 +594,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=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_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
-                            del 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=dtype)
+                            #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_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
-                            del 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):
@@ -615,7 +638,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:
@@ -748,14 +779,26 @@ 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)
-            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)
+            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)
+            if self.rank == 0:
+                f = h5py.File(self.cov_filename, 'r+')
+                cov_inv_dset = f['inv_cov']
+                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
@@ -771,15 +814,28 @@ 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)
-            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)
+            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)
+            if self.rank == 0:
+                f = h5py.File(self.cov_filename, 'r+')
+                cov_inv_dset = f['inv_cov']
+                attrs = cov_inv_dset.attrs
+                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', 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))
+                #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
@@ -871,7 +927,8 @@ def cross(set_list):
                 remaining.insert(0,cross_2)
                 return cross(remaining)
 
-def lock_and_write_buffer(obj, fname, offset, size):
+@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.
 
@@ -906,12 +963,30 @@ def lock_and_write_buffer(obj, fname, offset, size):
         print "Could not obtain lock"
 
     os.lseek(fd, offset, 0)
-
+    
+    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
+        #buf = buf[nb:]
+        buf = buffer(obj, 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: