writer.hpp

Go to the documentation of this file.

#ifndef __STAN__IO__WRITER_HPP__
#define __STAN__IO__WRITER_HPP__

#include <cstddef>
#include <stdexcept>
#include <vector>

#include <boost/multi_array.hpp>
#include <boost/throw_exception.hpp>

#include <stan/math/matrix.hpp>
#include <stan/math/special_functions.hpp>

#include <stan/prob/transform.hpp>

namespace stan {

  namespace io {

    template <typename T>
    class writer {
    private:
      std::vector<T> data_r_;
      std::vector<int> data_i_;
    public:

      typedef typename stan::math::EigenType<T>::matrix matrix_t;
      typedef typename stan::math::EigenType<T>::vector vector_t;
      typedef typename stan::math::EigenType<T>::row_vector row_vector_t;
      
      typedef Eigen::Array<T,Eigen::Dynamic,1> array_vec_t;

      const double CONSTRAINT_TOLERANCE;

      writer(std::vector<T>& data_r,
             std::vector<int>& data_i)
        : data_r_(data_r),
          data_i_(data_i),
          CONSTRAINT_TOLERANCE(1E-8) {
      }

      ~writer() { }

      std::vector<T>& data_r() {
        return data_r_;
      }


      std::vector<int>& data_i() {
        return data_i_;
      }

      void integer(int n) {
        data_i_.push_back(n);
      }

      void scalar_unconstrain(T& y) {
        data_r_.push_back(y);
      }

      void scalar_pos_unconstrain(T& y) {
        if (y < 0.0)
          BOOST_THROW_EXCEPTION(std::runtime_error ("y is negative"));
        data_r_.push_back(log(y));
      }

      void scalar_lb_unconstrain(double lb, T& y) {
        if (y < lb)
          BOOST_THROW_EXCEPTION(std::runtime_error ("y is lower than the lower bound"));
        data_r_.push_back(log(y - lb));
      }

      void scalar_ub_unconstrain(double ub, T& y) {
        if (y > ub)
          BOOST_THROW_EXCEPTION(std::runtime_error ("y is higher than the lower bound"));
        data_r_.push_back(log(ub - y));
      }

      void scalar_lub_unconstrain(double lb, double ub, T& y) {
        if (y < lb || y > ub)
          BOOST_THROW_EXCEPTION(std::runtime_error ("y is not between the lower and upper bounds"));
        data_r_.push_back(stan::math::logit((y - lb) / (ub - lb)));
      }

      void corr_unconstrain(T& y) {
        if (y > 1.0 || y < -1.0)
          BOOST_THROW_EXCEPTION(std::runtime_error ("y is not between -1.0 and 1.0"));
        data_r_.push_back(atanh(y));
      }

      void prob_unconstrain(T& y) {
        if (y > 1.0 || y < 0.0)
          BOOST_THROW_EXCEPTION(std::runtime_error ("y is not between 0.0 and 1.0"));
        data_r_.push_back(stan::math::logit(y));
      }

      void ordered_unconstrain(vector_t& y) {
        if (y.size() == 0) return;
        stan::math::check_ordered("stan::io::ordered_unconstrain(%1%)", y, "Vector");
        data_r_.push_back(y[0]);
        for (typename vector_t::size_type i = 1; i < y.size(); ++i) {
          data_r_.push_back(log(y[i] - y[i-1]));
        }
      }

      void positive_ordered_unconstrain(vector_t& y) {
        if (y.size() == 0) return;
        stan::math::check_positive_ordered("stan::io::positive_ordered_unconstrain(%1%)", y, "Vector");
        data_r_.push_back(log(y[0]));
        for (typename vector_t::size_type i = 1; i < y.size(); ++i) {
          data_r_.push_back(log(y[i] - y[i-1]));
        }
      }


      void vector_unconstrain(const vector_t& y) {
        for (typename vector_t::size_type i = 0; i < y.size(); ++i)
          data_r_.push_back(y[i]);
      }

      void row_vector_unconstrain(const vector_t& y) {
        for (typename vector_t::size_type i = 0; i < y.size(); ++i)
          data_r_.push_back(y[i]);
      }

      void matrix_unconstrain(const matrix_t& y) {
        for (typename matrix_t::size_type i = 0; i < y.rows(); ++i)
          for (typename matrix_t::size_type j = 0; j < y.cols(); ++j) 
            data_r_.push_back(y(i,j));
      }

