Contributing New Functions to Stan / Expose Function Signature to Stan Models

[Arseniy Tsipenyuk]

Hi everybody,

I am trying to implement my own function in C++ and expose it to STAN just as it is described here:
https://github.com/stan-dev/stan/wiki/Contributing-New-Functions-to-Stan

(The function is intended for private use only, I do not intend to test or submit it to stan git project.)
To start with something simple, I just want to copy the function 'abs' (only for the unary, not vector-valued arguments)
and name it 'my_func'. I did the following:

1. I copied stan/src/stan/math/functions/abs.hpp to stan/src/stan/math/functions/my_func.hpp
2. I made some mutatis mutandis changes in my_func.hpp (basically, I replaced 'abs' by 'my_func' everywhere - code
is included below).
3. I added "#include <stan/math/functions/my_func.hpp>" after line 64 (i.e., inserted this line below line 64) to stan/src/stan/math/functions.hpp
4. I added  "add_unary("my_func");" after line 586 to stan/src/stan/gm/function_signatures.h

Finally, I have written the following "model.stan" file:
"""
parameters {
  real y;
}
model {
  print(my_func(2.));
}
"""

Unfortunately, I still get the message
"""
SYNTAX ERROR, MESSAGE(S) FROM PARSER:
no matches for function name="my_func"
    arg 0 type=real
available function signatures for my_func :
no matches for function name="my_func"
    arg 0 type=real
available function signatures for resonances:
expression is ill formed

ERROR at line 5
  3:    }
  4:    model {
  5:      print(my_func(2.));
                              ^
  6:    }
"""
when I try to build the model. I have a feeling that I need to re-link the STAN library, but how do I do that?

N.B. I am still using STAN 2.5.0, I do not suppose that there is a difference in that issue, though.

Best regards,
Arseniy.


-----
my_func.hpp :
"""
#ifndef STAN__MATH__FUNCTIONS__MY_FUNC_HPP
#define STAN__MATH__FUNCTIONS__MY_FUNC_HPP

#include <cmath>

namespace stan {

  namespace math {

    /**
     * Return floating-point absolute value.
     *
     * Delegates to <code>fabs(double)</code> rather than
     * <code>std::abs(int)</code>.
     *
     * @param x scalar
     * @return absolute value of scalar
     */
    double my_func(double x) {
      return std::fabs(x);
    }

  }
}

#endif
"""

[Daniel Lee]

Hi Arseniy,

You're almost there. The stanc compiler needs to be rebuilt with your changes in it.

What interface are you using? In CmdStan, type 'make clean-all' then try to build your Stan program again.

Daniel

--
You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Arseniy Tsipenyuk]

Hi Daniel,

thanks a lot, that's exactly what I needed!

Best regards,
Arseniy.

[Arseniy Tsipenyuk]

Hi Daniel,

I have some more troubles with user-defined functions - maybe you could help me out? I am trying to
write a function A_1 (that I will use in STAN), and this function calls another function c_one that does
some computation with Eigen matrices. I am having the following issue:

 - the function A_1 takes (and returns) some vectors with base type stan::agrad::var;
 - the function c_one takes a double and returns an instance of Eigen::Matrix<double,2,1>.

I am having conversion troubles in between - see the example below.

Here is my MWE:
File containing the function "c_one" (the file is called "c_func.hpp"):
"""
#include <stan/math/matrix/Eigen.hpp>
#include <vector>
namespace c_func
{
    using Eigen::Matrix;
    using std::vector;

    /**
     * vector c_one(scalar)
     *
     * Complex number 1.
     * Takes a (dummy) double and returns the vector (1., 0.).
     */
    inline
    Eigen::Matrix<double,2,1>
        c_one(const double& y) {
        Eigen::Matrix<double,2,1>  res;
        res(0,0) = 1.0;
        res(1,0) = 0.0;
        return res;
    }
    }
}
"""


