ProSHADE_legendre.h

Go to the documentation of this file.

#ifndef __LEGENDRE_PROSHADE__
#define __LEGENDRE_PROSHADE__

//============================================ Standard Library
#include <iostream>
#include <cmath>
#include <vector>

//============================================ RVAPI
#include <rvapi_interface.h>

namespace ProSHADE_internal_legendre
{
    inline void getPolyAtZero ( unsigned int order,
                                double *polyValue,
                                double *deriValue )
    {
        //==================================== Initialise
        double hlpVal                         = 0.0;
        double prevPoly                       = 1.0;
        double prevPrevPoly                   = 0.0;
        double prevDeri                       = 0.0;
        double prevPrevDeri                   = 0.0;
        
        for ( unsigned int ordIt = 0; ordIt < order; ordIt++ )
        {
            hlpVal                            = static_cast<double> ( ordIt );
            *polyValue                        = -hlpVal * prevPrevPoly / ( hlpVal + 1.0 );
            *deriValue                        = ( ( 2.0 * hlpVal + 1.0 ) * prevPoly - hlpVal * prevPrevDeri ) / ( hlpVal + 1.0 );
            prevPrevPoly                      = prevPoly;
            prevPoly                          = *polyValue;
            prevPrevDeri                      = prevDeri;
            prevDeri                          = *deriValue;
        }
        
        //==================================== Done
        return ;
        
    }
    
    inline double advancePolyValue ( double from,
                                     double to,
                                     double valAtFrom,
                                     unsigned int noSteps )
    {
        //==================================== Initialise
        double hlpVal                         = 0.0;
        double stepSize                       = 0.0;
        double valChange                      = 0.0;
        double valSecChange                   = 0.0;
        double squareSteps                    = 0.0;
        double curVal                         = 0.0;
        unsigned int taylorSeriesCap          = 10;
        
        stepSize                              = ( to - from ) / static_cast<double> ( taylorSeriesCap );
        squareSteps                           = sqrt ( static_cast<double> ( noSteps * ( noSteps + 1 ) ) );
        curVal                                = from;
        
        for ( unsigned int iter = 0; iter < taylorSeriesCap; iter++ )
        {
            hlpVal                            = ( 1.0 - valAtFrom ) * ( 1.0 + valAtFrom );
            valChange                         = - stepSize * hlpVal / ( squareSteps * sqrt ( hlpVal ) - 0.5 * valAtFrom * sin ( 2.0 * curVal ) );
            valAtFrom                         = valAtFrom + valChange;
            
            curVal                            = curVal + stepSize;
            
            hlpVal                            = ( 1.0 - valAtFrom ) * ( 1.0 + valAtFrom );
            valSecChange                      = - stepSize * hlpVal / ( squareSteps * sqrt ( hlpVal ) - 0.5 * valAtFrom * sin ( 2.0 * curVal ) );
            valAtFrom                         = valAtFrom + 0.5 * ( valSecChange - valChange );
        }
        
        //==================================== Done
        return valAtFrom;
        
    }
    
    inline double evaluateSeries ( double *series,
                                   double target,
                                   unsigned int terms )
    {
        //==================================== Initalise
        double factorialValue                 = 1.0;
        double value                          = 0.0;
        
        for ( unsigned int iter = 1; iter <= terms; iter++ )
        {
            value                             = value + series[iter] * factorialValue;
            factorialValue                    = factorialValue * target;
        }
        
        //==================================== Done
        return ( value );
        
    }
    
