gsl_vector.h

Go to the documentation of this file.

/*
 * $Id: gsl_vector.h 17393 2014-08-23 20:19:14Z predoehl $
 */

#ifndef GSL_VECTOR_WRAP_H_KJBLIB_UARIZONAVISION
#define GSL_VECTOR_WRAP_H_KJBLIB_UARIZONAVISION

#include <l/l_sys_sys.h>
#include <m_cpp/m_cpp_incl.h>
#include "gsl_cpp/gsl_util.h"

#include <iterator>

#ifdef KJB_HAVE_GSL
#include "gsl/gsl_vector.h" /* no need for extern "C" */
#else

#warning "GNU GSL is absent, so GNU GSL wrapper will not work properly."

/* Suprisingly, it seems we must define a nominal fake vector here, because the
 * compiler apparently does a lot of inlining and it sees through most
 * transparent ruses.  So I'm basically giving it a working vector
 * implementation, although the class Gsl_Vector itself should not be able to
 * be instantiated (it should throw Missing_dependency in all ctors).
 *
 * (Later) I'm not sure this is really necessary, but if not it is no worse
 * than a harmless delusion.
 */
struct gsl_vector {
    double* p; size_t size;
    gsl_vector( size_t n ) : p( new double[ n ] ), size( n ) {}
    ~gsl_vector() { delete[] p; }
    double& at( size_t n ) { return p[n % size]; }
    const double& at( size_t n ) const { return p[n % size]; }
};
#define gsl_vector_alloc(x) (new gsl_vector(x)) /* don't actually use this! */
/*
 * It is sad but true:  because multimin foists upon us raw, unwrapped
 * gsl_vectors in its callback functions, it is simply easiest to use the
 * native access methods.  So it is also easiest to provide fake stand-ins when
 * the library is absent.  That is why the macros below are defined:  they
 * stand in when we compile without the library.
 * Although the code below looks plausible, remember IT WILL NEVER BE EXECUTED.
 * This is just elaborate scenery.
 */
#define gsl_vector_set(w,i,x) ((w)->at(i)=(x))
#define gsl_vector_get(w,i) ((w)->at(i))

#endif

#if ! defined( KJB_HAVE_ISNAN ) || ! defined( KJB_HAVE_FINITE )
#warning "This code tests for NAN and finite values but your system lacks"
#warning "this feature.  Method Gsl_Vector::is_normal will always return true."
#endif

#include <ostream>
#include <cmath>                /* required for isnormal(3) macro */

namespace kjb {

class Gsl_Vector {

    gsl_vector* m_vec;

    /*
     * this private method should be called shortly after each invocation of
     * gsl_vector_alloc() because we would like to maintain the invariant that
     * m_vec is never allowed to remain equal to NULL.
     */
    void throw_if_null()
    {
        ETX_2( 00 == m_vec, "Memory allocation failure for Gsl_Vector" );
    }

    // called by the at() methods
    void check_bounds( size_t index ) const
    {
        /*
         * GCC complains about this test:  "that could never happen!!"
        if ( index < 0 )
        {
            KJB_THROW( Index_out_of_bounds );
        }
        */
        if ( size() <= index ) 
        {
            KJB_THROW( Index_out_of_bounds );
        }
    }

public:

    Gsl_Vector( size_t dims )
    :   m_vec( gsl_vector_alloc( dims ) )
    {
#ifdef KJB_HAVE_GSL
        throw_if_null();
#else
        KJB_THROW_2( Missing_dependency, "GNU GSL" );
#endif
    }


    explicit Gsl_Vector( const kjb::Vector& kv )
    :   m_vec( gsl_vector_alloc( kv.size() ) )
    {
#ifdef KJB_HAVE_GSL
        throw_if_null();
        size_t KVSZ = std::max( kv.size(), 0 );
        for( size_t iii = 0; iii < KVSZ; ++iii )
        {
            at( iii ) = kv.at( iii );
        }
#else
        KJB_THROW_2( Missing_dependency, "GNU GSL" );
#endif
    }


    Gsl_Vector( const Gsl_Vector& src )
    :   m_vec( gsl_vector_alloc( src.size() ) )
    {
#ifdef KJB_HAVE_GSL
        throw_if_null();
        GSL_ETX( gsl_vector_memcpy( m_vec, src.m_vec ) );
#else
        KJB_THROW_2( Missing_dependency, "GNU GSL" );
#endif
    }


    Gsl_Vector( const gsl_vector& vec_to_be_copied )
    :   m_vec( gsl_vector_alloc( vec_to_be_copied.size ) )
    {
#ifdef KJB_HAVE_GSL
        throw_if_null();
        GSL_ETX( gsl_vector_memcpy( m_vec, &vec_to_be_copied ) );
#else
        KJB_THROW_2( Missing_dependency, "GNU GSL" );
#endif
    }