File containing the function "A_1" (called 'model.hpp'):
"""
#include <stan/math/matrix/Eigen.hpp>
#include <boost/math/tools/promotion.hpp>
#include <vector>
#include "c_func.hpp"

namespace stan {
  namespace math {
   
    using Eigen::Matrix;
    using std::vector;

    /**
     * vector A_1(vector)
     *
     * Takes a vector y as an argument, and returns c_one(y(0,0)),
     * i.e. the result is c_one evaluated at the first element of y.
     *
     * @tparam T0__ - Scalar type of the input vector
     */
    template <typename T0__>
    inline
    Eigen::Matrix<typename boost::math::tools::promote_args<T0__>::type, Eigen::Dynamic,1>
    A_1(const Eigen::Matrix<T0__, Eigen::Dynamic,1>& y) {

        typedef typename boost::math::tools::promote_args<T0__>::type fun_scalar_t__;
        typedef fun_scalar_t__ fun_return_scalar_t__;

        Eigen::Matrix<typename boost::math::tools::promote_args<T0__>::type, Eigen::Dynamic,1> res;
        res = stan::math::promote_scalar<fun_return_scalar_t__>(c_func::c_one(y(0,0).val()));

        return res;
    }
  }
}
"""

Here is the relevant part of the error message:
"""
stan/src/../../meson_deca/lib/c_lib/model.hpp:29:85: error: member reference
      base type 'const double' is not a structure or union
  ...y(0,0).val()));
     ~~~~~~^~~~
"""
However, if I remove '.val()', I get another error, since I am calling 'c_one' with argument
of the wrong type. I have poked around here and there, but I still can not force it to work.
Any ideas?

Thanks in advance,
Regards,
Arseniy.

[Arseniy Tsipenyuk]

N.B. Had some weird thing with copy-pasting - the file "c_func.hpp" actually looks like that
(no extra '}' in the second-to-last row):

[Daniel Lee]

HI Arseniy,

Your c_one function needs to be templated to allow for an arbitrary scalar type. So:

1. prefix your function with 'template <typename T>'

2. replace 'double' with 'T'

3. Replace y(0,0).val() with y(0,0)

That should get you further. I think you'll still have problems with the Eigen rows and columns, though.

Daniel

[Arseniy Tsipenyuk]

Hi Daniel,

thanks for the reply. Here is what I get:

#### If I use ####
"""
template <typename T>
    inline
    Eigen::Matrix<T,2,1>
        c_one(const T& y) {
        Eigen::Matrix<T,2,1>  res;

        res(0,0) = 1.0;
        res(1,0) = 0.0;
        return res;
    }
"""

