Checking whether lower/upper bounds are finite breaks backward compatibility

[andrjohns]

Description

The lb_constrain and ub_constrain functions were recently updated to check that the provided limits are finite, however this is breaking some stan_surv models in the rstanarm package.

The stan_surv function uses a conditional to determine whether some parameters should be lower-bounded at zero or at negative_infinity() (i.e., unrestricted):

  /**
  * Return the lower bound for the baseline hazard parameters
  *
  * @param type An integer indicating the type of baseline hazard
  * @return A real
  */
  real coefs_lb(int type) {
    real lb;
    if (type == 2) // B-splines, on log haz scale
      lb = negative_infinity();
    else if (type == 3) // piecewise constant, on log haz scale
      lb = negative_infinity();
    else
      lb = 0;
    return lb;
  }

The recent updates to the constrains is then causing these stan_surv models to fail: https://discourse.mc-stan.org/t/stan-surv-crashes-r-repeatedly/20715/11

Current Version:

v4.0.0

[SteveBronder]

Oh geez sorry about that. Let me see if I can make a simple patch for this otherwise we can revert the bounds pr

[bbbales2]

Oooooooof we changed that because of this: #2291

I think now to do it in Stan you'd need to do the optional parameters thing:

real myparam_unconstrained[type == 1];
real<lower = lb> myparam_constrained[type > 2];

And then code in transformed parameters to select which myparam is the one that is being used.

[rok-cesnovar]

So to make sure this isnt in 2.26.0?

[SteveBronder]

Yeah we might need to revert and represent if we are doing a quickfix release today. Otherwise I should have a fix up today

[bbbales2]

This is the pull: #2285 and it was after 2.26.

[bbbales2]

I think any reverting would be because StanHeaders is pulling from develop or something and we can't get a patch in reasonably quickly?

[bbbales2]

@SteveBronder I saved a copy of lb_constrain when I had it written this way cause I remember it was particularly annoying to write. Dunno if there's a way to bake it into prim since there's the conditional thing.

ub_constrain will be similar to write as this. lub_constrain just has more combinatorics:

#ifndef STAN_MATH_REV_FUN_LB_CONSTRAIN_HPP
#define STAN_MATH_REV_FUN_LB_CONSTRAIN_HPP

#include <stan/math/rev/core.hpp>
#include <stan/math/prim/meta.hpp>
#include <stan/math/prim/fun/lb_constrain.hpp>
#include <cmath>

