rbf.cpp

/*---------------------------------------------------------------------------*\
 *
 *  bitpit
 *
 *  Copyright (C) 2015-2021 OPTIMAD engineering Srl
 *
 *  -------------------------------------------------------------------------
 *  License
 *  This file is part of bitpit.
 *
 *  bitpit is free software: you can redistribute it and/or modify it
 *  under the terms of the GNU Lesser General Public License v3 (LGPL)
 *  as published by the Free Software Foundation.
 *
 *  bitpit 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 Lesser General Public
 *  License for more details.
 *
 *  You should have received a copy of the GNU Lesser General Public License
 *  along with bitpit. If not, see <http://www.gnu.org/licenses/>.
 *
\*---------------------------------------------------------------------------*/
 
#include <cassert>
#include <cmath>
#include <set>
 
#include "bitpit_private_lapacke.hpp"
 
#include "bitpit_common.hpp"
#include "bitpit_operators.hpp"
 
#include "rbf.hpp"
 
namespace bitpit {
 
RBFKernel::RBFKernel()
    : m_supportRadii(1, 1.)
{
    m_nodes         = 0;
    m_fields        = 0;
 
    m_mode = RBFMode::INTERP;
 
    m_maxFields = -1;
    m_value.clear();
    m_weight.clear();
    m_activeNodes.clear();
 
    setFunction( RBFBasisFunction::WENDLANDC2);
 
    m_polyEnabled = false;
    m_polyActiveBasis.clear();
    m_polynomial.clear();
}
 
RBFKernel::RBFKernel(const RBFKernel & other)
    : m_fields(other.m_fields), m_mode(other.m_mode),
      m_supportRadii(other.m_supportRadii), m_typef(other.m_typef),
      m_fPtr(other.m_fPtr), m_error(other.m_error), m_value(other.m_value),
      m_weight(other.m_weight), m_activeNodes(other.m_activeNodes),
      m_maxFields(other.m_maxFields), m_nodes(other.m_nodes),
      m_polyEnabled(other.m_polyEnabled), m_polyActiveBasis(other.m_polyActiveBasis), 
      m_polynomial(other.m_polynomial)
{
}
 
void RBFKernel::swap(RBFKernel &other) noexcept
{
   std::swap(m_fields, other.m_fields);
   std::swap(m_mode, other.m_mode);
   std::swap(m_supportRadii, other.m_supportRadii);
   std::swap(m_typef, other.m_typef);
   std::swap(m_fPtr, other.m_fPtr);
   std::swap(m_error, other.m_error);
   std::swap(m_value, other.m_value);
   std::swap(m_weight, other.m_weight);
   std::swap(m_activeNodes, other.m_activeNodes);
   std::swap(m_maxFields, other.m_maxFields);
   std::swap(m_nodes, other.m_nodes);
   std::swap(m_polyEnabled, other.m_polyEnabled);
   std::swap(m_polyActiveBasis, other.m_polyActiveBasis);
   std::swap(m_polynomial, other.m_polynomial);
}
 
void RBFKernel::setFunction( RBFBasisFunction bfunc )
{
    switch(bfunc){
 
    case( RBFBasisFunction::WENDLANDC2):
        setFunction( rbf::wendlandc2);
        m_typef = RBFBasisFunction::WENDLANDC2;
        break;
 
    case( RBFBasisFunction::LINEAR):
        setFunction( rbf::linear);
        m_typef = RBFBasisFunction::LINEAR;
        break;
 
    case( RBFBasisFunction::GAUSS90):
        setFunction( rbf::gauss90);
        m_typef = RBFBasisFunction::GAUSS90;
        break;
 
    case( RBFBasisFunction::GAUSS95):
        setFunction( rbf::gauss95);
        m_typef = RBFBasisFunction::GAUSS95;
        break;
 
    case( RBFBasisFunction::GAUSS99):
        setFunction( rbf::gauss99);
        m_typef = RBFBasisFunction::GAUSS99;
        break;
 
    case( RBFBasisFunction::C1C0):
        setFunction( rbf::c1c0);
        m_typef = RBFBasisFunction::C1C0;
        break;
 
    case( RBFBasisFunction::C2C0):
        setFunction( rbf::c2c0);
        m_typef = RBFBasisFunction::C2C0;
        break;
 
    case( RBFBasisFunction::C0C1):
        setFunction( rbf::c0c1);
        m_typef = RBFBasisFunction::C0C1;
        break;
 
    case( RBFBasisFunction::C1C1):
        setFunction( rbf::c1c1);
        m_typef = RBFBasisFunction::C1C1;
        break;
 
    case( RBFBasisFunction::C2C1):
        setFunction( rbf::c2c1);
        m_typef = RBFBasisFunction::C2C1;
        break;
 
    case( RBFBasisFunction::C0C2):
        setFunction( rbf::c0c2);
        m_typef = RBFBasisFunction::C0C2;
        break;
 
    case( RBFBasisFunction::C1C2):
        setFunction( rbf::c1c2);
        m_typef = RBFBasisFunction::C1C2;
        break;
 
    case( RBFBasisFunction::C2C2):
        setFunction( rbf::c2c2);
        m_typef = RBFBasisFunction::C2C2;
        break;
 
    case( RBFBasisFunction::COSINUS):
        setFunction( rbf::cosinus);
        m_typef = RBFBasisFunction::COSINUS;
        break;
 
    case (RBFBasisFunction::THINPLATE):
        setFunction(rbf::thinplate);
        m_typef = RBFBasisFunction::THINPLATE;
        break;
 
    default:
        setFunction( rbf::wendlandc2);
        m_typef = RBFBasisFunction::WENDLANDC2;
        break;
    }
}
 
void RBFKernel::setFunction( double (&bfunc)(double ) )
{
    m_fPtr = bfunc;
    m_typef = RBFBasisFunction::CUSTOM;
}
 
RBFBasisFunction RBFKernel::getFunctionType(  )
{
    return m_typef;
}
 
int RBFKernel::getDataCount(  )
{
    return m_fields;
}
 
int RBFKernel::getActiveCount(  )
{
    int nActive(0);
 
    for( bool active : m_activeNodes)
        nActive += (int) active;
 
    return nActive;
}
 
std::vector<int> RBFKernel::getActiveSet(  )
{
    int                 i(0);
    std::vector<int>    activeSet;
 
    activeSet.reserve( getActiveCount() );
 
    for( bool active : m_activeNodes) {
        if( active )
            activeSet.push_back(i);
        ++i;
    }
 
    return activeSet;
}
 
bool RBFKernel::isActive( int n )
{
    if(n<0 || n >= int(m_activeNodes.size())) {
        return false;
    }
 
    return m_activeNodes[n];
}
 
bool RBFKernel::activateNode(int  n)
{
    bool check = false;
    if(n>=0 && n<m_nodes) {
        m_activeNodes[n] = true;
        check = true;
    }
    return check;
}
 
bool RBFKernel::activateNode(const std::vector<int> & list)
{
    if(list.empty()) {
        return false;
    }
 
    bool check = true;
    for(int index : list) {
        check = check && activateNode(index);
    }
    return check;
}
 
void RBFKernel::activateAllNodes()
{
    std::fill(m_activeNodes.begin(), m_activeNodes.end(), true);
}
 
bool RBFKernel::deactivateNode(int  n )
{
    bool check = false;
    if(n>=0 && n<m_nodes) {
        m_activeNodes[n] = false;
        check=  true;
    }
    return check;
}
 
bool RBFKernel::deactivateNode(const std::vector<int> & list)
{
    if(list.empty()) {
        return false;
    }
 
    bool check = true;
    for(int index : list) {
        check = check && deactivateNode(index);
    }
    return check;
}
 
void RBFKernel::deactivateAllNodes()
{
    std::fill(m_activeNodes.begin(), m_activeNodes.end(), false);
}
 
void RBFKernel::setSupportRadius( double  radius )
{
    m_supportRadii.resize(1);
    m_supportRadii[0] = radius;
}
 
void RBFKernel::setSupportRadius( const std::vector<double> & radius )
{
    m_supportRadii = radius;
}
 
double RBFKernel::getSupportRadius()
{
    return m_supportRadii[0];
}
 
double RBFKernel::getSupportRadius(int i)
{
    bool variableSupportRadius = (m_supportRadii.size() != 1);
    if (!variableSupportRadius) {
        return m_supportRadii[0];
    }
 
    return m_supportRadii[i];
}
 
RBFMode RBFKernel::getMode()
{
    return m_mode;
}
 
void RBFKernel::setMode(RBFMode mode)
{
    m_mode = mode;
}
 
void RBFKernel::setDataToNode( int id, const std::vector<double> &value )
{
    if(id<0 || id >= m_fields) {
        return;
    }
 
    if((int)(value.size()) != m_fields) {
        log::cout() << "Mismatch dimension between value vector size and number of data attached to rbf.";
        log::cout() << "This may lead to nasty errors. Check it with getDataCount()!" << std::endl;
        log::cout() << "Data could not be set" << std::endl;
        return;
    }
 
    int i;
    if(m_mode != RBFMode::PARAM) {
        for( i=0; i<m_fields; ++i ) {
                m_value[i][id] = value[i];
        }
    }else{
        for( i=0; i<m_fields; ++i ) {
            m_weight[i][id] = value[i];
        }
    }
}
 
void RBFKernel::setDataToAllNodes( int id, const std::vector<double> &value )
{
    if(id<0 || id >= m_fields) {
        return;
    }
 
    int size = m_value[id].size();
 
    if((int)(value.size()) != size) {
        log::cout() << "Mismatch dimension between data vector and current data container. One or both does not match RBFKernel nodes count.";
        log::cout() << "This may lead to nasty errors. Use fitDataToNodes to reshape container or fit your data vector first!" << std::endl;
        log::cout() << "Data could not be set" << std::endl;
        return;
    }
    if(m_mode != RBFMode::PARAM) {
        m_value[id] = value;
    }else{
        m_weight[id] = value;
    }
}
 
int RBFKernel::addData( )
{
    if(m_fields == m_maxFields) {
        log::cout() << "Maximum number of data set reached" << std::endl;
        return -1;
    }
    m_fields++;
    fitDataToNodes(m_fields-1);
    return m_fields;
}
 
int RBFKernel::addData( const std::vector<double> & data )
{
    if(m_fields == m_maxFields) {
        log::cout() << "Maximum number of data set reached" << std::endl;
        return -1;
    }
    if(m_mode == RBFMode::INTERP)    m_value.push_back(data);
    else                            m_weight.push_back(data);
    m_fields++;
 
    return m_fields;
}
 
bool RBFKernel::removeData(int id)
{
    if(id<0 || id >=m_fields) {
        return false;
    }
 
    m_fields--;
    if(m_mode == RBFMode::INTERP)    m_value.erase(m_value.begin()+id);
    else                             m_weight.erase(m_weight.begin()+id);
 
    return(true);
}
 
bool RBFKernel::removeData(std::vector<int> & list)
{
    // List should be processed in reversed order because data
    // are removed using their position in the storages.
    std::sort(list.begin(), list.end(), std::greater<int>());
 
    int nInitialFields = getDataCount();
    for(int id : list) {
        removeData(id);
    }
    int nFinalFields = getDataCount();
 
    return ((nInitialFields - nFinalFields) == static_cast<int>(list.size()));
}
 
void RBFKernel::removeAllData()
{
    m_fields = 0;
    if(m_mode == RBFMode::INTERP)    m_value.clear();
    else                            m_weight.clear();
}
 
std::vector<double> RBFKernel::evalRBF( const std::array<double,3> &point)
{
    std::vector<double> values(m_fields, 0.);
    int                 i, j;
    double              dist, basis;
 
    for( i=0; i<m_nodes; ++i ){
        if( m_activeNodes[i] ) {
            dist = calcDist(point, i) / getSupportRadius(i);
            basis = evalBasis( dist );
 
            for( j=0; j<m_fields; ++j) {
                values[j] += basis * m_weight[j][i];
            }
        }
    }
 
    // If INTERP mode add the polynomial contribution
    if (m_mode == RBFMode::INTERP && m_polyEnabled) {
        for (j = 0; j < m_fields; ++j) {
            double *polynomialCoefficients = m_polynomial.getCoefficients(j);
            std::vector<double> basis      = evalPolynomialBasis(point);
            int z                          = 0;
            for (int iactive : m_polyActiveBasis) {
                values[j] += polynomialCoefficients[iactive] * basis[z];
                z++;
            }
        }
    }
 
    return values;
}
 
std::vector<double> RBFKernel::evalRBF(int jnode)
{
    std::vector<double> values(m_fields, 0.);
    int                 i, j;
    double              dist, basis;
 
    if(jnode<0 || jnode>= m_nodes ) {
        return values;
    }
 
    for( i=0; i<m_nodes; ++i ) {
        if( m_activeNodes[i] ) {
            dist = calcDist(jnode, i) / getSupportRadius(i);
            basis = evalBasis( dist );
 
            for( j=0; j<m_fields; ++j) {
                values[j] += basis * m_weight[j][i];
            }
        }
    }
 
    // If INTERP mode add the polynomial contribution
    if (m_mode == RBFMode::INTERP && m_polyEnabled) {
        for (j = 0; j < m_fields; ++j) {
            double *polynomialCoefficients = m_polynomial.getCoefficients(j);
            std::vector<double> basis      = evalPolynomialBasis(jnode);
            int z                          = 0;
            for (int iactive : m_polyActiveBasis) {
                values[j] += polynomialCoefficients[iactive] * basis[z];
                z++;
            }
        }
    }
 
    return values;
}
 
int RBFKernel::solve()
{
    if(m_mode == RBFMode::PARAM) {
        return -1;
    }
 
    if (m_polyEnabled) {
        // Initialize polynomial
        initializePolynomial();
 
        // Check which parameter of the polynomial has to be activated
        // in order to avoid undetermined system
        initializePolynomialActiveBasis();
    }
 
    double dist;
 
    std::vector<int> activeSet = getActiveSet();
    int nActive = activeSet.size();
 
    int nPoly   = m_polyEnabled ? m_polyActiveBasis.size() : 0;
    int nS      = nActive + nPoly;
    int nrhs    = getDataCount();
 
    int lda     = nS;
    int ldb     = nS;
    int info;
 
    std::vector<int> ipiv(nS);
 
 
    std::vector<double> A(lda * nS);
    std::vector<double> b(ldb * nrhs);
 
    for (int j = 0; j < nrhs; ++j) {
        for (const auto &i : activeSet) {
            int k = j * ldb + i;
            b[k]  = m_value[j][i];
        }
    }
 
    int k = 0;
    for (const auto &i : activeSet) {
        for (const auto &j : activeSet) {
            dist                = calcDist(j, i) / getSupportRadius(i); //order by column!
            int row             = k % nActive;
            int col             = k / nActive;
            A[(col * nS) + row] = evalBasis(dist);
            k++;
        }
    }
 
    // Filling terms given by the polynomial block.
    // Symmetric terms set in the loop.
    if (m_polyEnabled) {
        for (int i = 0; i < nActive; ++i) {
            std::vector<double> polynomialTerms = evalPolynomialBasis(activeSet[i]);
            int j                               = 0;
            for (const double &val : polynomialTerms) {
                A[(j + nActive) * nS + i] = val;
                A[i * nS + (j + nActive)] = val;
                j++;
            }
        }
    }
 
    info = LAPACKE_dgesv(LAPACK_COL_MAJOR, nS, nrhs, A.data(), lda, ipiv.data(), b.data(), ldb);
 
    if( info > 0 ) {
        printf( "The diagonal element of the triangular factor of the linear system matrix \n" );
        printf( "U(%i,%i) is zero, so that matrix is singular;\n", info, info );
        printf( "the solution could not be computed.\n" );
        return 1;
    }
 
    m_weight.resize(nrhs);
 
    for (int j = 0; j < nrhs; ++j) {
        m_weight[j].resize(m_nodes, 0);
        int k = 0;
        for (const auto &i : activeSet) {
            m_weight[j][i] = b[j * ldb + k];
            ++k;
        }
    }
 
    if (m_polyEnabled) {
        for (int j = 0; j < nrhs; ++j) {
            double *polynomialCoefficients = m_polynomial.getCoefficients(j);
            int i                          = 0;
            for (int iactive : m_polyActiveBasis) {
                polynomialCoefficients[iactive] = b[j * ldb + nActive + i];
                i++;
            }
        }
    }
    return 0;
}
 
int RBFKernel::greedy( double tolerance)
{
    if(m_mode == RBFMode::PARAM)    return -1;
 
    int                     i, j;
    double                  error(1.e18);
    std::vector<double>     local(m_fields);
 
    m_error.resize(m_nodes);
 
    deactivateAllNodes();
 
    for( i=0; i<m_nodes; ++i){
 
        for( j=0; j<m_fields; ++j) {
            local[j] = m_value[j][i];
        }
 
        m_error[i] = norm2(local);
    }
 
    std::ios::fmtflags streamFlags(log::cout().flags());
 
    int errorFlag = 0;
    while( error > tolerance) {
        i = addGreedyPoint();
 
        if( i != -1) {
            m_activeNodes[i] = true;
 
            //solve();
            solveLSQ();
 
            error = evalError();
 
            log::cout() << std::scientific;
            log::cout() << " error now " << error << " active nodes" << getActiveCount() << " / " << m_nodes << std::endl;
        } else {
            errorFlag = 1;
            break;
        }
    }
 
    log::cout().flags(streamFlags);
 
    return errorFlag;
}
 
const std::vector<std::vector<double>> & RBFKernel::getWeights() const
{
    return m_weight;
}
 
void RBFKernel::fitDataToNodes()
{
    for (int i=0;i<m_fields; ++i) {
        fitDataToNodes(i);
    }
}
 
void RBFKernel::fitDataToNodes(int id)
{
    if(m_mode != RBFMode::PARAM)    m_value[id].resize(m_nodes, 0.0);
    else                            m_weight[id].resize(m_nodes,0.0);
}
 
double RBFKernel::evalBasis( double dist )
{
    return (*m_fPtr)(dist);
}
 
double RBFKernel::evalBasisPair(const int i, const int j)
{
    double dist = calcDist(i, j) / getSupportRadius(j);
    double value = evalBasis(dist);
    return value;
}
 
int RBFKernel::addGreedyPoint( )
{
    if(m_mode == RBFMode::PARAM) {
        return -1;
    }
 
    int     i(0), index(-1);
    double  maxError(0.); 
 
    std::vector<int>     active( getActiveSet() );
 
    for( auto error : m_error ){
 
        if(!m_activeNodes[i] ){
 
            if( error > maxError ) {
                maxError = error;
                index = i;
            }
        }
 
        ++i;
    }
 
    return index;
}
 
double RBFKernel::evalError( )
{
    if(m_mode == RBFMode::PARAM) {
        return -1.0;
    }
 
    int                     i(0), j(0);
    //int                     index;
    double                  maxError(0), relError, realValue;
    std::vector<double>     reconValues;
 
    for( int iX=0; iX< m_nodes; ++iX ) {
        reconValues = evalRBF( iX );
 
        j=0;
        relError = 0.;
        for( auto &val : reconValues ) {
            realValue = m_value[j][i];
            relError += std::pow( (val - realValue), 2  );
 
            ++j;
        }
 
        relError = sqrt(relError);
        m_error[i] = relError;
 
        if( relError > maxError ) {
            maxError = relError;
        }
 
        ++i;
    }
 
    return maxError;
}
 
int RBFKernel::solveLSQ()
{
    if(m_mode == RBFMode::PARAM) {
        return -1;
    }
 
    if (m_polyEnabled) {
        // Initialize polynomial
        initializePolynomial();
 
        // Check which parameter of the polynomial has to be activated
        // in order to avoid undetermined system
        initializePolynomialActiveBasis();
    }
 
    double dist;
 
    std::vector<int> activeSet = getActiveSet();
    int nActive = activeSet.size();
 
    int nPoly   = m_polyEnabled ? getPolynomialWeightsCount() : 0;
    int nS      = nActive + nPoly;
    int nP      = m_nodes;
    int nrhs    = getDataCount();
 
    int     n ,m, lda, ldb, info, rank;
    double rcond = -1.0;
 
    m = nP;
    n = nS;
 
    lda = m + nPoly;
    ldb = m + nPoly;
 
    std::vector<double> A(lda * n);
    std::vector<double> b(ldb * nrhs, 0.0);
    std::vector<double> s(m);
 
    for (int j = 0; j < nrhs; ++j) {
        for (int i = 0; i < m; ++i) {
            int k = j * ldb + i;
            b[k]  = m_value[j][i];
        }
    }
 
    int k = 0;
    for (const auto &i : activeSet) {
        for (int j = 0; j < nP; ++j) {
            dist                 = calcDist(j, i) / getSupportRadius(i); //order by column!
            int row              = k % nP;
            int col              = k / nP;
            A[(col * ldb) + row] = evalBasis(dist);
            k++;
        }
    }
 
    // Filling column terms given by the polynomial block.
    if (m_polyEnabled) {
        // Loop on columns
        for (int i = 0; i < nActive; ++i) {
            std::vector<double> polynomialTerms = evalPolynomialBasis(activeSet[i]);
            int j                               = 0;
            for (const double &val : polynomialTerms) {
                A[lda * i + (nP + j)] = val;
                j++;
            }
        }
    }
 
    // Filling row terms given by the polynomial block.
    if (m_polyEnabled) {
        // Loop on rows
        for (int i = 0; i < nP; ++i) {
            std::vector<double> polynomialTerms = evalPolynomialBasis(i);
            int j                               = 0;
            for (const double &val : polynomialTerms) {
                A[(lda * (nActive + j) + i)] = val;
                j++;
            }
        }
    }
 
    info = LAPACKE_dgelsd(LAPACK_COL_MAJOR, lda, nS, nrhs, A.data(), lda, b.data(), ldb, s.data(), rcond, &rank);
 
    if (info > 0) {
        return 1;
    }
 
    m_weight.resize(nrhs);
 
    for (int j = 0; j < nrhs; ++j) {
        m_weight[j].clear();
        m_weight[j].resize(m_nodes,0);
 
        k=0;
        for( const auto &i : activeSet ) {
            m_weight[j][i] = b[j*ldb+k];
            ++k;
        }
    }
 
    if (m_polyEnabled) {
        for (int j = 0; j < nrhs; ++j) {
            double *polynomialCoefficients = m_polynomial.getCoefficients(j);
            int i                          = 0;
            for (int iactive : m_polyActiveBasis) {
                polynomialCoefficients[iactive] = b[j * ldb + nActive + i];
                i++;
            }
        }
    }
        return 0;
    }
 
void RBFKernel::enablePolynomial(bool enable)
{
    m_polyEnabled = enable;
}
 
int RBFKernel::getPolynomialWeightsCount()
{
    int nTerms = m_polyActiveBasis.size();
    return nTerms;
}
 
RBFKernel::LinearPolynomial::LinearPolynomial()
{
    m_dim    = 0;
    m_fields = 0;
}
 
void RBFKernel::LinearPolynomial::clear()
{
    ReconstructionPolynomial::clear(true);
    m_dim    = 0;
    m_fields = 0;
}
 
 
void RBFKernel::LinearPolynomial::initialize()
{
    std::array<double, 3> m_origin;
    m_origin.fill(0.);
    ReconstructionPolynomial::initialize(1, m_dim, m_origin, m_fields, true);
    double *coeffs = getCoefficients();
    for (auto i = 0; i < getCoefficientCount(); i++) {
        coeffs[i] = 0.;
    }
}
 
void RBFKernel::LinearPolynomial::setDimension(int dim)
{
    m_dim = (dim > 0) ? (dim < 4 ? dim : 3) : 0;
}
 
void RBFKernel::LinearPolynomial::setDataCount(int fields)
{
    m_fields = std::max(0, fields);
}
 
void RBFKernel::LinearPolynomial::evalBasis(const double *x, double *basis)
{
    // Set 0-th degree coefficients
    basis[0] = 1.;
 
    // Set 1-st degree coefficients
    for (int i = 0; i < m_dim; ++i) {
        basis[i + 1] = x[i];
    }
}
 
// RBF derived class implementation
 
RBF::~RBF()
{
}
 
// /*!
//  * Default constructor. RBFBasisFunction::WENDLANDC2 is default. RBFMode is
//  * INTERP, by default. Use setMode for changing it.
//  */
// RBF::RBF(): RBFKernel()
// {
//     m_node.clear();
// }
 
RBF::RBF( RBFBasisFunction bfunc)
{
    m_node.clear();
    setFunction( bfunc );
}
 
RBF::RBF(const RBF & other)
    : RBFKernel(other),
      m_node(other.m_node)
{
}
 
RBF & RBF::operator=(RBF other)
{
    swap(other);
 
    return *this;
}
 
void RBF::swap(RBF &other) noexcept
{
    RBFKernel::swap(other);
 
    std::swap(m_node, other.m_node);
}
 
int RBF::getTotalNodesCount(  )
{
    return m_nodes;
}
 
int RBF::addNode( const std::array<double,3> &node )
{
    m_node.push_back(node);
    m_activeNodes.push_back(true);
    m_nodes++;
 
    return m_nodes;
}
 
std::vector<int> RBF::addNode( const std::vector<std::array<double,3>> &node )
{
    int                 i( m_nodes );
    std::vector<int>    ids;
 
    ids.resize( node.size() );
 
    for( auto & id:ids ) {
        id = i;
        ++i;
    }
 
    m_node.insert( m_node.end(), node.begin(), node.end() );
    m_nodes += node.size();
 
    m_activeNodes.resize( m_nodes, true );
 
    return ids;
}
 
bool RBF::removeNode(int id)
{
    if(id < 0 || id >=m_nodes) {
        return false;
    }
 
    m_nodes--;
    m_node.erase(m_node.begin()+id);
    m_activeNodes.erase(m_activeNodes.begin()+id);
 
    return true;
}
 
bool RBF::removeNode(std::vector<int> & list)
{
    // List should be processed in reversed order because nodes
    // are removed using their position in the storages.
    std::sort(list.begin(), list.end(), std::greater<int>());
 
    int nInitialNodes = getTotalNodesCount();
    for(int id : list) {
        removeNode(id);
    }
    int nFinalNodes = getTotalNodesCount();
 
    return ((nInitialNodes - nFinalNodes) == static_cast<int>(list.size()));
}
 
void RBF::removeAllNodes()
{
    m_nodes = 0;
    m_node.clear();
    m_activeNodes.clear();
}
 
double RBF::calcDist(int i, int j)
{
    return norm2(m_node[i]-m_node[j]);
}
 
double RBF::calcDist(const std::array<double,3>& point, int j)
{
    return norm2(point - m_node[j]);
}
 
void RBF::initializePolynomialActiveBasis()
{
    std::vector<int> activeSet(RBFKernel::getActiveSet());
 
    // Initialize active terms for an only constant linear polynomial
    m_polyActiveBasis.clear();
    m_polyActiveBasis.push_back(0);
 
    // Initialize reference coordinates with the first node
    std::array<double, 3> coord(m_node[0]);
 
    // Check if at least one node has one of the coordinates
    // different from the reference ones and enable the related
    // polynomial term
    for (int d = 0; d < 2; ++d) {
        for (const auto &i : activeSet) {
            const std::array<double, 3> &point = m_node[i];
            if (!utils::DoubleFloatingEqual()(coord[d], point[d])) {
                utils::addToOrderedVector<int>(d + 1, m_polyActiveBasis);
                break;
            }
        }
    }
}
 
void RBF::initializePolynomial()
{
    m_polynomial.clear();
    m_polynomial.setDimension(3);
    m_polynomial.setDataCount(getDataCount());
    m_polynomial.initialize();
}
 
std::vector<double> RBF::evalPolynomialBasis(int i)
{
    int nPoly = m_polyActiveBasis.size();
    std::vector<double> result(nPoly);
 
    if (nPoly < 1)
        return result;
 
    const std::array<double, 3> &point = m_node[i];
 
    std::vector<double> completeBasis(m_polynomial.getCoefficientCount());
    m_polynomial.evalBasis(point.data(), completeBasis.data());
 
    int k = 0;
    for (int j : m_polyActiveBasis) {
        result[k] = completeBasis[j];
        k++;
    }
 
    return result;
}
 
std::vector<double> RBF::evalPolynomialBasis(const std::array<double,3> &point)
{
    int nPoly = m_polyActiveBasis.size();
    std::vector<double> result(nPoly);
 
    if (nPoly < 1)
        return result;
 
    std::vector<double> completeBasis(m_polynomial.getCoefficientCount());
    m_polynomial.evalBasis(point.data(), completeBasis.data());
 
    int k = 0;
    for (int j : m_polyActiveBasis) {
        result[k] = completeBasis[j];
        k++;
    }
 
    return result;
}
 
// RBF NAMESPACE UTILITIES
 
double rbf::wendlandc2( double dist )
{
    if( dist > 1) {
        return 0.;
    } else{
        return std::pow(1.-dist,4)*(4.*dist+1.);
    }
}
 
double rbf::linear( double dist )
{
    if( dist > 1) {
        return 0.;
    } else{
        return (1-dist);
    }
}
 
double rbf::gauss90( double dist )
{
    double eps = std::pow(-1.0*std::log(0.1),0.5);
 
    return std::exp(-1.0*std::pow(dist*eps,2));
}
 
double rbf::gauss95( double dist )
{
    double eps = std::pow(-1.0*std::log(0.05),0.5);
 
    return std::exp(-1.0*std::pow(dist*eps,2));
}
 
double rbf::gauss99( double dist )
{
    double eps = std::pow(-1.0*std::log(0.01),0.5);
 
    return std::exp(-1.0*std::pow(dist*eps,2));
}
 
double rbf::c1c0( double dist )
{
    if( dist > 1) {
        return 0.;
    } else{
        return (1.0-std::pow(dist,2));
    }
}
 
double rbf::c2c0( double dist )
{
    if( dist > 1) {
        return 0.;
    } else{
        return (1.0-std::pow(dist,3));
    }
}
 
double rbf::c0c1( double dist )
{
    if( dist > 1) {
        return 0.;
    } else{
        return (1.0- 2.0*dist + std::pow(dist,2));
    }
}
 
double rbf::c1c1( double dist )
{
    if( dist > 1) {
        return 0.;
    } else{
        return (1.0-3.0*std::pow(dist,2)+2.0*std::pow(dist,3));
    }
}
 
double rbf::c2c1( double dist )
{
    if( dist > 1) {
        return 0.;
    } else{
        return (1.0- 4.0*std::pow(dist,3) + 3.0*std::pow(dist,4));
    }
}
 
double rbf::c0c2( double dist )
{
    if( dist > 1) {
        return 0.;
    } else{
        return (1.0 -3.0*dist +3.0*std::pow(dist,2) - std::pow(dist,3));
    }
}
 
double rbf::c1c2( double dist )
{
    if( dist > 1) {
        return 0.;
    } else{
        return (1.0 -6.0*std::pow(dist,2) + 8.0*std::pow(dist,3) - 3.0*std::pow(dist,4));
    }
}
 
double rbf::c2c2( double dist )
{
    if( dist > 1) {
        return 0.;
    } else{
        return (1.0 -10.0*std::pow(dist,3) +15.0*std::pow(dist,4) -6.0*std::pow(dist,5));
    }
}
 
    double rbf::cosinus( double dist )
    {
        if( dist > 1) {
            return 0.;
        } else{
            return 0.5 + 0.5 * std::cos(dist*BITPIT_PI);
        }
    }
 
double rbf::thinplate(double dist)
{
    if (utils::DoubleFloatingGreater()(dist, 0.)) {
        return dist * dist * std::log(dist);
    } else {
        return 0.;
    }
}
 
} // namespace bitpit