common.h

Go to the documentation of this file.

/************************************************************************/
/*                                                                      */
/*    vspline - a set of generic tools for creation and evaluation      */
/*              of uniform b-splines                                    */
/*                                                                      */
/*            Copyright 2015 - 2023 by Kay F. Jahnke                    */
/*                                                                      */
/*    The git repository for this software is at                        */
/*                                                                      */
/*    https://bitbucket.org/kfj/vspline                                 */
/*                                                                      */
/*    Please direct questions, bug reports, and contributions to        */
/*                                                                      */
/*    kfjahnke+vspline@gmail.com                                        */
/*                                                                      */
/*    Permission is hereby granted, free of charge, to any person       */
/*    obtaining a copy of this software and associated documentation    */
/*    files (the "Software"), to deal in the Software without           */
/*    restriction, including without limitation the rights to use,      */
/*    copy, modify, merge, publish, distribute, sublicense, and/or      */
/*    sell copies of the Software, and to permit persons to whom the    */
/*    Software is furnished to do so, subject to the following          */
/*    conditions:                                                       */
/*                                                                      */
/*    The above copyright notice and this permission notice shall be    */
/*    included in all copies or substantial portions of the             */
/*    Software.                                                         */
/*                                                                      */
/*    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND    */
/*    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   */
/*    OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          */
/*    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT       */
/*    HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,      */
/*    WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING      */
/*    FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR     */
/*    OTHER DEALINGS IN THE SOFTWARE.                                   */
/*                                                                      */
/************************************************************************/
 
/*! \file common.h
 
    \brief definitions common to all files in this project, utility code
    
    This file contains
    
    - some common enums and strings
    
    - definition of a few utility types used throughout vspline
    
    - exceptions used throughout vspline
    
    It includes vspline/vector.h which defines vspline's use of
    vectorization (meaning SIMD operation) and associated types ans code. 
*/
 
#ifndef VSPLINE_COMMON
#define VSPLINE_COMMON
 
#include <vigra/multi_array.hxx>
#include <vigra/tinyvector.hxx>
 
