add output fn #17

hillalex · 2022-04-21T12:23:11Z

Ok, still not sure if this is right but see what you think. mod$run returns full state including output (constrained by index if present), and mod$update_state allows passing in state objects that have the full output size but in that case just truncates to the non-derivable values.

codecov · 2022-04-21T14:48:36Z

Codecov Report

Merging #17 (1ff17f7) into main (65aa8b9) will not change coverage.
The diff coverage is n/a.

❗ Current head 1ff17f7 differs from pull request most recent head 17602c7. Consider uploading reports for the commit 17602c7 to get more accurate results

@@           Coverage Diff           @@
##             main      #17   +/-   ##
=======================================
  Coverage   95.71%   95.71%           
=======================================
  Files           3        3           
  Lines          70       70           
=======================================
  Hits           67       67           
  Misses          3        3

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 65aa8b9...17602c7. Read the comment docs.

…op extra cols

hillalex · 2022-04-25T12:37:44Z

inst/include/mode/mode.hpp

-  size_t n_state_full() {
-    return solver_[0].size();
+  size_t n_output_full() {
+    return m_.size() + m_.n_output();


Have renamed these to distinguish between total output including derived values, and just the model size. Not sure about these names though, thoughts?

hillalex · 2022-04-25T12:39:53Z

inst/template/mode.R.template

    },

    update_state = function(pars = NULL, time = NULL,
                            state = NULL, index = NULL,
                            set_initial_state = NULL, reset_step_size = NULL) {
+      dim_state <- dim(state)


Here just dropping extra values/columns if a vector or matrix corresponds to the full output size, but leaving all other validation to the c++ layer as before. So also not sure whether this splitting of validation logic makes sense, but it keeps the changes here to a minimum.

I'd prefer this in the C++ (mostly because that's the pattern that we've followed in dust) but also happy to defer moving it there until later

richfitz

I've had a think about this and apologise for what looks like some false leads I've suggested. Below is a proposal that might help simplify this by tidying up the names and centralising into stepper the control over these things. We talked before about separation of concerns in terms of who knows about the model, and I am fairly sure that in the end we can end up with only the stepper holding a persistent copy of an initialised model object - this refactor should help with that.

Happy to discuss further tomorrow as this has raised way ambiguous bits than I was expecting

I think that we might refer to a few things

deterministic - the things that enter into the ODEs
stochastic - the things that enter into the stochastic update
variables - the set of deterministic + stochastic
output - quantities derived from the above at the point of reporting
state_full - variables + output
state_run - filtered state from run

We always pack our variables so that they are deterministic followed by stochastic, and so we can exclude the stochastic variables from the error calculation eventually. Or we could leave this distinction as a detail.

Our model then provides methods:

size_t n_variables() const;

The number of deterministic and stochastic variables

size_t n_output() const;

Returns the number of outputs

void rhs(double t,
         const std::vector<double>& y,
         std::vector<double>& dydt);

y and dydt are length n_variables

void output(double t, 
            const std::vector<double>& y,
            std::vector<double>& output);

y is a vector of length n_variables and output is a vector of length n_outputs (sorry, see below for getting this working).

I think that the stepper can lose the state_full_, but pick up a small vector output_ which is of length n_outputs()

We'll need to reflect these changes out into the solver object (replacing size_ with n_variables_ and n_output_.

More complicated on the container object though because this has the default index set which controls what is returned when we call run()

So that then ends up with methods

n_variables():, i.e., what does the model need to start up
n_outputs(): The number of additional outputs
n_state_full: n_variables() + n_outputs()
n_state_run: Filtered state, on [0, n_state_full()]

Then our interface run method looks like:

template <typename T>
cpp11::sexp mode_run(SEXP ptr, double end_time) {
  T *obj = cpp11::as_cpp<cpp11::external_pointer<T>>(ptr).get();
  auto time = obj->time();
  if (end_time < time) {
    cpp11::stop("'end_time' (%f) must be greater than current time (%f)",
                end_time, time);
  }
  obj->run(end_time);

  // these 3 change minimally:
  std::vector<double> dat(obj->n_state_run() * obj->n_particles());
  obj->state_run(dat);
  return mode::r::state_array(dat, obj->n_state_run(), obj->n_particles());
}

The solver::state() method needs changing then

  std::vector<double>::iterator
  state(const std::vector<size_t>& index,
             std::vector<double>::iterator end_state) {
    return stepper_.state(it, index);
  }

And the stepper (which is the thing that actually holds the model object - and I think we can reduce to just this holding it)

  void state(std::vector<double>::iterator end_state) {
    const n_variables = m_.n_variables();
    bool have_run_output = false;
    for (size_t i = 0; i < n; ++i, ++end_state) {
      if (i < n_variables) {
        *end_state = y[i];
      } else {
        if (!have_run_output) {
          m_.output(t_, y, output);
          have_run_output = true;
        }
        *end_state = output[i - n_variables];
      }
    }
  }

richfitz · 2022-04-25T13:31:24Z

inst/examples/logistic.cpp

    const double N1 = y[0];
    const double N2 = y[1];
    dydt[0] = shared->r1 * N1 * (1 - N1 / shared->K1);
    dydt[1] = shared->r2 * N2 * (1 - N2 / shared->K2);
  }

+  std::vector<double>::iterator
+  output(double t,


This looks like its setting the variables as well as the output - I think that we want the mode function doing the first part

richfitz · 2022-04-25T13:34:52Z

inst/include/mode/mode.hpp

  }

-  size_t n_state() const {
+  size_t n_output() const {


I think that n_state_index or n_state_run or something might capture this one best - it's the number of state variables that we will return when running the model due to the index and is capped as the total of state + output. If we come up with something good here we should backport to dust

Possibly if we introduce n_variables as a term we can make this work

richfitz · 2022-04-25T13:51:06Z

inst/template/mode.R.template

    },

    update_state = function(pars = NULL, time = NULL,
                            state = NULL, index = NULL,
                            set_initial_state = NULL, reset_step_size = NULL) {
+      dim_state <- dim(state)


I'd prefer this in the C++ (mostly because that's the pattern that we've followed in dust) but also happy to defer moving it there until later

richfitz · 2022-04-25T16:34:39Z

inst/include/mode/solver.hpp

-             std::vector<double>::iterator end_state) const {
-    auto y = stepper_.output();
+             std::vector<double>::iterator end_state) {
+    auto it = state_full_.begin();


There's definitely more churn of data here than we need

richfitz · 2022-04-25T16:48:39Z

inst/include/mode/solver.hpp

+             std::vector<double>::iterator end_state) {
+    auto it = state_full_.begin();
+    stepper_.state(it);
+    m_.output(t_, it);


We can move to only the stepper holding m_ which would be preferable - so this definitely needs moving deeper

richfitz

This looks good - you said

[11:24] Hill, Alexandra E
I'll also move that logic about dropping output values on a state update from the R to the c++

did you want to do that in this PR or a following one?

richfitz · 2022-04-26T15:49:49Z

inst/include/mode/solver.hpp

@@ -180,31 +179,22 @@ class solver {
  }

  void set_model(Model m) {
-    m_ = m;


hillalex · 2022-04-27T08:46:18Z

did you want to do that in this PR or a following one?

I'll add here

hillalex · 2022-04-27T12:49:54Z

this code is confusing and fiddly. is there a better way?!

hillalex added 2 commits April 21, 2022 13:22

add output fn

0b85857

tests passing

d78db0f

hillalex added 3 commits April 25, 2022 13:35

return full output from run, accept updating state with output and dr…

b2f7e94

…op extra cols

Merge branch 'main' into mrc-3150

d0d76d0

update merged test

d1c70c8

hillalex marked this pull request as ready for review April 25, 2022 12:46

hillalex requested a review from richfitz April 25, 2022 12:47

hillalex commented Apr 25, 2022

View reviewed changes

richfitz added 2 commits April 25, 2022 14:29

Merge branch 'main' into mrc-3150

483b4c7

Update stochastic model

8eb4828

richfitz requested changes Apr 25, 2022

View reviewed changes

try resolve ambiguities and inefficiencies

84499fb

hillalex requested a review from richfitz April 26, 2022 10:23

richfitz reviewed Apr 26, 2022

View reviewed changes

inst/include/mode/solver.hpp

@@ -180,31 +179,22 @@ class solver {

}

void set_model(Model m) {

m_ = m;

Copy link

Member

richfitz Apr 26, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

nice

move logic to c++

17602c7

hillalex requested a review from richfitz April 27, 2022 12:49

richfitz approved these changes May 12, 2022

View reviewed changes

richfitz merged commit 7907b59 into main May 12, 2022

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

add output fn #17

add output fn #17

hillalex commented Apr 21, 2022 •

edited

Loading

codecov bot commented Apr 21, 2022 •

edited

Loading

hillalex Apr 25, 2022

hillalex Apr 25, 2022

richfitz Apr 25, 2022

richfitz left a comment

richfitz Apr 25, 2022

richfitz Apr 25, 2022

richfitz Apr 25, 2022

richfitz Apr 25, 2022

richfitz Apr 25, 2022

richfitz Apr 25, 2022

richfitz left a comment

richfitz Apr 26, 2022

hillalex commented Apr 27, 2022

hillalex commented Apr 27, 2022

add output fn #17

add output fn #17

Conversation

hillalex commented Apr 21, 2022 • edited Loading

codecov bot commented Apr 21, 2022 • edited Loading

Codecov Report

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

richfitz left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

richfitz left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

hillalex commented Apr 27, 2022

hillalex commented Apr 27, 2022

hillalex commented Apr 21, 2022 •

edited

Loading

codecov bot commented Apr 21, 2022 •

edited

Loading