Calib.cc

Go to the documentation of this file.

// -*- lsst-c++ -*-

/* 
 * LSST Data Management System
 * Copyright 2008, 2009, 2010 LSST Corporation.
 * 
 * This product includes software developed by the
 * LSST Project (http://www.lsst.org/).
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 * 
 * You should have received a copy of the LSST License Statement and 
 * the GNU General Public License along with this program.  If not, 
 * see <http://www.lsstcorp.org/LegalNotices/>.
 */
 
#include <cmath>

#include "boost/format.hpp"
#include "boost/lexical_cast.hpp"
#include "boost/algorithm/string/trim.hpp"

#include "ndarray.h"
#include "lsst/pex/exceptions.h"
#include "lsst/daf/base/PropertySet.h"
#include "lsst/afw/image/Calib.h"
#include "lsst/afw/cameraGeom/Detector.h"
#include "lsst/afw/table/io/OutputArchive.h"
#include "lsst/afw/table/io/InputArchive.h"
#include "lsst/afw/table/io/CatalogVector.h"

namespace lsst { namespace afw { namespace image {
Calib::Calib() : _midTime(), _exptime(0.0), _fluxMag0(0.0), _fluxMag0Sigma(0.0) {}
Calib::Calib(std::vector<CONST_PTR(Calib)> const& calibs 
            ) :
    _midTime(), _exptime(0.0), _fluxMag0(0.0), _fluxMag0Sigma(0.0)
{
    if (calibs.empty()) {
        throw LSST_EXCEPT(lsst::pex::exceptions::InvalidParameterError,
                          "You must provide at least one input Calib");
    }

    double const fluxMag00 = calibs[0]->_fluxMag0;
    double const fluxMag0Sigma0 = calibs[0]->_fluxMag0Sigma;

    double midTimeSum = 0.0;            // sum(time*expTime)
    for (std::vector<CONST_PTR(Calib)>::const_iterator ptr = calibs.begin(); ptr != calibs.end(); ++ptr) {
        Calib const& calib = **ptr;

        if (::fabs(fluxMag00 - calib._fluxMag0) > std::numeric_limits<double>::epsilon() ||
            ::fabs(fluxMag0Sigma0 - calib._fluxMag0Sigma) > std::numeric_limits<double>::epsilon()) {
            throw LSST_EXCEPT(lsst::pex::exceptions::InvalidParameterError,
                              (boost::format("You may only combine calibs with the same fluxMag0: "
                                             "%g +- %g v %g +- %g")
                               % calib.getFluxMag0().first % calib.getFluxMag0().second
                               % calibs[0]->getFluxMag0().first % calibs[0]->getFluxMag0().second
                              ).str());
        }

        double const exptime = calib._exptime;

        midTimeSum += calib._midTime.get()*exptime;
        _exptime += exptime;
    }

    daf::base::DateTime tmp(midTimeSum/_exptime); // there's no way to set the double value directly
    using std::swap;
    swap(_midTime, tmp);
}

Calib::Calib(CONST_PTR(lsst::daf::base::PropertySet) metadata) {
    lsst::daf::base::DateTime midTime;
    double exptime = 0.0;
    double fluxMag0 = 0.0, fluxMag0Sigma = 0.0;

    std::string key = "TIME-MID";
    if (metadata->exists(key)) {
        midTime  = lsst::daf::base::DateTime(boost::algorithm::trim_right_copy(metadata->getAsString(key)));
    }

    key = "EXPTIME";
    if (metadata->exists(key)) {
        try {
            exptime = metadata->getAsDouble(key);
        } catch (lsst::pex::exceptions::TypeError & err) {
            std::string exptimeStr = metadata->getAsString(key);
            exptime = boost::lexical_cast<double>(exptimeStr);
        }
    }

    key = "FLUXMAG0";
    if (metadata->exists(key)) {
        fluxMag0 = metadata->getAsDouble(key);
        
        key = "FLUXMAG0ERR";
        if (metadata->exists(key)) {
            fluxMag0Sigma = metadata->getAsDouble(key);
        }
    }
        
    _midTime = midTime;
    _exptime = exptime;
    _fluxMag0 = fluxMag0;
    _fluxMag0Sigma = fluxMag0Sigma;
}
bool Calib::_throwOnNegativeFlux = true;
void
Calib::setThrowOnNegativeFlux(bool raiseException 
                             )
{
    _throwOnNegativeFlux = raiseException;
}

bool
Calib::getThrowOnNegativeFlux()
{
    return _throwOnNegativeFlux;
}

namespace detail {
int stripCalibKeywords(PTR(lsst::daf::base::PropertySet) metadata 
                      )
{
    int nstripped = 0;

    std::string key = "TIME-MID";
    if (metadata->exists(key)) {
        metadata->remove(key);
        nstripped++;
    }

    key = "EXPTIME";
    if (metadata->exists(key)) {
        metadata->remove(key);
        nstripped++;
    }

    key = "FLUXMAG0";
    if (metadata->exists(key)) {
        metadata->remove(key);
        nstripped++;
    }
        
    key = "FLUXMAG0ERR";
    if (metadata->exists(key)) {
        metadata->remove(key);
        nstripped++;
    }

    return nstripped;
}
}

bool Calib::operator==(Calib const& rhs) const {
    return
        ::fabs(_midTime.get() - rhs._midTime.get()) < std::numeric_limits<double>::epsilon() &&
        _exptime == rhs._exptime &&
        _fluxMag0 == rhs._fluxMag0 &&
        _fluxMag0Sigma == rhs._fluxMag0Sigma;
}

void Calib::setMidTime(lsst::daf::base::DateTime const& midTime 
                      )
{
    _midTime = midTime;
}

lsst::daf::base::DateTime Calib::getMidTime () const
{
    return _midTime;
}

lsst::daf::base::DateTime Calib::getMidTime(
        boost::shared_ptr<const lsst::afw::cameraGeom::Detector>, 
        lsst::afw::geom::Point2I const&            
                                           ) const
{
    return _midTime;
}

void Calib::setExptime(double exptime      
                       ) {
    _exptime = exptime;
}

double Calib::getExptime() const
{
    return _exptime;
}

void Calib::setFluxMag0(double fluxMag0,      
                        double fluxMag0Sigma  
                       )
{
    _fluxMag0 = fluxMag0;
    _fluxMag0Sigma = fluxMag0Sigma;
}
void Calib::setFluxMag0(std::pair<double, double> fluxMag0AndSigma 
                       )
{
    _fluxMag0 = fluxMag0AndSigma.first;
    _fluxMag0Sigma = fluxMag0AndSigma.second;
}

std::pair<double, double> Calib::getFluxMag0() const
{
    return std::make_pair(_fluxMag0, _fluxMag0Sigma);
}

namespace {
inline void checkNegativeFlux0(double fluxMag0) {
    if (fluxMag0 <= 0) {
        throw LSST_EXCEPT(lsst::pex::exceptions::DomainError,
                          (boost::format("Flux of 0-mag object must be >= 0: saw %g") % fluxMag0).str());
    }
}
inline bool isNegativeFlux(double flux, bool doThrow)
{
    if (flux <= 0) {
        if (doThrow) {
            throw LSST_EXCEPT(lsst::pex::exceptions::DomainError,
                              (boost::format("Flux must be >= 0: saw %g") % flux).str());
        }
        return true;
    }
    return false;
}
inline double convertToFlux(double fluxMag0, double mag) {
    return fluxMag0 * ::pow(10.0, -0.4 * mag);
}
inline double convertToFluxErr(double fluxMag0InvSNR, double flux, double magErr) {
    // Want to:
    //     return flux * hypot(_fluxMag0Sigma/_fluxMag0, 0.4*std::log(10)*magSigma/mag);
    // But hypot is not standard C++ so use <http://en.wikipedia.org/wiki/Hypot#Implementation>
    double a = fluxMag0InvSNR;
    double b = 0.4 * std::log(10.0) * magErr;
    if (std::abs(a) < std::abs(b)) {
        std::swap(a, b);
    }
    return flux * std::abs(a) * std::sqrt(1 + std::pow(b / a, 2));
}
inline double convertToMag(double fluxMag0, double flux) {
    return -2.5*::log10(flux/fluxMag0);
}
inline void convertToMagWithErr(double *mag, double *magErr, double fluxMag0, double fluxMag0InvSNR,
                                double flux, double fluxErr)
{
    double const rat = flux/fluxMag0;
    double const ratErr = ::sqrt((::pow(fluxErr, 2) + ::pow(flux*fluxMag0InvSNR, 2))/::pow(fluxMag0, 2));
    
    *mag = -2.5*::log10(rat);
    *magErr = 2.5/::log(10.0)*ratErr/rat;
}

} // anonymous namespace

double Calib::getFlux(double const mag 
                        ) const {
    checkNegativeFlux0(_fluxMag0);
    return convertToFlux(_fluxMag0, mag);
}
ndarray::Array<double,1> Calib::getFlux(ndarray::Array<double const,1> const & mag) const {
    checkNegativeFlux0(_fluxMag0);
    ndarray::Array<double,1> flux = ndarray::allocate(mag.size());
    ndarray::Array<double const,1>::Iterator inIter = mag.begin();
    ndarray::Array<double,1>::Iterator outIter = flux.begin();
    for (; inIter != mag.end(); ++inIter, ++outIter) {
        *outIter = convertToFlux(_fluxMag0, *inIter);
    }
    return flux;
}

std::pair<double, double> Calib::getFlux(
        double const mag,       
        double const magSigma   
    ) const
{
    checkNegativeFlux0(_fluxMag0);
    double const flux = convertToFlux(_fluxMag0, mag);
    double const fluxErr = convertToFluxErr(_fluxMag0Sigma/_fluxMag0, flux, magSigma);
    return std::make_pair(flux, fluxErr);
}

std::pair<ndarray::Array<double,1>, ndarray::Array<double,1> > Calib::getFlux(
    ndarray::Array<double const,1> const & mag,
    ndarray::Array<double const,1> const & magErr
) const {
    checkNegativeFlux0(_fluxMag0);
    if (mag.size() != magErr.size()) {
        throw LSST_EXCEPT(lsst::pex::exceptions::LengthError,
                          (boost::format("Size of mag (%d) and magErr (%d) don't match") % 
                           mag.size() % magErr.size()).str());
    }

    ndarray::Array<double,1> flux = ndarray::allocate(mag.size());
    ndarray::Array<double,1> fluxErr = ndarray::allocate(mag.size());
    ndarray::Array<double const,1>::Iterator magIter = mag.begin();
    ndarray::Array<double const,1>::Iterator magErrIter = magErr.begin();
    ndarray::Array<double,1>::Iterator fluxIter = flux.begin();
    ndarray::Array<double,1>::Iterator fluxErrIter = fluxErr.begin();

    double fluxMag0InvSNR = _fluxMag0Sigma/_fluxMag0;
    for (; magIter != mag.end(); ++magIter, ++magErrIter, ++fluxIter, ++fluxErrIter) {
        *fluxIter = convertToFlux(_fluxMag0, *magIter);
        *fluxErrIter = convertToFluxErr(fluxMag0InvSNR, *fluxIter, *magErrIter);
    }

    return std::make_pair(flux, fluxErr);
}

double Calib::getMagnitude(double const flux 
                         ) const
{
    checkNegativeFlux0(_fluxMag0);
    if (isNegativeFlux(flux, Calib::getThrowOnNegativeFlux())) {
        return std::numeric_limits<double>::quiet_NaN();
    }
    return convertToMag(_fluxMag0, flux);
}

std::pair<double, double> Calib::getMagnitude(double const flux, 
                                              double const fluxErr 
                                              ) const
{
    checkNegativeFlux0(_fluxMag0);
    if (isNegativeFlux(flux, Calib::getThrowOnNegativeFlux())) {
        double const NaN = std::numeric_limits<double>::quiet_NaN();
        return std::make_pair(NaN, NaN);
    }

    double mag, magErr;
    convertToMagWithErr(&mag, &magErr, _fluxMag0, _fluxMag0Sigma/_fluxMag0, flux, fluxErr);
    return std::make_pair(mag, magErr);
}

ndarray::Array<double,1> Calib::getMagnitude(ndarray::Array<double const,1> const & flux) const {
    checkNegativeFlux0(_fluxMag0);
    ndarray::Array<double,1> mag = ndarray::allocate(flux.size());
    ndarray::Array<double const,1>::Iterator fluxIter = flux.begin();
    ndarray::Array<double,1>::Iterator magIter = mag.begin();
    int nonPositive = 0;
    for (; fluxIter != flux.end(); ++fluxIter, ++magIter) {
        if (isNegativeFlux(*fluxIter, false)) {
            ++nonPositive;
            *magIter = std::numeric_limits<double>::quiet_NaN();
            continue;
        }
        *magIter = convertToMag(_fluxMag0, *fluxIter);
    }
    if (nonPositive && Calib::getThrowOnNegativeFlux()) {
        throw LSST_EXCEPT(lsst::pex::exceptions::DomainError,
                          (boost::format("Flux must be >= 0: %d non-positive seen") % nonPositive).str());
    }
    return mag;
}

std::pair<ndarray::Array<double,1>, ndarray::Array<double,1> > Calib::getMagnitude(
    ndarray::Array<double const,1> const & flux,
    ndarray::Array<double const,1> const & fluxErr
) const {
    checkNegativeFlux0(_fluxMag0);
    if (flux.size() != fluxErr.size()) {
        throw LSST_EXCEPT(lsst::pex::exceptions::LengthError,
                          (boost::format("Size of flux (%d) and fluxErr (%d) don't match") % 
                           flux.size() % fluxErr.size()).str());
    }

    ndarray::Array<double,1> mag = ndarray::allocate(flux.size());
    ndarray::Array<double,1> magErr = ndarray::allocate(flux.size());
    ndarray::Array<double const,1>::Iterator fluxIter = flux.begin();
    ndarray::Array<double const,1>::Iterator fluxErrIter = fluxErr.begin();
    ndarray::Array<double,1>::Iterator magIter = mag.begin();
    ndarray::Array<double,1>::Iterator magErrIter = magErr.begin();
    int nonPositive = 0;
    double fluxMag0InvSNR = _fluxMag0Sigma/_fluxMag0;
    for (; fluxIter != flux.end(); ++fluxIter, ++fluxErrIter, ++magIter, ++magErrIter) {
        if (isNegativeFlux(*fluxIter, false)) {
            ++nonPositive;
            double const NaN = std::numeric_limits<double>::quiet_NaN();
            *magIter = NaN;
            *magErrIter = NaN;
            continue;
        }
        double f, df;
        convertToMagWithErr(&f, &df, _fluxMag0, fluxMag0InvSNR, *fluxIter, *fluxErrIter);
        *magIter = f;
        *magErrIter = df;
    }
    if (nonPositive && Calib::getThrowOnNegativeFlux()) {
        throw LSST_EXCEPT(lsst::pex::exceptions::DomainError,
                          (boost::format("Flux must be >= 0: %d non-positive seen") % nonPositive).str());
    }
    return std::make_pair(mag, magErr);
}
   
namespace {

class CalibSchema : private boost::noncopyable {
public:
    table::Schema schema;
    table::Key<boost::int64_t> midTime;
    table::Key<double> expTime;
    table::Key<double> fluxMag0;
    table::Key<double> fluxMag0Sigma;
    
