Try solving merge conflicts

Michael Mayer
2 parents 45a80b60 b7dd76b3
Showing 8 changed files with 212 additions and 28 deletions
cscripts/cstsmapmerge.py
cscripts/cstsmapsplit.par
src/cttsmap/cttsmap.cpp
src/cttsmap/cttsmap.hpp
src/cttsmap/cttsmap.par
test/data/tsmap_0.fits
test/data/tsmap_1.fits
test/test_cttsmap.py
@@ -46,7 +46,7 @@ class cstsmapmerge(ctools.cscript):
         self._files        = None
         self._in_filename  = ''
         self._tsmap        = gammalib.GSkyMap() # Empty sky map
-        self._statusmap    = gammalib.GSkyMap() # Empty sky map
+        self._statusmap  = gammalib.GSkyMap() # Empty sky map
         self._maps         = []
         self._mapnames     = []
         self._merged_files = []
@@ -118,7 +118,11 @@ class cstsmapmerge(ctools.cscript):
         if self._read_ahead():
             self['outmap'].filename()
  
+<<<<<<< HEAD
         #  Write input parameters into logger
+=======
+        # Write input parameters into logger
+>>>>>>> b7dd76b35e6d2fd74ff010258ffa7606e072dbf5
         self._log_parameters(gammalib.TERSE)
  
         # Return
@@ -153,7 +157,11 @@ class cstsmapmerge(ctools.cscript):
         for hdu in fits:
  
             # Leave out primary and status extension
+<<<<<<< HEAD
             if hdu.extname() != 'IMAGE' and hdu.extname() != 'STATUS MAP':
+=======
+            if hdu.extname() != "IMAGE" and hdu.extname() != 'STATUS MAP':
+>>>>>>> b7dd76b35e6d2fd74ff010258ffa7606e072dbf5
  
                 # Add present maps
                 skymap = gammalib.GSkyMap()
@@ -192,7 +200,11 @@ class cstsmapmerge(ctools.cscript):
         for hdu in fits:
  
             # Leave out primary and status extension
+<<<<<<< HEAD
             if hdu.extname() != 'IMAGE' and hdu.extname() != 'STATUS MAP':
+=======
+            if hdu.extname() != "IMAGE" and hdu.extname() != 'STATUS MAP':
+>>>>>>> b7dd76b35e6d2fd74ff010258ffa7606e072dbf5
  
                 # Add present maps
                 skymap = gammalib.GSkyMap()
@@ -213,14 +225,22 @@ class cstsmapmerge(ctools.cscript):
         for i in range(self._tsmap.npix()):
  
             # Consider only bins that have been computed
-            if add_statusmap[i] > 0.5:
+            if add_statusmap[i] > -0.5:
  
                 # Raise exception if this bin has already been computed
+<<<<<<< HEAD
                 if self._statusmap[i] > 0.5 and not self._overwrite:
                     msg = ('Attempt to merge bin which apparently has already '
                            'been merged. File "%s" contains already merged '
                            'bins. Set hidden parameter "overwrite=yes" to '
                            'avoid this error.' % fitsfile)
+=======
+                if self._statusmap[i] > -0.5 and not self._overwrite:
+                    msg = "Attempt to merge bin which apparently has "+\
+                          "already been merged. File \""+fitsfile+"\" "+\
+                          "contains already merged bins. Set hidden "+\
+                          "parameter \"overwrite=yes\" to avoid this error."
+>>>>>>> b7dd76b35e6d2fd74ff010258ffa7606e072dbf5
                     raise RuntimeError(msg)
  
                 # Copy TS values
@@ -256,7 +276,7 @@ class cstsmapmerge(ctools.cscript):
         # Count number of pixels with TS status set
         count = 0
         for pix in status:
-            if pix > 0.5:
+            if pix > -0.5:
                 count += 1
  
         # Return number of pixels
@@ -290,8 +310,16 @@ class cstsmapmerge(ctools.cscript):
  
             # Skip file if it's not a FITS file
             if not gammalib.GFilename(fitsfile).is_fits():
+<<<<<<< HEAD
                 self._log_value(gammalib.EXPLICIT, 'Skip file',
                                 '%s (not a FITS file)' % fitsfile)