I get
"""
In file included from <built-in>:167:
In file included from <command line>:5:
In file included from model.cpp:3:
In file included from stan/src/stan/model/model_header.hpp:16:
In file included from stan/src/stan/agrad/rev.hpp:7:
In file included from stan/src/stan/agrad/rev/var.hpp:5:
In file included from stan/src/stan/math/matrix/Eigen.hpp:4:
In file included from stan/lib/eigen_3.2.0/Eigen/Dense:1:
In file included from stan/lib/eigen_3.2.0/Eigen/Core:284:
stan/lib/eigen_3.2.0/Eigen/src/Core/Assign.h:493:32: error: no member named
      'YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY'
      in 'Eigen::internal::static_assertion<false>'
  EIGEN_STATIC_ASSERT(SameType,YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NE...

stan/lib/eigen_3.2.0/Eigen/src/Core/util/StaticAssert.h:112:65: note: expanded
      from macro 'EIGEN_STATIC_ASSERT'
        if (Eigen::internal::static_assertion<bool(CONDITION)>::MSG) {}
                                                                ^
stan/lib/eigen_3.2.0/Eigen/src/Core/PlainObjectBase.h:411:20: note: in
      instantiation of function template specialization
      'Eigen::DenseBase<Eigen::Matrix<stan::agrad::var, -1, 1, 0, -1, 1>
      >::lazyAssign<Eigen::Matrix<double, 2, 1, 0, 2, 1> >' requested here
      return Base::lazyAssign(other.derived());
                   ^
stan/lib/eigen_3.2.0/Eigen/src/Core/Assign.h:520:97: note: in instantiation of
      function template specialization
      'Eigen::PlainObjectBase<Eigen::Matrix<stan::agrad::var, -1, 1, 0, -1, 1>
      >::lazyAssign<Eigen::Matrix<double, 2, 1, 0, 2, 1> >' requested here
  ...run(Derived& dst, const OtherDerived& other) { return dst.lazyAssign(oth...
                                                               ^
stan/lib/eigen_3.2.0/Eigen/src/Core/PlainObjectBase.h:618:69: note: in
      instantiation of member function
      'Eigen::internal::assign_selector<Eigen::Matrix<stan::agrad::var, -1, 1,
      0, -1, 1>, Eigen::Matrix<double, 2, 1, 0, 2, 1>, false, false>::run'
      requested here
      return internal::assign_selector<Derived,OtherDerived,false>::run(...
                                                                    ^
stan/lib/eigen_3.2.0/Eigen/src/Core/PlainObjectBase.h:603:102: note: in
      instantiation of function template specialization
      'Eigen::PlainObjectBase<Eigen::Matrix<stan::agrad::var, -1, 1, 0, -1, 1>
      >::_set_noalias<Eigen::Matrix<double, 2, 1, 0, 2, 1> >' requested here
  ...OtherDerived& other, const internal::false_type&) { _set_noalias(other); }
                                                         ^
stan/lib/eigen_3.2.0/Eigen/src/Core/PlainObjectBase.h:595:7: note: in
      instantiation of function template specialization
      'Eigen::PlainObjectBase<Eigen::Matrix<stan::agrad::var, -1, 1, 0, -1, 1>
      >::_set_selector<Eigen::Matrix<double, 2, 1, 0, 2, 1> >' requested here
      _set_selector(other.derived(), typename internal::conditional<stat...
      ^
stan/lib/eigen_3.2.0/Eigen/src/Core/Matrix.h:172:20: note: (skipping 2 contexts
      in backtrace; use -ftemplate-backtrace-limit=0 to see all)
      return Base::_set(other);
                   ^
model.cpp:154:43: note: in instantiation of function template specialization
      'stan::math::A_1<stan::agrad::var>' requested here
                stan_print(pstream__,dims(A_1(v)));
                                          ^
stan/src/stan/model/util.hpp:99:21: note: in instantiation of function template
      specialization 'model_model_namespace::model_model::log_prob<true, true,
      stan::agrad::var>' requested here
          .template log_prob<propto,
                    ^
stan/src/stan/model/util.hpp:277:24: note: in instantiation of function template
      specialization 'stan::model::log_prob_grad<true, true,
      model_model_namespace::model_model>' requested here
        = stan::model::log_prob_grad<propto,
                       ^
stan/src/stan/common/command.hpp:208:28: note: in instantiation of function
      template specialization 'stan::model::test_gradients<true, true,
      model_model_namespace::model_model>' requested here
            = stan::model::test_gradients<true,true>(model,cont_vector, ...
                           ^
src/cmdstan/main.cpp:7:26: note: in instantiation of function template
      specialization 'stan::common::command<model_model_namespace::model_model>'
      requested here
    return stan::common::command<stan_model>(argc,argv);
                         ^
1 error generated.
make: *** [model] Error 1

"""

####

I think, the problem is that I am trying to define an Eigen matrix of type
stan::agrad::var, and I do not suppose that Eigen is fine with it.

I can try to use another approach - technically, I can implement all the functions
I need in STAN language, but I still need those functions in some sort of callable form.
I could write a script that extracts the necessary functions from the 'model.cpp' file
produced by STAN, but then I do not know what to do with the pstream__ variable...
Or would it be simpler to write a script that makes C-functions from STAN?

[Arseniy Tsipenyuk]

Oh wait, it seems to work! I will check everything carefully and then post the working code
for the sake of thread completeness.

Best regards,
Arseniy.

[Bob Carpenter]

Boost doesn't allow mixing types. We've had to
work around that extensively in Stan. All of
the Stan functions accept mixed type inputs, but nothing
in straight-up Boost does.

- Bob