    Gsl_Vector& operator=( const Gsl_Vector& src )
    {
#ifdef KJB_HAVE_GSL
        if ( m_vec != src.m_vec )
        {
            if ( size() == src.size() )
            {
                GSL_ETX( gsl_vector_memcpy( m_vec, src.m_vec ) );
            }
            else
            {
                Gsl_Vector clone( src );    // copy ctor makes a duplicate
                swap( clone );              // then, become that duplicate
            }
        }
#endif
        return *this;
    }

    double* begin()
    {
#ifdef KJB_HAVE_GSL
        return gsl_vector_ptr( m_vec, 0 );
#else
        return 00;
#endif
    }

    const double* begin() const
    {
#ifdef KJB_HAVE_GSL
        return gsl_vector_const_ptr( m_vec, 0 );
#else
        return 00;
#endif
    }

    double* end()
    {
#ifdef KJB_HAVE_GSL
        return gsl_vector_ptr( m_vec, size() );
#else
        return 00;
#endif
    }

    const double* end() const
    {
#ifdef KJB_HAVE_GSL
        return gsl_vector_const_ptr( m_vec, size() );
#else
        return 00;
#endif
    }


    void swap( Gsl_Vector& v )
    {
        using std::swap;
        KJB(ASSERT( m_vec ));
        KJB(ASSERT( v.m_vec ));
        swap( m_vec, v.m_vec );
    }


    ~Gsl_Vector()
    {
        KJB(ASSERT( m_vec ));
#ifdef KJB_HAVE_GSL
        gsl_vector_free( m_vec );
#else
        delete m_vec;
#endif
    }


    template < class Iterator >
    Gsl_Vector( Iterator begin, Iterator end )
    :   m_vec( gsl_vector_alloc( std::distance( begin, end ) ) )
    {
#ifdef KJB_HAVE_GSL
        throw_if_null();
        size_t iii = 0;
        for( Iterator cur = begin; cur != end; ++cur )
        {
            at( iii++ ) = *cur;
        }
#else
        KJB_THROW_2( Missing_dependency, "GNU GSL" );
#endif
    }


    size_t size() const
    {
        KJB(ASSERT( m_vec ));
        return m_vec -> size;
    }

    
    /*
     * Because of the transparent convert-to-pointer operator, we lose the
     * ability to set the vector using square brackets, because square brackets
     * are ambiguous -- are they an implicit pointer dereference, or do they
     * refer to operator[] applied to the object?
     * But you can still use at().
     *
    double operator[]( size_t index ) const
    {
        KJB(ASSERT( m_vec ));
        return gsl_vector_get( m_vec, index );
    }

    double& operator[]( size_t index )
    {
        KJB(ASSERT( m_vec );
        return *gsl_vector_ptr( m_vec, index );
    }
    */


    double at( size_t index ) const
    {
        KJB(ASSERT( m_vec ));
#ifdef KJB_HAVE_GSL
        check_bounds( index );
        return gsl_vector_get( m_vec, index );
#else
        return m_vec -> at( index );
#endif
    }


    double& at( size_t index )
    {
        KJB(ASSERT( m_vec ));
#ifdef KJB_HAVE_GSL
        check_bounds( index );
        return *gsl_vector_ptr( m_vec, index );
#else
        return m_vec -> at( index );
#endif
    }


    double& front()
    {
        KJB(ASSERT( m_vec ));
        return at( 0 );
    }


    double front() const
    {
        KJB(ASSERT( m_vec ));
        return at( 0 );
    }


    double& back()
    {
        KJB(ASSERT( m_vec ));
        return at( size() - 1 );
    }


    double back() const
    {
        KJB(ASSERT( m_vec ));
        return at( size() - 1 );
    }


    Gsl_Vector& set_all( double val )
    {
        KJB(ASSERT( m_vec ));
#ifdef KJB_HAVE_GSL
        gsl_vector_set_all( m_vec, val );
#endif
        return *this;
    }


    Gsl_Vector& set_zero()
    {
        KJB(ASSERT( m_vec ));
#ifdef KJB_HAVE_GSL
        gsl_vector_set_zero( m_vec );
#endif
        return *this;
    }


    Gsl_Vector& set_basis( size_t index )
    {
        KJB(ASSERT( m_vec ));
#ifdef KJB_HAVE_GSL
        GSL_ETX( gsl_vector_set_basis( m_vec, index ) );
#endif
        return *this;
    }

    
    operator gsl_vector*()
    {
        return m_vec;
    }


    kjb::Vector vec() const
    {
        kjb::Vector kv( size() );
        for( size_t iii = 0; iii < size(); ++iii )
        {
            kv.at( iii ) = at( iii );
        }
        return kv;
    }


    Gsl_Vector ew_multiply( const Gsl_Vector& that ) const
    {
        if ( size() != that.size() )
        {
            KJB_THROW( Dimension_mismatch );
        }
        Gsl_Vector product( *this );    // copy this vector
#ifdef KJB_HAVE_GSL
        GSL_ETX( gsl_vector_mul( product.m_vec, that.m_vec ) );
#endif
        return product;
    }