    static CalibSchema const & get() {
        static CalibSchema instance;
        return instance;
    }

private:
    CalibSchema() :
        schema(),
        midTime(schema.addField<boost::int64_t>(
                    "midtime", "middle of the time of the exposure relative to Unix epoch", "nanoseconds"
                )),
        expTime(schema.addField<double>("exptime", "exposure time", "seconds")),
        fluxMag0(schema.addField<double>("fluxmag0", "flux of a zero-magnitude object", "dn")),
        fluxMag0Sigma(schema.addField<double>("fluxmag0.err", "1-sigma error on fluxmag0", "dn"))
    {
        schema.getCitizen().markPersistent();
    }
};

class CalibFactory : public table::io::PersistableFactory {
public:

    virtual PTR(table::io::Persistable)
    read(InputArchive const & archive, CatalogVector const & catalogs) const {
        CalibSchema const & keys = CalibSchema::get();
        LSST_ARCHIVE_ASSERT(catalogs.size() == 1u);
        LSST_ARCHIVE_ASSERT(catalogs.front().size() == 1u);
        LSST_ARCHIVE_ASSERT(catalogs.front().getSchema() == keys.schema);
        table::BaseRecord const & record = catalogs.front().front();
        PTR(Calib) result(new Calib());
        result->setMidTime(daf::base::DateTime(static_cast<long long>(record.get(keys.midTime))));
        result->setExptime(record.get(keys.expTime));
        result->setFluxMag0(record.get(keys.fluxMag0), record.get(keys.fluxMag0Sigma));
        return result;
    }

    explicit CalibFactory(std::string const & name) : table::io::PersistableFactory(name) {}

};

std::string getCalibPersistenceName() { return "Calib"; }

CalibFactory registration(getCalibPersistenceName());

} // anonymous

std::string Calib::getPersistenceName() const { return getCalibPersistenceName(); }

void Calib::write(OutputArchiveHandle & handle) const {
    CalibSchema const & keys = CalibSchema::get();
    table::BaseCatalog cat = handle.makeCatalog(keys.schema);
    PTR(table::BaseRecord) record = cat.addNew();
    record->set(keys.midTime, getMidTime().nsecs());
    record->set(keys.expTime, getExptime());
    std::pair<double,double> fluxMag0 = getFluxMag0();
    record->set(keys.fluxMag0, fluxMag0.first);
    record->set(keys.fluxMag0Sigma, fluxMag0.second);
    handle.saveCatalog(cat);
}

}}}  // lsst::afw::image