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 dump(std::ostream &stream = std::cout) 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 dump(std::ostream &stream = std::cout) const override {
        stream << "index: " << index << " fixed: " << fixed << " transform parameters: ";
        _transform->dump(stream);
    }

    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 dump(std::ostream &stream = std::cout) const override {
        stream << "index: " << index << " chipMapping: ";
        _chipMapping->dump(stream);
        stream << "visitMapping: ";
        _visitMapping->dump(stream);
    }

    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