moved subtraction/datalike to pwr operations; working autopower pipel…

…ine; cleaning up transfer pipeline

input/ers/autopower/GBT_15hr_map_oldcal.ini

@@ -16,13 +16,12 @@ import os
 #basemap = 'GBT_15hr_map_oldcalpol'
 basemap = 'GBT_15hr_map_oldcal'
 # and an identifer for this run (any keywords)
-#output_tag = basemap + "_" + "order1_ext"
 output_tag = basemap + "_" + "order1"
 
 # control the main operations
-do_cleaning = False
+do_cleaning = True
 do_ext_cleaning = False
-do_power = False
+do_power = True
 do_power_compile = True
 do_analysis = True
 
@@ -63,8 +62,8 @@ else:
 pwrout_analysis_plots = "./pwrspec_plots/" + output_tag + "_analysis/"
 
 trans_root = "/mnt/raid-project/gmrt/eswitzer/GBT/pwrspec/"
-trans_beam = trans_root + "GBT_15hr_map_oldcal_transfer_function_beamtransfer.hd5"
-trans_mode = trans_root + "GBT_15hr_map_oldcal_transfer_function_modetransfer.hd5"
+trans_beam = trans_root + basemap + "_transfer_function_beamtransfer.hd5"
+trans_mode = trans_root + basemap + "_transfer_function_modetransfer.hd5"
 #trans_root = "/mnt/raid-project/gmrt/eswitzer/GBT/bulksim/"
 #trans_beam = trans_root + "GBT_15hr_map_oldcal_blackman_order1_one-sided_beamtransfer.hd5"
 #trans_mode = trans_root + "GBT_15hr_map_oldcal_blackman_order1_one-sided_modetransfer.hd5"
@@ -128,7 +127,7 @@ xs1_order = pwr_order
 xs1_bins = pwr_bins
 xs1_freq_list = freq_list
 xs1_map_key = basemap + "_cleaned"
-xs1_ncpu = 6
+xs1_ncpu = 24
 xs1_outfile = pwrout_base + basemap + ".shelve"
 
 if do_power:
@@ -146,7 +145,7 @@ ns1_order = pwr_order
 ns1_bins = pwr_bins
 ns1_freq_list = freq_list
 ns1_map_key = basemap + "_cleaned"
-ns1_ncpu = 6
+ns1_ncpu = 24
 ns1_outfile = pwrout_base + basemap + "_noise.shelve"
 
 #-----------------------------------------------------------------------------
@@ -200,8 +199,8 @@ if not os.path.isdir(pwrout_plots):
 
 if do_analysis:
     pipe_modules.append(auto_pwrspec_analyze.AnalyzeAutopower)
-analyzeautopower_sim_auto_summary = pwrout_root + "GBT_15hr_map_oldcal_noise_simulations_auto.shelve"
-analyzeautopower_sim_xspec_summary = pwrout_root + "GBT_15hr_map_oldcal_noise_simulations_xspec.shelve"
+analyzeautopower_sim_auto_summary = pwrout_root + "GBT_15hr_map_oldcal_noise_simulations_stat_auto.hd5"
+analyzeautopower_sim_xspec_summary = pwrout_root + "GBT_15hr_map_oldcal_noise_simulations_stat_xspec.hd5"
 analyzeautopower_data_auto_summary = pwrout_base + basemap + "_self.hd5"
 analyzeautopower_data_xspec_summary = pwrout_base + basemap + "_signal.hd5"
 analyzeautopower_outdir = pwrout_analysis_plots

input/ers/autopower/GBT_15hr_map_oldcal_errors.ini

@@ -20,7 +20,7 @@ fix_weight_treatment = None
 outplotdir = "/cita/h/home-2/eswitzer/code/analysis_IM/pwrspec_plots/"
 #outplotdir = "/home/r/rbond/eswitzer/code/analysis_IM/pwrspec_plots/"
 
-do_compile = False
+do_compile = True
 do_statistics = True
 
 # the rest is fairly automated and should not need to change

input/ers/autopower/GBT_15hr_map_oldcal_transfer.ini

@@ -4,7 +4,9 @@ pipe_processes = 1
 from quadratic_products import aggregate_bulksim
 import os
 
-basemap = "GBT_15hr_map_oldcal"
+basemap = "GBT_15hr_map_oldcalpol"
+#basemap = "GBT_15hr_map_oldcal"
+wbeam_basemap = "GBT_15hr_map_oldcal"
 pwrout_root = '/mnt/raid-project/gmrt/eswitzer/GBT/pwrspec/'
 output_tag = basemap + "_" + "transfer_function"
 
@@ -19,8 +21,8 @@ outplotdir = "/cita/h/home-2/eswitzer/code/analysis_IM/pwrspec_plots/"
 
 do_compile_phys = False
 do_compile_wbeam = False
