vector_ops.cpp

//***************************************************************************//
//* libtrivia: various little routines for C++                              *//
//* vector_ops: vector operations                                           *//
//* (C) Joachim Wuttke 2001, 2012-                                          *//
//* http://apps.jcns.fz-juelich.de                                          *//
//***************************************************************************//

#include <cmath>
#include <string>
#include <vector>
#include <algorithm>
#include <numeric>

using namespace std;

#include "vector_ops.hpp"


//***************************************************************************//
//*  vector operations                                                      *//
//***************************************************************************//

//! Insert val into vector *V, presuming the latter is sorted.

void triv::merge( vector<double> *V, double val )
{
    for (size_t i=0; i<V->size(); ++i) {
        if (val>(*V)[i]) continue;
        if (val==(*V)[i]) return;
        V->insert(V->begin()+i, val);
        return;
    }
    V->push_back(val);
}

//! Eliminate double entries from vector V.

void triv::unique( vector<double> *V, double tolabs, double tolrel )
{
    // tolrel z.Zt. nicht benutzt
    vector<double> aux;
    double diff;
    if (!V->size()) return;
    aux.push_back((*V)[0]);
    for (size_t ii=1; ii<V->size(); ++ii)
        if((diff=fabs((*V)[ii]-aux.back()))>tolabs)
            aux.push_back((*V)[ii]);
    *V = aux;
}

//! True unless V[i]<=V[i-1] for some i.

bool triv::is_ascending( const vector<double>& V )
{
    for( size_t i=2; i<V.size(); ++i )
        if ( V[i] <= V[i-1] )
            return false;
    return true;
}

//! Detect equidistant grid, set *step.

bool triv::is_equidist( double *step, const vector<double>& V )
{
    size_t n = V.size();
    if( n < 2 )
        return false;  // vector too short to determine equidistance
    *step = ( V[n-1] - V[0] ) / (n-1);
    for ( size_t i=1; i<n-1; ++i )
        if( fabs( V[i] - (V[0]+i*(*step)) ) > 1e-12*fabs(*step)+1e-200 )
            return false;
    return true;
}

//! Compute array of sorted indices.

vector<size_t> triv::sorted_indices(vector<double> const& V)
{
    vector<size_t> I(V.size());
    // C++11, not yet supported by gcc-4.4:
    // iota(begin(I), end(I), 0);
    // workaround:
    for( size_t i=0; i<I.size(); ++i )
	I[i] = i;
    // end workaround.
    sort( begin(I), end(I), [&](size_t a, size_t b) { return V[a] < V[b]; } );
    return I;
}