Autodetect file format from filename (#1863)

Jürgen Knödlseder
1 parent e724faa8
Showing 11 changed files with 385 additions and 139 deletions
.gitignore
ChangeLog
NEWS
include/GFilename.hpp
include/GObservations.hpp
pyext/GFilename.i
src/obs/GObservations_likelihood.cpp
src/support/GFilename.cpp
test/test_GObservation.cpp
test/test_GObservation.hpp
test/test_GSupport.cpp
@@ -113,6 +113,7 @@ test/*.fits
 test/*.xml
 test/*.trs
 test/*.html
+test/*.csv
 test/test_csv*.dat
 test/test_url.dat
 test/reports/*.xml
-2017-04-05
+2017-04-10
         * Version 1.3.0 released
         =========================
+        Add GFilename::type() method (#1863)
         Copy shared Python libraries instead of creating a symbolic link (#1987)
         Added FITS image conversion operators (#1986)
         If code originates from Git, always compile the Python bindings (#1985)
 New Features and Important Changes in GammaLib 1.3.0.dev1
-5 April 2017
+10 April 2017
 1. Introduction
@@ -28,6 +28,7 @@ The following methods have been added:
 - GCTAEdisp2D::fetch
 - GEbounds::set
 - GEnergies::set
+- GFilename::type()
 The following methods have been removed:
 - ...
@@ -117,7 +118,8 @@ The GWcs base class now also supports CD keywords instead of CDELT keywords
 13. Support module
 ------------------
-None
+A GFilename::type() method has been added that returns the type of the file
+based on the file type extension (e.g. ".fits") (#1863).
 14. XML module
 /***************************************************************************
  *                      GFilename.hpp - Filename class                     *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2015-2016 by Juergen Knoedlseder                         *
+ *  copyright (C) 2015-2017 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -62,8 +62,10 @@
 class GFilename : public GBase {
     // Friend functions
-    friend std::string operator+(const GFilename& filename, const std::string& string);
-    friend std::string operator+(const std::string& string, const GFilename& filename);
+    friend std::string operator+(const GFilename&   filename,
+                                 const std::string& string);
+    friend std::string operator+(const std::string& string,
+                                 const GFilename&   filename);
     friend bool        operator==(const GFilename &a, const GFilename &b);
     friend bool        operator!=(const GFilename &a, const GFilename &b);
@@ -89,6 +91,7 @@ public:
     std::string        protocol(void) const;
     std::string        path(void) const;
     std::string        file(void) const;
+    std::string        type(void) const;
     bool               exists(void) const;
     bool               is_fits(void) const;
     void               remove(void) const;
 /***************************************************************************
  *              GObservations.hpp - Observation container class            *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2009-2016 by Juergen Knoedlseder                         *
+ *  copyright (C) 2009-2017 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -158,13 +158,16 @@ public:
         double        npred(void) const;
         GMatrixSparse hessian(const GOptimizerPars& pars);
         GMatrixSparse covariance(void) const;
-        void          save(const GFilename& filename, bool csv = true) const;
+        void          save(const GFilename& filename) const;
     protected:
         // Protected methods
-        void           init_members(void);
-        void           copy_members(const likelihood& fct);
-        void           free_members(void);
+        void                     init_members(void);
+        void                     copy_members(const likelihood& fct);
+        void                     free_members(void);
+        void                     save_csv(const GFilename& filename) const;
+        void                     save_fits(const GFilename& filename) const;
+        std::vector<std::string> covariance_names(void) const;
         // Protected data members
         double         m_value;     //!< Function value
 /***************************************************************************
  *                       GFilename.i - Filename class                      *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2015-2016 by Juergen Knoedlseder                         *
+ *  copyright (C) 2015-2017 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -54,6 +54,7 @@ public:
     std::string        protocol(void) const;
     std::string        path(void) const;
     std::string        file(void) const;
+    std::string        type(void) const;
     bool               exists(void) const;
     bool               is_fits(void) const;
     void               remove(void) const;
 /***************************************************************************
  *         GObservations_likelihood.cpp - Likelihood function class        *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2009-2016 by Juergen Knoedlseder                         *
+ *  copyright (C) 2009-2017 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -621,88 +621,27 @@ GMatrixSparse GObservations::likelihood::covariance(void) const
 /***********************************************************************//**
- * @brief Save likelihood fit results into CSV or FITS file.
+ * @brief Save likelihood fit results into a CSV or FITS file.
  *
- * @param[in] filename CSV/FITS filename.
- * @param[in] csv save as CSV? (default: true)
+ * @param[in] filename CSV or FITS filename.
  *
- * Saves the likelihood fit results into a CSV/FITS file. For the moment the
- * method only writes the covariance matrix.
+ * Saves the likelihood fit results into a CSV or FITS file. The result
+ * format depends on the filename extension. If the extension is `.fits` or
+ * `.fit` the file is written into a FITS file, otherwise it is written into
+ * a CSV file.
  ***************************************************************************/
-void GObservations::likelihood::save(const GFilename& filename, bool csv) const
+void GObservations::likelihood::save(const GFilename& filename) const
 {
-    // Get covariance matrix
-    GMatrixSparse covmat = covariance();
+    // Get file type
+    std::string filetype = filename.type();
-    // Create covariance matrix entry names
-    std::vector<std::string> covmat_entries;
-    for (int i = 0 ; i < m_this->m_models.size(); ++i) {
-	    for (int j = 0; j < m_this->m_models[i]->size(); ++j) {
-            covmat_entries.push_back(m_this->m_models[i]->at(j).name() + "(" +
-                                     m_this->m_models[i]->name() + ")");
-        }
+    // If file type is a FITS file then write save covariance matrix into a
+    // FITS file, otherwise save the covariance matrix into a CSV file
+    if (filetype == "fits") {
+        save_fits(filename);
     }
-
-    // Create binary table and columns
-    int size = covmat_entries.size();
-
-    if (csv) {
-        GCsv table(size + 1, size);
-
-        // Fill CSV header
-        for (int i = 0; i < size; ++i) {
-            table.string(0, i, covmat_entries[i]);
-        }
-
-        // Fill CSV table
-        for (int i = 0; i < size; ++i) {
-            for (int j = 0; j < size; j ++) {
-                table.real(i+1, j, covmat(i,j));
-            }
-        }
-
-        // Save CSV table
-        table.save(filename, " ", true);
-    } else {
-
-        GFitsBinTable       covmat_table;
-        GFitsTableStringCol par("Parameters", 1, 50, size);
-        GFitsTableDoubleCol cov("Covariance", 1, size*size);
-
-        // Fill tables
-        int counter = 0;
-        for (int i = 0; i < size; ++i) {
-            par(0, i) = covmat_entries[i];
-            for (int j = 0; j < size; ++j) {
-                cov(0, counter) = covmat(i,j);
-                ++counter;
-            }
-        }
-
-        // Set dimension for covariance matrix column
-        std::vector<int> dim;
-        dim.push_back(size);
-        dim.push_back(size);
-        cov.dim(dim);
-
-        // Append columns to table
-        covmat_table.append(par);
-        covmat_table.append(cov);
-
-        // Set extension name
-        covmat_table.extname("Covariance Matrix");
-
-        // Allocate FITS object
-        GFits fits;
-
-        // Append covariance matrix table to FITS object
-        fits.append(covmat_table);
-
-        // Save FITS file to disk
-        fits.saveto(filename, true);
-
-        // Close FITS object
-        fits.close();
+    else {
+        save_csv(filename);
     }
     // Return
@@ -776,3 +715,134 @@ void GObservations::likelihood::free_members(void)
     // Return
     return;
 }
+
+
+/***********************************************************************//**
+ * @brief Save likelihood fit results into CSV file.
+ *
+ * @param[in] filename CSV filename.
+ *
+ * Saves the likelihood fit results into a CSV file. For the moment the
+ * method only writes the covariance matrix.
+ ***************************************************************************/
+void GObservations::likelihood::save_csv(const GFilename& filename) const
+{
+    // Get covariance matrix
+    GMatrixSparse covmat = covariance();
+
+    // Get covariance matrix row and column names
+    std::vector<std::string> names = covariance_names();
+
+    // Create binary table and columns
+    int size = names.size();
+
+    // Create binary table and columns
+    GCsv table(size + 1, size);
+
+    // Fill CSV header
+    for (int i = 0; i < size; ++i) {
+        table.string(0, i, names[i]);
+    }
+
+    // Fill CSV table
+    for (int i = 0; i < size; ++i) {
+        for (int j = 0; j < size; j ++) {
+            table.real(i+1, j, covmat(i,j));
+        }
+    }
+
+    // Save CSV table
+    table.save(filename, " ", true);
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Save likelihood fit results into FITS file.
+ *
+ * @param[in] filename FITS filename.
+ *
+ * Saves the likelihood fit results into a FITS file. For the moment the
+ * method only writes the covariance matrix.
+ ***************************************************************************/
+void GObservations::likelihood::save_fits(const GFilename& filename) const
+{
+    // Get covariance matrix
+    GMatrixSparse covmat = covariance();
+
+    // Get covariance matrix row and column names
+    std::vector<std::string> names = covariance_names();
+
+    // Create binary table and columns
+    int size = names.size();
+
+    // Create binary table and columns
+    GFitsBinTable       covmat_table;
+    GFitsTableStringCol par("Parameters", 1, 50, size);
+    GFitsTableDoubleCol cov("Covariance", 1, size*size);
+
+    // Fill tables
+    int counter = 0;
+    for (int i = 0; i < size; ++i) {
+        par(0, i) = names[i];
+        for (int j = 0; j < size; ++j) {
+            cov(0, counter) = covmat(i,j);
+            ++counter;
+        }
+    }
+
+    // Set dimension for covariance matrix column
+    std::vector<int> dim;
+    dim.push_back(size);
+    dim.push_back(size);
+    cov.dim(dim);
+
+    // Append columns to table
+    covmat_table.append(par);
+    covmat_table.append(cov);
+
+    // Set extension name
+    covmat_table.extname("Covariance Matrix");
+
+    // Allocate FITS object
+    GFits fits;
+
+    // Append covariance matrix table to FITS object
+    fits.append(covmat_table);
+
+    // Save FITS file to disk
+    fits.saveto(filename, true);
+
+    // Close FITS object
+    fits.close();
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Return covariance matrix row and column names
+ *
+ * @return Covariance matrix row and column names.
+ *
+ * Returns the row and column names of the covariance matrix.
+ ***************************************************************************/
+std::vector<std::string> GObservations::likelihood::covariance_names(void) const
+{
+    // Initialise covariance matrix row and column names
+    std::vector<std::string> names;
+
+    // Create covariance matrix row and column names
+    for (int i = 0 ; i < m_this->m_models.size(); ++i) {
+	    for (int j = 0; j < m_this->m_models[i]->size(); ++j) {
+            names.push_back(m_this->m_models[i]->at(j).name() + "(" +
+                            m_this->m_models[i]->name() + ")");
+        }
+    }
+
+    // Return covariance matrix row and column names
+    return names;
+}
 /***************************************************************************
  *                      GFilename.cpp - Filename class                     *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2015-2016 by Juergen Knoedlseder                         *
+ *  copyright (C) 2015-2017 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -248,6 +248,61 @@ bool GFilename::exists(void) const
 /***********************************************************************//**
+ * @brief Return file type
+ *
+ * @return File type.
+ *
+ * Returns the file type based on the file extension name. Possible file
+ * types are
+ *
+ *        Type   Extension names
+ *        fits   .fits / .fit
+ *
+ * For any other file types an empty string will be returned. Any compression
+ * extensions such as `.gz`, `.Z`, `.z` or `.zip` will be stripped.
+ ***************************************************************************/
+std::string GFilename::type(void) const
+{
+    // Initialise empty file type
+    std::string type = "";
+
+    // Initialise type with the name of the file
+    std::string file = m_file;
+
+    // Strip compression extensions
+    size_t start = file.rfind(".gz");
+    if (start != std::string::npos) {
+        file = file.substr(0, start);
+    }
+    start = file.rfind(".Z");
+    if (start != std::string::npos) {
+        file = file.substr(0, start);
+    }
+    start = file.rfind(".z");
+    if (start != std::string::npos) {
+        file = file.substr(0, start);
+    }
+    start = file.rfind(".zip");
+    if (start != std::string::npos) {
+        file = file.substr(0, start);
+    }
+
+    // Extract file type
+    start = file.rfind(".fits");
+    if (start != std::string::npos) {
+        type = "fits";
+    }
+    start = file.rfind(".fit");
+    if (start != std::string::npos) {
+        type = "fits";
+    }
+
+    // Return file type
+    return type;
+}
+
+
+/***********************************************************************//**
  * @brief Checks whether file is a FITS file
  *
  * @return True if file is a FITS file.
@@ -55,14 +55,24 @@ void TestGObservation::set(void)
     name("Observation module");
     // Append tests
-    append(static_cast<pfunction>(&TestGObservation::test_time_reference), "Test GTimeReference class");
-    append(static_cast<pfunction>(&TestGObservation::test_time), "Test GTime class");
-    append(static_cast<pfunction>(&TestGObservation::test_times), "Test GTimes class");
-    append(static_cast<pfunction>(&TestGObservation::test_gti), "Test GGti class");
-    append(static_cast<pfunction>(&TestGObservation::test_energy), "Test GEnergy class");
-    append(static_cast<pfunction>(&TestGObservation::test_energies), "Test GEnergies class");
-    append(static_cast<pfunction>(&TestGObservation::test_ebounds), "Test GEbounds class");
-    append(static_cast<pfunction>(&TestGObservation::test_photons), "Test GPhotons class");
+    append(static_cast<pfunction>(&TestGObservation::test_time_reference),
+           "Test GTimeReference class");
+    append(static_cast<pfunction>(&TestGObservation::test_time),
+           "Test GTime class");
+    append(static_cast<pfunction>(&TestGObservation::test_times),
+           "Test GTimes class");
+    append(static_cast<pfunction>(&TestGObservation::test_gti),
+           "Test GGti class");
+    append(static_cast<pfunction>(&TestGObservation::test_energy),
+           "Test GEnergy class");
+    append(static_cast<pfunction>(&TestGObservation::test_energies),
+           "Test GEnergies class");
+    append(static_cast<pfunction>(&TestGObservation::test_ebounds),
+           "Test GEbounds class");
+    append(static_cast<pfunction>(&TestGObservation::test_photons),
+           "Test GPhotons class");
+    append(static_cast<pfunction>(&TestGObservation::test_observations_optimizer),
+           "Test GObservations::likelihood class");
     // Return
     return;
@@ -1262,6 +1272,112 @@ void TestGObservation::test_photons(void)
 }
+/***********************************************************************//**
+ * @brief Test GObservations::likelihood
+ ***************************************************************************/
+void TestGObservation::test_observations_optimizer(void)
+{
+    // Create Test Model
+    GTestModelData model;
+
+    // Create Models conteners
+    GModels models;
+    models.append(model);
+
+    // Time iterval
+    GTime tmin(0.0);
+    GTime tmax(1800.0);
+
+    // Rate : events/sec
+    double rate = RATE;
+
+    // Create observations
+    GObservations obs;
+
+    // Add observations to container
+    for (int i = 0; i < 6; ++i) {
+
+        // Set random Generator
+        GRan ran(i);
+
+        // Create an event list
+        GEvents *events = model.generateList(rate, tmin, tmax, ran);
+
+        // Create a observation
+        GTestObservation observation;
+        observation.id(gammalib::str(i));
+
+        // Add events to the observation
+        observation.events(*events);
+        observation.ontime(tmax.secs()-tmin.secs());
+
+        // Append observation to container
+        obs.append(observation);
+
+        // Free the event list
+        delete events;
+        
+    } // endfor: added observations to container
+
+    // Add the model to the observation
+    obs.models(models);
+
+    // Create a GLog for show the interations of optimizer.
+    GLog log;
+
+    // Create an optimizer.
+    GOptimizerLM opt(&log);
+
+    // Set number of stalls
+    opt.max_stalls(50);
+
+    // Optimize
+    obs.optimize(opt);
+
+    // Compute errors
+    obs.errors(opt);
+
+    // Get the result
+    GModelPar result = (*(obs.models()[0]))[0];
+
+    // Check if converged
+    test_value(opt.status(), 0, "Check if converged",
+                                "Optimizer did not converge");
+
+    // Check if value is correct
+    test_value(result.factor_value(), RATE, result.factor_error()*3.0);
+
+    // Save covariance matrix as FITS file
+    obs.function().save("test_observations_optimizer.fits");
+
+    // Try loading covariance matrix as FITS file
+    test_try("Load covariance matrix as FITS file");
+    try {
+        GFits file("test_observations_optimizer.fits");
+        test_try_success();
+    }
+    catch (std::exception &e) {
+        test_try_failure(e);
+    }
+
+    // Save covariance matrix as CSV file
+    obs.function().save("test_observations_optimizer.csv");
+
+    // Try loading covariance matrix as CSV file
+    test_try("Load covariance matrix as CSV file");
+    try {
+        GCsv file("test_observations_optimizer.csv");
+        test_try_success();
+    }
+    catch (std::exception &e) {
+        test_try_failure(e);
+    }
+
+    // Return
+    return;
+}
+
+
 #ifdef _OPENMP
 /***********************************************************************//**
 * @brief Set tests
@@ -1363,6 +1479,7 @@ void TestOpenMP::test_observations_optimizer(const int&amp; mode)
     // Create an optimizer.
     GOptimizerLM opt(&log);
+    // Set number of stalls
     opt.max_stalls(50);
     // Optimize
 /***************************************************************************
  *              test_GObservation.hpp - Test observation module            *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2012-2014 by Jean-Baptiste Cayrou                        *
+ *  copyright (C) 2012-2017 by Jean-Baptiste Cayrou                        *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -55,6 +55,7 @@ public:
     void                      test_times(void);
     void                      test_energy(void);
     void                      test_energies(void);
+    void                      test_observations_optimizer(void);
 };
 /***************************************************************************
  *                 test_GSupport.cpp - test support module                 *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2010-2016 by Juergen Knoedlseder                         *
+ *  copyright (C) 2010-2017 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -710,31 +710,21 @@ void TestGSupport::test_filename(void)
     // Set filename without extension
     GFilename filename = "myfile.fits";
-    test_assert(filename.url() == "myfile.fits",
-                "Expected \"myfile.fits\" filename, "
-                "found \""+filename.url()+"\"");
+    test_value(filename.url(), "myfile.fits");
     test_assert(!filename.has_extname(),
                 "Expected that extension name is not set.");
     // Set filename with extension name
     filename = "myfile.fits[EVENTS]";
-    test_assert(filename.url() == "myfile.fits",
-                "Expected \"myfile.fits\" filename, "
-                "found \""+filename.url()+"\"");
-    test_assert(filename.extname() == "EVENTS",
-                "Expected \"EVENTS\" extension name, "
-                "found \""+filename.extname()+"\"");
+    test_value(filename.url(), "myfile.fits");
+    test_value(filename.extname(), "EVENTS");
     test_assert(filename.has_extname(),
                 "Expected that extension name is set.");
     // Set filename with extension name and version
     filename = "myfile.fits[EVENTS,2]";
-    test_assert(filename.url() == "myfile.fits",
-                "Expected \"myfile.fits\" filename, "
-                "found \""+filename.url()+"\"");
-    test_assert(filename.extname() == "EVENTS",
-                "Expected \"EVENTS\" extension name, "
-                "found \""+filename.extname()+"\"");
+    test_value(filename.url(), "myfile.fits");
+    test_value(filename.extname(), "EVENTS");
     test_assert(filename.has_extname(),
                 "Expected that extension name is set.");
     test_value(filename.extver(), 2);
@@ -743,21 +733,15 @@ void TestGSupport::test_filename(void)
     // Set filename with lower case extension name
     filename = "myfile.fits[events]";
-    test_assert(filename.extname() == "EVENTS",
-                "Expected \"EVENTS\" extension name, "
-                "found \""+filename.extname()+"\"");
+    test_value(filename.extname(), "EVENTS");
     // Set filename with number in extension name
     filename = "myfile.fits[EVENTS_2D]";
-    test_assert(filename.extname() == "EVENTS_2D",
-                "Expected \"EVENTS_2D\" extension name, "
-                "found \""+filename.extname()+"\"");
+    test_value(filename.extname(), "EVENTS_2D");
     // Set filename with extension number and version
     filename = "myfile.fits[1,2]";
-    test_assert(filename.url() == "myfile.fits",
-                "Expected \"myfile.fits\" filename, "
-                "found \""+filename.url()+"\"");
+    test_value(filename.url(), "myfile.fits");
     test_value(filename.extno(), 1);
     test_assert(filename.has_extno(),
                 "Expected that extension number is set.");
@@ -767,15 +751,9 @@ void TestGSupport::test_filename(void)
     // Set filename with extension name and expression
     filename = "myfile.fits[EVENTS][ENERGY>0.1]";
-    test_assert(filename.url() == "myfile.fits",
-                "Expected \"myfile.fits\" filename, "
-                "found \""+filename.url()+"\"");
-    test_assert(filename.extname() == "EVENTS",
-                 "Expected \"EVENTS\" extension name, "
-                 "found \""+filename.extname()+"\"");
-    test_assert(filename.expression() == "ENERGY>0.1",
-                 "Expected \"EVENTS>0.1\" expression, "
-                  "found \""+filename.expression()+"\"");
+    test_value(filename.url(), "myfile.fits");
+    test_value(filename.extname(), "EVENTS");
+    test_value(filename.expression(), "ENERGY>0.1");
     // Test missing closing symbol
     test_try("Missing ] symbol");
@@ -854,13 +832,9 @@ void TestGSupport::test_filename(void)
     // Test default extension name
     filename = "myfile.fits";
-    test_assert(filename.extname("EVENTS") == "EVENTS",
-                "Expected \"EVENTS\" extension name, "
-                "found \""+filename.extname()+"\"");
+    test_value(filename.extname("EVENTS"), "EVENTS");
     filename = "myfile.fits[LUCKY]";
-    test_assert(filename.extname("EVENTS") == "LUCKY",
-                "Expected \"LUCKY\" extension name, "
-                "found \""+filename.extname()+"\"");
+    test_value(filename.extname("EVENTS"), "LUCKY");
     // Test default extension number
     filename = "myfile.fits";
@@ -874,6 +848,24 @@ void TestGSupport::test_filename(void)
     filename = "myfile.fits[3,3]";
     test_value(filename.extver(1), 3);
+    // Test type method
+    filename = "myfile.fits";
+    test_value(filename.type(), "fits");
+    filename = "myfile.fit";
+    test_value(filename.type(), "fits");
+    filename = "myfile.fits.gz";
+    test_value(filename.type(), "fits");
+    filename = "myfile.fits.Z";
+    test_value(filename.type(), "fits");
+    filename = "myfile.fits.z";
+    test_value(filename.type(), "fits");
+    filename = "myfile.fits.zip";
+    test_value(filename.type(), "fits");
+    filename = "myfile.csv";
+    test_value(filename.type(), "");
+    filename = "myfile.csv.gz";
+    test_value(filename.type(), "");
+
     // Return
     return;
 }