-do_compile_mode = False
-do_statistics = False
+do_compile_mode = True
+do_statistics = True
 do_transfer = True
 
 # the rest is fairly automated and should not need to change
@@ -45,8 +47,8 @@ as2_apply_2d_beamtransfer = None
 as2_apply_2d_modetransfer = None
 as2_noiseweights_2dto1d = None
 as2_fix_weight_treatment = None
-as2_directory = pwrout_root + basemap + "_" + "noise_simulations_xspec"
-as2_basefile = "%s_sim_" % (basemap + "_" + "noise_simulations")
+as2_directory = pwrout_root + wbeam_basemap + "_" + "noise_simulations_xspec"
+as2_basefile = "%s_sim_" % (wbeam_basemap + "_" + "noise_simulations")
 as2_outfile = "%s/%s_xspec_wbeam.hd5" % (pwrout_root, output_tag)
 
 if do_compile_mode:
@@ -59,6 +61,17 @@ as3_directory = "%s_modeclean_plussim" % modeout_base
 as3_basefile = "%s_sim_" % modeoutput_tag
 as3_outfile = "%s/%s_plussim.hd5" % (pwrout_root, output_tag)
 
+if do_compile_mode:
+    pipe_modules.append((aggregate_bulksim.AggregateSummary, ('as4_', 'as_')))
+as4_apply_2d_beamtransfer = None
+as4_apply_2d_modetransfer = None
+as4_noiseweights_2dto1d = None
+as4_fix_weight_treatment = None
+as4_subtract_pwrspec = sigpwr_base + basemap + ".shelve"
+as4_directory = "%s_modeclean_plussim" % modeout_base
+as4_basefile = "%s_sim_" % modeoutput_tag
+as4_outfile = "%s/%s_modesim.hd5" % (pwrout_root, output_tag)
+
 #----------------------------------------------------------------------------
 # then find some statistics on the above sets
 #----------------------------------------------------------------------------
@@ -83,15 +96,15 @@ ast3_aggfile_in = as3_outfile
 ast3_statfile_out = "%s/%s_stat_plussim.hd5" % (pwrout_root, output_tag)
 ast3_outplotdir = "%s/%s/sim_plussim" % (outplotdir, output_tag)
 
+if do_statistics:
+    pipe_modules.append((aggregate_bulksim.AggregateStatistics, ('ast4_', 'ast_')))
+ast4_aggfile_in = as4_outfile
+ast4_statfile_out = "%s/%s_stat_modesim.hd5" % (pwrout_root, output_tag)
+ast4_outplotdir = "%s/%s/sim_modesim" % (outplotdir, output_tag)
+
 #----------------------------------------------------------------------------
 # find the transfer functions
 #----------------------------------------------------------------------------
-if do_transfer:
-    pipe_modules.append((aggregate_bulksim.SubtractMap, ('sub1_', 'sub_')))
-sub1_plussim_file = ast3_statfile_out
-sub1_mappower_file = sigpwr_base + basemap + ".shelve"
-sub1_output_file = "%s/%s_stat_modes.hd5" % (pwrout_root, output_tag)
-
 if do_transfer:
     pipe_modules.append((aggregate_bulksim.CalculateTransfer, ('atr1_', 'atr_')))
 atr1_powerfile_in = ast1_statfile_out
@@ -102,6 +115,7 @@ atr1_outplotdir = "%s/%s/transfer" % (outplotdir, output_tag)
 if do_transfer:
     pipe_modules.append((aggregate_bulksim.CalculateTransfer, ('atr2_', 'atr_')))
 atr2_powerfile_in = ast2_statfile_out
-atr2_powerfile_out = sub1_output_file
+#atr2_powerfile_out = sub1_output_file
+atr2_powerfile_out = ast4_statfile_out
 atr2_transferfile = "%s/%s_modetransfer.hd5" % (pwrout_root, output_tag)
 atr2_outplotdir = "%s/%s/modetransfer" % (outplotdir, output_tag)

quadratic_products/aggregate_bulksim.py

@@ -22,6 +22,7 @@
         "fix_weight_treatment": None,
         "apply_2d_beamtransfer": None,
         "apply_2d_modetransfer": None,
+        "subtract_pwrspec": None,
         "outfile": "file"
     }
 
@@ -57,6 +58,11 @@ def produce_summary(self, filelist, outfile, debug=False):
         num_sim = len(filelist)
         print "number of simulations to aggregate: %d" % num_sim
 
