Functions.h

//
//  Copyright (C) 2017, Alexander Stukowski and Paul Erhart
//
//  This file is part of atomicrex.
//
//  Atomicrex 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.
//
//  Atomicrex 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 <http://www.gnu.org/licenses/>.
//

#pragma once

#include "../../job/FitObject.h"
#include "../../util/NumericalIntegration.h"
#include "../../dof/DegreeOfFreedom.h"

namespace atomicrex {

class FunctionBase : public FitObject, public NumericalIntegration
{
public:

    enum NormalizationMode {
      NONE,
      SPHERICAL_INTEGRAL,
      UNITY_ORIGIN,
    };

public:

    FunctionBase(const FPString& id, FitJob* job, const FPString& tag = FPString())
        : FitObject(id, job, tag), _normalizationMode(NONE), _prefactor(0), _cutoff(FLOATTYPE_MAX) {}

    bool hasCutoff() const { return cutoff() != FLOATTYPE_MAX; }

    double cutoff() const { return _cutoff; }

    void setCutoff(double cutoff) { _cutoff = cutoff; }

    double prefactor() const { return _prefactor; }

    double evaluate(double r);

    double evaluate(double r, double& deriv);

    void writeTabulated(std::ostream& out, double rmin, double rmax, double dr);

    NormalizationMode normalizationMode() const { return _normalizationMode; }

    void setNormalizationMode(NormalizationMode mode) { _normalizationMode = mode; }

    void computePrefactor();

    virtual void dofValueChanged(DegreeOfFreedom& dof) override {
        FitObject::dofValueChanged(dof);
        _prefactor = 0;
    }

    const std::shared_ptr<FunctionBase>& screeningFunction() const { return _screeningFunction; }

    void setScreeningFunction(const std::shared_ptr<FunctionBase>& func) {
        BOOST_ASSERT(!_screeningFunction);
        _screeningFunction = func;
        registerSubObject(func.get());
    }

    void verifyDerivative();

    void checkDerivativeNumerically(double r);

    double centralDerivative(double r, double h, double& abserr_trunc, double& abserr_round);

public:

    virtual void parse(XML::Element functionElement) override;

    static std::shared_ptr<FunctionBase> createAndParse(XML::Element functionElement, const FPString& id, FitJob* job);

    virtual XML::OElement generateXMLDefinition(const FPString& elementName);

protected:

    virtual double evaluateInternal(double r) override = 0;

    virtual double evaluateInternal(double r, double& deriv) = 0;

    virtual std::vector<DegreeOfFreedom*> DOFToExport() { return DOF(); }

protected:

    double _cutoff;

    double _prefactor;

    NormalizationMode _normalizationMode;

    std::shared_ptr<FunctionBase> _screeningFunction;
};


/************************************************************
 * exp[1/(r-rc)]
 ************************************************************/
class FunctionExpA : public FunctionBase
{
public:
    FunctionExpA(const FPString& id, FitJob* job)
        : FunctionBase(id, job, "exp") {}

protected:
    double evaluateInternal(double r) override;
    double evaluateInternal(double r, double& deriv) override;
};


/************************************************************
 * exp[sgn(n)*alpha/(1-(r/rc)^n)]
 ************************************************************/
class FunctionExpB : public FunctionBase
{
public:
    FunctionExpB(const FPString& id, FitJob* job);

    virtual void parse(XML::Element functionElement) override;

    virtual XML::OElement generateXMLDefinition(const FPString& elementName) override;

protected:
    double evaluateInternal(double r) override;
    double evaluateInternal(double r, double& deriv) override;

private:
    ScalarDOF _alpha, _exponent, _rc;
};


/************************************************************
 * ExpB * gaussian
 ************************************************************/
class FunctionExpGaussian : public FunctionBase
{
public:
    FunctionExpGaussian(const FPString& id, FitJob* job);

    virtual void parse(XML::Element functionElement) override;

    virtual XML::OElement generateXMLDefinition(const FPString& elementName) override;

protected:
    double evaluateInternal(double r) override;
    double evaluateInternal(double r, double& deriv) override;

private:
    ScalarDOF _stddev, _alpha, _exponent;
};

/************************************************************
 * Standard Morse potential (type A).
 *
 * D0 * [exp(-2*alpha*(r-r0)) - 2*exp(-alpha*(r-r0))]
 ************************************************************/
class FunctionMorseA : public FunctionBase
{
public:

    FunctionMorseA(const FPString& id, FitJob* job);

    virtual void parse(XML::Element functionElement) override;

    virtual XML::OElement generateXMLDefinition(const FPString& elementName) override;

protected:
    double evaluateInternal(double r) override;
    double evaluateInternal(double r, double& deriv) override;

private:
    ScalarDOF _D0, _alpha, _r0;
};

/************************************************************
 * Modified Morse potential - type B.
 *
 * [D0 / (S-1) * exp(-beta sqrt(2 S) (r - r0)] - [D0 S / (S-1) * exp(-beta sqrt(2 S) (r - r0))] + delta
 ************************************************************/
class FunctionMorseB : public FunctionBase
{
public:
    FunctionMorseB(const FPString& id, FitJob* job);

    virtual void parse(XML::Element functionElement) override;

    virtual XML::OElement generateXMLDefinition(const FPString& elementName) override;

protected:
    double evaluateInternal(double r) override;
    double evaluateInternal(double r, double& deriv) override;

private:
    ScalarDOF _D0, _r0, _S, _beta, _delta;
};

/************************************************************
 * Modified Morse potential - type C.
 ************************************************************/
class FunctionMorseC : public FunctionBase
{
public:
    FunctionMorseC(const FPString& id, FitJob* job);

