extrapolation.h

Go to the documentation of this file.

#pragma once
 
// #include <random>
#include "blas.h"
#include "topology/operator.h"
 
namespace dg
{
 
template<class ContainerType0, class ContainerType1, class ContainerType2>
std::vector<double> least_squares( const std::vector<ContainerType0>& bs, const ContainerType1 & b, const ContainerType2& weights)
{
    // it would be interesting to see how this algorithm fares against
    // Gram-Schmidt/QR-factorization
    // This implementation should have as many scalar dots as Gram-Schmidt
    // namely (size^2+size)/2
    // Solve B^T B a = B^T b
    unsigned size = bs.size();
    // B^T B is the "Gram matrix"
    dg::SquareMatrix<double> op( size, 0.); // B^T B
    std::vector<double> rhs( size, 0.);
    std::vector<double> a(size,0.);
    for( unsigned i=0; i<size; i++)
    {
        for( unsigned j=i; j<size; j++)
            op(i,j) = dg::blas2::dot( bs[i], weights, bs[j]);
        for( unsigned j=0; j<i; j++)
            op(i,j) = op(j,i);
        rhs[i] = dg::blas2::dot( bs[i], weights, b);
    }
    // possibly replace with Cholesky factorization?
    std::vector<unsigned> p;
    dg::create::lu_pivot( op, p);
    dg::lu_solve<double>( op, p, rhs);
    return rhs;
}
 
template<class ContainerType0, class ContainerType1>
std::vector<double> least_squares( const std::vector<ContainerType0>& bs, const ContainerType1 & b)
{
    return least_squares( bs, b, 1.);
}
template<class ContainerType0, class ContainerType1>
struct LeastSquaresExtrapolation
{
    using value_type = get_value_type<ContainerType0>;
    using container_type = ContainerType0;
    LeastSquaresExtrapolation( ){ m_counter = 0; }
    LeastSquaresExtrapolation( unsigned max, const ContainerType0& copyable0, const ContainerType1& copyable1) {
        set_max(max, copyable0, copyable1);
    }
    void set_max( unsigned max, const ContainerType0& copyable0,
                    const ContainerType1& copyable1)
    {
        m_counter = 0;
        m_x.assign( max, copyable0);
        m_y.assign( max, copyable1);
        m_max = max;
    }
    unsigned get_max( ) const{
        return m_counter;
    }
 
    void extrapolate( double alpha, const ContainerType0& x, double beta,
            ContainerType1& y) const{
        unsigned size = m_counter;
        thrust::host_vector<double> rhs( size, 0.), a(rhs), opIi(rhs); // B^T b
        for( unsigned i=0; i<size; i++)
            rhs[i] = dg::blas1::dot( m_x[i], x);
        // a = op_inv * rhs
        dg::blas1::scal( y, beta);
        for( unsigned i=0; i<size; i++)
        {
            for( unsigned j=0; j<size; j++)
                opIi[j] = m_op_inv(i,j);
            a[i] = dg::blas1::dot( rhs, opIi) ;
            dg::blas1::axpby( alpha*a[i], m_y[i], 1., y);
        }
    }
    void extrapolate( const ContainerType0& x, ContainerType1& y) const{
        extrapolate( 1., x, 0., y);
    }
 
