(** Discrete-time node *)
type ('i, 'o, 'r, 's) dnode =
  DNode of {
    state : 's;                        (** current state                      *)
    step  : 's -> 'i -> 's * 'o;       (** step function                      *)
    reset : 's -> 'r -> 's;            (** reset function                     *)
  }

(** Run a discrete node on a list of inputs *)
val drun : ('i, 'o, 'r, 's) dnode -> 'i list -> 'o list

type time =
  float                                (** [≥ 0.0]                            *)

(** Interval-defined functions *)
type 'a dense =
  { h : time;                          (** horizon                            *)
    f : time -> 'a }                   (** [f : [0, h] -> α]                  *)

(** Continuous-time signal *)
type 'a signal =
  'a dense option

(** Initial value problem (IVP) *)
type ('y, 'yder) ivp =
  { y0   : 'y;                         (** initial position                   *)
    fder : time -> 'y -> 'yder;        (** derivative function on [[0.0, h]]  *)
    h    : time; }                     (** maximal horizon                    *) 

(** ODE solver *)
type ('y, 'yder, 's) csolver =
  (time,                               (** requested horizon                  *) 
   'y dense,                           (** solution approximation             *) 
   ('y, 'yder) ivp,                    (** initial value problem              *) 
   's) dnode

(** Zero-crossing problem (ZCP) *)
type ('y, 'zin) zcp =
  { y0   : 'y;                         (** initial position                   *)
    fzer : time -> 'y -> 'zin;         (** zero-crossing function             *)
    h    : time; }                     (** maximal horizon                    *)

(** Zero-crossing solver *)
type ('y, 'zin, 'zout, 's) zsolver =
  ('y dense,                           (** input value                        *)
   time * 'zout,                       (** horizon and zero-crossing events   *)
   ('y, 'zin) zcp,                     (** zero-crossing problem              *)
   's) dnode

(** Full solver (composition of an ODE and zero-crossing solver) *)
type ('y, 'yder, 'zin, 'zout, 's) solver =
  (time,                               (** requested horizon                  *)
   'y dense * 'zout,                   (** output and zero-crossing events    *)
   ('y, 'yder) ivp * ('y, 'zin) zcp,   (** (re)initialization parameters      *)
   's) dnode

(** Compose an ODE solver and a zero-crossing solver *)
val build_solver : ('y, 'yder, 'cs) csolver ->
                   ('y, 'zin, 'zout, 'zs) zsolver ->
                   ('y, 'yder, 'zin, 'zout, 'cs * 'zs) solver

(** Hybrid (discrete-time and continuous-time) node *)
type ('i, 'o, 'r, 's, 'y, 'yder, 'zin, 'zout) hnode =
  HNode of {
    state : 's;                        (** current state                      *)
    step  : 's -> 'i -> 's * 'o;       (** discrete step function             *)
    reset : 's -> 'r -> 's;            (** reset function                     *)
    fder  : 's -> 'i -> 'y -> 'yder;   (** derivative function                *)
    fzer  : 's -> 'i -> 'y -> 'zin;    (** zero-crossing function             *)
    fout  : 's -> 'i -> 'y -> 'o;      (** continuous output function         *)
    cget  : 's -> 'y;                  (** continuous state getter            *)
    cset  : 's -> 'y -> 's;            (** continuous state setter            *)
    zset  : 's -> 'zout -> 's;         (** zero-crossing information setter   *)
    jump  : 's -> bool;                (** discrete go-again indicator        *)
  }

(** Simulation mode (either discrete or continuous) *)
type mode = D | C

(** Simulation state *)
type ('i, 'o, 'r, 'y, 'yder, 'zin, 'zout, 'ms, 'ss) state =
  State of {
    solver :                           (** solver state                       *)
      ('y, 'yder, 'zin, 'zout, 'ss) solver;
    model :                            (** model state                        *)
      ('i, 'o, 'r, 'ms, 'y, 'yder, 'zin, 'zout) hnode;
    input : 'i signal;                 (** current input                      *)
    time : time;                       (** current time                       *)
    mode : mode;                       (** current step mode                  *)
  }

(** Discrete simulation step *)
val dstep : ('i, 'o, 'r, 'y, 'yder, 'zin, 'zout, 'ms, 'ss) state ->
            ('i, 'o, 'r, 'y, 'yder, 'zin, 'zout, 'ms, 'ss) state * 'o signal

(** Continuous simulation step *)
val cstep : ('i, 'o, 'r, 'y, 'yder, 'zin, 'zout, 'ms, 'ss) state ->
            ('i, 'o, 'r, 'y, 'yder, 'zin, 'zout, 'ms, 'ss) state * 'o signal

(** Simulate a hybrid model with a solver *)
val hsim : ('i, 'o, 'r, 'ms, 'y, 'yder, 'zin, 'zout) hnode ->
           ('y, 'yder, 'zin, 'zout, 'ss) solver ->
           ('i signal, 'o signal, 'r, ('i, 'o, 'r, 'y, 'yder, 'zin, 'zout, 'ms, 'ss) state) dnode

(** Run a simulation on a list of inputs *)
val hrun : ('i, 'o, 'r, 'ms, 'y, 'yder, 'zin, 'zout) hnode ->
           ('y, 'yder, 'zin, 'zout, 'ss) solver ->
           'i dense list ->
           'o dense list