    inline void getFirstEvenRoot ( double polyAtZero,
                                   unsigned int order,
                                   double *abscAtZero,
                                   double *weighAtZero )
    {
        //==================================== Initialise
        *abscAtZero                           = advancePolyValue ( 0.0, -M_PI / 2.0, 0.0, order );
        double hlp                            = 0.0;
        double hlpVal                         = static_cast<double> ( order );
        unsigned int taylorSeriesCap          = 30;
        double *abscSteps;
        double *weightSteps;

        //==================================== Allocate memory
        abscSteps                             = new double [taylorSeriesCap+2];
        weightSteps                           = new double [taylorSeriesCap+1];
        
        //==================================== Pre-set values
        abscSteps[0]                          = 0.0;
        abscSteps[1]                          = polyAtZero;
        weightSteps[0]                        = 0.0;
        
        //==================================== Fill in abscissa and weight steps
        for ( unsigned int iter = 0; iter <= taylorSeriesCap - 2; iter = iter + 2 )
        {
            hlp                               = static_cast<double> ( iter );
            
            abscSteps[iter+2]                 = 0.0;
            abscSteps[iter+3]                 = ( hlp * ( hlp + 1.0 ) - hlpVal * ( hlpVal + 1.0 ) ) * abscSteps[iter+1] / (hlp + 1.0) / (hlp + 2.0 );
            
            weightSteps[iter+1]               = 0.0;
            weightSteps[iter+2]               = ( hlp + 2.0 ) * abscSteps[iter+3];
        }
        
        //==================================== Find abscissa and weights
        for ( unsigned int iter = 0; iter < 5; iter++ )
        {
            *abscAtZero                       = *abscAtZero - evaluateSeries ( abscSteps, *abscAtZero, taylorSeriesCap ) / evaluateSeries ( weightSteps, *abscAtZero, taylorSeriesCap-1 );
        }
        *weighAtZero                          = evaluateSeries ( weightSteps,
                                                                *abscAtZero,
                                                                 taylorSeriesCap-1 );
        
        //==================================== Free memory
        delete abscSteps;
        delete weightSteps;
        
        //==================================== Done
        return ;
        
    }
    
    inline void completeLegendreSeries ( unsigned int order,
                                         std::vector<double>* abscissa,
                                         std::vector<double>* weights )
    {
        //==================================== Initialise
        double hlpTaylorVal                   = 0.0;
        double hlpOrderVal                    = static_cast<double> ( order );
        double abscValueChange                = 0.0;
        double prevAbsc                       = 0.0;
        double *hlpAbscSeries;
        double *hlpWeightSeries;
        unsigned int taylorSeriesCap          = 30;
        unsigned int noSeriesElems            = 0;
        unsigned int oddEvenSwitch            = 0;
        
        //==================================== Pre-set values
        if ( order % 2 == 1 )
        {
            noSeriesElems                     = ( order - 1 ) / 2 - 1;
            oddEvenSwitch                     = 1;
        }
        else
        {
            noSeriesElems                     = order / 2 - 1;
            oddEvenSwitch                     = 0;
        }
        
        //==================================== Allocate memory
        hlpAbscSeries                         = new double[taylorSeriesCap+2];
        hlpWeightSeries                       = new double[taylorSeriesCap+1];
        
        //==================================== For each series element
        for ( unsigned int serIt = noSeriesElems + 1; serIt < order - 1; serIt++ )
        {
            //================================ Init loop
            prevAbsc                          = abscissa->at(serIt);
            abscValueChange                   = advancePolyValue ( M_PI/2.0, -M_PI/2.0, prevAbsc, order ) - prevAbsc;
            
            //================================ Init abscissas
            hlpAbscSeries[0]                  = 0.0;
            hlpAbscSeries[1]                  = 0.0;
            hlpAbscSeries[2]                  = weights->at(serIt);
            
            //================================ Init weights
            hlpWeightSeries[0]                = 0.0;
            hlpWeightSeries[1]                = hlpAbscSeries[2];
            
            //================================ Taylor expansion
            for ( unsigned int tayIt = 0; tayIt <= taylorSeriesCap - 2; tayIt++ )
            {
                hlpTaylorVal                  = static_cast<double> ( tayIt );
                
                hlpAbscSeries[tayIt+3]        = ( 2.0 * prevAbsc * ( hlpTaylorVal + 1.0 ) * hlpAbscSeries[tayIt+2]
                                                  + ( hlpTaylorVal * ( hlpTaylorVal + 1.0 ) - hlpOrderVal * ( hlpOrderVal + 1.0 ) ) * hlpAbscSeries[tayIt+1] /
                                                    ( hlpTaylorVal + 1.0 ) ) / ( 1.0 - prevAbsc ) / ( 1.0 + prevAbsc ) / ( hlpTaylorVal + 2.0 );
                
                hlpWeightSeries[tayIt+2]      = ( hlpTaylorVal + 2.0 ) * hlpAbscSeries[tayIt+3];
            }
            
            //================================ Sum over results
            for ( unsigned int iter = 0; iter < 5; iter++ )
            { 
                abscValueChange               = abscValueChange - evaluateSeries ( hlpAbscSeries,   abscValueChange, taylorSeriesCap   ) /
                                                                  evaluateSeries ( hlpWeightSeries, abscValueChange, taylorSeriesCap-1 );
            }
            
            //================================ Save results
            abscissa->at(serIt+1)             = prevAbsc + abscValueChange;
            weights->at(serIt+1)              = evaluateSeries ( hlpWeightSeries, abscValueChange, taylorSeriesCap - 1 );
        }
        
        for ( unsigned int serIt = 0; serIt <= noSeriesElems + oddEvenSwitch; serIt++ )
        {
            abscissa->at(serIt)               = - abscissa->at(order-serIt-1);
            weights->at(serIt)                = weights->at(order-serIt-1);
        }
        
        //==================================== Free memory
        delete hlpAbscSeries;
        delete hlpWeightSeries;
        
        //==================================== Done
        return ;
        
    }
    
