Function.h

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/>.
 */

#ifndef LSST_AFW_MATH_FUNCTION_H
#define LSST_AFW_MATH_FUNCTION_H
/*
 * Define the basic Function classes.
 */
#include <cmath>
#include <stdexcept>
#include <sstream>
#include <vector>

#include "boost/format.hpp"

#include "lsst/pex/exceptions.h"

#include "lsst/afw/table/io/Persistable.h"

namespace lsst {
namespace afw {
namespace math {

template <typename T>
bool constexpr IS_NOTHROW_INIT = noexcept(static_cast<T>(1.0));

template <typename ReturnT>
class Function : public afw::table::io::PersistableFacade<Function<ReturnT>>,
                 public afw::table::io::Persistable {
public:
    explicit Function(unsigned int nParams)  
            : _params(nParams), _isCacheValid(false) {}

    explicit Function(std::vector<double> const& params)  
            : _params(params), _isCacheValid(false) {}

    Function(Function const&) = default;
    Function(Function&&) = default;
    Function& operator=(Function const&) = default;
    Function& operator=(Function&&) = default;

    ~Function() noexcept override = default;

    unsigned int getNParameters() const noexcept { return _params.size(); }

    virtual double getParameter(unsigned int ind)  
            const {
        return _params[ind];
    }

    std::vector<double> const& getParameters() const noexcept { return _params; }

    virtual bool isLinearCombination() const noexcept { return false; }

    void setParameter(unsigned int ind,  
                      double newValue)   
    {
        _isCacheValid = false;
        _params[ind] = newValue;
    }

    void setParameters(std::vector<double> const& params)  
    {
        if (_params.size() != params.size()) {
            throw LSST_EXCEPT(
                    pexExcept::InvalidParameterError,
                    (boost::format("params has %d entries instead of %d") % params.size() % _params.size())
                            .str());
        }
        _isCacheValid = false;
        _params = params;
    }

    virtual std::string toString(std::string const&) const {
        std::stringstream os;
        os << "parameters: [ ";
        for (std::vector<double>::const_iterator i = _params.begin(); i != _params.end(); ++i) {
            if (i != _params.begin()) os << ", ";
            os << *i;
        }
        os << " ]";
        return os.str();
    }

protected:
    std::vector<double> _params;
    mutable bool _isCacheValid;

    std::string getPythonModule() const override { return "lsst.afw.math"; }

    /* Default constructor: intended only for serialization */
    explicit Function() : _params(0), _isCacheValid(false) {}
};

template <typename ReturnT>
class Function1 : public afw::table::io::PersistableFacade<Function1<ReturnT>>, public Function<ReturnT> {
public:
    explicit Function1(unsigned int nParams)  
            : Function<ReturnT>(nParams) {}

    explicit Function1(std::vector<double> const& params)  
            : Function<ReturnT>(params) {}

    Function1(Function1 const&) = default;
    Function1(Function1&&) = default;
    Function1& operator=(Function1 const&) = default;
    Function1& operator=(Function1&&) = default;

    ~Function1() noexcept override = default;

    virtual std::shared_ptr<Function1<ReturnT>> clone() const = 0;

    virtual ReturnT operator()(double x) const = 0;

    std::string toString(std::string const& prefix = "") const override {
        return std::string("Function1: ") + Function<ReturnT>::toString(prefix);
    }

    virtual void computeCache(int const n) {}

protected:
    /* Default constructor: intended only for serialization */
    explicit Function1() : Function<ReturnT>() {}
};

template <typename ReturnT>
class Function2 : public afw::table::io::PersistableFacade<Function2<ReturnT>>, public Function<ReturnT> {
public:
    explicit Function2(unsigned int nParams)  
            : Function<ReturnT>(nParams) {}

    explicit Function2(std::vector<double> const& params)  
            : Function<ReturnT>(params) {}

    Function2(Function2 const&) = default;
    Function2(Function2&&) = default;
    Function2& operator=(Function2 const&) = default;
    Function2& operator=(Function2&&) = default;

    ~Function2() noexcept override = default;

    virtual std::shared_ptr<Function2<ReturnT>> clone() const = 0;

    virtual ReturnT operator()(double x, double y) const = 0;