    void update( const ContainerType0& x_new, const ContainerType1& y_new){
        if( m_max == 0) return;
        unsigned size = m_counter < m_max ? m_counter + 1 : m_max;
        dg::SquareMatrix<double> op( size, 0.), op_inv( size, 0.); // B^T B
        //i = 0
        op(0,0) = dg::blas1::dot( x_new, x_new);
        for( unsigned j=1; j<size; j++)
            op(0,j) = op( j, 0) = dg::blas1::dot( x_new, m_x[j-1]);
        // recursively fill in previous values
        for( unsigned i=1; i<size; i++)
            for( unsigned j=1; j<size; j++)
                op(i,j) = m_op(i-1,j-1);
        // test if new value is linearly independent or zero
        // maybe one can get a better control (with a tolerance value) over
        // this test
        try{
            op_inv = dg::create::inverse( op);
        }
        catch ( std::runtime_error & e){
            return;
        }
        m_op_inv = op_inv, m_op = op;
        if( m_counter < m_max)
            m_counter++;
        //push out last value (keep track of what is oldest value
        std::rotate( m_x.rbegin(), m_x.rbegin()+1, m_x.rend());
        std::rotate( m_y.rbegin(), m_y.rbegin()+1, m_y.rend());
        blas1::copy( x_new, m_x[0]);
        blas1::copy( y_new, m_y[0]);
    }
 
    private:
    unsigned m_max, m_counter;
    std::vector<ContainerType0> m_x;
    std::vector<ContainerType1> m_y;
    dg::SquareMatrix<double> m_op, m_op_inv;
};
 
template<class ContainerType, class real_type = dg::get_value_type<ContainerType>> // The real_type is if ContainerType is complex
struct Extrapolation
{
    using value_type = real_type;
    using container_type = ContainerType;
    Extrapolation( ){ m_counter = 0; }
    Extrapolation( unsigned max, const ContainerType& copyable) {
        set_max(max, copyable);
    }
    void set_max( unsigned max, const ContainerType& copyable)
    {
        m_counter = 0;
        m_x.assign( max, copyable);
        m_t.assign( max, 0);
        m_max = max;
    }
    unsigned get_max( ) const{
        return m_counter;
    }
 
 
    bool exists( value_type t)const{
        if( m_max == 0) return false;
        for( unsigned i=0; i<m_counter; i++)
            if( fabs(t - m_t[i]) <1e-14)
                return true;
        return false;
    }
 
    template<class ContainerType0>
    void extrapolate( value_type t, ContainerType0& new_x) const{
        switch(m_counter)
        {
            case(0): dg::blas1::copy( 0, new_x);
                     break;
            case(1): dg::blas1::copy( m_x[0], new_x);
                     break;
            case(3): {
                value_type f0 = (t-m_t[1])*(t-m_t[2])/(m_t[0]-m_t[1])/(m_t[0]-m_t[2]);
                value_type f1 =-(t-m_t[0])*(t-m_t[2])/(m_t[0]-m_t[1])/(m_t[1]-m_t[2]);
                value_type f2 = (t-m_t[0])*(t-m_t[1])/(m_t[2]-m_t[0])/(m_t[2]-m_t[1]);
                dg::blas1::evaluate( new_x, dg::equals(), dg::PairSum(),
                        f0, m_x[0], f1, m_x[1], f2, m_x[2]);
                 break;
            }
            default: {
                value_type f0 = (t-m_t[1])/(m_t[0]-m_t[1]);
                value_type f1 = (t-m_t[0])/(m_t[1]-m_t[0]);
                dg::blas1::axpby( f0, m_x[0], f1, m_x[1], new_x);
            }
        }
    }
 
    template<class ContainerType0>
    void derive( value_type t, ContainerType0& dot_x) const{
        switch(m_counter)
        {
            case(0): dg::blas1::copy( 0, dot_x);
                     break;
            case(1): dg::blas1::copy( 0, dot_x);
                     break;
            case(3): {
                value_type f0 =-(-2.*t+m_t[1]+m_t[2])/(m_t[0]-m_t[1])/(m_t[0]-m_t[2]);
                value_type f1 = (-2.*t+m_t[0]+m_t[2])/(m_t[0]-m_t[1])/(m_t[1]-m_t[2]);
                value_type f2 =-(-2.*t+m_t[0]+m_t[1])/(m_t[2]-m_t[0])/(m_t[2]-m_t[1]);
                dg::blas1::evaluate( dot_x, dg::equals(), dg::PairSum(),
                        f0, m_x[0], f1, m_x[1], f2, m_x[2]);
                break;
            }
            default: {
                value_type f0 = 1./(m_t[0]-m_t[1]);
                value_type f1 = 1./(m_t[1]-m_t[0]);
                dg::blas1::axpby( f0, m_x[0], f1, m_x[1], dot_x);
            }
        }
    }
 