    inline void getLegendreAbscAndWeights ( unsigned int order,
                                            std::vector<double>* abscissa,
                                            std::vector<double>* weights,
                                            ProSHADE::ProSHADE_settings* settings )
    {
        //==================================== Sanity check
        if ( order < 1 )
        {
            std::cerr << "!!! ProSHADE ERROR !!! Trying to compute legendre quadrature with order " << order << " which seems too low!" << std::endl;
            
            if ( settings->htmlReport )
            {
                std::stringstream hlpSS;
                hlpSS << "<font color=\"red\">" << "The automatic determination of the Gauss-Legendre integration order has determined the optimal order to be less than 1. This may happen for very small maps, if this is the case, please manually set the integration value to 10-ish and re-run. If the map is reasonably large, please report this as a bug." << "</font>";
                rvapi_set_text                    ( hlpSS.str().c_str(),
                                                   "ProgressSection",
                                                   settings->htmlReportLineProgress,
                                                   1,
                                                   1,
                                                   1 );
                settings->htmlReportLineProgress += 1;
                rvapi_flush                       ( );
            }
            
            exit ( -1 );
        }
        
        //==================================== Initialise
        double polyValue                      = 0.0;
        double deriValue                      = 0.0;
        double weightSum                      = 0.0;
        
        //==================================== Find the polynomial and derivative values at 0
        getPolyAtZero                         ( order,
                                               &polyValue,
                                               &deriValue );
        
        //==================================== If the order is odd, then 0 is a root ...
        if ( order % 2 == 1 )
        {
            abscissa->at((order-1)/2)         = polyValue;
            weights->at((order-1)/2)          = deriValue;
        }
        else
        {
            // ... and if order is even, find the first root
            getFirstEvenRoot                  ( polyValue, order, &abscissa->at(order/2), &weights->at(order/2) );
        }
        
        //==================================== Now, having computed the first roots, complete the series
        completeLegendreSeries                ( order,
                                                abscissa,
                                                weights );
        
        //==================================== Correct weights by anscissa values
        for ( unsigned int iter = 0; iter < order; iter++ )
        {
            weights->at(iter)                 = 2.0 / ( 1.0 - abscissa->at(iter) ) / ( 1.0 + abscissa->at(iter) ) / weights->at(iter) / weights->at(iter);
            weightSum                         = weightSum + weights->at(iter);
        }
        
        //==================================== Normalise weights
        for ( unsigned int iter = 0; iter < order; iter++ )
        {
            weights->at(iter)                 = 2.0 * weights->at(iter) / weightSum;
        }
        
        //==================================== Done
        return ;
        
    }
}

#endif