+=======
+                
+                # Log message
+                self._log_value(gammalib.EXPLICIT, "Skip file", fitsfile +" (not a FITS file)")
+
+                # Continue
+>>>>>>> b7dd76b35e6d2fd74ff010258ffa7606e072dbf5
                 continue
  
             # Open FITS file
@@ -299,18 +327,35 @@ class cstsmapmerge(ctools.cscript):
  
             # If file contains a status map then use it
             if fits.contains('STATUS MAP'):
+<<<<<<< HEAD
+=======
+                
+                # Close FITS file
+>>>>>>> b7dd76b35e6d2fd74ff010258ffa7606e072dbf5
                 fits.close()
+                
+                # Set filename to initial filename
                 file0 = fitsfile
+<<<<<<< HEAD
                 if self._logTerse():
                     self._log(gammalib.parformat('Initial TS map file'))
                     self._log(fitsfile)
                     count = self._get_number_of_ts_pixels(fitsfile)
                     self._log(' (%d TS pixels computed)' % count)
                     self._log('\n')
+=======
+                
+                # Log message
+                pix_info = " (%d TS pixels computed)" % self._get_number_of_ts_pixels(fitsfile)
+                self._log_value(gammalib.TERSE, "Initial TS map file", fitsfile + pix_info)
+                
+                # Leave loop
+>>>>>>> b7dd76b35e6d2fd74ff010258ffa7606e072dbf5
                 break
  
             # ... otherwise signal that file is useless
             else:
+<<<<<<< HEAD
                 fits.close()
                 self._log_value(gammalib.EXPLICIT, 'Skip file',
                                 '%s (no "STATUS MAP" extension)' % fitsfile)
@@ -322,6 +367,24 @@ class cstsmapmerge(ctools.cscript):
                 self._log('None of the provided files seems to be a sliced '
                           'TS map file (none has a "STATUS MAP" extension).\n')
  
+=======
+                
+                # Close fits file
+                fits.close()  
+                
+                # Info message             
+                self._log_value(gammalib.EXPLICIT, "Skip file", fitsfile +" (no \"STATUS MAP\" extension)")
+                
+                # Continue
+                continue
+
+        # Signal if no suitable file was found
+        if file0 == "":
+            self._log_string(gammalib.TERSE, "None of the provided files seems to be a sliced "
+                          "TS map file (none has a \"STATUS MAP\" "
+                          "extension).\n")
+
+>>>>>>> b7dd76b35e6d2fd74ff010258ffa7606e072dbf5
         # ... otherwise merge files
         else:
  
@@ -342,9 +405,14 @@ class cstsmapmerge(ctools.cscript):
  
                 # Skip if file is not FITS
                 if not gammalib.GFilename(fitsfile).is_fits():
+<<<<<<< HEAD
                     self._log_value(gammalib.EXPLICIT, 'Skip file',
                                     '%s (not a FITS file)' % fitsfile)
                     continue
+=======
+                    self._log_value(gammalib.EXPLICIT, "Skip file", fitsfile +" (not a FITS file)")
+                    continue   
+>>>>>>> b7dd76b35e6d2fd74ff010258ffa7606e072dbf5
  
                 # Open FITS file
                 fits = gammalib.GFits(fitsfile)
@@ -352,20 +420,29 @@ class cstsmapmerge(ctools.cscript):
                 # Skip if file does not contain status map
                 if not fits.contains('STATUS MAP'):
                     fits.close()
+<<<<<<< HEAD
                     self._log_value(gammalib.EXPLICIT, 'Skip file',
                                     '%s (no "STATUS MAP" extension)' % fitsfile)
+=======
+                    self._log_value(gammalib.EXPLICIT, "Skip file", fitsfile +" (no \"STATUS MAP\" extension)")
+>>>>>>> b7dd76b35e6d2fd74ff010258ffa7606e072dbf5
                     continue
  
                 # Close FITS file
                 fits.close()
  
                 # Logging
