#pragma once

#include <stan/math/rev/mat/functor/ode_system.hpp>

// define an analytic jacobian for the hornberg for Stan's CVODES
// solver by partial template specialisation


namespace stan {
  namespace math {

    // todo: idealy, we derive the specialised ode_model from the
    // general one which gives possibly syntax problem which is why we
    // chose to use a typedef, maybe we need another helper... not
    // sure.
    template<>
    struct ode_system<pbpkautojac_model_namespace::pbpk_functor__> {
      typedef pbpkautojac_model_namespace::pbpk_functor__ F;
      const F& f_;
      const std::vector<double> theta_;
      const std::vector<double>& x_;
      const std::vector<int>& x_int_;
      std::ostream* msgs_;

      ode_system(const F& f,
                 const std::vector<double> theta,
                 const std::vector<double>& x,
                 const std::vector<int>& x_int,
                 std::ostream* msgs)
	: f_(f),
	  theta_(theta),
	  x_(x),
	  x_int_(x_int),
	  msgs_(msgs)
      {}

      inline void operator()(const double t,
                             const std::vector<double>& y,
//                             std::vector<double>& dy_dt) const {
//	dy_dt = f_(t, y, theta_, x_, x_int_, msgs_);
// the previous two lines for 2.14.0
// the next two lines for 2.17.1
                             Eigen::Map<Eigen::Matrix<double, -1, 1>>& dy_dt) const {
 	std::vector<double> dy_dt_ = f_(t, y, theta_, x_, x_int_, msgs_);
 	for(size_t i = 0; i < dy_dt_.size(); i++)
 	  dy_dt(i) = dy_dt_[i];      }

      template <typename Derived1, typename Derived2>
      inline void
      jacobian(const double t,
               const std::vector<double>& y,
               Eigen::MatrixBase<Derived1>& fy,
               Eigen::MatrixBase<Derived2>& Jy
               ) const {
	using Eigen::VectorXd;
	using std::vector;
	using std::pow;

	const vector<double> f = f_(t, y, theta_, x_, x_int_, msgs_);
	fy = VectorXd::Map(&f[0], f.size());

	Jy.setZero();
	Eigen::Map<Eigen::VectorXd> J = Eigen::Map<Eigen::VectorXd>(&Jy(0,0), Jy.cols()*Jy.rows());
// dy/dt=f(y,theta)
// df/dy	
		#include "states.txt"
      }

      template <typename Derived1, typename Derived2>
      inline void
      jacobian(const double t,
               const std::vector<double>& y,
               Eigen::MatrixBase<Derived1>& fy,
               Eigen::MatrixBase<Derived2>& Jy,
               Eigen::MatrixBase<Derived2>& Jtheta
               ) const {
	using Eigen::VectorXd;
	using std::vector;
	using std::pow;

	jacobian(t, y, fy, Jy);

	//const vector<double>& parms = theta_;

	Jtheta.setZero();

	Eigen::Map<Eigen::VectorXd> J = Eigen::Map<Eigen::VectorXd>(&Jtheta(0,0), Jtheta.cols()*Jtheta.rows());
// df/dtheta
		#include "params.txt"
      }

    };

  }
}