PhotometryMapping.h

Go to the documentation of this file.

// -*- LSST-C++ -*-
/*
 * This file is part of jointcal.
 *
 * Developed for the LSST Data Management System.
 * This product includes software developed by the LSST Project
 * (https://www.lsst.org).
 * See the COPYRIGHT file at the top-level directory of this distribution
 * for details of code ownership.
 *
 * 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 GNU General Public License
 * along with this program.  If not, see <https://www.gnu.org/licenses/>.
 */
 
#ifndef LSST_JOINTCAL_PHOTOMETRY_MAPPING_H
#define LSST_JOINTCAL_PHOTOMETRY_MAPPING_H
 
#include <memory>
 
#include "lsst/afw/image/PhotoCalib.h"
 
#include "lsst/jointcal/CcdImage.h"
#include "lsst/jointcal/Eigenstuff.h"
#include "lsst/jointcal/MeasuredStar.h"
#include "lsst/jointcal/Point.h"
#include "lsst/jointcal/PhotometryTransform.h"
 
namespace lsst {
namespace jointcal {
 
class PhotometryMappingBase {
public:
    PhotometryMappingBase() : index(-1), fixed(false) {}
    virtual ~PhotometryMappingBase(){};
 
    PhotometryMappingBase(PhotometryMappingBase const &) = delete;
    PhotometryMappingBase(PhotometryMappingBase &&) = delete;
    PhotometryMappingBase &operator=(PhotometryMappingBase const &) = delete;
    PhotometryMappingBase &operator=(PhotometryMappingBase &&) = delete;
 
    virtual std::size_t getNpar() const = 0;
 
    virtual double transform(MeasuredStar const &measuredStar, double value) const = 0;
 
    virtual double transformError(MeasuredStar const &measuredStar, double value, double valueErr) const = 0;
 
    virtual void freezeErrorTransform() = 0;
 
    virtual void computeParameterDerivatives(MeasuredStar const &measuredStar, double value,
                                             Eigen::Ref<Eigen::VectorXd> derivatives) const = 0;
 
    void setFixed(bool _fixed) { fixed = _fixed; }
    bool isFixed() { return fixed; }
 
    virtual Eigen::VectorXd getParameters() = 0;
 
    virtual void getMappingIndices(IndexVector &indices) const = 0;
 
    virtual void print(std::ostream &out) const = 0;
 
    Eigen::Index getIndex() { return index; }
 
    void setIndex(Eigen::Index i) { index = i; }
 
protected:
    // Start index of this mapping in the "grand" fit
    Eigen::Index index;
    // Should this mapping be varied during fitting?
    bool fixed;
};
 
class PhotometryMapping : public PhotometryMappingBase {
public:
    explicit PhotometryMapping(std::shared_ptr<PhotometryTransform> transform)
            : PhotometryMappingBase(), _transform(std::move(transform)), _transformErrors(_transform) {}
 
    std::size_t getNpar() const override {
        if (fixed) {
            return 0;
        } else {
            return _transform->getNpar();
        }
    }
 
    double transform(MeasuredStar const &measuredStar, double value) const override {
        return _transform->transform(measuredStar.x, measuredStar.y, value);
    }
 
    double transformError(MeasuredStar const &measuredStar, double value, double valueErr) const override {
        return _transformErrors->transformError(measuredStar.x, measuredStar.y, value, valueErr);
    }
 
    void freezeErrorTransform() override {
        _transformErrors = std::shared_ptr<PhotometryTransform>(_transform->clone());
    }
 
    void computeParameterDerivatives(MeasuredStar const &measuredStar, double value,
                                     Eigen::Ref<Eigen::VectorXd> derivatives) const override {
        if (fixed) {
            return;
        } else {
            _transform->computeParameterDerivatives(measuredStar.x, measuredStar.y, value, derivatives);
        }
    }
 
    void offsetParams(Eigen::VectorXd const &delta) { _transform->offsetParams(delta); }
 
    Eigen::VectorXd getParameters() override { return _transform->getParameters(); }
 
    void getMappingIndices(IndexVector &indices) const override {
        if (indices.size() < getNpar()) indices.resize(getNpar());
        for (std::size_t k = 0; k < getNpar(); ++k) {
            indices[k] = index + k;
        }
    }
 
