[SciPy-dev] Which matrix library in C++ for scipy

[Neal Becker]

Matthieu Brucher wrote:

>>
>> I've got loads of examples.  Good is subjective.  What kind of example
>> would you like?
>>
> 
> For instance a simple example with an array in input and another one in
> ouput :)
> 
> Matthieu

boost::python can expose class as well as functions, but since you are only asking
for a (lowly) function, this is a simple example.

First, generic c++ code for computing a difference of 1 sample:

template<typename scalar_t>
inline std::complex<scalar_t> diff_demod (std::complex<scalar_t> z1, std::complex<scalar_t> z2, scalar_t epsilon) {
  return z2 * std::conj (Limit2<std::complex<scalar_t> > (z1, epsilon));
}

Now code which takes 2 input containers and applies the above, returning an
output container:

template<typename out_t, typename in1_t, typename in2_t, typename scalar_t>
inline out_t diff_demod2 (in1_t const& in1, in2_t const& in2, scalar_t epsilon) {
  if (boost::size (in1) != boost::size (in2))
    throw std::runtime_error ("diff_demod size mismatch");

  out_t out (boost::size (in1));
  
  typename boost::range_const_iterator<in1_t>::type i1 = boost::begin(in1);
  typename boost::range_const_iterator<in2_t>::type i2 = boost::begin(in2);
  typename boost::range_iterator<out_t>::type o = boost::begin(out);
  for (; i1 != boost::end (in1); ++i1, ++i2, ++o)
    *o = diff_demod (*i1, *i2, epsilon);

  return out;
}

Another variant:

template<typename out_t, typename in_t, typename scalar_t>
inline out_t diff_demod1 (in_t const& in, scalar_t epsilon) {
  out_t out (boost::size (in));
  
  typename boost::range_const_iterator<in_t>::type i = boost::begin(in);
  typename boost::range_iterator<out_t>::type o = boost::begin(out);
  typename boost::range_value<in_t>::type prev = 0;
  for (; i != boost::end (in); ++i, ++o) {
    *o = diff_demod (prev, *i, epsilon);
    prev = *i;
  }

  return out;
}

 
Now expose this to python:

BOOST_PYTHON_MODULE(limit)
{

  def ("Limit", &Compute<ublas::vector<Complex>, ublas::vector<Complex>, double >,
       (arg ("in"),
        arg ("epsilon")=1e-6));


  def ("DiffDemod2", &diff_demod2<ublas::vector<Complex>,ublas::vector<Complex>,ublas::vector<Complex>,double>,
       (arg ("in1"), arg ("in2"), arg ("epsilon")=1e-6))
       ;

  def ("DiffDemod1", &diff_demod1<ublas::vector<Complex>,ublas::vector<Complex>,double>,
       (arg ("in"), arg ("epsilon")=1e-6))
       ;
}