namespace stan {
namespace math {

/**                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                    
 * Return the lower-bounded value for the specified unconstrained input                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                
 * and specified lower bound.                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                          
 *                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                     
 * <p>The transform applied is                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                         
 *                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                     
 * <p>\f$f(x) = \exp(x) + L\f$                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                         
 *                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                     
 * <p>where \f$L\f$ is the constant lower bound.                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                       
 *                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                     
 * <p>If the lower bound is negative infinity, this function                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                           
 * reduces to <code>identity_constrain(x)</code>.                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                      
 *                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                     
 * @tparam T type of Matrix                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                            
 * @tparam L type of lower bound                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                       
 * @param[in] x Unconstrained Matrix input                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                             
 * @param[in] lb lower bound on constrained output                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                     
 * @return lower bound constrained value corresponding to inputs                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                       
 */
template <typename T, typename L,
          require_all_matrix_t<T, L>* = nullptr,
          require_any_st_var<T, L>* = nullptr>
inline auto lb_constrain(const T& x, const L& lb) {
  using std::exp;
  using ret_type = return_var_matrix_t<T, T, L>;

  if(!is_constant<T>::value && !is_constant<L>::value) {
    arena_t<promote_scalar_t<var, T>> arena_x = x;
    arena_t<promote_scalar_t<var, L>> arena_lb = lb;

    arena_t<ret_type> ret = lb_constrain(arena_x.val(), arena_lb.val());

    reverse_pass_callback([arena_x, arena_lb, ret]() mutable {
      for(size_t j = 0; j < arena_x.cols(); ++j) {
        for(size_t i = 0; i < arena_x.rows(); ++i) {
          double x_val = arena_x.val().coeff(i, j);
          double lb_val = arena_lb.val().coeff(i, j);
          double ret_adj = ret.adj().coeff(i, j);

          if (unlikely(lb_val == NEGATIVE_INFTY)) {
            arena_x.adj().coeffRef(i, j) += ret_adj;
          } else {
            arena_x.adj().coeffRef(i, j) += ret_adj * exp(x_val);
            arena_lb.adj().coeffRef(i, j) += ret_adj;
          }
        }
      }
      });

    return ret_type(ret);
  } else if(!is_constant<T>::value) {
    arena_t<promote_scalar_t<var, T>> arena_x = x;
    arena_t<promote_scalar_t<double, L>> arena_lb = value_of(lb);

    arena_t<ret_type> ret = lb_constrain(arena_x.val(), arena_lb.val());

    reverse_pass_callback([arena_x, arena_lb, ret]() mutable {
      for(size_t j = 0; j < arena_x.cols(); ++j) {
        for(size_t i = 0; i < arena_x.rows(); ++i) {
          double x_val = arena_x.val().coeff(i, j);
          double lb_val = arena_lb.val().coeff(i, j);
          double ret_adj = ret.adj().coeff(i, j);

          if (unlikely(lb_val == NEGATIVE_INFTY)) {
            arena_x.adj().coeffRef(i, j) += ret_adj;
          } else {
            arena_x.adj().coeffRef(i, j) += ret_adj * exp(x_val);
          }
        }
      }
    });

    return ret_type(ret);
  } else {
    arena_t<promote_scalar_t<var, L>> arena_lb = lb;

    arena_t<ret_type> ret = lb_constrain(value_of(x), arena_lb.val());

    reverse_pass_callback([arena_lb, ret]() mutable {
      for(size_t j = 0; j < arena_lb.cols(); ++j) {
        for(size_t i = 0; i < arena_lb.rows(); ++i) {
          double lb_val = arena_lb.val().coeff(i, j);
          double ret_adj = ret.adj().coeff(i, j);

          if (lb_val != NEGATIVE_INFTY) {
            arena_lb.adj().coeffRef(i, j) += ret_adj;
          }
        }
      }
    });

    return ret_type(ret);
  }
}

/**                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                    
 * Return the lower-bounded value for the specified unconstrained input                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                
 * and specified lower bound.                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                          
 *                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                     
 * <p>The transform applied is                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                         
 *                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                     
 * <p>\f$f(x) = \exp(x) + L\f$                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                         
 *                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                     
 * <p>where \f$L\f$ is the constant lower bound.                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                       
 *                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                     
 * <p>If the lower bound is negative infinity, this function                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                           
 * reduces to <code>identity_constrain(x)</code>.                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                      
 *                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                     
 * @tparam T type of Matrix                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                            
 * @tparam L type of lower bound                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                       
 * @param[in] x Unconstrained Matrix input                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                             
 * @param[in] lb lower bound on constrained output                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                     
 * @return lower bound constrained value corresponding to inputs                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                       
 */
template <typename T, typename L,
          require_all_matrix_t<T, L>* = nullptr,
          require_any_st_var<T, L>* = nullptr>
inline auto lb_constrain(const T& x, const L& lb, return_type_t<T, L>& lp) {
  using std::exp;
  using ret_type = return_var_matrix_t<T, T, L>;

  if(!is_constant<T>::value && !is_constant<L>::value) {
    arena_t<promote_scalar_t<var, T>> arena_x = x;
    arena_t<promote_scalar_t<var, L>> arena_lb = lb;

    double lp_val = 0.0;
    arena_t<ret_type> ret = lb_constrain(arena_x.val(), arena_lb.val(), lp_val);
    lp += lp_val;

    reverse_pass_callback([arena_x, arena_lb, ret, lp]() mutable {
      double lp_adj = lp.adj();
      for(size_t j = 0; j < arena_x.cols(); ++j) {
        for(size_t i = 0; i < arena_x.rows(); ++i) {
          double x_val = arena_x.val().coeff(i, j);
          double lb_val = arena_lb.val().coeff(i, j);
          double ret_adj = ret.adj().coeff(i, j);

          if (unlikely(lb_val == NEGATIVE_INFTY)) {
            arena_x.adj().coeffRef(i, j) += ret_adj;
          } else {
            arena_x.adj().coeffRef(i, j) += ret_adj * exp(x_val) + lp_adj;
            arena_lb.adj().coeffRef(i, j) += ret_adj;
          }
        }
      }
      });

    return ret_type(ret);
  } else if(!is_constant<T>::value) {
    arena_t<promote_scalar_t<var, T>> arena_x = x;
    arena_t<promote_scalar_t<double, L>> arena_lb = value_of(lb);

    double lp_val = 0.0;
    arena_t<ret_type> ret = lb_constrain(arena_x.val(), arena_lb.val(), lp_val);
    lp += lp_val;

    reverse_pass_callback([arena_x, arena_lb, ret, lp]() mutable {
      double lp_adj = lp.adj();

      for(size_t j = 0; j < arena_x.cols(); ++j) {
        for(size_t i = 0; i < arena_x.rows(); ++i) {
          double x_val = arena_x.val().coeff(i, j);
          double lb_val = arena_lb.val().coeff(i, j);
          double ret_adj = ret.adj().coeff(i, j);

          if (unlikely(lb_val == NEGATIVE_INFTY)) {
            arena_x.adj().coeffRef(i, j) += ret_adj;
          } else {
            arena_x.adj().coeffRef(i, j) += ret_adj * exp(x_val) + lp_adj;
          }
        }
      }
    });

    return ret_type(ret);
  } else {
    arena_t<promote_scalar_t<var, L>> arena_lb = lb;

    double lp_val = 0.0;
    arena_t<ret_type> ret = lb_constrain(value_of(x), arena_lb.val(), lp_val);
    lp += lp_val;

    reverse_pass_callback([arena_lb, ret]() mutable {
      for(size_t j = 0; j < arena_lb.cols(); ++j) {
        for(size_t i = 0; i < arena_lb.rows(); ++i) {
          double lb_val = arena_lb.val().coeff(i, j);
          double ret_adj = ret.adj().coeff(i, j);

          if (lb_val != NEGATIVE_INFTY) {
            arena_lb.adj().coeffRef(i, j) += ret_adj;
          }
        }
      }
    });

    return ret_type(ret);
  }
}

}  // namespace math                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                   
}  // namespace stan                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                   

#endif

[rok-cesnovar]

I think any reverting would be because StanHeaders is pulling from develop or something and we can't get a patch in reasonably quickly?

the PR was done after 2.26.0 so no need to do this now. The development version of stanheaders needs to not use develop.

[hsbadr]

The development version of stanheaders needs to not use develop.

I've to disagree. I think the development version should track develop branches (for incremental updates and testing, to find and fix bugs like this) while releases track specific branches at specific tags. For example, StanHeaders_2.26 branches in all modules (math v4.0.0, stan v2.26.0, and StanHeaders v2.26) are being patched for backward compatibility.

[rok-cesnovar]

@bbbales2 This is closed right?

[bbbales2]

Hopefully this in math and this in Stan did it.