    template<class ContainerType0>
    void extrapolate( ContainerType0& new_x) const{
        value_type t = m_t[0] +1.;
        extrapolate( t, new_x);
    }
 
 
    template<class MatrixType0, class ContainerType0, class ContainerType1>
    void extrapolate_least_squares( MatrixType0&& A, const ContainerType0& b,
            ContainerType0& new_x, const ContainerType1& weights)
    {
        if( m_counter < m_max)
        {
            extrapolate( new_x);
            return;
        }
        // allocate m_b if not yet done
        if( m_b.empty())
            m_b.assign( m_max, b);
        std::vector<const ContainerType*> x_ptrs = dg::asPointers( m_x);
        // An attempt at the algorithm in https://arxiv.org/abs/2309.02156
        //if( m_b.empty())
        //    m_b.assign( mm, b);
        // if( m_zx.empty())
        //     m_zx.assign( mm, b);
        // // compress the subspace via random linear combinations
        // if( mm < m_max)
        // {
        //     std::random_device rd{};
        //     std::mt19937 gen{rd()};
        //     std::normal_distribution<double> dist{0.0, 1.0};
        //     std::vector<double> z(m_max);
        //     for( unsigned i=0; i<mm; i++)
        //     {
        //         for( unsigned k=0; k<m_max; k++)
        //             z[k] = dist(gen);
        //         dg::blas2::gemv( 1., dg::asDenseMatrix(x_ptrs), z, 0., m_zx[i]);
        //     }
        //     x_ptrs = dg::asPointers( m_zx);
        // }
        // First compute bs
        for( unsigned u=0; u<m_max; u++)
            dg::apply( A, *x_ptrs[u], m_b[u]);
 
        try{
            std::vector<double> a = least_squares( m_b, b, weights);
            dg::blas2::gemv( 1., dg::asDenseMatrix(x_ptrs), a, 0., new_x);
        }
        catch( std::runtime_error& err)
        {
            return extrapolate( new_x);
        }
 
    }
 
 
    template<class ContainerType0>
    void derive( ContainerType0& dot_x) const{
        derive( m_t[0], dot_x);
    }
 
    template<class ContainerType0>
    void update( value_type t_new, const ContainerType0& new_entry){
        if( m_max == 0) return;
        //check if entry is already there to avoid division by zero errors
        for( unsigned i=0; i<m_counter; i++)
            if( fabs(t_new - m_t[i]) <1e-14)
            {
                blas1::copy( new_entry, m_x[i]);
                return;
            }
        if( m_counter < m_max) //don't update counter if Time entry was rejected
            m_counter++;
        //push out last value (keep track of what is oldest value
        std::rotate( m_x.rbegin(), m_x.rbegin()+1, m_x.rend());
        std::rotate( m_t.rbegin(), m_t.rbegin()+1, m_t.rend());
        m_t[0] = t_new;
        blas1::copy( new_entry, m_x[0]);
    }
    template<class ContainerType0>
    void update( const ContainerType0& new_entry){
        value_type t_new = m_t[0] + 1;
        update( t_new, new_entry);
    }
 
    const ContainerType& head()const{
        return m_x[0];
    }
    ContainerType& tail(){
        return m_x[m_max-1];
    }
    const ContainerType& tail()const{
        return m_x[m_max-1];
    }
 
    private:
    unsigned m_max, m_counter;
    std::vector<value_type> m_t;
    std::vector<ContainerType> m_x;
    //only allocated if least squares extrapolate is used
    std::vector<ContainerType> m_b;
};
 
 
}//namespace dg