add prototype calendar timeline #47

src/Module.cpp

@@ -75,6 +75,22 @@ PYBIND11_MODULE(neworder, m)
     .def("at_end", &no::Timeline::at_end, timeline_at_end_docstr)
     .def("__repr__", &no::Timeline::repr, timeline_repr_docstr);
 
+  py::class_<no::CalendarTimeline>(m, "CalendarTimeline", "Timestepping functionality")
+    .def(py::init<std::chrono::system_clock::time_point, std::chrono::system_clock::time_point>(), "start"_a, "end"_a)
+    // TODO remove next once tested
+    .def("next", &no::CalendarTimeline::next, timeline_start_docstr)
+    .def("start", &no::CalendarTimeline::start, timeline_start_docstr)
+    .def("end", &no::CalendarTimeline::end, timeline_end_docstr)
+    .def("index", &no::CalendarTimeline::index, timeline_index_docstr)
+    .def("time", &no::CalendarTimeline::time, timeline_time_docstr)
+    .def("dt", &no::CalendarTimeline::dt, timeline_dt_docstr)
+    .def("nsteps", &no::CalendarTimeline::nsteps, timeline_nsteps_docstr)
+    //.def("next", &no::Timeline::next) not exposed
+    //.def("at_checkpoint", &no::CalendarTimeline::at_checkpoint, timeline_at_checkpoint_docstr)
+    .def("at_end", &no::CalendarTimeline::at_end, timeline_at_end_docstr)
+    .def("dow", &no::CalendarTimeline::dow)
+    .def("__repr__", &no::CalendarTimeline::repr, timeline_repr_docstr);
+
   // Microsimulation (or ABM) model class
   py::class_<no::Model>(m, "Model", "The base model class from which all neworder models should be subclassed")
     .def(py::init<no::Timeline&, const py::function&>(),

src/NPArray.cpp

@@ -10,11 +10,11 @@ py::array_t<bool> no::time::isnever_a(const py::array_t<double>& x)
 
 py::array_t<double> no::logistic(const py::array_t<double>& x, double x0, double k)
 {
-  return no::unary_op<double, double>(x, [&](double x) { return 1.0 / (1.0 + exp(-k*(x-x0))); });
+  return no::unary_op<double, double>(x, [&](double x) { return 1.0 / (1.0 + std::exp(-k*(x-x0))); });
 }
 
 py::array_t<double> no::logit(const py::array_t<double>& x)
 {
   // NB no check for x not in [0,1)
-  return no::unary_op<double, double>(x, [](double x) { return log(x/(1.0 - x)); });
+  return no::unary_op<double, double>(x, [](double x) { return std::log(x/(1.0 - x)); });
 }

src/Timeline.cpp

@@ -102,6 +102,127 @@ std::string no::Timeline::repr() const
           % start() % end() % checkpoints() % index();
 }
 