+<<<<<<< HEAD
                 if self._logTerse():
                     self._log(gammalib.parformat('Merge TS map file'))
                     self._log(fitsfile)
                     count = self._get_number_of_ts_pixels(fitsfile)
                     self._log(' (%d TS pixels computed)' % count)
                     self._log('\n')
+=======
+                pix_info = " (%d TS pixels computed)" % self._get_number_of_ts_pixels(fitsfile)
+                self._log_value(gammalib.TERSE, "Merge TS map file", fitsfile + pix_info)
+>>>>>>> b7dd76b35e6d2fd74ff010258ffa7606e072dbf5
  
                 # Merge TS map
                 self._merge_ts_map(fitsfile)
@@ -380,6 +457,7 @@ class cstsmapmerge(ctools.cscript):
         """ 
         Save TS map and remove slices if requested
         """
+<<<<<<< HEAD
         # Write header
         self._log_header1(gammalib.TERSE, 'Save TS map')
  
@@ -388,6 +466,11 @@ class cstsmapmerge(ctools.cscript):
  
         # Log filename
         self._log_value(gammalib.TERSE, 'TS map file', outmap.url())
+=======
+
+        # Get output filename in case it was not read ahead
+        outmap = self["outmap"].filename()
+>>>>>>> b7dd76b35e6d2fd74ff010258ffa7606e072dbf5
  
         # Create FITS file
         fits = gammalib.GFits()
@@ -403,26 +486,67 @@ class cstsmapmerge(ctools.cscript):
         for i in range(len(self._mapnames)):   
             fits[i+1].extname(self._mapnames[i])
  
+        # Treat status map separately 
+        self._statusmap.write(fits)
+        fits[fits.size()-1].extname('STATUS MAP')
+        
         # Check if map is fully done
-        done = True
+        bins_merged = 0
         for pix in self._statusmap:
-            if pix < 0.5:
-                done = False
-                break
+            if pix > -0.5:
+                bins_merged += 1
+
+        # Set boolean if map is fully done
+        done = (bins_merged == self._tsmaps.npix())
  
+        # Log summary header
+        self._log_header2(gammalib.TERSE, "Merging Summary")
+                
         # Write status map if we are not done yet
         if not done:
+<<<<<<< HEAD
             self._statusmap.write(fits)
             fits[fits.size()-1].extname('STATUS MAP')
+=======
+            
+            # Log bins that were computed and those that were missing
+            self._log_value(gammalib.TERSE, "TS map bins", self._tsmap.npix())
+            self._log_value(gammalib.TERSE, "Bins computed and merged", bins_merged )
+            self._log_value(gammalib.TERSE, "Bins missing", self._tsmap.npix() - bins_merged)
+        
+        else:
+            
+            # Write success message
+            self._log_string(gammalib.TERSE, "TS map was fully computed and successfully merged") 
+        
+        # Write header
+        self._log_header1(gammalib.TERSE, "Save TS map")
+        
+        # Log filename
+        self._log_value(gammalib.TERSE, "TS map file", outmap.url())
+>>>>>>> b7dd76b35e6d2fd74ff010258ffa7606e072dbf5
  
         # Save FITS file
-        fits.saveto(outmap, self._clobber())
+        fits.saveto(outmap, self._clobber()) 
  
-        # Delete TS input maps if requested
-        if self._delete:
+        # Delete TS input maps if requested and map is fully done
+        if self._delete and done:
+            
+            # Log header
+            self._log_header2(gammalib.TERSE, "Delete slices TS map files")
+            
+            # Loop over sliced files
             for filename in self._merged_files:
+                
+                # Remove file
                 os.remove(filename)
+<<<<<<< HEAD
                 self._log_value(gammalib.TERSE, 'Deleted input file', filename)
+=======
+                
+                # Log message about file
+                self._log_value(gammalib.TERSE,"Deleted input file", filename)
+>>>>>>> b7dd76b35e6d2fd74ff010258ffa7606e072dbf5
  
         # Return
         return