      void vector_lb_unconstrain(double lb, vector_t& y) {
        for (int i = 0; i < y.size(); ++i)
          scalar_lb_unconstrain(lb,y(i));
      }
      void row_vector_lb_unconstrain(double lb, row_vector_t& y) {
        for (int i = 0; i < y.size(); ++i)
          scalar_lb_unconstrain(lb,y(i));
      }
      void matrix_lb_unconstrain(double lb, matrix_t& y) {
        for (int i = 0; i < y.rows(); ++i)
          for (int j = 0; j < y.cols(); ++j)
            scalar_lb_unconstrain(lb,y(i,j));
      }

      void vector_ub_unconstrain(double ub, vector_t& y) {
        for (int i = 0; i < y.size(); ++i)
          scalar_ub_unconstrain(ub,y(i));
      }
      void row_vector_ub_unconstrain(double ub, row_vector_t& y) {
        for (int i = 0; i < y.size(); ++i)
          scalar_ub_unconstrain(ub,y(i));
      }
      void matrix_ub_unconstrain(double ub, matrix_t& y) {
        for (int i = 0; i < y.rows(); ++i)
          for (int j = 0; j < y.cols(); ++j)
            scalar_ub_unconstrain(ub,y(i,j));
      }


      void vector_lub_unconstrain(double lb, double ub, vector_t& y) {
        for (int i = 0; i < y.size(); ++i)
          scalar_lub_unconstrain(lb,ub,y(i));
      }
      void row_vector_lub_unconstrain(double lb, double ub, row_vector_t& y) {
        for (int i = 0; i < y.size(); ++i)
          scalar_lub_unconstrain(lb,ub,y(i));
      }
      void matrix_lub_unconstrain(double lb, double ub, matrix_t& y) {
        for (int i = 0; i < y.rows(); ++i)
          for (int j = 0; j < y.cols(); ++j)
            scalar_lub_unconstrain(lb,ub,y(i,j));
      }

      

      void simplex_unconstrain(vector_t& y) {
        stan::math::check_simplex("stan::io::simplex_unconstrain(%1%)", y, "Vector");
        typename vector_t::size_type k_minus_1 = y.size() - 1;
        double log_y_k = log(y[k_minus_1]);
        for (typename vector_t::size_type i = 0; i < k_minus_1; ++i) {
          data_r_.push_back(log(y[i]) - log_y_k);
        }
      }

      void corr_matrix_unconstrain(matrix_t& y) {
        stan::math::check_corr_matrix("stan::io::corr_matrix_unconstrain(%1%)", y, "Matrix");
        size_t k = y.rows();
        size_t k_choose_2 = (k * (k-1)) / 2;
        array_vec_t cpcs(k_choose_2);
        array_vec_t sds(k);
        bool successful = stan::prob::factor_cov_matrix(cpcs,sds,y);
        if (!successful)
          BOOST_THROW_EXCEPTION(std::runtime_error ("y cannot be factorized by factor_cov_matrix"));
        for (size_t i = 0; i < k; ++i) {
          // sds on log scale unconstrained
          if (fabs(sds[i] - 0.0) >= CONSTRAINT_TOLERANCE)
            BOOST_THROW_EXCEPTION(std::runtime_error ("sds on log scale are unconstrained"));
        }
        for (size_t i = 0; i < k_choose_2; ++i)
          data_r_.push_back(cpcs[i]);
      }

      void cov_matrix_unconstrain(matrix_t& y) {
        typename matrix_t::size_type k = y.rows();
        if (k == 0 || y.cols() != k)
          BOOST_THROW_EXCEPTION(
              std::runtime_error ("y must have elements and y must be a square matrix"));
        typename matrix_t::size_type k_choose_2 = (k * (k-1)) / 2;
        array_vec_t cpcs(k_choose_2);
        array_vec_t sds(k);
        bool successful = stan::prob::factor_cov_matrix(cpcs,sds,y);
        if(!successful)
          BOOST_THROW_EXCEPTION(std::runtime_error ("factor_cov_matrix failed"));
        for (typename matrix_t::size_type i = 0; i < k_choose_2; ++i)
          data_r_.push_back(cpcs[i]);
        for (typename matrix_t::size_type i = 0; i < k; ++i)
          data_r_.push_back(sds[i]);
      }
    };
  }

}

#endif