namespace vspline
{
/// This enumeration is used for codes connected to boundary conditions.
/// There are two aspects to boundary conditions: During prefiltering,
/// the initial causal and anticausal coefficients have to be calculated
/// in a way specific to the chosen boundary conditions. Bracing also needs
/// these codes to pick the appropriate extrapolation code to extrapolate
/// the coefficients beyond the core array.
 
typedef enum { 
  MIRROR ,    // mirror on the bounds, so that f(-x) == f(x)
  PERIODIC,   // periodic boundary conditions
  REFLECT ,   // reflect, so  that f(-1) == f(0) (mirror between bounds)
  NATURAL,    // natural boundary conditions, f(-x) + f(x) == 2 * f(0)
  CONSTANT ,  // clamp. used for framing, with explicit prefilter scheme
  ZEROPAD ,   // used for boundary condition, bracing
  GUESS ,     // used instead of ZEROPAD to keep margin errors lower
  INVALID
} bc_code;
 
/// bc_name is for diagnostic output of bc codes
 
const std::string bc_name[] =
{
  "MIRROR   " ,
  "PERIODIC ",
  "REFLECT  " ,
  "NATURAL  ",
  "CONSTANT " ,
  "ZEROPAD  " ,
  "GUESS    "
} ;
 
// why 'xlf'? it's for eXtra Large Float'. I use this for the 'most precise
// floating point type which can be used with standard maths' - currently
// this is long double, but in the future this might be quads or even better.
// I use this type where I have very precise constants (like in poles.h)
// which I try and preserve in as much precision as feasible before they
// end up being cast down to some lesser type.
 
typedef long double xlf_type ;
 
/// using definition for the 'elementary type' of a type via vigra's
/// ExpandElementResult mechanism. Since this is a frequently used idiom
/// in vspline but quite a mouthful, here's an abbreviation:
 
template < typename T >
using ET =
typename vigra::ExpandElementResult < T > :: type ;
 
/// produce a std::integral_constant from the size obtained from
/// vigra's ExpandElementResult mechanism
 
template < typename T >
using EN =
typename
std::integral_constant
  < int ,
    vigra::ExpandElementResult < T > :: size
  > :: type ;
 
/// is_element_expandable tests if a type T is known to vigra's
/// ExpandElementResult mechanism. If this is so, the type is
/// considered 'element-expandable'.
  
template < class T >
using is_element_expandable = typename
  std::integral_constant
  < bool ,
    ! std::is_same
      < typename vigra::ExpandElementResult < T > :: type ,
        vigra::UnsuitableTypeForExpandElements
      > :: value
  > :: type ;
 
// for mathematics, we use a few using declarations. First we have
// promote_type, which yields a type an arithmetic operation between
// a T1 and a T2 should yield. We use vigra's PromoteTraits to obtain
// this type. Note how, by starting out with the formation of the
// Promote type, we silently enforce that T1 and T2 have equal channel
// count or are both fundamental types.
 
template < class T1 , class T2 >
using promote_type
      = typename vigra::PromoteTraits < T1 , T2 > :: Promote ;
 
// using vigra's RealPromote on the promote_type gives us a real type
// suitable for real arithmetics involving T1 and T2. We'll use this
// type in template argument lists where T1 and T2 are given already,
// as default value for 'math_type', which is often user-selectable in
// vspline. With this 'sensible' default, which is what one would
// 'usually' pick 'anyway', oftentimes the whole template argument
// list of a specific vspline function can be fixed by ATD, making it
// easy for the users: they merely have to pass, say, an input and output
// array (fixing T1 and T2), and the arithmetics will be done in a suitable
// type derived from these types. This allows users to stray from the
// 'normal' path (by, for example, stating they want the code to perform
// arithemtics in double precision even if T1 and T2 are merely float),
// but makes 'normal' use of the code concise and legible, requiring
// no explicit template arguments.
 
template < class T1 , class T2 >
using common_math_type
      = typename vigra::NumericTraits
                 < promote_type < T1 , T2 > > :: RealPromote ;
 
// while most of the T1, T2 we'll be dealing with will be small aggregates
// of fundametal types - like TinyVectors of float - we also want a
// convenient way to get the fundamental type involved - or the 'ele_type'
// in vspline parlance. Again we use vigra to obtain this value. The
// elementary type can be obtained by using vigra's ExpandElementResult
// mechanism, which, as a traits class, can easily be extended to any
// type holding an aggregate of equally typed fundamentals. Note that we might
// start out with the elementary types of T1 and T2 and take their Promote,
// yielding the same resultant type. If promote_type<T1,T2> is fundamental,
// the ExpandElementResult mechanism passes it through unchanged.
 
template < class T1 , class T2 >
using promote_ele_type
      = typename vigra::ExpandElementResult
                 < promote_type < T1 , T2 > > :: type ;
 
// when doing arithmetic in vspline, typically we use real fundametal types
// which 'suit' the data types we want to put into the calculations. Again we
// use vigra to determine an apprpriate fundamental type by obtaining the
// RealPromote of the promote_ele_type. 
 
template < class T1 , class T2 >
using common_math_ele_type
      = typename vigra::NumericTraits
                 < promote_ele_type < T1 , T2 > > :: RealPromote ;
 
/// produce the 'canonical type' for a given type T. If T is
/// single-channel, this will be the elementary type itself,
/// otherwise it's a TinyVector of the elementary type.
/// optionally takes the number of elements the resulting
/// type should have, to allow construction from a fundamental
/// and a number of channels.
  
// vspline 'respects the singular'. this is to mean that throughout
// vspline I avoid using fixed-size containers holding a single value
// and use the single value itself. While this requires some type
// artistry in places, it makes using the code more natural. Only
// when it comes to providing generic code to handle data which may
// or may not be aggregates, vspline moves to using 'synthetic' types,
// which are always TinyVectors, possibly with only one element.
// user code can pass data as canonical or synthetic types. All
// higher-level operations will produce the same type as output
// as they receive as input. The move to 'synthetic' types is only
// done internally to ease the writing of generic code.
// both the 'canonical' and the 'synthetic' type have lost all
// 'special' meaning of their components (like their meaning as
// the real and imaginary part of a complex number). vspline will
// process all types in it's scope as such 'neutral' aggregates of
// several identically-typed elementary values.
 
template < typename T , int N = EN < T > :: value >
using canonical_type =
typename
  std::conditional
  < N == 1 ,
    ET < T > ,
    vigra::TinyVector < ET < T > , N >
  > :: type ;
 
template < typename T , int N = EN < T > :: value >
using synthetic_type =
vigra::TinyVector < ET < T > , N > ;
 
template < class T , size_t sz = 1 >
struct invalid_scalar
{
  static_assert ( sz == 1 , "scalar values can only take size 1" ) ;
} ;
 
/// definition of a scalar with the same template argument list as
/// a simdized type, to use 'scalar' in the same syntactic slot
 
template < class T , size_t sz = 1 >
using scalar =
  typename std::conditional
    < sz == 1 ,
      T ,
      invalid_scalar < T , sz >
    > :: type ;
 
// KFJ 2019-02-21 I want to switch to using size_t for all size-like
// arguments. vigra::TinyVector takes int as it's size argument, so I
// introduce 'tiny_vector' which takes size_t instead and static-casts
// the size_t template argument to int. As long as there are no sizes
// beyond an int's maximal size this is safe.
    
template < typename T , size_t SZ >
using tiny_vector = typename
vigra::TinyVector < T , static_cast < int > ( SZ ) > ;
 
/// vspline creates vigra::MultiArrays of vectorized types. As long as
/// the vectorized types are Vc::SimdArray or vspline::simd_type, using
/// std::allocator is fine, but when using other types, using a specific
/// allocator may be necessary. Currently this is never the case, but I
/// have the lookup of allocator type from this traits class in place if
/// it should become necessary.
 
template < typename T >
struct allocator_traits
{
  typedef std::allocator < T > type ;
} ;
 
// TODO My use of exceptions is a bit sketchy...
 
/// for interfaces which need specific implementations we use:
 
struct not_implemented
: std::invalid_argument
{
  not_implemented ( const char * msg )
  : std::invalid_argument ( msg ) { }  ;
} ;
 
/// dimension-mismatch is thrown if two arrays have different dimensions
/// which should have the same dimensions.
 
struct dimension_mismatch
: std::invalid_argument
{
  dimension_mismatch ( const char * msg )
  : std::invalid_argument ( msg ) { }  ;
} ;
 
/// shape mismatch is the exception which is thrown if the shapes of
/// an input array and an output array do not match.
 
struct shape_mismatch
: std::invalid_argument
{
  shape_mismatch  ( const char * msg )
  : std::invalid_argument ( msg ) { }  ;
} ;
 
/// exception which is thrown if an opertion is requested which vspline
/// does not support
 
struct not_supported
: std::invalid_argument
{
  not_supported  ( const char * msg )
  : std::invalid_argument ( msg ) { }  ;
} ;
 
/// out_of_bounds is thrown by mapping mode REJECT for out-of-bounds coordinates
/// this exception is left without a message, it only has a very specific application,
/// and there it may be thrown often, so we don't want anything slowing it down.
 
struct out_of_bounds
{
} ;
 
/// exception which is thrown when assigning an rvalue which is larger than
/// what the lvalue can hold
 
struct numeric_overflow
: std::invalid_argument
{
  numeric_overflow  ( const char * msg )
  : std::invalid_argument ( msg ) { }  ;
} ;
 
} ; // end of namespace vspline
 