    virtual void parse(XML::Element functionElement) override;

    virtual XML::OElement generateXMLDefinition(const FPString& elementName) override;

protected:
    double evaluateInternal(double r) override;
    double evaluateInternal(double r, double& deriv) override;

private:
    ScalarDOF _A, _B, _mu, _lambda, _delta;
};

/************************************************************
 * Polynomial
 ************************************************************/
class FunctionPoly : public FunctionBase
{
public:

    FunctionPoly(const FPString& id, FitJob* job)
            : FunctionBase(id, job, "poly"), _isPrepared(false) {}

    double integrateSpherical(double a, double b) override;

    virtual void parse(XML::Element functionElement) override;

    virtual XML::OElement generateXMLDefinition(const FPString& elementName) override;

    virtual void dofValueChanged(DegreeOfFreedom& dof) override {
        FunctionBase::dofValueChanged(dof);
        if(&dof != &_nodes.back().coeff) {
            _isPrepared = false;
        }
    }

protected:

    double evaluateInternal(double r) override;
    double evaluateInternal(double r, double& deriv) override;

    void imposeBoundaryConditions();

    virtual std::vector<DegreeOfFreedom*> DOFToExport() override { return std::vector<DegreeOfFreedom*>(); }

private:

    struct PolyNode {
        ScalarDOF coeff;
        int exponent;
    };

    bool _isPrepared;

    std::vector<PolyNode> _nodes;

    int _degree;
};


/************************************************************
 * Cubic spline
 ************************************************************/
class FunctionSpline : public FunctionBase
{
public:

    FunctionSpline(const FPString& id, FitJob* job)
            : FunctionBase(id, job, "spline"), _isPrepared(false) {}

    virtual void parse(XML::Element functionElement) override;

    virtual XML::OElement generateXMLDefinition(const FPString& elementName) override;

    virtual void dofValueChanged(DegreeOfFreedom& dof) override {
        FunctionBase::dofValueChanged(dof);
        _isPrepared = false;
    }

protected:

    double evaluateInternal(double r) override;
    double evaluateInternal(double r, double& deriv) override;
    void prepareSpline();

    virtual std::vector<DegreeOfFreedom*> DOFToExport() override { return std::vector<DegreeOfFreedom*>(); }

private:

    struct SplineNode {
        double x, y2;
        ScalarDOF y;

        bool operator<(const SplineNode& other) const { return x < other.x; }       // Used for node sorting.
    };

    bool _isPrepared;

    // Vector of spline nodes.
    std::vector<SplineNode> _nodes;

    // Derivatives at left and right most nodes.
    ScalarDOF _deriv0, _derivN;
};

/************************************************************
 * Function of function of 1-x
 ************************************************************/
class FunctionInverseArgument : public FunctionBase
{
public:

    FunctionInverseArgument(const FPString& id, FitJob* job, const std::shared_ptr<FunctionBase>& func)
            : FunctionBase(id, job, "functionInverseArgument"), _innerFunction(func) {}

    const std::shared_ptr<FunctionBase>& innerFunction() const { return _innerFunction; }

protected:

    double evaluateInternal(double r) override {
        return _innerFunction->evaluate(1.0 - r);
    }

    double evaluateInternal(double r, double& deriv) override {
        double h = _innerFunction->evaluate(1.0 - r, deriv);
        deriv = -deriv;
        return h;
    }

private:
    std::shared_ptr<FunctionBase> _innerFunction;
};

/************************************************************
 * Sum of two or more functions.
 ************************************************************/
class FunctionSum : public FunctionBase
{
public:

    FunctionSum(const FPString& id, FitJob* job) : FunctionBase(id, job, "functionSum") {}

    virtual void parse(XML::Element functionElement) override;

protected:

    double evaluateInternal(double r) override {
        double s = 0;
        for(const auto& f : _functions)
            s += f->evaluate(r);
        return s;
    }

    double evaluateInternal(double r, double& deriv) override {
        double s = 0;
        deriv = 0;
        for(const auto& f : _functions) {
            double d;
            s += f->evaluate(r, d);
            deriv += d;
        }
        return s;
    }

private:

    std::vector<std::shared_ptr<FunctionBase>> _functions;
};

/************************************************************
 * A constant function, with the constant being a degree of freedom.
 ************************************************************/
class FunctionConstant : public FunctionBase
{
public:

    FunctionConstant(const FPString& id, FitJob* job);

    virtual void parse(XML::Element functionElement) override;

    virtual XML::OElement generateXMLDefinition(const FPString& elementName) override;

protected:

    double evaluateInternal(double r) override;
    double evaluateInternal(double r, double& deriv) override;

private:

    ScalarDOF _const;
};

/************************************************************
 * Product of two functions.
 ************************************************************/
class FunctionProduct : public FunctionBase
{
public:

    FunctionProduct(const FPString& id, FitJob* job) : FunctionBase(id, job, "functionProduct") {}

    virtual void parse(XML::Element functionElement) override;

protected:

    double evaluateInternal(double r) override {
        return _function1->evaluate(r) * _function2->evaluate(r);
    }

    double evaluateInternal(double r, double& deriv) override {
        double p1, p2;
        double d1, d2;
        p1 = _function1->evaluate(r, d1);
        p2 = _function2->evaluate(r, d2);
        deriv = p1 * d2 + p2 * d1;
        return p1 * p2;
    }

private:

    std::shared_ptr<FunctionBase> _function1;
    std::shared_ptr<FunctionBase> _function2;
};

} // End of namespace