interpolation.h

Go to the documentation of this file.

#pragma once
//#include <iomanip>
 
#include "dg/backend/typedefs.h"
#include "dg/backend/view.h"
#include "grid.h"
#include "evaluation.h"
#include "functions.h"
#include "operator_tensor.h"
#include "xspacelib.h"
 
namespace dg{
 
namespace create{
namespace detail{
template<class real_type>
void shift( bool& negative, real_type & x, dg::bc bc, real_type x0, real_type x1)
{
    // This is now a free function because not all Topologies may have it
    if( bc == dg::PER)
    {
        real_type N0 = floor((x-x0)/(x1-x0)); // ... -2[ -1[ 0[ 1[ 2[ ...
        x = x == x1 ? x0 : x - N0*(x1-x0); //shift
    }
    //mirror along boundary as often as necessary
    while( (x<x0) || (x>x1) )
    {
        if( x < x0){
            x = 2.*x0 - x;
            //every mirror swaps the sign if Dirichlet
            if( bc == dg::DIR || bc == dg::DIR_NEU)
                negative = !negative;//swap sign
        }
        if( x >= x1){
            x = 2.*x1 - x;
            if( bc == dg::DIR || bc == dg::NEU_DIR) //notice the different boundary NEU_DIR to the above DIR_NEU !
                negative = !negative; //swap sign
        }
    }
}
 
template<class real_type>
std::vector<real_type> legendre( real_type xn, unsigned n)
{
    assert( xn <= 1. && xn >= -1.);
    std::vector<real_type> px(n);
    if( xn == -1)
    {
        for( unsigned u=0; u<n; u++)
            px[u] = (u%2 == 0) ? +1. : -1.;
    }
    else if( xn == 1)
    {
        for( unsigned i=0; i<n; i++)
            px[i] = 1.;
    }
    else
    {
        px[0] = 1.;
        if( n > 1)
        {
            px[1] = xn;
            for( unsigned i=1; i<n-1; i++)
                px[i+1] = ((real_type)(2*i+1)*xn*px[i]-(real_type)i*px[i-1])/(real_type)(i+1);
        }
    }
    return px;
}
 
// evaluate n base polynomials for n given abscissas
template<class real_type>
std::vector<real_type> lagrange( real_type x, const std::vector<real_type>& xi)
{
    unsigned n = xi.size();
    std::vector<real_type> l( n , 1.);
    for( unsigned i=0; i<n; i++)
        for( unsigned k=0; k<n; k++)
        {
            if ( k != i)
                l[i] *= (x-xi[k])/(xi[i]-xi[k]);
        }
    return l;
}
 
template<class real_type>
void choose_1d_abscissas( real_type X,
        unsigned points_per_line, double lx,// const RealGrid1d<real_type>& g,
        const thrust::host_vector<real_type>& abs, dg::bc bcx,
        thrust::host_vector<int>& cols,
        thrust::host_vector<real_type>& px
        )
{
    // Select points for nearest, linear or cubic interpolation
    // (Optimisation 9.9.25: the fact that abs are equidistant does not gain a speedup)
    // lx is needed for PER bondary conditions
    //determine which cell (X) lies in
    // abs must be sorted for std::lower_bound to work
    // std::lower_bound finds the first abs that is >= X
    auto it = std::lower_bound( abs.begin(), abs.end(), X);
    // Test if point already exists since then no interpolation is needed
    int idxX = -1;
    if( fabs( X - *it) < 1e-13)
        idxX = it - abs.begin();
    if( it != abs.begin() and fabs(  X - *(it-1)) < 1e-13)
        idxX = (it - abs.begin())-1;
    if( idxX >= 0)
    {
        cols.assign(1, idxX );
        px.assign( 1, 1.0);
        return;
    }
 
    std::vector<real_type> xs( points_per_line, 0);
    cols.resize( points_per_line, 0);
    switch( points_per_line)
    {
        case 1: xs[0] = 1.;
                if( it == abs.begin())
                    cols[0] = 0;
                else if( it == abs.end())
                    cols[0] = it - abs.begin() - 1;
                else
                {
                    if ( fabs(X - *it) < fabs( X - *(it-1)))
                        cols[0] = it - abs.begin();
                    else
                        cols[0] = it - abs.begin() -1;
                }
                break;
        case 2: if( it == abs.begin())
                {
                    if( bcx == dg::PER)
                    {
                        xs[0] = *it;
                        xs[1] = *(abs.end() -1)-lx;;
                        cols[0] = 0, cols[1] = abs.end()-abs.begin()-1;
                    }
                    else
                    {
                        // Nearest Neigbor: This makes it consistent with fem_t
                        xs.resize(1);
                        cols.resize(1);
                        xs[0] = *it;
                        cols[0] = 0;
                    }
                }
                else if( it == abs.end())
                {
                    if( bcx == dg::PER)
                    {
                        xs[0] = *(abs.begin())+lx;
                        xs[1] = *(it-1);
                        cols[0] = 0, cols[1] = it-abs.begin()-1;
                    }
                    else
                    {
                        // Nearest neighbor: this makes it consistent with fem_t
                        xs.resize(1);
                        cols.resize(1);
                        xs[0] = *(it-1);
                        cols[0] = it-abs.begin()-1;
                    }
                }
                else
                {
                    xs[0] = *(it-1);
                    xs[1] = *it;
                    cols[0] = it - abs.begin() - 1;
                    cols[1] = cols[0]+1;
                }
                break;
        case 4: if( it <= abs.begin() +1)
                {
                    if( bcx == dg::PER)
                    {
                        xs[0] = *abs.begin(), cols[0] = 0;
                        xs[1] = *(abs.begin()+1), cols[1] = 1;
                        xs[2] = it == abs.begin() ? *(abs.end() -2) -lx: *(abs.begin()+2);
                        cols[2] = it == abs.begin() ? abs.end()-abs.begin() -2 : 2;
                        xs[3] = *(abs.end() -1) -lx;
                        cols[3] = abs.end()-abs.begin() -1;
                    }
                    else
                    {   // use boundary polynomial
                        it = abs.begin();
                        xs[0] = *it,     xs[1] = *(it+1);
                        xs[2] = *(it+2), xs[3] = *(it+3);
                        cols[0] = 0, cols[1] = 1;
                        cols[2] = 2, cols[3] = 3;
                    }
                }
                else if( it >= abs.end() -1)
                {
                    if( bcx == dg::PER)
                    {
                        xs[0] = *abs.begin()+lx, cols[0] = 0;
                        xs[1] = it == abs.end() ? *(abs.begin()+1) + lx : *(abs.end() -3) ;
                        cols[1] = it == abs.end() ? 1 :  abs.end()-abs.begin()-3 ;
                        xs[2] = *(abs.end() - 2), cols[2] = abs.end()-abs.begin()-2;
                        xs[3] = *(abs.end() - 1), cols[3] = abs.end()-abs.begin()-1;
                    }
                    else
                    {   // use boundary polynomial
                        it = abs.end();
                        xs[0] = *(it-4), xs[1] = *(it-3);
                        xs[2] = *(it-2), xs[3] = *(it-1);
                        cols[0] = it - abs.begin() - 4;
                        cols[1] = cols[0]+1;
                        cols[2] = cols[1]+1;
                        cols[3] = cols[2]+1;
                    }
                }
                else
                {
                    xs[0] = *(it-2), xs[1] = *(it-1);
                    xs[2] = *(it  ), xs[3] = *(it+1);
                    cols[0] = it - abs.begin() - 2;
                    cols[1] = cols[0]+1;
                    cols[2] = cols[1]+1;
                    cols[3] = cols[2]+1;
                }
                break;
    }
    px = detail::lagrange( X, xs);
}
 
template<class real_type>
void interpolation_row( dg::space sp, real_type X,
    const RealGrid1d<real_type>& g, dg::bc bcx,
    thrust::host_vector<int>& cols,
    thrust::host_vector<real_type>& vals)
{
    bool negative = false;
    detail::shift( negative, X, bcx, g.x0(), g.x1());
 
    //determine which cell (x) lies in
    real_type xnn = (X-g.x0())/g.h();
    unsigned nn = (unsigned)floor(xnn);
    //determine normalized coordinates
    real_type xn = 2.*xnn - (real_type)(2*nn+1);
    //intervall correction
    if (nn==g.N()) { //only happens fo X == x1
        nn-=1;
        xn = 1.;
    }
    //Test if the point is a Gauss point since then no interpolation is needed
    int idxX =-1;
    std::vector<real_type> gauss_nodesx = dg::DLT<real_type>::abscissas(g.n());
    for( unsigned k=0; k<g.nx(); k++)
    {
        // Due to rounding errors xn will not perfectly be gauss_nodex[k] ( errors ~ 1e-14 possible)
        // even if x are the abscissas of the grid
        if( fabs( xn - gauss_nodesx[k]) < 1e-13)
            idxX = nn*g.nx() + k; //determine which grid column it is
    }
    if( idxX < 0 or sp == dg::lspace) //there is no corresponding point
    {
        //evaluate 1d Legendre polynomials at (xn, yn)...
        std::vector<real_type> px = detail::legendre( xn, g.n());
        //...these are the matrix coefficients with which to multiply
        if( sp == dg::xspace)
        {
            std::vector<real_type> pxF(px.size(), 0);
            std::vector<real_type> forward =
                dg::DLT<real_type>::forward(g.n());
            for( unsigned l=0; l<g.n(); l++)
                for( unsigned k=0; k<g.n(); k++)
                    pxF[l]+= px[k]*forward[k*g.n() + l];
            px.swap( pxF);
        }
        for ( unsigned l=0; l<g.n(); l++)
        {
            cols.push_back( nn*g.n() + l);
            vals.push_back( negative ? -px[l] : px[l]);
        }
    }
    else //the point already exists
    {
        cols.push_back( idxX);
        vals.push_back( negative ? -1. : 1.);
    }
}
 
 
// dG interpolation
template<class host_vector, class real_type>
dg::SparseMatrix<int, real_type, thrust::host_vector> interpolation1d(
    dg::space sp,
    const host_vector& x, // can be a view...
    const RealGrid1d<real_type>& g,
    dg::bc bcx )
{
    thrust::host_vector<real_type> values;
    thrust::host_vector<int> row_offsets;
    thrust::host_vector<int> column_indices;
    auto ptr = x.begin();
    row_offsets.push_back(0);
    for( unsigned i=0; i<x.size(); i++)
    {
        interpolation_row( sp, *ptr, g, bcx, column_indices, values);
        row_offsets.push_back( values.size());
        ptr++;
    }
    return {x.size(), g.size(), row_offsets, column_indices, values};
}
 
// nearest, linear or cubic interpolation
template<class host_vector1, class host_vector2 >
dg::SparseMatrix<int, dg::get_value_type<host_vector2>, thrust::host_vector> interpolation1d(
        const host_vector1& x,
        const host_vector2& abs, // must be sorted
        dg::bc bcx, dg::get_value_type<host_vector2> x0, dg::get_value_type<host_vector2> x1,
        std::string method )
{
    using real_type = dg::get_value_type<host_vector2>;
    // boundary condidions for dg::Box likely won't work | Box is now removed
    // from library ...
    thrust::host_vector<real_type> values;
    thrust::host_vector<int> row_offsets;
    thrust::host_vector<int> column_indices;
    unsigned points_per_line = 1;
    if( method == "nearest")
        points_per_line = 1;
    else if( method == "linear")
        points_per_line = 2;
    else if( method == "cubic")
        points_per_line = 4;
    else
        throw std::runtime_error( "Interpolation method "+method+" not recognized!\n");
    if( abs.size() < points_per_line)
        throw std::runtime_error( "There must be more points to interpolate\n");
    auto ptr = x.begin();
    row_offsets.push_back(0);
    thrust::host_vector<int> cols;
    thrust::host_vector<real_type> px;
    for( unsigned i=0; i<x.size(); i++)
    {
        row_offsets.push_back(row_offsets[i]);
        real_type X = *ptr;
        ptr++;
        bool negative = false;
        detail::shift( negative, X, bcx, x0, x1);
        // THIS IS A VERY BAD IDEA PERFORMANCE WISE
        //for( unsigned u=0; u<abs.size(); u++)
        //    if( fabs( X - abs[u]) <1e-13)
        //        idxX = u;
        detail::choose_1d_abscissas( X,
            points_per_line, x1-x0, abs, bcx, cols, px);
 
        // px may have size != points_per_line (at boundary)
        for ( unsigned l=0; l<px.size(); l++)
        {
            row_offsets[i+1]++;
            column_indices.push_back( cols[l]);
            values.push_back(negative ? -px[l] : px[l]);
        }
    }
    return {x.size(), abs.size(), row_offsets, column_indices, values};
}
 
}//namespace detail
 
template<class RecursiveHostVector, class real_type, size_t Nd>
dg::SparseMatrix<int, real_type, thrust::host_vector> interpolation(
        const RecursiveHostVector& x,
        const aRealTopology<real_type, Nd>& g,
        std::array<dg::bc, Nd> bcx,
        std::string method = "dg")
{
 
    std::array<dg::SparseMatrix<int,real_type,thrust::host_vector>,Nd> axes;
    for( unsigned u=0; u<Nd; u++)
    {
        if( x[u].size() != x[0].size())
            throw dg::Error( dg::Message(_ping_)<<"All coordinate lists must have same size "<<x[0].size());
        if( method == "dg")
        {
            axes[u] = detail::interpolation1d( dg::xspace, x[u], g.grid(u), bcx[u]);
        }
        else
        {
            axes[u] = detail::interpolation1d( x[u], g.abscissas(u), bcx[u], g.p(u), g.q(u), method);
        }
    }
    for( unsigned u=1; u<Nd; u++)
    {
        axes[0] = dg::tensorproduct_cols( axes[u], axes[0]);
    }
    return axes[0];
}
 
template<class host_vector, class real_type, typename = std::enable_if_t<dg::is_vector_v<host_vector>>>
dg::SparseMatrix<int, real_type, thrust::host_vector> interpolation(
        const host_vector& x,
        const RealGrid1d<real_type>& g,
        dg::bc bcx = dg::NEU,
        std::string method = "dg")
{
    dg::View<const host_vector> vx( x.data(), x.size());
    return interpolation(
        std::vector{vx}, g,
        std::array<dg::bc,1>{bcx}, method);
}
 
template<class host_vector, class real_type>
dg::SparseMatrix<int, real_type, thrust::host_vector> interpolation(
        const host_vector& x,
        const host_vector& y,
        const aRealTopology2d<real_type>& g,
        dg::bc bcx = dg::NEU, dg::bc bcy = dg::NEU,
        std::string method = "dg")
{
    dg::View<const host_vector> vx( x.data(), x.size());
    dg::View<const host_vector> vy( y.data(), y.size());
    return interpolation( std::vector{vx,vy}, g, {bcx,bcy}, method);
}
 
 
 
template<class host_vector, class real_type>
dg::SparseMatrix<int, real_type, thrust::host_vector> interpolation(
        const host_vector& x,
        const host_vector& y,
        const host_vector& z,
        const aRealTopology3d<real_type>& g,
        dg::bc bcx = dg::NEU, dg::bc bcy = dg::NEU, dg::bc bcz = dg::PER,
        std::string method = "dg")
{
    dg::View<const host_vector> vx( x.data(), x.size());
    dg::View<const host_vector> vy( y.data(), y.size());
    dg::View<const host_vector> vz( z.data(), z.size());
    return interpolation( std::vector{vx,vy,vz}, g, {bcx,bcy,bcz}, method);
}
template<class real_type, size_t Nd>
dg::SparseMatrix<int, real_type, thrust::host_vector> interpolation(
    const aRealTopology<real_type,Nd>& g_new,
    const aRealTopology<real_type,Nd>& g_old, std::string method = "dg")
{
    //assert both grids are on the same box
    for( unsigned u=0; u<Nd; u++)
    {
        if( g_new.p(u) < g_old.p(u))
            std::cerr << "ERROR: New grid boundary number "<<u<<" with value p "<<g_new.p(u)<<" lies outside old grid "<<g_old.p(u)<<" "<<g_old.p(u)-g_new.p(u)<<"\n";
        assert( g_new.p(u) >= g_old.p(u));
        if( g_new.q(u) > g_old.q(u))
            std::cerr << "ERROR: New grid boundary number "<<u<<" with value q "<<g_new.q(u)<<" lies outside old grid "<<g_old.q(u)<<" "<<g_old.q(u)-g_new.q(u)<<"\n";
        assert( g_new.q(u) <= g_old.q(u));
    }
    std::array<dg::SparseMatrix<int,real_type,thrust::host_vector>,Nd> axes;
    for( unsigned u=0; u<Nd; u++)
    {
        if( method == "dg")
        {
            axes[u] = detail::interpolation1d( dg::xspace, g_new.abscissas(u),
                g_old.grid(u), g_old.bc(u));
        }
        else
        {
            axes[u] = detail::interpolation1d( g_new.abscissas(u),
                g_old.abscissas(u), g_old.bc(u), g_old.p(u), g_old.q(u),
                method);
        }
    }
    for( unsigned u=1; u<Nd; u++)
        axes[0] = dg::tensorproduct( axes[u], axes[0]);
    return axes[0];
 
}
 
 
}//namespace create
 
 
template<class host_vector, class real_type>
real_type interpolate(
    dg::space sp,
    const host_vector& v,
    real_type x,
    const RealGrid1d<real_type>& g,
    dg::bc bcx = dg::NEU)
{
    assert( v.size() == g.size());
    thrust::host_vector<real_type> vals;
    thrust::host_vector<int> cols;
    create::detail::interpolation_row( sp, x, g, bcx, cols, vals);
    //multiply x
    real_type value = 0;
    for( unsigned j=0; j<vals.size(); j++)
        value += v[cols[j]]*vals[j];
    return value;
}
 
template<class host_vector, class real_type>
real_type interpolate(
    dg::space sp,
    const host_vector& v,
    real_type x, real_type y,
    const aRealTopology<real_type, 2>& g,
    dg::bc bcx = dg::NEU, dg::bc bcy = dg::NEU )
{
    assert( v.size() == g.size());
    thrust::host_vector<real_type> valsx;
    thrust::host_vector<int> colsx;
    create::detail::interpolation_row( sp, x, g.gx(), bcx, colsx, valsx);
    thrust::host_vector<real_type> valsy;
    thrust::host_vector<int> colsy;
    create::detail::interpolation_row( sp, y, g.gy(), bcy, colsy, valsy);
    //multiply x
    real_type value = 0;
    for( unsigned i=0; i<valsy.size(); i++)
        for( unsigned j=0; j<valsx.size(); j++)
            value += v[colsy[i]*g.shape(0) + colsx[j]]*valsx[j]*valsy[i];
    return value;
}
 
 
} //namespace dg