    void print(std::ostream &out) const override {
        out << "index: " << index << " fixed: " << fixed << " ";
        _transform->print(out);
    }
 
    std::shared_ptr<PhotometryTransform> getTransform() const { return _transform; }
 
    std::shared_ptr<PhotometryTransform> getTransformErrors() const { return _transformErrors; }
 
private:
    // the actual transformation to be fit
    std::shared_ptr<PhotometryTransform> _transform;
    // the transformation used for errors
    std::shared_ptr<PhotometryTransform> _transformErrors;
};
 
class ChipVisitPhotometryMapping : public PhotometryMappingBase {
public:
    ChipVisitPhotometryMapping(std::shared_ptr<PhotometryMapping> chipMapping,
                               std::shared_ptr<PhotometryMapping> visitMapping)
            : PhotometryMappingBase(),
              _nParChip(0),
              _nParVisit(0),
              _chipMapping(std::move(chipMapping)),
              _visitMapping(std::move(visitMapping)) {}
 
    std::size_t getNpar() const override { return _nParChip + _nParVisit; }
 
    double transform(MeasuredStar const &measuredStar, double value) const override {
        double temp = _chipMapping->getTransform()->transform(measuredStar.x, measuredStar.y, value);
        return _visitMapping->getTransform()->transform(measuredStar.getXFocal(), measuredStar.getYFocal(),
                                                        temp);
    }
 
    void freezeErrorTransform() override {
        _chipMapping->freezeErrorTransform();
        _visitMapping->freezeErrorTransform();
    }
 
    Eigen::VectorXd getParameters() override {
        Eigen::VectorXd joined(getNpar());
        joined << _chipMapping->getParameters(), _visitMapping->getParameters();
        return joined;
    }
 
    void getMappingIndices(IndexVector &indices) const override;
 
    void setWhatToFit(bool const fittingChips, bool const fittingVisits);
 
    void print(std::ostream &out) const override {
        out << "index: " << index << std::endl;
        out << "chip mapping: ";
        _chipMapping->print(out);
        out << std::endl << "visit mapping: ";
        _visitMapping->print(out);
    }
 
    std::shared_ptr<PhotometryMapping> getChipMapping() const { return _chipMapping; }
    std::shared_ptr<PhotometryMapping> getVisitMapping() const { return _visitMapping; }
 
    std::size_t getNParChip() const { return _nParChip; }
    std::size_t getNParVisit() const { return _nParVisit; }
 
protected:
    // These are either transform.getNpar() or 0, depending on whether we are fitting that component or not.
    std::size_t _nParChip, _nParVisit;
 
    // the actual transformation to be fit
    std::shared_ptr<PhotometryMapping> _chipMapping;
    std::shared_ptr<PhotometryMapping> _visitMapping;
};
 
class ChipVisitFluxMapping : public ChipVisitPhotometryMapping {
public:
    ChipVisitFluxMapping(std::shared_ptr<PhotometryMapping> chipMapping,
                         std::shared_ptr<PhotometryMapping> visitMapping)
            : ChipVisitPhotometryMapping(chipMapping, visitMapping) {}
 
    double transformError(MeasuredStar const &measuredStar, double value, double valueErr) const override;
 
    void computeParameterDerivatives(MeasuredStar const &measuredStar, double value,
                                     Eigen::Ref<Eigen::VectorXd> derivatives) const override;
};
 
class ChipVisitMagnitudeMapping : public ChipVisitPhotometryMapping {
public:
    ChipVisitMagnitudeMapping(std::shared_ptr<PhotometryMapping> chipMapping,
                              std::shared_ptr<PhotometryMapping> visitMapping)
            : ChipVisitPhotometryMapping(chipMapping, visitMapping) {}
 
    double transformError(MeasuredStar const &measuredStar, double value, double valueErr) const override;
 
    void computeParameterDerivatives(MeasuredStar const &measuredStar, double value,
                                     Eigen::Ref<Eigen::VectorXd> derivatives) const override;
};
 
}  // namespace jointcal
}  // namespace lsst
 
#endif  // LSST_JOINTCAL_PHOTOMETRY_MAPPING_H