+        data_subtract = None
+        if self.params['subtract_pwrspec'] is not None:
+            print "agg WARNING: you are subtracting a power spectrum"
+            data_subtract = ps.PowerSpectrum(self.params['subtract_pwrspec'])
+
         sim_toread = ps.PowerSpectrum(filelist[0])
         # run this here just to get the 2D->1D bin centers
         sim_toread.convert_2d_to_1d()
@@ -121,6 +127,18 @@ def produce_summary(self, filelist, outfile, debug=False):
         for (simfile, index) in zip(filelist, range(num_sim)):
             print "processing ", simfile, index
             sim_toread = ps.PowerSpectrum(simfile)
+
+            # optionally subtract some reference spectrum
+            if data_subtract is not None:
+                print "agg WARNING: you are subtracting a power spectrum"
+                # assuming these both have the same treatments
+                for pwrspec_case in sim_toread.pwrspec_1d:
+                    sim_toread.pwrspec_1d[pwrspec_case] -= \
+                        data_subtract.pwrspec_1d[pwrspec_case]
+
+                    sim_toread.pwrspec_2d[pwrspec_case] -= \
+                        data_subtract.pwrspec_2d[pwrspec_case]
+
             sim_toread.apply_2d_trans_by_treatment(transfer_dict)
             sim_toread.convert_2d_to_1d(weights_2d=weights_2d)
 
@@ -322,6 +340,18 @@ def calc_stat_2d(self, trial_array, calc_corr=False):
         old_settings = np.seterr(invalid="warn", under="warn")
         mtrial_array = ma.masked_invalid(trial_array)
 
+        hcheck = []
+        for index in range(mtrial_array.shape[0]):
+            hcheck.append(np.ma.sum(mtrial_array[index, :, :]))
+
+        hcheck = np.array(hcheck)
+        meanval = np.mean(hcheck)
+        hcheck = (hcheck - meanval) / meanval
+        if not np.all(hcheck == 0.):
+            print hcheck
+
+        #print np.apply_over_axes(np.ma.sum, mtrial_array, axes=[1,2])
+
         stat_2d['mean'] = np.ma.filled(np.ma.mean(mtrial_array, axis=0),
                                        fill_value=np.nan)
 
@@ -384,52 +414,6 @@ def aggregate_2d_statistics(self, treatment):
         return (stat_2d, counts_2d)
 
 