@@ -27,6 +27,7 @@ bkgcube,       f, a, NONE,,, &quot;Input background cube file (only needed for stacke
 caldb,         s, a, prod2,,, "Calibration database"
 irf,           s, a, South_0.5h,,, "Instrument response function"
 edisp,         b, h, no,,, "Apply energy dispersion?"
+errors,        b, h, no,,, "Compute/store parameter errors?"
 outmap,        f, a, tsmap.fits,,, "Output Test Statistic map file"
 bins_per_job,  i, a, 5,,, "Number of TS map bins per task"
 compute_null,  b, a, yes,,,"Compute null hypothesis first?"
@@ -215,6 +215,12 @@ void cttsmap::run(void)
     // Get task parameters
     get_parameters();
  
+<<<<<<< HEAD
+=======
+    // Write parameters into logger
+    log_parameters(TERSE);
+
+>>>>>>> b7dd76b35e6d2fd74ff010258ffa7606e072dbf5
     // Set energy dispersion flag for all CTA observations and save old
     // values in save_edisp vector
     std::vector<bool> save_edisp;
@@ -315,6 +321,14 @@ void cttsmap::run(void)
     	// Optimize observation container
     	m_obs.optimize(*opt);
  
+    	// Compute errors if necessary
+    	if (m_errors) {
+    		m_obs.errors(*opt);
+    	}
+
+    	// Get status of optimization
+    	int status = opt->status();
+
     	// Retrieve the Likelihood value
     	double logL1 = -(opt->value());
  
@@ -326,6 +340,7 @@ void cttsmap::run(void)
     		log << " TS value ....: ";
             log << ts << std::endl;
     	}
+
         else if (logTerse()) {
     		log << "TS for bin number ";
             log << i;
@@ -345,12 +360,23 @@ void cttsmap::run(void)
     	// Extract fitted test source parameters
         if (m_testsource != NULL) {
             for (int j = 0; j < m_mapnames.size(); ++j) {
-                m_maps[j](i) = (*testsource)[m_mapnames[j]].value();
+
+            	// Set parameter value and error maps
+            	if (m_mapnames[j].substr(0,2) == "e_") {
+
+            		// Get Parameter name by removing the error prefix
+            		std::string parname = m_mapnames[j].substr(2, m_mapnames[j].size());
+            		m_maps[j](i) = (*testsource)[parname].error();
+
+            	}
+            	else {
+            		m_maps[j](i) = (*testsource)[m_mapnames[j]].value();
+            	}
             }
     	}
  
-    	// Set status of bin to true
-    	m_statusmap(i) = 1.0;
+        // Set Fit status of the bin
+        m_statusmap(i) = status;
  
     	// Remove model from container
     	models.remove(m_srcname);
@@ -417,11 +443,9 @@ void cttsmap::save(void)
             fits[i+1]->extname(m_mapnames[i]);
         }
  
-        // Add computation log if not all bins are computed
-        if (m_binmin != -1 || m_binmax != -1) {
-            m_statusmap.write(fits);
-            fits[m_mapnames.size()+1]->extname("STATUS MAP");
-        }
+        // Write Fit status map
+        m_statusmap.write(fits);
+        fits[fits.size()-1]->extname("FITSTATUS");
  
         // Save FITS file
         fits.saveto(m_outmap, clobber());
@@ -482,6 +506,7 @@ void cttsmap::init_members(void)
     m_outmap.clear();
     m_apply_edisp = false;
     m_publish     = false;
+    m_errors      = false;
  
     // Initialise protected members
     m_obs.clear();
@@ -511,16 +536,17 @@ void cttsmap::copy_members(const cttsmap&amp; app)
     m_outmap      = app.m_outmap;
     m_apply_edisp = app.m_apply_edisp;
     m_publish     = app.m_publish;
+    m_errors      = app.m_errors;
  
     // Copy protected members
-    m_binmin    = app.m_binmin;
-    m_binmax    = app.m_binmax;
-    m_logL0     = app.m_logL0;
-    m_obs       = app.m_obs;
-    m_tsmap     = app.m_tsmap;
-    m_mapnames  = app.m_mapnames;
-    m_maps      = app.m_maps;
-    m_statusmap = app.m_statusmap;
+    m_binmin      = app.m_binmin;
+    m_binmax      = app.m_binmax;
+    m_logL0       = app.m_logL0;
+    m_obs         = app.m_obs;
+    m_tsmap       = app.m_tsmap;
+    m_mapnames    = app.m_mapnames;
+    m_maps        = app.m_maps;
+    m_statusmap   = app.m_statusmap;
  
     // Clone protected members
     m_testsource = (app.m_testsource != NULL) ? app.m_testsource->clone() : NULL;
