protocol DSL

[martinjrobins]

Might be worth implementing a protocol DSL based on diffsl?

field of hybrid automata looks interesting:

• Alur, R., & Dill, D. L. (1994). A theory of timed automata. Theoretical computer science, 126(2), 183-235.
• https://www.researchgate.net/profile/Jean-Francois-Raskin/publication/227222476_An_Introduction_to_Hybrid_Automata/links/56efac0508ae4b8b5e75764f/An-Introduction-to-Hybrid-Automata.pdf
• Henzinger, T. A. (1996, July). The theory of hybrid automata. In Proceedings 11th Annual IEEE Symposium on Logic in Computer Science (pp. 278-292). IEEE.
• Raskin, J. F. (2005). An introduction to hybrid automata. In Handbook of networked and embedded control systems (pp. 491-517). Boston, MA: Birkhäuser Boston. https://link.springer.com/chapter/10.1007/0-8176-4404-0_21
• Alur, R., Courcoubetis, C., Henzinger, T. A., & Ho, P. H. (1991, June). Hybrid automata: An algorithmic approach to the specification and verification of hybrid systems. In International hybrid systems workshop (pp. 209-229). Berlin, Heidelberg: Springer Berlin Heidelberg. https://ecommons.cornell.edu/server/api/core/bitstreams/4acee44b-8ef7-4ff4-8c08-8a3e3e6c2301/content

e.g.

// Example: single-compartment with drug pulse
model PumpedCell {
  parameters:
    k = 0.1
    kdeg = 0.05
    pulse_amp = 10
    pulse_dur = 5

  variables:
    x: real  // species level
    t_clock: clock

  locations:
    idle:
      invariant: true
      flow:
        dx/dt = -kdeg * x
        dt_clock/dt = 1
      transitions:
        -> pulse
          guard: t_clock >= 10
          reset: t_clock = 0

    pulse:
      invariant: t_clock <= pulse_dur
      flow:
        dx/dt = pulse_amp - kdeg * x
        dt_clock/dt = 1
      transitions:
        -> idle
          guard: t_clock >= pulse_dur
          reset: t_clock = 0

  init:
    location: idle
    x = 0
    t_clock = 0

  outputs:
    sample: sample_every(0.1, x)
}

[martinjrobins]

in diffsl:

// Example: single-compartment with drug pulse
in_i {
    k = 0.1
    kdeg = 0.05
    pulse_amp = 10
    pulse_dur = 5
}

u_i [idle] {
    x = 0  // species level
    t_clock = 0,
}

F_i [idle] { -kdeg * x, 1 }

stop_i [idle->pulse] { 10 - t_clock }
reset_i [idle->pulse] {  x, 0 }
 
F_i [pulse] { pulse_amp - kdeg * x, 1 }
stop_i [pulse->idle] { pulse_dur - t_clock }
reset_i [pulse->idle] { x, 0}

out_i { x }

[martinjrobins]

would this codegen to multiple rhs, stops, resets?

need to keep track of what intermediate tensors are relevent to each rhs/stop/resets

fn rhs0()
fn rhs1()
fn stop0() // this would have num elem equal to possible transitions
fn stop1()
fn reset0_0() //num elem equal to num tstates
fn reset1_0() 

[martinjrobins]

or would this be syntactical sugar that generated n modules?

question would be what to do about resets? generate a separate module for these? Could this be a separate DSL to encode reset functions

u_i { x = 1, y=1}
reset_i [idle->pulse] { 2 * x, y + 10}
reset_i [pulse->idle] { 2 * x, y + 10}

how about if you have a bunch of doses? t = [0, 1, 2], d = [2.0, 3.0. 1.1]?

u_i { x = 1, y=1}
d { 2.0, 3.0, 1.1 }
reset_i [idle->idle] { 2 * x, y + d[n]}

[martinjrobins]

Maybe I could introduce a couple of reserved integer variables indicating the current model index N, and a model type T, and an augentation of an expression into expr | expr | ... to indicate different expressions for different model types T.

the length of stop is the number of different model types and indicates a transition to model T. The index N increments every time you transition to model T (different N's for different Ts)

Could start with just N (implication that T = 1), can do a lot with that. Having T would be more flexible, but can leave that till later.

// Example: single-compartment with drug pulse
in_i {
    k = 0.1
    kdeg = 0.05
    pulse_amp = 10
    pulse_dur = 5
}

amp { 0, pulse_amp }
dur { 10, pulse_dur }
u_i {
    x = 0  // species level
    t_clock = 0,
}

F_i { amp[N % 2] - kdeg * x , 1 }

reset_i {  x, 0 }
stop_i {dur[N % 2] }

out_i { x }

eeep, need integer operations like %.....

[martinjrobins]

semantics for reset: if reset exists, stop_i[0] is th trigger for incrementing N, running the reset and continuing. stop_i[1:] are normal stopping conditions. If reset does not exist, then stop_i[:] are all normal stopping conditions

[martinjrobins]

this N index is more like switched systems in control theory, with resets and with N only being able to increment at each transition

[martinjrobins]

would be nice to be able to control how N is set during transitions (rather than just increment)

switch [N -> N + 1] { 
   t_clock - dur[N % 2], 
}

maybe have both switch and stop?

could start with just switch with no control over N, implication is N -> N + 1

switch {  t_clock - dur[N % 2], }
reset_i { x, 0 }
stop { 100 - N }

update: leave this for now. I don't want to introduce yet another tensor type, lets just do reset for now

[martinjrobins]

Perhaps the stop index gives the next value for N? It can be inefficient in some cases, but the semantics are good. If the value of N gives an out-of-bounds error then we naturally stop, so we can have both transitions and pure stop conditions

setting piecewise dose levels, stop gives the next N. Disadvantage is stop is wasteful since < N will never trigger

start_i  {0, 1, 2}
level_i  {1, 2, 3}
u { x }
F { x + level[N] }
stop_i {  t - start_i }

instantaneous dosing

start_i  {0, 1, 2}
level_i  {1, 2, 3}
u { x }
F { x }
reset_i { x + level_i[N] }
stop_i {  t - start_i }

bouncing ball, this at least is efficient

g { 9.81 } h { 10.0 }
u_i {x = h, v = 0, }
F_i { v, -g, }
stop { x, }
reset {max(x, 0), -e * v }

[martinjrobins]

perhaps this?? truncate the stop vector to those elements we know can be triggered?

stop_i {  (N:N+1): t - start_i[N], (3): 0 }

So stop_i is a sparse vector of length 3 but with 1 nnz, the way this tensor is evaluated is exactly the same as

stop_i {  t - start_i[N] }

but the former will correctly set N

[martinjrobins]

I'm pretty sure I assume everywhere that stop is dense :( this will have to be refactored.

And the api is slightly more complicated.

fn calc_stop(t, x, data, stop)

the length of stop should be the nnz of stop, no longer its length. Before I'd always assumed that the nnz and length of x, stop, F etc were equal. This is probably built in a few places, so will need to be fixed

[martinjrobins]

note that this:

stop_i {  (N:3): t - start_i[N:3], (3): 0 }

should error because the nnz of stop is not constant, it varies with N. Using N in the element indice is only allowed if the nnz is constant. And we'd need to restrict it to vectors