Basis.h

#include <cmath>
#include "json/json.h"
#include "../Grids/Grid.h"
 
namespace SSAGES
{
 
    struct BasisLUT
    {
        std::vector<double> values;
 
        std::vector<double> derivs;
 
 
        BasisLUT(const std::vector<double>& values,
            const std::vector<double>& derivs) :
            values(values), derivs(derivs)
        {}
    };
 
 
    struct Map
    {
        double value;
 
        std::vector<int> map;
 
 
        Map(const std::vector<int>& map,
            double value) :
            value(value), map(map)
        {}
    };
 
    class BasisFunction
    {
    protected:
        unsigned int polyOrd_; 
        unsigned int nbins_; 
        bool isFinite_; 
        bool zeroOrder_; 
        double boundLow_; 
        double boundUp_; 
 
    public:
 
        BasisFunction(unsigned int polyOrd,
                      unsigned int nbins,
                      bool isFinite,
                      bool zeroOrder,
                      double boundLow,
                      double boundUp) :
        polyOrd_(polyOrd), nbins_(nbins), isFinite_(isFinite), zeroOrder_(zeroOrder),
        boundLow_(boundLow), boundUp_(boundUp)
        {
        }
 
 
        unsigned int GetOrder() {return polyOrd_;}
 
 
        unsigned int GetBins() {return nbins_;}
 
 
        bool GetZeroOrder() {return zeroOrder_;}
 
 
        double GetLower() {return boundLow_;}
 
 
        double GetUpper() {return boundUp_;}
 
 
        double GetRange() 
        {
            if(isFinite_) 
                return boundUp_ - boundLow_;
            // No infinitely bounded basis functions are included currently so this is going to return nothing for right now
            else
                return 0.0;
        }
 
 
        virtual double GetNorm(int /* order */) {
            return 0;
        }
 
 
        virtual double Evaluate(double /* val */, int /* order */) {
            return 0;
        }
 
 
        virtual double EvalGrad(double /* val */, int /* order */) {
            return 0;
        }
 
 
        virtual double Weight(double /* val */) {
            return 1;
        }
 
 
        static BasisFunction* Build(const Json::Value& json, const std::string& path, unsigned int nbins);
 
        virtual ~BasisFunction() {}
    };
 
    class Chebyshev : public BasisFunction
    {
    private:  
    protected:
    public:
 
        Chebyshev(unsigned int polyOrd, double boundLow, double boundUp, unsigned int nbins) :
        BasisFunction(polyOrd, nbins, true, false, boundLow, boundUp)
        {
        }
 
 
        double ChangeVariable(double x)
        {
            return (x-boundLow_)*2.0/(boundUp_ - boundLow_)-1.0;
        }
 
        double Evaluate(double x, int n)
        {
            double xMod = ChangeVariable(x);
            return n == 0 ? 1.0 : 
                    n == 1 ? xMod :
                    2.0*xMod*Evaluate(x,n-1) - Evaluate(x,n-2);
        }
 
        double EvalGrad(double x, int n)
        {
            double xMod = ChangeVariable(x);
            return n == 0 ? 0.0 :
                    n == 1 ? 1.0 :
                    2.0*(Evaluate(x,n-1) + xMod * EvalGrad(x,n-1)) - EvalGrad(x,n-2);
        }
 
        double Weight(double x)
        {
            double xMod = ChangeVariable(x);
            return 1.0/sqrt(1.0-xMod*xMod);
        }
 
        double GetNorm(int /* n */)
        {
            return 2.0/M_PI;
        }
        
        static Chebyshev* Build(const Json::Value& json, const std::string& path, unsigned int nbins);
 
    };
 
    class Legendre : public BasisFunction
    {
    private:  
    protected:
    public:
 
        Legendre(unsigned int polyOrd, unsigned int nbins) :
        BasisFunction(polyOrd, nbins, true, false, -1.0, 1.0)
        {
        }
 
        double Evaluate(double x, int n)
        {
            return n == 0 ? 1.0 : 
                    n == 1 ? x :
                    (2.0*n-1.0)/(double)n*x*Evaluate(x,n-1) - (n-1.0)/(double)n*Evaluate(x,n-2);
        }
 
        double EvalGrad(double x, int n)
        {
            return n == 0 ? 0.0 :
                    n == 1 ? 1.0 :
                    (2.0*n-1.0)/(double)n*(Evaluate(x,n-1) + x * EvalGrad(x,n-1)) - (n-1.0)/(double)n*EvalGrad(x,n-2);
        }
 
        double GetNorm(int n)
        {
            return n + 0.5;
        }
 
        static Legendre* Build(const Json::Value& json, const std::string& path, unsigned int nbins);
    };
 
    class Fourier : public BasisFunction
    {
    private:
    protected:
    public:
 
        Fourier(unsigned int polyOrd, double boundLow, double boundUp, unsigned int nbins) :
        BasisFunction(polyOrd, nbins, true, true, boundLow, boundUp)
        {
        }
 
        //For the Fourier series, I am choosing to split the coefficients in two
        //The first odd coefficients are the sine series and the even are cosine
        //However when evaluated, they will be treated as n, just stored differently
        double Evaluate(double x, int n)
        {
            return n == 0 ? 0.5 : 
                    n % 2 == 1 ? std::sin(2.0*floor(double(n)*0.5)*M_PI*x/(boundUp_-boundLow_)) :
                      std::cos(2.0*floor(double(n)*0.5)*M_PI*x/(boundUp_-boundLow_));
        }
 
        double EvalGrad(double x, int n)
        {
            return n == 0 ? 0.0 :
                    n % 2 == 0 ? (-2.0*floor(double(n)*0.5)*M_PI/(boundUp_-boundLow_))*(std::sin(2.0*floor(double(n)*0.5)*M_PI*x/(boundUp_-boundLow_))) : 
                     (2.0*floor(double(n)*0.5)*M_PI/(boundUp_-boundLow_))*std::cos(2.0*floor(double(n)*0.5)*M_PI*x/(boundUp_-boundLow_));
        }
 
        double GetNorm(int /* n */)
        {
            //Assumes the period of the function is the range of the CV
            return 2.0/(boundUp_-boundLow_);
        }
 
        static Fourier* Build(const Json::Value& json, const std::string& path, unsigned int nbins);
    };
    
    class BasisEvaluator
    {
    private:
        std::vector<Map> coeff_; 
        std::vector<BasisFunction*> functions_; 
        std::vector<BasisLUT> lookup_; 
 
    public:
 
        BasisEvaluator(const std::vector<BasisFunction*>& functions) : 
            functions_(functions)
        {
            CoeffInit();
            BasisInit();
        }
    
        void CoeffInit(); 
        void BasisInit(); 
 
 
        void UpdateBias(Grid<double> *bias, Grid<std::vector<double>> *grad);
 
 
        double UpdateCoeff(const std::vector<double> &array, Grid<unsigned int> *hist);
 
 
        std::vector<double> GetCoeff(void)
        {
            std::vector<double> coeffArray (coeff_.size(),0);
            for (size_t i=0; i<coeff_.size(); i++)
            {
                coeffArray[i] = coeff_[i].value;
            }
            return coeffArray; 
        }
 
 
        void SetCoeff(const std::vector<double>& coeff)
        {
            size_t ii = 0;
            for(auto& val : coeff_) 
            {    
                val.value = coeff[ii];
                ii++;
            }
        }
 
        ~BasisEvaluator() {}
    };
}