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

#include <stdexcept>
#include <sstream>
#include <vector>

#include "boost/format.hpp"
#include "boost/serialization/nvp.hpp"
#include "boost/serialization/vector.hpp"
#include "boost/serialization/void_cast.hpp"
#include "boost/serialization/export.hpp"

#include "lsst/daf/base/Citizen.h"
#include "lsst/pex/exceptions.h"

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

namespace lsst {
namespace afw {
namespace math {

#ifndef SWIG
using boost::serialization::make_nvp;
#endif

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

        explicit Function(
            std::vector<double> const &params)   
        :
            lsst::daf::base::Citizen(typeid(this)),
            _params(params),
            _isCacheValid(false)
        {}

        virtual ~Function() {}

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

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

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

        virtual bool isLinearCombination() const { 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;

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

        /* Default constructor: intended only for serialization */
        explicit Function() : lsst::daf::base::Citizen(typeid(this)), _params(0), _isCacheValid(false) {}

    private: // serialization support
        friend class boost::serialization::access;
        template <class Archive>
        void serialize(Archive& ar, unsigned int const version) {
            ar & make_nvp("params", _params);
        }
    };


    template<typename ReturnT>
    class Function1 : public afw::table::io::PersistableFacade< Function1<ReturnT> >,
                      public Function<ReturnT>
    {
    public:
        typedef std::shared_ptr<Function1<ReturnT> > Ptr;

        explicit Function1(
            unsigned int nParams)   
        :
            Function<ReturnT>(nParams)
        {}

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

        virtual ~Function1() {}

        virtual Ptr clone() const = 0;

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

        virtual std::string toString(std::string const& prefix="") const {
            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>() {}

    private: // serialization
        friend class boost::serialization::access;
#ifndef SWIG // SWIG doesn't like base_object
        template <class Archive>
        void serialize(Archive& ar, unsigned const int version) {
            ar & make_nvp("fn", boost::serialization::base_object<Function<ReturnT> >(*this));
        }
#endif
    };

    template<typename ReturnT>
    class Function2 : public afw::table::io::PersistableFacade< Function2<ReturnT> >,
                      public Function<ReturnT>
    {
    public:
        typedef std::shared_ptr<Function2<ReturnT> > Ptr;

        explicit Function2(
            unsigned int nParams)   
        :
            Function<ReturnT>(nParams)
        {}

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

        virtual ~Function2() {}

        virtual Ptr clone() const = 0;

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

        virtual std::string toString(std::string const& prefix="") const {
            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>() {}

    private:
        friend class boost::serialization::access;
#ifndef SWIG // SWIG doesn't like base_object
        template <class Archive>
        void serialize(Archive& ar, unsigned const int version) {
            ar & make_nvp("fn", boost::serialization::base_object<Function<ReturnT> >(*this));
        }
#endif
    };


    template<typename ReturnT>
    class BasePolynomialFunction2: public Function2<ReturnT> {
    public:
        typedef typename Function2<ReturnT>::Ptr Function2Ptr;

        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())))
        {}

        virtual ~BasePolynomialFunction2() {}

        int getOrder() const { return _order; }

        virtual bool isLinearCombination() const { 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;
        }

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

            Function2Ptr 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) {}

    private:
        friend class boost::serialization::access;
#ifndef SWIG // SWIG doesn't like base_object
        template <class Archive>
        void serialize(Archive& ar, unsigned const int version) {
            ar & make_nvp("fn2", boost::serialization::base_object<Function2<ReturnT> >(*this));
            ar & make_nvp("order", _order);
        }
#endif
    };


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

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

    private:
        friend class boost::serialization::access;
#ifndef SWIG // SWIG doesn't like base_object
        template <class Archive>
        void serialize(Archive& ar, unsigned int const version) {
            ar & make_nvp("fn1", boost::serialization::base_object<Function1<ReturnT> >(*this));
        }
#endif
    };

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

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

    private:
        friend class boost::serialization::access;
#ifndef SWIG // SWIG doesn't like base_object
        template <class Archive>
        void serialize(Archive& ar, unsigned int const version) {
            ar & make_nvp("fn2", boost::serialization::base_object<Function2<ReturnT> >(*this));
        }
#endif
    };


}}}   // lsst::afw::math

#endif // #ifndef LSST_AFW_MATH_FUNCTION_H