removed: t/k arguments in NonLinModel signature

Can be implemented with a new measured disturbance. Seems like a bloat feature.

Author: franckgaga

src/estimator/execute.jl

@@ -40,7 +40,7 @@ function f̂!(x̂next0, û0, estim::StateEstimator, model::SimModel, x̂0, u0,
     @views x̂d_next, x̂s_next = x̂next0[1:model.nx], x̂next0[model.nx+1:end]
     mul!(û0, estim.Cs_u, x̂s)
     û0 .+= u0
-    f!(x̂d_next, model, x̂d, û0, d0)
+    f!(x̂d_next, model, x̂d, û0, d0, model.p)
     mul!(x̂s_next, estim.As, x̂s)
     return nothing
 end
@@ -65,7 +65,7 @@ Mutating output function ``\mathbf{ĥ}`` of the augmented model, see [`f̂!`](@
 function ĥ!(ŷ0, estim::StateEstimator, model::SimModel, x̂0, d0)
     # `@views` macro avoid copies with matrix slice operator e.g. [a:b]
     @views x̂d, x̂s = x̂0[1:model.nx], x̂0[model.nx+1:end]
-    h!(ŷ0, model, x̂d, d0)
+    h!(ŷ0, model, x̂d, d0, model.p)
     mul!(ŷ0, estim.Cs_y, x̂s, 1, 1)
     return nothing
 end

src/estimator/internal_model.jl

@@ -171,16 +171,16 @@ end
 
 State function ``\mathbf{f̂}`` of [`InternalModel`](@ref) for [`NonLinModel`](@ref).
 
-It calls `model.f!(x̂next0, x̂0, u0 ,d0)` since this estimator does not augment the states.
+It calls `model.f!(x̂next0, x̂0, u0 ,d0, model.p)` since this estimator does not augment the states.
 """
-f̂!(x̂next0, _, ::InternalModel, model::NonLinModel, x̂0, u0, d0) = model.f!(x̂next0, x̂0, u0, d0)
+f̂!(x̂next0, _, ::InternalModel, model::NonLinModel, x̂0, u0, d0) = model.f!(x̂next0, x̂0, u0, d0, model.p)
 
 @doc raw"""
     ĥ!(ŷ0, estim::InternalModel, model::NonLinModel, x̂0, d0)
 
 Output function ``\mathbf{ĥ}`` of [`InternalModel`](@ref), it calls `model.h!`.
 """
-ĥ!(x̂next0, ::InternalModel, model::NonLinModel, x̂0, d0) = model.h!(x̂next0, x̂0, d0)
+ĥ!(x̂next0, ::InternalModel, model::NonLinModel, x̂0, d0) = model.h!(x̂next0, x̂0, d0, model.p)
 
 
 @doc raw"""
@@ -246,7 +246,7 @@ Compute the current stochastic output estimation `ŷs` for [`InternalModel`](@r
 """
 function correct_estimate!(estim::InternalModel, y0m, d0)
     ŷ0d = estim.buffer.ŷ
-    h!(ŷ0d, estim.model, estim.x̂d, d0)
+    h!(ŷ0d, estim.model, estim.x̂d, d0, estim.model.p)
     ŷs = estim.ŷs
     for j in eachindex(ŷs) # broadcasting was allocating unexpectedly, so we use a loop
         if j in estim.i_ym
@@ -279,7 +279,7 @@ function update_estimate!(estim::InternalModel, _ , d0, u0)
     x̂d, x̂s, ŷs = estim.x̂d, estim.x̂s, estim.ŷs
     # -------------- deterministic model ---------------------
     x̂dnext = estim.buffer.x̂
-    f!(x̂dnext, model, x̂d, u0, d0) 
+    f!(x̂dnext, model, x̂d, u0, d0, model.p) 
     x̂d .= x̂dnext # this also updates estim.x̂0 (they are the same object)
     # --------------- stochastic model -----------------------
     x̂snext = estim.x̂snext
@@ -317,7 +317,7 @@ function init_estimate!(estim::InternalModel, model::LinModel{NT}, y0m, d0, u0)
     # TODO: use estim.buffer.x̂ to reduce the allocation:
     x̂d .= (I - model.A)\(model.Bu*u0 + model.Bd*d0 + model.fop - model.xop)
     ŷ0d = estim.buffer.ŷ
-    h!(ŷ0d, model, x̂d, d0)
+    h!(ŷ0d, model, x̂d, d0, model.p)
     ŷs = ŷ0d
     ŷs[estim.i_ym] .= @views y0m .- ŷ0d[estim.i_ym]
     # ŷs=0 for unmeasured outputs :

src/estimator/kalman.jl

@@ -582,9 +582,10 @@ See [`SteadyKalmanFilter`](@ref) for details on ``\mathbf{v}(k), \mathbf{w}(k)``
 \text{diag}(Q, Q_{int_u}, Q_{int_{ym}})}`` and ``\mathbf{R̂ = R}``. The functions
 ``\mathbf{f̂, ĥ}`` are `model` state-space functions augmented with the stochastic model of
 the unmeasured disturbances, which is specified by the numbers of integrator `nint_u` and
-`nint_ym` (see Extended Help). The ``\mathbf{ĥ^m}`` function represents the measured outputs
-of ``\mathbf{ĥ}`` function (and unmeasured ones, for ``\mathbf{ĥ^u}``). The matrix 
-``\mathbf{P̂}`` is the estimation error covariance of `model` state augmented with the 
+`nint_ym` (see Extended Help). Model parameters ``\mathbf{p}`` are not an argument of
+``\mathbf{f̂, ĥ}`` functions for conciseness. The ``\mathbf{ĥ^m}`` function represents the
+measured outputs of ``\mathbf{ĥ}`` function (and unmeasured ones, for ``\mathbf{ĥ^u}``). The
+matrix ``\mathbf{P̂}`` is the estimation error covariance of `model` state augmented with the 
 stochastic ones. Three keyword arguments specify its initial value with ``\mathbf{P̂}_{-1}(0) = 
 \mathrm{diag}\{ \mathbf{P}(0), \mathbf{P_{int_{u}}}(0), \mathbf{P_{int_{ym}}}(0) \}``. The 
 initial state estimate ``\mathbf{x̂}_{-1}(0)`` can be manually specified with [`setstate!`](@ref).

src/model/linearization.jl

@@ -121,22 +121,22 @@ function linearize!(
     A::Matrix{NT}, Bu::Matrix{NT}, Bd::Matrix{NT}  = linmodel.A, linmodel.Bu, linmodel.Bd
     C::Matrix{NT}, Dd::Matrix{NT} = linmodel.C, linmodel.Dd
     xnext0::Vector{NT}, y0::Vector{NT} = linmodel.buffer.x, linmodel.buffer.y
-    myf_x0!(xnext0, x0) = f!(xnext0, nonlinmodel, x0, u0, d0)
-    myf_u0!(xnext0, u0) = f!(xnext0, nonlinmodel, x0, u0, d0)
-    myf_d0!(xnext0, d0) = f!(xnext0, nonlinmodel, x0, u0, d0)
-    myh_x0!(y0, x0) = h!(y0, nonlinmodel, x0, d0)
-    myh_d0!(y0, d0) = h!(y0, nonlinmodel, x0, d0)
+    myf_x0!(xnext0, x0) = f!(xnext0, nonlinmodel, x0, u0, d0, model.p)
+    myf_u0!(xnext0, u0) = f!(xnext0, nonlinmodel, x0, u0, d0, model.p)
+    myf_d0!(xnext0, d0) = f!(xnext0, nonlinmodel, x0, u0, d0, model.p)
+    myh_x0!(y0, x0) = h!(y0, nonlinmodel, x0, d0, model.p)
+    myh_d0!(y0, d0) = h!(y0, nonlinmodel, x0, d0, model.p)
     ForwardDiff.jacobian!(A,  myf_x0!, xnext0, x0)
     ForwardDiff.jacobian!(Bu, myf_u0!, xnext0, u0)
     ForwardDiff.jacobian!(Bd, myf_d0!, xnext0, d0)
     ForwardDiff.jacobian!(C,  myh_x0!, y0, x0)
     ForwardDiff.jacobian!(Dd, myh_d0!, y0, d0)
     # --- compute the nonlinear model output at operating points ---
-    h!(y0, nonlinmodel, x0, d0)
+    h!(y0, nonlinmodel, x0, d0, model.p)
     y  = y0
     y .= y0 .+ nonlinmodel.yop
     # --- compute the nonlinear model next state at operating points ---
-    f!(xnext0, nonlinmodel, x0, u0, d0)
+    f!(xnext0, nonlinmodel, x0, u0, d0, model.p)
     xnext  = xnext0
     xnext .= xnext0 .+ nonlinmodel.fop .- nonlinmodel.xop
     # --- modify the linear model operating points ---

src/model/linmodel.jl

@@ -5,7 +5,7 @@ struct LinModel{NT<:Real} <: SimModel{NT}
     Bd  ::Matrix{NT}
     Dd  ::Matrix{NT}
     x0::Vector{NT}
-    k::Vector{Int}
+    p ::Nothing
     Ts::NT
     t::Vector{NT}
     nu::Int
@@ -36,8 +36,8 @@ struct LinModel{NT<:Real} <: SimModel{NT}
         size(C)  == (ny,nx) || error("C size must be $((ny,nx))")
         size(Bd) == (nx,nd) || error("Bd size must be $((nx,nd))")
         size(Dd) == (ny,nd) || error("Dd size must be $((ny,nd))")
+        p = nothing # the parameter field p is not used for LinModel
         Ts > 0 || error("Sampling time Ts must be positive")
-        k = [0]
         uop = zeros(NT, nu)
         yop = zeros(NT, ny)
         dop = zeros(NT, nd)
@@ -52,8 +52,9 @@ struct LinModel{NT<:Real} <: SimModel{NT}
         buffer = SimModelBuffer{NT}(nu, nx, ny, nd)
         return new{NT}(
             A, Bu, C, Bd, Dd, 
-            x0, 
-            k, Ts, t,
+            x0,
+            p,
+            Ts, t,
             nu, nx, ny, nd, 
             uop, yop, dop, xop, fop,
             uname, yname, dname, xname,
@@ -256,11 +257,11 @@ function steadystate!(model::LinModel, u0, d0)
 end
 
 """
-    f!(xnext0, model::LinModel, x0, u0, d0, _ , _ ) -> nothing
+    f!(xnext0, model::LinModel, x0, u0, d0, p) -> nothing
 
 Evaluate `xnext0 = A*x0 + Bu*u0 + Bd*d0` in-place when `model` is a [`LinModel`](@ref).
 """
-function f!(xnext0, model::LinModel, x0, u0, d0, _ , _ )
+function f!(xnext0, model::LinModel, x0, u0, d0, _ )
     mul!(xnext0, model.A,  x0)
     mul!(xnext0, model.Bu, u0, 1, 1)
     mul!(xnext0, model.Bd, d0, 1, 1)
@@ -269,11 +270,11 @@ end
 
 
 """
-    h!(y0, model::LinModel, x0, d0, _ , _ ) -> nothing
+    h!(y0, model::LinModel, x0, d0, p) -> nothing
 
 Evaluate `y0 = C*x0 + Dd*d0` in-place when `model` is a [`LinModel`](@ref).
 """
-function h!(y0, model::LinModel, x0, d0, _ , _ )
+function h!(y0, model::LinModel, x0, d0, _ )
     mul!(y0, model.C,  x0)
     mul!(y0, model.Dd, d0, 1, 1)
     return nothing

src/model/nonlinmodel.jl

@@ -6,7 +6,6 @@ struct NonLinModel{
     h!::H
     p::P
     solver::DS
-    k::Vector{Int}
     Ts::NT
     t::Vector{NT}
     nu::Int
@@ -27,7 +26,6 @@ struct NonLinModel{
         f!::F, h!::H, Ts, nu, nx, ny, nd, p::P, solver::DS
     ) where {NT<:Real, F<:Function, H<:Function, P<:Any, DS<:DiffSolver}
         Ts > 0 || error("Sampling time Ts must be positive")
-        k = [0]
         uop = zeros(NT, nu)
         yop = zeros(NT, ny)
         dop = zeros(NT, nd)
@@ -45,7 +43,7 @@ struct NonLinModel{
             f!, h!,
             p,
             solver, 
-            k, Ts, t,
+            Ts, t,
             nu, nx, ny, nd, 
             uop, yop, dop, xop, fop,
             uname, yname, dname, xname,
@@ -65,18 +63,18 @@ continuous dynamics. Use `solver=nothing` for the discrete case (see Extended He
 functions are defined as:
 ```math
 \begin{aligned}
-    \mathbf{ẋ}(t) &= \mathbf{f}\Big( \mathbf{x}(t), \mathbf{u}(t), \mathbf{d}(t), \mathbf{p}(t), t \Big) \\
-    \mathbf{y}(t) &= \mathbf{h}\Big( \mathbf{x}(t), \mathbf{d}(t), \mathbf{p}(t), t \Big)
+    \mathbf{ẋ}(t) &= \mathbf{f}\Big( \mathbf{x}(t), \mathbf{u}(t), \mathbf{d}(t), \mathbf{p} \Big) \\
+    \mathbf{y}(t) &= \mathbf{h}\Big( \mathbf{x}(t), \mathbf{d}(t), \mathbf{p} \Big)
 \end{aligned}
 ```
 where ``\mathbf{x}``, ``\mathbf{y}, ``\mathbf{u}``, ``\mathbf{d}`` and ``\mathbf{p}`` are
 respectively the state, output, manipulated input, measured disturbance and parameter vectors,
 and ``t`` the time in second. The functions can be implemented in two possible ways:
 
-1. **Non-mutating functions** (out-of-place): define them as `f(x, u, d, p, t) -> ẋ` and
-   `h(x, d, p, t) -> y`. This syntax is simple and intuitive but it allocates more memory.
-2. **Mutating functions** (in-place): define them as `f!(ẋ, x, u, d, p, t) -> nothing` and
-   `h!(y, x, d, p, t) -> nothing`. This syntax reduces the allocations and potentially the 
+1. **Non-mutating functions** (out-of-place): define them as `f(x, u, d, p) -> ẋ` and
+   `h(x, d, p) -> y`. This syntax is simple and intuitive but it allocates more memory.
+2. **Mutating functions** (in-place): define them as `f!(ẋ, x, u, d, p) -> nothing` and
+   `h!(y, x, d, p) -> nothing`. This syntax reduces the allocations and potentially the 
    computational burden as well.
 
 `Ts` is the sampling time in second. `nu`, `nx`, `ny` and `nd` are the respective number of 
@@ -85,7 +83,7 @@ is the parameters of the model passed to the two functions. It can be of any Jul
 but use a mutable type if you want to change them later e.g.: a vector.
 
 !!! tip
-    Replace the `d`, `p` or `t`  argument with `_` if not needed (see Examples below).
+    Replace the `d` or `p` argument with `_` in your functions if not needed (see Examples below).
     
 A 4th order [`RungeKutta`](@ref) solver discretizes the differential equations by default. 
 The rest of the documentation assumes discrete dynamics since all models end up in this 
@@ -100,9 +98,9 @@ See also [`LinModel`](@ref).
 
 # Examples
 ```jldoctest
-julia> f!(ẋ, x, u, _ ) = (ẋ .= -0.2x .+ u; nothing);
+julia> f!(ẋ, x, u, _ , _ ) = (ẋ .= -0.2x .+ u; nothing);
 
-julia> h!(y, x, _ ) = (y .= 0.1x; nothing);
+julia> h!(y, x, _ , _ ) = (y .= 0.1x; nothing);
 
 julia> model1 = NonLinModel(f!, h!, 5.0, 1, 1, 1)               # continuous dynamics
 NonLinModel with a sample time Ts = 5.0 s, RungeKutta solver and:
@@ -111,9 +109,9 @@ NonLinModel with a sample time Ts = 5.0 s, RungeKutta solver and:
  1 outputs y
  0 measured disturbances d
 
-julia> f(x, u, _ ) = 0.1x + u;
+julia> f(x, u, _ , _ ) = 0.1x + u;
 
-julia> h(x, _ ) = 2x;
+julia> h(x, _ , _ ) = 2x;
 
 julia> model2 = NonLinModel(f, h, 2.0, 1, 1, 1, solver=nothing) # discrete dynamics
 NonLinModel with a sample time Ts = 2.0 s, empty solver and:
@@ -128,16 +126,16 @@ NonLinModel with a sample time Ts = 2.0 s, empty solver and:
     State-space functions are similar for discrete dynamics:
     ```math
     \begin{aligned}
-        \mathbf{x}(k+1) &= \mathbf{f}\Big( \mathbf{x}(k), \mathbf{u}(k), \mathbf{d}(k), \mathbf{p}(k), k \Big) \\
-        \mathbf{y}(k)   &= \mathbf{h}\Big( \mathbf{x}(k), \mathbf{d}(k), \mathbf{p}(k), k \Big)
+        \mathbf{x}(k+1) &= \mathbf{f}\Big( \mathbf{x}(k), \mathbf{u}(k), \mathbf{d}(k), \mathbf{p}(k) \Big) \\
+        \mathbf{y}(k)   &= \mathbf{h}\Big( \mathbf{x}(k), \mathbf{d}(k), \mathbf{p}(k) \Big)
     \end{aligned}
     ```
     with two possible implementations as well:
 
-    1. **Non-mutating functions**: define them as `f(x, u, d, p, k) -> xnext` and 
-       `h(x, d, p, k) -> y`.
-    2. **Mutating functions**: define them as `f!(xnext, x, u, d, p, k) -> nothing` and
-       `h!(y, x, d, p, k) -> nothing`.
+    1. **Non-mutating functions**: define them as `f(x, u, d, p) -> xnext` and 
+       `h(x, d, p) -> y`.
+    2. **Mutating functions**: define them as `f!(xnext, x, u, d, p) -> nothing` and
+       `h!(y, x, d, p) -> nothing`.
 """
 function NonLinModel{NT}(
     f::Function, h::Function, Ts::Real, nu::Int, nx::Int, ny::Int, nd::Int=0;
@@ -146,8 +144,8 @@ function NonLinModel{NT}(
     isnothing(solver) && (solver=EmptySolver())
     ismutating_f = validate_f(NT, f)
     ismutating_h = validate_h(NT, h)
-    f! = ismutating_f ? f : f!(xnext, x, u, d, p, k) = (xnext .= f(x, u, d, p, k); nothing)
-    h! = ismutating_h ? h : h!(y, x, d, p, k) = (y .= h(x, d, p, k); nothing)
+    f! = ismutating_f ? f : f!(xnext, x, u, d, p) = (xnext .= f(x, u, d, p); nothing)
+    h! = ismutating_h ? h : h!(y, x, d, p) = (y .= h(x, d, p); nothing)
     f!, h! = get_solver_functions(NT, solver, f!, h!, Ts, nu, nx, ny, nd)
     F, H, P, DS = get_types(f!, h!, p, solver)
     return NonLinModel{NT, F, H, P, DS}(f!, h!, Ts, nu, nx, ny, nd, p, solver)
@@ -173,13 +171,13 @@ end
 Validate `f` function argument signature and return `true` if it is mutating.
 """
 function validate_f(NT, f)
-    ismutating = hasmethod(f, Tuple{Vector{NT}, Vector{NT}, Vector{NT}, Vector{NT}, Any, NT})
-    if !(ismutating || hasmethod(f, Tuple{Vector{NT}, Vector{NT}, Vector{NT}, Any, NT}))
+    ismutating = hasmethod(f, Tuple{Vector{NT}, Vector{NT}, Vector{NT}, Vector{NT}, Any})
+    if !(ismutating || hasmethod(f, Tuple{Vector{NT}, Vector{NT}, Vector{NT}, Any}))
         error(
             "the state function has no method with type signature "*
-            "f(x::Vector{$(NT)}, u::Vector{$(NT)}, d::Vector{$(NT)}, p::Any, t::$(NT)) or "*
+            "f(x::Vector{$(NT)}, u::Vector{$(NT)}, d::Vector{$(NT)}, p::Any) or "*
             "mutating form f!(xnext::Vector{$(NT)}, x::Vector{$(NT)}, u::Vector{$(NT)}, "*
-                                                   "d::Vector{$(NT)}, p::Any, t::$(NT))"
+                                                   "d::Vector{$(NT)}, p::Any)"
         )
     end
     return ismutating
@@ -191,12 +189,12 @@ end
 Validate `h` function argument signature and return `true` if it is mutating.
 """
 function validate_h(NT, h)
-    ismutating = hasmethod(h, Tuple{Vector{NT}, Vector{NT}, Vector{NT}, Any, NT})
-    if !(ismutating || hasmethod(h, Tuple{Vector{NT}, Vector{NT}, Any, NT}))
+    ismutating = hasmethod(h, Tuple{Vector{NT}, Vector{NT}, Vector{NT}, Any})
+    if !(ismutating || hasmethod(h, Tuple{Vector{NT}, Vector{NT}, Any}))
         error(
             "the output function has no method with type signature "*
-            "h(x::Vector{$(NT)}, d::Vector{$(NT)}, p::Any, t::$(NT)) or mutating form "*
-            "h!(y::Vector{$(NT)}, x::Vector{$(NT)}, d::Vector{$(NT)}, p::Any, t::$(NT))"
+            "h(x::Vector{$(NT)}, d::Vector{$(NT)}, p::Any) or mutating form "*
+            "h!(y::Vector{$(NT)}, x::Vector{$(NT)}, d::Vector{$(NT)}, p::Any)"
         )
     end
     return ismutating
@@ -205,11 +203,11 @@ end
 "Do nothing if `model` is a [`NonLinModel`](@ref)."
 steadystate!(::SimModel, _ , _ ) = nothing
 
-"Call `f!(xnext0, x0, u0, d0, p, k)` with `model.f!` method for [`NonLinModel`](@ref)."
-f!(xnext0, model::NonLinModel, x0, u0, d0, p, k) = model.f!(xnext0, x0, u0, d0, p, k)
+"Call `model.f!(xnext0, x0, u0, d0, p)` for [`NonLinModel`](@ref)."
+f!(xnext0, model::NonLinModel, x0, u0, d0, p) = model.f!(xnext0, x0, u0, d0, p)
 
-"Call `h!(y0, x0, d0, p, k)` with `model.h` method for [`NonLinModel`](@ref)."
-h!(y0, model::NonLinModel, x0, d0, p, k) = model.h!(y0, x0, d0, p, k)
+"Call `model.h!(y0, x0, d0, p)` for [`NonLinModel`](@ref)."
+h!(y0, model::NonLinModel, x0, d0, p) = model.h!(y0, x0, d0, p)
 
 detailstr(model::NonLinModel) = ", $(typeof(model.solver).name.name) solver"
 detailstr(::NonLinModel{<:Real, <:Function, <:Function, <:Any, <:EmptySolver}) = ", empty solver"

src/model/solver.jl

@@ -45,7 +45,7 @@ function get_solver_functions(NT::DataType, solver::RungeKutta, fc!, hc!, Ts, _
     k2_cache::DiffCache{Vector{NT}, Vector{NT}}   = DiffCache(zeros(NT, nx), Nc)
     k3_cache::DiffCache{Vector{NT}, Vector{NT}}   = DiffCache(zeros(NT, nx), Nc)
     k4_cache::DiffCache{Vector{NT}, Vector{NT}}   = DiffCache(zeros(NT, nx), Nc)
-    f! = function inner_solver_f!(xnext, x, u, d, p, k)
+    f! = function inner_solver_f!(xnext, x, u, d, p)
         CT = promote_type(eltype(x), eltype(u), eltype(d))
         # dummy variable for get_tmp, necessary for PreallocationTools + Julia 1.6 :
         var::CT = 0
@@ -55,23 +55,22 @@ function get_solver_functions(NT::DataType, solver::RungeKutta, fc!, hc!, Ts, _
         k3   = get_tmp(k3_cache, var)
         k4   = get_tmp(k4_cache, var)
         @. xcur = x
-        for k_inner=0:solver.supersample-1
+        for i=1:solver.supersample
             xterm = xnext # TODO: move this out of the loop, just above (to test) ?
-            t = k*Ts + k_inner*Ts_inner
             @. xterm = xcur
-            fc!(k1, xterm, u, d, p, t)
+            fc!(k1, xterm, u, d, p)
             @. xterm = xcur + k1 * Ts_inner/2
-            fc!(k2, xterm, u, d, p, t)
+            fc!(k2, xterm, u, d, p)
             @. xterm = xcur + k2 * Ts_inner/2
-            fc!(k3, xterm, u, d, p, t)
+            fc!(k3, xterm, u, d, p)
             @. xterm = xcur + k3 * Ts_inner
-            fc!(k4, xterm, u, d, p, t)
+            fc!(k4, xterm, u, d, p)
             @. xcur = xcur + (k1 + 2k2 + 2k3 + k4)*Ts_inner/6
         end
         @. xnext = xcur
         return nothing
     end
-    h!(y, x, d, p, k) = (hc!(y, x, d, p, k*Ts); nothing)
+    h! = hc!
     return f!, h!
 end