// with these common definitions done, we now include 'vector.h', which
// defines vectorization methods used throughout vspline.
 
#include "vector.h"
 
#ifdef USE_FMA
 
// 2019-02-13 tentative use of fma in vspline
// we need definitions of fma for fundamentals, vigra::TinyVectors and
// vspline::simd_type. The least spcific template delegates to std::fma,
// then we have two more specific templates for the container types.
 
template < typename arg_t >
arg_t fma ( const arg_t & arg1 , const arg_t & arg2 , const arg_t & arg3 )
{
  return std::fma ( arg1 , arg2 , arg3 ) ;
}
 
template < typename T , int SZ >
vigra::TinyVector < T , SZ > fma ( const vigra::TinyVector < T , SZ > & arg1 ,
                                   const vigra::TinyVector < T , SZ > & arg2 ,
                                   const vigra::TinyVector < T , SZ > & arg3 )
{
  vigra::TinyVector < T , SZ > result ;
  for ( size_t i = 0 ; i < SZ ; i++ )
    result[i] = fma ( arg1[i] , arg2[i] , arg3[i] ) ;
  return result ;
}
 
template < typename T , size_t SZ >
vspline::simd_type < T , SZ > fma ( const vspline::simd_type < T , SZ > & arg1 ,
                                  const vspline::simd_type < T , SZ > & arg2 ,
                                  const vspline::simd_type < T , SZ > & arg3 )
{
  vspline::simd_type < T , SZ > result ;
  for ( size_t i = 0 ; i < SZ ; i++ )
    result[i] = fma ( arg1[i] , arg2[i] , arg3[i] ) ;
  return result ;
}
 
#endif // USE_FMA
 
#endif // VSPLINE_COMMON