-subtractmap_init = {
-        "plussim_file": "file",
-        "mappower_file": "file",
-        "output_file": "file"
-    }
-subtractmap_prefix = 'sub_'
-
-
-class SubtractMap(object):
-    """Input: clean(map+sim)xclean(map+sim), clean(map)xclean(map)
-    subtract these for the purpose of finding the transfer funcion
-    """
-
-    def __init__(self, parameter_file=None, params_dict=None, feedback=0):
-        self.params = params_dict
-        np.seterr(invalid='raise')
-
-        if parameter_file:
-            self.params = parse_ini.parse(parameter_file,
-                                          subtractmap_init,
-                                          prefix=subtractmap_prefix)
-
-        print self.params["plussim_file"], "-", self.params["mappower_file"]
-        print "writing to ", self.params["output_file"]
-
-    def execute(self, processes):
-        shutil.copyfile(self.params["plussim_file"], self.params["output_file"])
-        outfile = h5py.File(self.params["output_file"], "r+")
-        pwr_map = ps.PowerSpectrum(self.params["mappower_file"])
-
-        signal2d_agg = pwr_map.agg_stat_2d_pwrspec()
-        print signal2d_agg.keys()
-
-        for treatment in outfile:
-            plussim_power = \
-                copy.deepcopy(outfile[treatment]['pk_2d_stat']['mean'].value)
-
-            map_power = signal2d_agg[treatment]["mean"]
-            del outfile[treatment]['pk_2d_stat']['mean']
-
-            outfile[treatment]['pk_2d_stat']['mean'] = plussim_power - \
-                                                       map_power
-
-        outfile.close()
-
-
 calculatetransfer_init = {
         "powerfile_in": "file",
         "powerfile_out": "file",
@@ -507,74 +491,3 @@ def execute(self, processes):
         self.stats_out.close()
 
 
-calculatedatalike_init = {
-        "powerfile": "file",
-        "powerdatalike_out": "file",
-    }
-calculatedatalike_prefix = 'cdl_'
-
-
-class CalculateDatalike(object):
-    """Reformat aggregated simulations to look like data (vestigial?)
-    """
-
-    def __init__(self, parameter_file=None, params_dict=None, feedback=0):
-        self.params = params_dict
-        np.seterr(invalid='raise')
-
-        if parameter_file:
-            self.params = parse_ini.parse(parameter_file,
-                                          calculatedatalike_init,
-                                          prefix=calculatedatalike_prefix)
-
-        print self.params["powerfile_in"], "->", self.params["powerdatalike_out"]
-
-        self.stats_in = h5py.File(self.params["powerfile"], "r")
-        self.treatments_in = self.stats_in["results"].keys()
-
-        # maybe make this hd5?
-        self.stats_dataout = shelve.open(self.params["powerdatalike_out"], "w")
-        #file_tools.convert_numpytree_hdf5(out_dicttree, outfile)
-
-    def execute(self, processes):
-        """produce a list of files to combine and run"""
-        for treatment in self.treatments_in:
-            print treatment
-            print self.stats_in["results"][treatment].keys()
-
-            stat_in = self.stats_in["stats"]
-
-            stats_1d_in = stat_in["0modes"]["pk_1d_stat"]["mean"]
-            stats_2d_in = stat_in["0modes"]["pk_2d_stat"]["mean"]
-            counts_1d_in = stat_in["0modes"]["pk_1d_counts"]["mean"]
-            counts_2d_in = stat_in["0modes"]["pk_2d_counts"]["mean"]
-
-            k_1d = self.stats_in["k_1d"]
-            k_1d_from_2d = self.stats_in["k_1d_from_2d"]
-            kx_2d = self.stats_in["kx_2d"]
-            ky_2d = self.stats_in["ky_2d"]
-
-            # make a package that looks like the data
-            pwrspec2d_product = {}
-            pwrspec2d_product["bin_x_left"] = kx_2d["left"]
-            pwrspec2d_product["bin_x_center"] = kx_2d["center"]
-            pwrspec2d_product["bin_x_right"] = kx_2d["right"]
-            pwrspec2d_product["bin_y_left"] = ky_2d["left"]
-            pwrspec2d_product["bin_y_center"] = ky_2d["center"]
-            pwrspec2d_product["bin_y_right"] = ky_2d["right"]
-            pwrspec2d_product["counts_histo"] = counts_2d_in
-            pwrspec2d_product["binavg"] = stats_2d_in
-
-            pwrspec1d_product = {}
-            pwrspec1d_product["bin_left"] = k_1d["left"]
-            pwrspec1d_product["bin_center"] = k_1d["center"]
-            pwrspec1d_product["bin_right"] = k_1d["right"]
-            pwrspec1d_product["counts_histo"] = counts_1d_in
-            pwrspec1d_product["binavg"] = stats_1d_in
-
-            datakey = "data:%s" % treatment
-            self.stats_dataout[datakey] = (0, (pwrspec2d_product,
-                                               pwrspec1d_product))
-
-        self.stats_in.close()
-        self.stats_dataout.close()

quadratic_products/auto_pwrspec_analyze.py

@@ -32,18 +32,21 @@ def __init__(self, parameter_file=None, params_dict=None, feedback=0):
         print self.params
         self.data_auto = h5py.File(self.params["data_auto_summary"], "r")
         self.data_xspec = h5py.File(self.params["data_xspec_summary"], "r")
-        self.sim_auto = shelve.open(self.params["sim_auto_summary"], "r")
-        self.sim_xspec = shelve.open(self.params["sim_xspec_summary"], "r")
+        self.sim_auto = h5py.File(self.params["sim_auto_summary"], "r")
+        self.sim_xspec = h5py.File(self.params["sim_xspec_summary"], "r")
 
     def execute(self, processes):
-        sim_auto_stats = self.sim_auto['stats']['0modes']['pk_1d_from_2d_stat']
-        sim_xspec_stats = self.sim_xspec['stats']['0modes']['pk_1d_from_2d_stat']
+        #sim_auto_stats = self.sim_auto['stats']['0modes']['pk_1d_from_2d_stat']
+        #sim_xspec_stats = self.sim_xspec['stats']['0modes']['pk_1d_from_2d_stat']
+        sim_auto_stats = self.sim_auto['0modes']['pk_1d_from_2d_stat']
+        sim_xspec_stats = self.sim_xspec['0modes']['pk_1d_from_2d_stat']
 
-        sim_xspec_mean = sim_xspec_stats['mean']
-        sim_xspec_cov = sim_xspec_stats['cov']
-        sim_auto_mean = sim_auto_stats['mean']
+        sim_xspec_mean = sim_xspec_stats['mean'].value
+        sim_xspec_cov = sim_xspec_stats['cov'].value
+        sim_auto_mean = sim_auto_stats['mean'].value
 
-        k_vec_sim = self.sim_auto['k_1d_from_2d']['center']
+        #k_vec_sim = self.sim_auto['k_1d_from_2d']['center']
+        k_vec_sim = self.sim_auto['0modes']['k_1d_from_2d'].value
         k_vec = self.data_xspec['0modes']['k_vec'].value
 
         assert np.array_equal(k_vec_sim, k_vec), "simulation and data unaligned kvecs"