Heaviside-type forcing functions

[dcnorris]

To apply 'instantaneous' impulses to a discrete-time rprocess, I have defined Heaviside-like cumulative time-series via covariate_table, and am differentiating them to obtain 'Dirac deltas':

             dM <- lookup(ct,t+delta.t)$M - lookup(ct,t)$M

Firstly, I would wonder if there might be some more idiomatic way to do this with pomp.
Secondly, if the above is acceptable practice, then how could one do such lookups in a C snippet?

[kingaa]

Help me try to understand what you're up to. Firstly, you say you have a discrete-time process, yet you simultaneously describe a desire to incorporate impulsive forcing, which typically requires a continuous-time setting to make sense. Perhaps I misunderstand.

You should also know that the lookup is handled for you, so that there is no need to call lookup yourself. Each variable mentioned in your covariate table (except the time) will be available to you in each basic model component. In particular, in a C snippet, you can simply refer to that variable.

[dcnorris]

My mental model is ultimately continuous-time, but I'm using the discrete_time process as an approximation — somewhat reassured in so doing by your arguments here BTW. If I could go straight to an idiomatic, continuous-time pomp treatment of these impulses (injections into pharmacological compartments) I would gladly side-step this detour through discretization.

While I apprecate that variables at the present time, t are available automatically in the snippets, this alone does not let me perform the differentiation. Maybe my question is as simple as wondering how to access lagged values of covariates?

[kingaa]

This is a really good question.

Before I address it, I recommend that, since your model is a continuous-time model, you use the euler plug-in rather than the discrete_time one. The main difference between these is that discrete_time assumes that the trajectories of your state process change exactly and only at multiples of the timestep delta.t, while euler makes no such assumption.

On to the main question. There are several ways of implementing the impulsive forcing you describe. I think the conceptually cleanest one is as follows.

First, include a covariate consisting of the sum of your Heaviside forcing functions. In other words, the covariate, call it f, is the cumulative sum of the delta functions. We think of f as the forcing to be applied.

Next, include an additional state variable, say F, which represents the cumulative forcing that has been applied. Initialize this to zero.

Finally, in the rprocess component, include a statement like

if (f != F) {
  X += f - F;
  F = f;
}

Here, I've written it as though the effect of the forcing is to increment the state variable X. In your model, presumably, it will enter in a different way.

This is the C snippet form. In an R function, it would be something like

function (X, f, F, ...) {
  if (f != F) {
    X <- X + f - F;
    F <- f;
  }
}

[dcnorris]

This seems like an elegant idiom, one that respects the continuous-time concepts nicely!

[dcnorris]

Interestingly, I find the following behavior when applying this idiom. Is this expected?

I'm using covariate V_ as the strictly increasing cumulative count of a handful of discrete events, and I accumulate this into state variable V. These events trigger the updating of state variable W. The following code (with strict inequality >) works just fine:

if (V_ > V) {
  V = V_;
  W = pow(1450/F, 0.3);
}

But... when I use the test (V_ != V) as in your example, this block of code recalculates continuously during one of the time intervals between events! I would have thought executing V = V_ once would ensure that (V_ != V) is false upon the next iteration and thereafter until V_ increments with the next occurrence of the event:

I can sanitize my code to post a full repro, but wanted to check first that there isn't an explanation straightaway.

[kingaa]

Some questions:

1. When you provide the lookup table using covariate_table, do you specify order="linear"?
2. What is F?
3. Can you examine the difference V - V_? Is it a question of these things differing by a very small amount?

[dcnorris]

1. Switching to order="constant" eliminates the problem; is this what I should have been using?
2. F is just another state variable, which is always positive and indeed bounded well away from zero.
3. Using the following debugging technique, I get all zeros printed. But with order="constant" I get just one printout per rising step of the forcing function, whereas with order="linear" I get many more printouts than there are rising steps.

           if (V_ != V) {
             V = V_;
             Rprintf(\"diff = %g\\n\", 1e17 * (V - V_));
             W = pow(1450/F, 0.3);
           }

[kingaa]

Yes, sorry, I should have mentioned that at the outset. Internally, pomp interpolates the lookup table at the current time. The default is to use linear interpolation. Setting order="constant" switches to order-0 interpolation (right-continuous with left limits).

[dcnorris]

Many thanks for clarification; so glad to know the fix is this simple!