+
+namespace {
+
+// for incrementing time in months preserving day of month
+// e.g. passing 2020,1 will return 29 (leap)
+int daysInFollowingMonth(int year, int month)
+{
+  month +=1;
+  if (month == 12)
+  {
+    year += 1;
+    month -= 12;
+  }
+  static const int days[]{31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31};
+  int d = days[month];
+  if (month == 1 && ((year % 100 != 0) ^ (year % 400 == 0)) && (year % 4 == 0)) // February of a leap year
+    ++d;
+  return d;
+}
+
+}
+
+
+no::CalendarTimeline::CalendarTimeline(time_point start, time_point end) : m_index(0)
+{
+  //size_t i = 0;
+  m_times.push_back(start);
+  time_point time = start;
+
+  for(;;)
+  {
+    std::time_t t = std::chrono::system_clock::to_time_t(time);
+    tm* local_tm = std::localtime(&t);
+    // track whether we cross a DST change
+    int dst_prev = local_tm->tm_isdst;
+    local_tm->tm_mday += daysInFollowingMonth(local_tm->tm_year, local_tm->tm_mon);
+    std::mktime(local_tm);
+    //no::log("h: %% dst: %% prev: %%"s % local_tm->tm_hour % local_tm->tm_isdst % dst_prev);
+
+    // adjust so that hour of day is preserved across DST changes
+    if (local_tm->tm_isdst == 0 && dst_prev == 1)
+    {
+      local_tm->tm_hour += 1;
+    }
+    else if (local_tm->tm_isdst == 1 && dst_prev == 0)
+    {
+      local_tm->tm_hour -= 1;
+    }
+
+    t = std::mktime(local_tm);
+
+    time = std::chrono::system_clock::from_time_t(t);
+    if (time >= end)
+      break;
+    m_times.push_back(time);
+  }
+  m_times.push_back(end);
+}
+
+no::CalendarTimeline::time_point no::CalendarTimeline::start() const
+{
+  return m_times.front();
+}
+
+no::CalendarTimeline::time_point no::CalendarTimeline::end() const
+{
+  return m_times.back();
+}
+
+
+size_t no::CalendarTimeline::index() const
+{
+  return m_index;
+}
+
+
+bool no::CalendarTimeline::at_end() const
+{
+  return m_index >= m_times.size();
+}
+
+no::CalendarTimeline::time_point no::CalendarTimeline::time() const
+{
+  return m_times[m_index];
+}
+
+void no::CalendarTimeline::next()
+{
+  ++m_index;
+}
+
+double no::CalendarTimeline::dt() const
+{
+  if (m_index < m_times.size() - 1)
+    return std::chrono::duration_cast<std::chrono::seconds>(m_times[m_index+1] - m_times[m_index]).count();
+  return 0.0;
+}
+
+size_t no::CalendarTimeline::nsteps() const
+{
+  return m_times.size();
+}
+
+
+// Sun=0, Mon=1 etc
+int no::CalendarTimeline::dow() const
+{
+  std::time_t t = std::chrono::system_clock::to_time_t(m_times[m_index]);
+  tm* local_tm = std::localtime(&t);
+  return local_tm->tm_wday;
+}
+
+std::string no::CalendarTimeline::repr() const
+{
+  std::time_t t = std::chrono::system_clock::to_time_t(m_times[m_index]);
+  std::string buf(64, 0);
+  std::strftime(buf.data(), buf.size(), "%F %T %Z", std::localtime(&t));
+  return std::string(buf);
+}
+
+
 // returns a floating point number that compares less than any other number
 double no::time::distant_past()
 {

src/Timeline.h

@@ -2,13 +2,47 @@
 
 #include "NewOrder.h"
 
+#include "Log.h"
+
+#include <pybind11/chrono.h>
+
 #include <vector>
+#include <chrono>
 #include <cstddef>
 
 
 namespace no {
 
-class NEWORDER_EXPORT Timeline final
+// // Abstract base class for timelines
+// template<typename T>
+// class Timeline
+// {
+// public:
+//   virtual ~Timeline() = default;
+
+//   T time() const;
+//   size_t index() const;
+
+//   T start() const;
+//   T end() const;
+//   size_t nsteps() const;
+
+//   double dt() const;
+//   const std::vector<size_t>& checkpoints() const;
+
+//   void next();
+
+//   bool at_checkpoint() const;
+
+//   bool at_end() const;
+
+//   // used by python __repr__
+//   std::string repr() const;
+
+
+// }
+
+class NEWORDER_EXPORT Timeline final //: public Timeline<double>
 {
 public:
 
@@ -17,7 +51,7 @@ class NEWORDER_EXPORT Timeline final
 
   Timeline(double start, double end, const std::vector<size_t>& checkpoints);
 
-  ~Timeline() = default;
+  virtual ~Timeline() = default;
 
   Timeline(const Timeline&) = default;
   Timeline& operator=(const Timeline&) = default;
@@ -37,7 +71,6 @@ class NEWORDER_EXPORT Timeline final
   void next();
 
   bool at_checkpoint() const;
-
   bool at_end() const;
 
   // used by python __repr__
@@ -52,6 +85,41 @@ class NEWORDER_EXPORT Timeline final
   std::vector<size_t> m_checkpoints;
 };
 
+
+class NEWORDER_EXPORT CalendarTimeline
+{
+public:
+  using time_point = std::chrono::system_clock::time_point;
+
+  // TODO specify time increment in days, months or years
+  // TODO specify checkpoints (as multiple of steps)
+  CalendarTimeline(time_point start, time_point end);
+
+  time_point time() const;
+  size_t index() const;
+
+  time_point start() const;
+  time_point end() const;
+  size_t nsteps() const;
+
+  double dt() const;
+  //const std::vector<size_t>& checkpoints() const;
+
+  void next();
+
+  //bool at_checkpoint() const;
+  bool at_end() const;
+
+  // weekday of current time point. Sun=0, Mon=1 etc
+  int dow() const;
+
+  std::string repr() const;
+
+private:
+  size_t m_index;
+  std::vector<time_point> m_times;
+};
+
 namespace time {
 
   // returns a floating point number that compares unequal (and unordered) to any other number