@@ -612,6 +638,9 @@ void cttsmap::get_parameters(void)
     // Create sky map based on task parameters
     GSkyMap map = create_map(m_obs);
  
+    // Check if errors should be computed
+    m_errors  = (*this)["errors"].boolean();
+
     // Initialise maps from user parameters
     init_maps(map);
  
@@ -654,9 +683,13 @@ void cttsmap::init_maps(const GSkyMap&amp; map)
     // Initialise map information
 	m_statusmap.clear();
  
-	// Create status map
+	// Create status map and progress map
 	m_statusmap = GSkyMap(map);
  
+	// Set status map to invalid value -1 to signal
+	// that bin wasn't computed yet
+	m_statusmap = -1.0;
+
 	// Initialise maps of free parameters
     if (m_testsource != NULL) {
  
@@ -674,6 +707,16 @@ void cttsmap::init_maps(const GSkyMap&amp; map)
  
             // Store parameter name
             m_mapnames.push_back((*m_testsource)[i].name());
+
+            // Add parameter errors if requested
+            if (m_errors) {
+
+				// Add sky map for fit error of the free parameter
+				m_maps.push_back(map);
+
+				// Store name of parameter error
+				m_mapnames.push_back("e_" + (*m_testsource)[i].name());
+            }
  
         } // endfor: looped over all model parameters
  
@@ -89,6 +89,7 @@ protected:
     GFilename                m_outmap;      //!< Output TS map file name
     bool                     m_apply_edisp; //!< Apply energy dispersion?
     bool                     m_publish;     //!< Publish TS map?
+    bool                     m_errors;      //!< Compute and store parameter errors?
  
     // Parameters to control speed and job splitting
     int                      m_binmin;      //!< Map bin number from which computation should start
@@ -98,7 +99,7 @@ protected:
     // Protected members
     GObservations            m_obs;         //!< Observation container
     GSkyMap                  m_tsmap;       //!< TS map
-    GSkyMap                  m_statusmap;   //!< Map of computed bins
+    GSkyMap                  m_statusmap;   //!< Map of optimizer fit status
     std::vector<std::string> m_mapnames;    //!< Names of free parameters
     std::vector<GSkyMap>     m_maps;        //!< Sky maps for each free parameter
     GModel*                  m_testsource;  //!< Pointer to test source for TS computation
@@ -27,6 +27,7 @@ bkgcube,   f, a, &quot;NONE&quot;,,, &quot;Input background cube file (only needed for stacked 
 caldb,     s, a, "prod2",,, "Calibration database"
 irf,       s, a, "South_0.5h",,, "Instrument response function"
 edisp,     b, h, no,,, "Apply energy dispersion?"
+errors,    b, h, no,,, "Compute/store parameter errors?"
 outmap,    f, a, "tsmap.fits",,, "Output Test Statistic map file"
  
 #
@@ -96,6 +96,20 @@ class Test(test):
         # Check if execution failed
         self.test_assert(self._execute(cmd) != 0,
              'Check invalid input file when executed from command line')
+        
+        # Setup cttsmap command
+        cmd = cttsmap+' inobs="'+self._events+'"'+ \
+                      ' inmodel="'+self._model+'" srcname="Crab"'+ \
+                      ' caldb="'+self._caldb+'" irf="'+self._irf+'"'+ \
+                      ' errors=yes'+ \
+                      ' outmap="cttsmap_cmd1.fits"'+ \
+                      ' nxpix=5 nypix=5 binsz=0.02'+ \
+                      ' coordsys="CEL" proj="CAR" xref=83.63 yref=22.01'+ \
+                      ' logfile="cttsmap_cmd1.log" chatter=1'
+
+        # Check if execution was successful
+        self.test_assert(self._execute(cmd) == 0,
+             'Check successful execution from command line')
  
         # Return
         return