#ifndef ARRAY_CONVOLUTION_H
#define ARRAY_CONVOLUTION_H

// For compatibility with boost <1.31 when apply_if replaced with eval_if
#include <boost/version.hpp>
#if BOOST_VERSION<=103100
#include <boost/mpl/apply_if.hpp>
#ifndef eval_if
#define eval_if apply_if
#endif
#else
#include <boost/mpl/eval_if.hpp>
#endif

#include <boost/array.hpp>
#include <boost/multi_array.hpp>
#include <map>
#include <vector>
#include "linalg.h"

namespace mimas {

/** @addtogroup arrayOp
    @{ */
/** Correlation of two multi-arrays.
    This method provides correlation of two n-dimensional arrays. The
    resulting array will have the same size as the input-array \c x.
    Elements outside of the array-boundaries are assumed to be zero.

    The algorithm is intented to be used for convoluting an array with
    a small filter. If the filter is very big, it may be more efficient, to
    perform the correlation in fourier-space.

    For a small array \c y this algorithm is <b>very fast</b>, because before
    doing the correlation the <b>loops are reordered for maximum
    performance</b>.

    @param x First array (the big one).
    @param y Second array (the small one).
    @return Result of correlation.
    @see fourierTransforms
    @author Jan Wedekind (jan@wedesoft.de)
    @date Fri Jul 15 21:37:32 UTC 2005 */
template< class _Array1, class _Array2 >
boost::multi_array< typename _Array1::element, _Array1::dimensionality >
  correlate( const _Array1 &x, const _Array2 &y );

/** Perform a two-dimensional separable correlation.
    If you have compiled Mimas with LAPACK-wrappers, you can also
    let Mimas separate the filter for you (if it is separable).

    @param x 2D array to be filtered.
    @param f Two 2D arrays (one with a single column and one with a single
    row)
    @see correlate_separable */
template< typename T >
boost::multi_array< T, 2 > correlate_separable
  ( const boost::const_multi_array_ref< T, 2 > &x,
    const std::vector< boost::multi_array< T, 2 > > &f )
{
  assert( f.size() == 2 );
  return correlate( correlate( x, f[0] ), f[1] );
}

/** Separate a 2D filter.
    Separate the 2D filter using singular value decomposition. If the filter
    is not separable, you will only get an approximate solution.

    The 2D filter is copied to a matrix. In an ideal case only the first
    singular value \f$\sigma_1\f$ is unequal to zero and the
    corresponding vectors \f$\vec{u_1}\f$ and \f$\vec{v_1}\f$ are holding
    the elements of the separable filter:
    \f[A=\vec{u_1}\sigma_1\vec{v_1}^\top\f] */
template< typename T >
std::vector< boost::multi_array< T, 2 > > separate
   ( const boost::const_multi_array_ref< T, 2 > &array );

/** Fast 2D correlation with a separable 2D filter.
    @see separate */
template< typename T >
boost::multi_array< T, 2 > correlate_separable
  ( const boost::const_multi_array_ref< T, 2 > &x,
    const boost::const_multi_array_ref< T, 2 > &f )
{
  boost::multi_array< T, 2 > result( boost::extents[ x.shape()[0] ]
                                                   [ x.shape()[1] ] );
  std::vector< boost::multi_array< T, 2 > > components( separate( f ) );
  return correlate_separable< T >( x, components );
}

///@}

};

#include "multi_array_conv.tcc"

#endif