    std::string toString(std::string const& prefix = "") const override {
        return std::string("Function2: ") + Function<ReturnT>::toString(prefix);
    }
    virtual std::vector<double> getDFuncDParameters(double, double) const {
        throw LSST_EXCEPT(lsst::pex::exceptions::NotFoundError,
                          "getDFuncDParameters is not implemented for this class");
    }

protected:
    /* Default constructor: intended only for serialization */
    explicit Function2() : Function<ReturnT>() {}
};

template <typename ReturnT>
class BasePolynomialFunction2 : public Function2<ReturnT> {
public:
    explicit BasePolynomialFunction2(unsigned int order)  
            : Function2<ReturnT>(BasePolynomialFunction2::nParametersFromOrder(order)), _order(order) {}

    explicit BasePolynomialFunction2(std::vector<double> params)  
            : Function2<ReturnT>(params),
              _order(BasePolynomialFunction2::orderFromNParameters(static_cast<int>(params.size()))) {}

    BasePolynomialFunction2(BasePolynomialFunction2 const&) = default;
    BasePolynomialFunction2(BasePolynomialFunction2&&) = default;
    BasePolynomialFunction2& operator=(BasePolynomialFunction2 const&) = default;
    BasePolynomialFunction2& operator=(BasePolynomialFunction2&&) = default;

    ~BasePolynomialFunction2() noexcept override = default;

    int getOrder() const noexcept { return _order; }

    bool isLinearCombination() const noexcept override { return true; }

    static int nParametersFromOrder(int order) {
        if (order < 0) {
            std::ostringstream os;
            os << "order=" << order << " invalid: must be >= 0";
            throw LSST_EXCEPT(lsst::pex::exceptions::InvalidParameterError, os.str());
        }
        return (order + 1) * (order + 2) / 2;
    }

    static int orderFromNParameters(int nParameters) {
        int order = static_cast<int>(
                0.5 + ((-3.0 + (std::sqrt(1.0 + (8.0 * static_cast<double>(nParameters))))) / 2.0));
        if (nParameters != BasePolynomialFunction2::nParametersFromOrder(order)) {
            std::ostringstream os;
            os << "nParameters=" << nParameters << " invalid: order is not an integer";
            throw LSST_EXCEPT(lsst::pex::exceptions::InvalidParameterError, os.str());
        }
        return order;
    }

    std::vector<double> getDFuncDParameters(double x, double y) const override {
        unsigned int const numParams = this->getNParameters();  // Number of parameters
        std::vector<double> deriv(numParams);                   // Derivatives, to return

        std::shared_ptr<Function2<ReturnT>> dummy =
                this->clone();  // Dummy function to evaluate for derivatives
        for (unsigned int i = 0; i < numParams; ++i) {
            dummy->setParameter(i, 0.0);
        }

        for (unsigned int i = 0; i < numParams; ++i) {
            dummy->setParameter(i, 1.0);
            deriv[i] = (*dummy)(x, y);
            dummy->setParameter(i, 0.0);
        }

        return deriv;
    }

protected:
    int _order;  

    /* Default constructor: intended only for serialization */
    explicit BasePolynomialFunction2() : Function2<ReturnT>(1), _order(0) {}
};

template <typename ReturnT>
class NullFunction1 : public Function1<ReturnT> {
public:
    explicit NullFunction1() : Function1<ReturnT>(0) {}
    std::shared_ptr<Function1<ReturnT>> clone() const override {
        return std::shared_ptr<Function1<ReturnT>>(new NullFunction1());
    }

private:
    ReturnT operator()(double) const noexcept(IS_NOTHROW_INIT<ReturnT>) override {
        return static_cast<ReturnT>(0);
    }
};

template <typename ReturnT>
class NullFunction2 : public Function2<ReturnT> {
public:
    explicit NullFunction2() : Function2<ReturnT>(0) {}
    std::shared_ptr<Function2<ReturnT>> clone() const override {
        return std::shared_ptr<Function2<ReturnT>>(new NullFunction2());
    }

private:
    ReturnT operator()(double, double) const noexcept(IS_NOTHROW_INIT<ReturnT>) override {
        return static_cast<ReturnT>(0);
    }
};
}  // namespace math
}  // namespace afw
}  // namespace lsst

#endif  // #ifndef LSST_AFW_MATH_FUNCTION_H