    double sum() const
    {
        double sum_vec = 0;
        for( size_t iii = 0; iii < size(); ++iii )
        {
            sum_vec += at( iii );
        }
        return sum_vec;
    }


    double dot( const Gsl_Vector& that ) const
    {
        Gsl_Vector this_x_that( ew_multiply( that ) );
        return this_x_that.sum();
    }


    double l2_norm() const
    {
        using std::sqrt;
        return sqrt( dot( *this ) );
    }


    bool is_normal() const
    {
        if ( 0 == m_vec )
        {
            return true;
        }
        for( size_t iii = 0; iii < size(); ++iii )
        {
            double x = at( iii );
            // The following uses KJB-lib wraps on the lower level facilities.
            // On systems lacking these capabilities, all values seem normal.
            if ( isnand( x ) || ! finite( x ) )
            {
                return false;
            }
        }
        return true;
    }


    Gsl_Vector& operator+=( const Gsl_Vector& that )
    {
#ifdef KJB_HAVE_GSL
        if ( that.size() != size() )
        {
            KJB_THROW( Dimension_mismatch );
        }
        GSL_ETX( gsl_vector_add( m_vec, that.m_vec ) );
#endif
        return *this;
    }


    Gsl_Vector operator+( const Gsl_Vector& that ) const
    {
        Gsl_Vector vsum( *this );
        return vsum += that;
    }


    Gsl_Vector& operator-=( const Gsl_Vector& that )
    {
#ifdef KJB_HAVE_GSL
        if ( that.size() != size() )
        {
            KJB_THROW( Dimension_mismatch );
        }
        GSL_ETX( gsl_vector_sub( m_vec, that.m_vec ) );
#endif
        return *this;
    }


    Gsl_Vector operator-( const Gsl_Vector& that ) const
    {
        Gsl_Vector vdif( *this );
        return vdif -= that;
    }


    Gsl_Vector& ow_multiply_vector_by_scalar( double scalar )
    {
#ifdef KJB_HAVE_GSL
        GSL_ETX( gsl_vector_scale( m_vec, scalar ) );
#endif
        return *this;
    }


    Gsl_Vector& operator*=( double scalar )
    {
        return ow_multiply_vector_by_scalar( scalar );
    }


    Gsl_Vector& ow_reverse()
    {
#ifdef KJB_HAVE_GSL
        GSL_ETX( gsl_vector_reverse( m_vec ) );
#endif
        return *this;
    }


    Gsl_Vector operator*( double scalar ) const
    {
        Gsl_Vector product( *this );
        return product *= scalar;
    }


    Gsl_Vector& operator+=( double scalar )
    {
#ifdef KJB_HAVE_GSL
        GSL_ETX( gsl_vector_add_constant( m_vec, scalar ) );
#endif
        return *this;
    }


    Gsl_Vector operator+( double scalar ) const
    {
        Gsl_Vector vsum( *this );
        return vsum += scalar;
    }

    Gsl_Vector& operator-=( double scalar )
    {
        return operator+=( -scalar );
    }


    Gsl_Vector operator-( double scalar ) const
    {
        Gsl_Vector vsum( *this );
        return vsum -= scalar;
    }


    Gsl_Vector operator-() const
    {
        return *this*-1.0;
    }


    double max() const
    {
#ifdef KJB_HAVE_GSL
        KJB(ASSERT( m_vec ));
        return gsl_vector_max( m_vec );
#endif
    }


    double min() const
    {
#ifdef KJB_HAVE_GSL
        KJB(ASSERT( m_vec ));
        return gsl_vector_min( m_vec );
#endif
    }


    /*
     * There are more GSL vector operations that are not wrapped:
     * fprintf, fscanf, fread, fwrite, additional math operations.
     */
};


inline
Gsl_Vector operator*( double scalar, const Gsl_Vector& vector )
{
    return vector * scalar;
}


} // namespace kjb


inline
std::ostream& operator<<( std::ostream& os, const kjb::Gsl_Vector& v )
{
    for( size_t iii = 0; iii < v.size(); ++iii )
    {
        os << v.at( iii ) << '\n';
    }
    return os;
}

inline
std::istream& operator>>( std::istream& is, kjb::Gsl_Vector& v )
{
    // To repeat:  v must already know what size it is supposed to be!
    for( size_t iii = 0; iii < v.size(); ++iii )
    {
        is >> v.at( iii );
    }
    return is;
}


namespace std {

    template<>
    inline void swap( kjb::Gsl_Vector& v1, kjb::Gsl_Vector& v2 )
    {
        v1.swap( v2 );
    }
}


#endif /* GSL_VECTOR_WRAP_H_KJBLIB_UARIZONAVISION */