API — ICoreBlocks/ICoreMath/Optimization
The public contract of 20 header(s) under src/ICoreBlocks/ICoreMath/Optimization — 20 class/struct definition(s), 33 declaration(s). Each section shows the header's banner and its public (and protected-virtual) surface exactly as the file writes it.
ICoreConstrainedNonlinearOptimization.h#
src/ICoreBlocks/ICoreMath/Optimization/ICoreConstrainedNonlinearOptimization.h
ICoreConstrainedNonlinearOptimization#
ICoreConstrainedNonlinearOptimization.h:58 · class · nested Program, Result · 1 declaration(s)
MATLAB's fmincon -- minimise a function HANDLE subject to linear constraints, bounds and NONLINEAR constraints (the toolbox board's T5.1).
class ICoreConstrainedNonlinearOptimization {
public:
// f(x1..xn) -> fx. False means the call itself failed (the caller's
// objective is a console function handle whose body can be in error, C2.8).
using Objective = std::function<bool(const std::vector<double>& x, double& fx)>;
// f(x1..xn) -> [fx, gradient], the two-output objective MATLAB reads
// under 'SpecifyObjectiveGradient'. `gradient` is null when only the
// value is wanted.
using ObjectiveWithGradient =
std::function<bool(const std::vector<double>& x, double& fx, std::vector<double>* gradient)>;
// nonlcon(x) -> [c, ceq], MATLAB's TWO-output constraint handle: `c <= 0`
// and `ceq == 0`. Either may be empty, and an empty one is written `[]`
// there and answers an empty vector here.
using Constraints = std::function<bool(const std::vector<double>& x, std::vector<double>& c,
std::vector<double>& ceq)>;
// The linear half of the problem, in MATLAB's own argument order. Every
// one may be empty, which is how MATLAB spells "no constraint of this
// kind" too -- `fmincon(fun, x0, [], [], [], [], lb, ub)`.
struct Program {
ICoreMatrix* inequality = nullptr; // A
ICoreMatrix* inequalityBound = nullptr; // b
ICoreMatrix* equality = nullptr; // Aeq
ICoreMatrix* equalityBound = nullptr; // beq
ICoreMatrix* lower = nullptr; // lb
ICoreMatrix* upper = nullptr; // ub
};
// What `fmincon` answers, in MATLAB's own vocabulary. `exitflag` carries
// MATLAB's integers unchanged: 1 first-order optimality met, 2 the step
// fell below StepTolerance at a feasible point, 5 the merit function could
// not be reduced further, 0 the iteration or evaluation budget was spent,
// -2 no feasible point was found. They are an answer about the problem,
// not an error in the call, which is why they are returned rather than
// reported as a failure.
struct Result {
std::vector<double> x; // the minimiser, in x0's own order
double fval = 0.0; // the objective there
int exitflag = 0;
int iterations = 0;
int funcCount = 0;
double constraintViolation = 0.0; // output.constrviolation
double stepsize = 0.0; // output.stepsize
double firstOrderOptimality = 0.0; // output.firstorderopt
std::vector<double> gradient; // GRAD, n entries
// LAMBDA, in MATLAB's own blocks. Each is empty when the problem
// carries no constraint of that kind.
std::vector<double> lambdaInequalityLinear; // lambda.ineqlin
std::vector<double> lambdaEqualityLinear; // lambda.eqlin
std::vector<double> lambdaLower; // lambda.lower
std::vector<double> lambdaUpper; // lambda.upper
std::vector<double> lambdaInequalityNonlinear;// lambda.ineqnonlin
std::vector<double> lambdaEqualityNonlinear; // lambda.eqnonlin
};
// fmincon(fun, x0, A, b, Aeq, beq, lb, ub, nonlcon, options) with the
// defaults MATLAB's own option set carries: OptimalityTolerance 1e-6,
// ConstraintTolerance 1e-6, StepTolerance 1e-10 (the 'sqp' default, which
// is NOT interior-point's 1e-10 by accident -- both are 1e-10 and
// 'active-set' alone uses 1e-6), MaxIterations 400,
// MaxFunctionEvaluations 100*n.
//
// `constraints` is null for a problem with no `nonlcon`. `options` is a
// resolved `optimoptions('fmincon', ...)` set (T5.11); null, or one whose
// fields are all unset, runs those defaults. `analytic` is the
// 'SpecifyObjectiveGradient' objective, used INSTEAD of the
// finite-difference gradient at every point -- which changes `funcCount`
// on both sides, as it does for `fminunc` (T5.2).
//
// ⚠ `whyNot` and the two option parameters are LAST for the reason
// `ICoreNonlinearOptimization` records: every existing caller passes
// `whyNot` positionally, so a parameter inserted before it would change
// the meaning of a call silently rather than failing to compile.
static bool minimum(const Objective& f, const std::vector<double>& x0,
const Program& program, const Constraints* constraints, Result& out,
std::string* whyNot = nullptr,
const ICoreOptimizationOptions::Values* options = nullptr,
const ObjectiveWithGradient* analytic = nullptr);
};
};ICoreConstrainedOptimization.h#
src/ICoreBlocks/ICoreMath/Optimization/ICoreConstrainedOptimization.h
ICoreConstrainedOptimization#
ICoreConstrainedOptimization.h:32 · class · nested Multipliers, Statistics · 0 declaration(s)
The constrained solvers of the toolbox board's T5.6 to T5.9 -- MATLAB's linprog, quadprog, lsqlin and intlinprog -- over one feasibility engine.
class ICoreConstrainedOptimization {
public:
// What a solver decided, in MATLAB's own vocabulary: `exitflag` 1 is a
// solution, -2 infeasible, -3 unbounded. They are returned rather than
// reported as failures because MATLAB returns them too -- an infeasible
// LP is an answer about the problem, not an error in the call.
enum class Outcome { Solved, Infeasible, Unbounded, Failed };
// The Lagrange multipliers at a quadratic programme's solution, in
// MATLAB's own blocks (`lambda.ineqlin`, `lambda.eqlin`, `lambda.lower`,
// `lambda.upper`). Every one is non-negative except `equality`, whose
// sign is the constraint's own.
//
// ⚠ THEY ARE RECOVERED AFTER THE FACT AND NOT READ OUT OF THE ITERATION.
// The active set carries multipliers for its WORKING SET only, and it
// leaves the loop by three different doors; recomputing them once at the
// answer -- the least-squares solution of `C' z = -(H x + f)` over the
// active constraints -- is one place instead of three, and it is the same
// number. A constraint that is not active gets zero, which is what
// complementarity says it must be.
struct Multipliers {
std::vector<double> inequality; // one per row of A
std::vector<double> equality; // one per row of Aeq
std::vector<double> lower; // one per variable
std::vector<double> upper; // one per variable
};
// ==== mb337 T5.14 begin ====
// What MATLAB's `output` record reports for the four programmes here (the
// toolbox board's T5.14). Every field was already being computed and
// thrown away: the simplex counts its own pivots, the active set its own
// iterations, and the branch and bound already knows its node and
// incumbent counts and its gaps -- they simply had nowhere to go.
//
// ⚠ `iterations` IS THIS CONSOLE'S OWN COUNT AND IS NOT MATLAB'S, and the
// difference is the algorithm rather than a defect: a two-phase simplex
// does not take the same number of steps as MATLAB's dual-simplex-highs,
// and a primal active set does not take an interior point's. The number is
// honest about what ran HERE, which is what `output` means; it compares as
// a PROPERTY and never digit-for-digit.
//
// `constraintViolation` and `firstOrderOptimality` are the two fields that
// do NOT depend on the path taken -- both are read off the ANSWER -- so
// they are the digit-comparable half of the record.
struct Statistics {
int iterations = 0; // output.iterations
double constraintViolation = 0.0; // output.constrviolation, at the answer
double firstOrderOptimality = 0.0; // output.firstorderopt, at the answer
int nodes = 0; // intlinprog: output.numnodes
int feasiblePoints = 0; // intlinprog: output.numfeaspoints
double absoluteGap = 0.0; // intlinprog: output.absolutegap
double relativeGap = 0.0; // intlinprog: output.relativegap
};
// ==== mb337 T5.14 end ====
// min f'x subject to A x <= b, Aeq x = beq, lb <= x <= ub.
//
// Any of A/Aeq/lb/ub may be empty. A bound of +/-infinity is absent, which
// is how MATLAB spells "no bound" too. The default when `lb` is empty is
// FREE variables -- linprog does not assume x >= 0, and a solver that did
// would answer a different (and often smaller) optimum without saying so.
//
// Two-phase simplex over the standard form, with Bland's rule on the
// second phase so a degenerate vertex cannot cycle.
static Outcome linearProgram(const ICoreMatrix& objective, const ICoreMatrix& inequality,
const ICoreMatrix& inequalityBound, const ICoreMatrix& equality,
const ICoreMatrix& equalityBound, const ICoreMatrix& lower,
const ICoreMatrix& upper, ICoreMatrix& solution, double& value,
std::string* whyNot = nullptr,
Multipliers* multipliers = nullptr,
Statistics* statistics = nullptr);
// min 0.5 x'H x + f'x subject to the same constraint set.
//
// H is symmetrised as `(H + H')/2`, which is what MATLAB does with a
// warning, and a non-convex H (any negative eigenvalue) is refused with
// MATLAB's own message rather than solved to a stationary point that is
// not a minimum.
//
// Primal active-set: start from a feasible vertex found by the simplex
// above, solve the equality-constrained problem on the working set through
// its KKT system, add the first constraint the step would cross, and drop
// a constraint whose multiplier has gone negative. It terminates on a
// finite set of working sets, and where it stops IS the minimiser.
static Outcome quadraticProgram(const ICoreMatrix& quadratic, const ICoreMatrix& objective,
const ICoreMatrix& inequality,
const ICoreMatrix& inequalityBound,
const ICoreMatrix& equality, const ICoreMatrix& equalityBound,
const ICoreMatrix& lower, const ICoreMatrix& upper,
ICoreMatrix& solution, double& value,
std::string* whyNot = nullptr,
Multipliers* multipliers = nullptr,
Statistics* statistics = nullptr);
// min 0.5 * ||C x - d||^2 over the same constraint set -- MATLAB's
// `lsqlin` (T5.6). It is `quadraticProgram` above on H = C'C and
// f = -C'd, which is exactly the reduction `lsqlin.m` makes for its
// 'active-set' algorithm, so nothing new is solved here: the row is the
// reduction and the two extra outputs.
//
// `residualNorm` is MATLAB's `resnorm` and it is NOT the halved objective:
// resnorm = ||C x - d||^2 while the quadratic being minimised is half of
// it. `residual` is `C x - d`, the signed vector, not its norm.
//
// ⚠ WITH NO CONSTRAINTS THIS DOES NOT GO THROUGH THE QP AT ALL. The
// unconstrained minimiser is the ordinary least-squares solution and is
// taken from a rank-revealing QR of C directly, because forming C'C
// SQUARES the condition number -- on a problem MATLAB answers to full
// precision, a normal-equations route would answer half the digits and
// the parity band would have to be widened to hide it. Once a constraint
// is present there is no such route and C'C is formed, which is what
// MATLAB does there too.
static Outcome constrainedLeastSquares(const ICoreMatrix& design, const ICoreMatrix& target,
const ICoreMatrix& inequality,
const ICoreMatrix& inequalityBound,
const ICoreMatrix& equality,
const ICoreMatrix& equalityBound,
const ICoreMatrix& lower, const ICoreMatrix& upper,
ICoreMatrix& solution, double& residualNorm,
ICoreMatrix& residual, std::string* whyNot = nullptr,
Multipliers* multipliers = nullptr,
Statistics* statistics = nullptr);
// min f'x over the same constraint set, with the variables named by
// `integerIndices` (MATLAB's `intcon`, ONE-BASED) required to be integers
// -- MATLAB's `intlinprog` (T5.9).
//
// Branch and bound over `linearProgram` above: solve the relaxation, pick
// a variable that came back fractional, and split the problem in two by
// TIGHTENING THAT VARIABLE'S BOUNDS -- x_i <= floor(v) on one side and
// x_i >= ceil(v) on the other. A node whose relaxation is already worse
// than the best integer point found so far cannot contain a better one
// and is dropped unopened; that bound is the whole method and it is why
// the search terminates on problems whose tree is astronomically large.
//
// The optimum VALUE is again the problem's and not the method's, which is
// what lets this be compared to MATLAB's own branch and bound (the same
// argument the class comment makes for the LP). x itself is unique only
// when the optimal integer point is.
static Outcome integerLinearProgram(const ICoreMatrix& objective,
const ICoreMatrix& integerIndices,
const ICoreMatrix& inequality,
const ICoreMatrix& inequalityBound,
const ICoreMatrix& equality,
const ICoreMatrix& equalityBound, const ICoreMatrix& lower,
const ICoreMatrix& upper, ICoreMatrix& solution,
double& value, std::string* whyNot = nullptr,
Statistics* statistics = nullptr);
};
};ICoreContinuousTimeIdentification.h#
src/ICoreBlocks/ICoreMath/Optimization/ICoreContinuousTimeIdentification.h
ICoreContinuousTimeIdentification#
ICoreContinuousTimeIdentification.h:56 · class · nested Model · 0 declaration(s)
MATLAB's tfest(data, np, nz[, ioDelay]) with no 'Ts' -- a CONTINUOUS-time transfer function fitted to a DISCRETE-time record (toolbox row T6.2).
class ICoreContinuousTimeIdentification {
public:
struct Model {
// Descending powers of s. `num` is nz + 1 long, `den` np + 1 and monic
// -- the shapes `m.Numerator` and `m.Denominator` read on the idtf.
std::vector<double> num;
std::vector<double> den;
// The io delay as it was asked for, in seconds, echoed back the way
// `m.IODelay` does. The samples are shifted by it before anything else
// happens and it is never estimated (T6.2: MATLAB disables iterative
// delay estimation for a continuous model over discrete data).
double ioDelay = 0.0;
// V at the answer -- e'e/N, `Report.Fit.LossFcn`.
double loss = 0.0;
size_t iterations = 0;
size_t functionCount = 0;
std::string whyStop;
};
// `input` and `output` are the record's samples, `sampleTime` its own Ts
// (> 0). `np` poles, `nz` zeros, `np >= nz` -- MATLAB refuses `nz > np` for
// a continuous estimate in as many words ("the number of poles must be
// greater than or equal to the number of zeros") and so does this.
//
// `ioDelay` is in SECONDS and must be a non-negative whole number of
// samples; a fractional one is refused with the reason. Measured on
// R2026a: an integer io delay answers EXACTLY the estimate over the record
// whose input has been shifted by that many samples with zeros prepended,
// bit for bit over 36 order sets, which is what this does. A fractional one
// does not reduce to a shift -- it becomes a fractional input delay inside
// the sampled simulation -- and is the one input format of this row that is
// refused rather than answered.
static bool estimate(const ICoreMatrix& input, const ICoreMatrix& output,
const double& sampleTime, const size_t& np, const size_t& nz,
const double& ioDelay, Model& out, std::string* whyNot = nullptr);
// The initialisation on its own -- `inival_time.m`'s `'iv'` mode. Exposed
// for the same reason the polynomial one is: an initialisation that is one
// filter out lands the search in a different minimum, and from the outside
// that failure looks exactly like a bad search.
static bool initialize(const ICoreMatrix& input, const ICoreMatrix& output,
const double& sampleTime, const size_t& np, const size_t& nz,
const double& ioDelay, Model& out, std::string* whyNot = nullptr);
};
};ICoreExcitationSignals.h#
src/ICoreBlocks/ICoreMath/Optimization/ICoreExcitationSignals.h
ICoreExcitationSignals#
ICoreExcitationSignals.h:47 · class · 1 declaration(s)
MATLAB's idinput (toolbox row T6.10) -- the four signals the System Identification Toolbox generates to EXCITE a plant before anything is identified from it.
class ICoreExcitationSignals {
public:
// MATLAB's four type words. `rs` is its own alias for `rgs`.
enum class Excitation { Prbs, RandomBinary, RandomGaussian, Sine };
static bool excitationOf(const std::string& word, Excitation& type);
// `[u, freqs] = idinput([P nu M], type, band, levels, sinedata)`.
//
// `period` is P, `inputs` is nu, `periods` is M (the whole signal repeated).
// `band` is the pair `[low high]` in units of the Nyquist frequency,
// `levels` the pair `[min max]`, and `sineData` MATLAB's
// `[sinusoids trials odd]` (defaulting to 10, 10 and 1).
//
// `frequencies` comes back empty for every type but `sine`, which is what
// MATLAB's second output does.
static bool generate(const Excitation& type, const size_t& period, const size_t& inputs,
const size_t& periods, const ICoreMatrix& band,
const ICoreMatrix& levels, const ICoreMatrix& sineData,
ICoreMatrix& signal, ICoreMatrix& frequencies,
std::string* whyNot = nullptr);
};
};ICoreIdMinimizer.h#
src/ICoreBlocks/ICoreMath/Optimization/ICoreIdMinimizer.h
ICoreIdMinimizer#
ICoreIdMinimizer.h:29 · class · nested Result · 0 declaration(s)
idminimizer.m with SearchMethod 'auto' -- the search EVERY iterative estimator of the System Identification Toolbox stops in, and therefore the half of those estimators' answers that is not the c...
class ICoreIdMinimizer {
public:
// What the search did, for the caller's own record and for a refusal that
// has to say why nothing better came back. `whyStop` is MATLAB's sentence.
struct Result {
double loss = 0.0;
size_t iterations = 0;
size_t functionCount = 0;
std::string whyStop;
};
// The criterion at a point. `loss` is MATLAB's `resnorm` -- e'e/N for the
// SISO det-criterion. When `jacobian` and `residual` are both non-null the
// criterion also answers the QR factors the step is taken from: `jacobian`
// is MATLAB's R and `residual` its `-Re`, i.e. `getErrorAndJacobian`'s
// final `e = -e` already applied. Answering false means the point is
// infeasible (an unstable model, a divergent filter) and the line search
// reads that as an infinite loss and bisects back.
using Criterion = std::function<bool(const std::vector<double>& x, double& loss,
ICoreMatrix* jacobian, ICoreMatrix* residual)>;
// `samples` is `sum(ut.DataSize)`, the N the expected-improvement test
// divides by.
//
// ⚠ THE NORMALISATION FACTOR `F` IS THE LOSS FOR BOTH CALLERS, and it gets
// there by two different routes -- which is why it is not a parameter.
// `idminimizer.m` sets `F = resnorm` when `isPoly || isTraceCrit`.
// `isPoly` is what makes it so for an idpoly (T6.4). For the structured
// idss of T6.2 `isPoly` is false, but `isTraceCrit` is
// `hasOutputWeight(option) && isnumeric(option.OutputWeight)` and
// `tfestOptions`' OutputWeight is `[]`, which IS numeric -- so it is true
// there too. Measured: with F = 1 the continuous `tfest` stops one to
// three iterations early on five of nine order sets; with F = resnorm
// every iteration and function count matches R2026a exactly.
static bool minimize(const Criterion& criterion, const size_t& samples,
std::vector<double>& x, Result& out, std::string* whyNot = nullptr);
};
};ICoreMinimaxOptimization.h#
src/ICoreBlocks/ICoreMath/Optimization/ICoreMinimaxOptimization.h
ICoreMinimaxOptimization#
ICoreMinimaxOptimization.h:53 · class · nested Result · 1 declaration(s)
The two EPIGRAPH solvers of the Optimization Toolbox -- MATLAB's fminimax and fgoalattain (toolbox board row T5.10).
class ICoreMinimaxOptimization {
public:
// F(x) -> the vector of objectives. False means the call itself failed --
// the caller's F is a console handle whose body can be in error (C2.8).
using VectorObjective =
std::function<bool(const std::vector<double>& x, std::vector<double>& F)>;
// What either solver answers. `attainment` is `fminimax`'s `maxfval` (the
// largest objective at the solution) and `fgoalattain`'s `attainfactor`
// (how far the goals were over- or under-attained); they are the same
// number in the same place -- the epigraph variable -- which is why one
// struct serves both.
struct Result {
std::vector<double> x;
std::vector<double> values; // F at the solution, in F's own order
double attainment = 0.0;
int iterations = 0;
int functionCalls = 0;
};
// MATLAB's exitflag vocabulary is NOT reproduced, and the reason is that
// its integers say which stopping TEST fired inside its SQP (1 first-order
// optimality, 4 the search direction was small, 5 the directional
// derivative was) -- a fact about MATLAB's iteration and not about the
// problem. These four are about the PROBLEM, the way `linprog`'s -2 and -3
// are, and the row refuses the exitflag output rather than answering a
// number under MATLAB's name that means something else. `Failed` also
// carries the nonlinear refusal, which is a fact about this console.
enum class Outcome { Solved, Infeasible, Unbounded, Failed };
// min max_i F_i(x) subject to A x <= b, Aeq x = beq, lb <= x <= ub.
// Any constraint matrix may be empty. `x0` is where F is read and where
// its Jacobian is measured; for an affine F the answer does not depend on
// it, which is worth knowing when a result looks surprising.
static Outcome minimax(const VectorObjective& objective, const std::vector<double>& start,
const ICoreMatrix& inequality, const ICoreMatrix& inequalityBound,
const ICoreMatrix& equality, const ICoreMatrix& equalityBound,
const ICoreMatrix& lower, const ICoreMatrix& upper, Result& result,
std::string* whyNot = nullptr);
// min gamma subject to F_i(x) - w_i * gamma <= goal_i, and the same
// constraint set.
//
// ⚠ A ZERO WEIGHT IS NOT "IGNORE THIS OBJECTIVE": it makes the goal a HARD
// constraint (F_i(x) <= goal_i, with no gamma to relax it), which is what
// MATLAB's own documentation means by a weight of zero and is the one
// thing about this solver that a reader guesses wrong.
static Outcome goalAttainment(const VectorObjective& objective,
const std::vector<double>& start, const ICoreMatrix& goal,
const ICoreMatrix& weight, const ICoreMatrix& inequality,
const ICoreMatrix& inequalityBound, const ICoreMatrix& equality,
const ICoreMatrix& equalityBound, const ICoreMatrix& lower,
const ICoreMatrix& upper, Result& result,
std::string* whyNot = nullptr);
};
};ICoreModelPrediction.h#
src/ICoreBlocks/ICoreMath/Optimization/ICoreModelPrediction.h
ICoreModelPrediction#
ICoreModelPrediction.h:67 · class · 0 declaration(s)
What an identified model DOES once it has been fitted -- MATLAB's sim, predict, forecast and resid (toolbox row T6.9).
class ICoreModelPrediction {
public:
// `sim(sys, u)` -- the model's forced response from a ZERO state, which is
// `filter(B, A, u)`. The console's own `simtf` computes the same thing and
// stays as a marked extra (X7).
//
// `input` is a vector written either way round; `output` comes back the
// same length and always a COLUMN, which is MATLAB's shape (see below).
static bool simulate(const ICoreTransferFunction& model, const ICoreMatrix& input,
ICoreMatrix& output, std::string* whyNot = nullptr);
// `predict(sys, u, y, Ts, k)` -- the console's numeric form of MATLAB's
// `predict(sys, iddata(y, u, Ts), k)` with `'InitialCondition', 'z'`.
//
// `horizon` is k. It is a whole number at least 1, or infinity for the
// pure simulation MATLAB spells `Inf`; anything else is refused by name
// rather than rounded, because a fractional horizon has no meaning and
// silently taking its floor is the one answer that cannot be right.
static bool predictHorizon(const ICoreTransferFunction& model, const ICoreMatrix& input,
const ICoreMatrix& output, const double& horizon,
ICoreMatrix& predicted, std::string* whyNot = nullptr);
// `forecast(sys, past, K)` -- K samples BEYOND the record, continuing the
// model's own recursion with the noise set to its mean of zero.
//
// `futureInput` may be empty, which is MATLAB's own default of a zero
// future input; when it is not empty it holds exactly `steps` samples and
// is refused otherwise. `forecast` answers only the K new samples, not the
// record with them appended -- MATLAB's shape.
static bool forecast(const ICoreTransferFunction& model, const ICoreMatrix& input,
const ICoreMatrix& output, const size_t& steps,
const ICoreMatrix& futureInput, ICoreMatrix& forecastOut,
std::string* whyNot = nullptr);
// `resid(u, y, Ts, sys)` -- the one-step prediction errors `e = y - yhat`,
// which is what MATLAB's first output is (verified digit for digit against
// `y - predict(m, d, 1)` on this row's fixture).
//
// ⚠ MATLAB'S SECOND OUTPUT IS NOT AVAILABLE HERE AND IS NOT APPROXIMATED.
// `[e, r] = resid(d, m)` answers an `r` measured at 26 x 2 x 2 on R2026a
// -- the residual autocorrelation and the residual/input cross-correlation
// stacked in a THIRD dimension -- and an N-D result is off this board by
// Z10. The console answers `e` and refuses the second output by name; it
// does not flatten the array to a shape MATLAB has no spelling for, which
// is C0.9's rule that a row states why rather than inventing a value.
static bool residuals(const ICoreTransferFunction& model, const ICoreMatrix& input,
const ICoreMatrix& output, ICoreMatrix& errors,
std::string* whyNot = nullptr);
};
};ICoreModelValidation.h#
src/ICoreBlocks/ICoreMath/Optimization/ICoreModelValidation.h
ICoreModelValidation#
ICoreModelValidation.h:25 · class · 1 declaration(s)
compare and goodnessOfFit -- toolbox row T6.8, the pair every other System Identification row is scored by.
class ICoreModelValidation {
public:
// MATLAB's three measures. `Mse` answers one number for the whole matrix
// (`trace(e' e) / Ns`, summed over every channel and divided by the sample
// count only); the other two answer one number per COLUMN.
enum class FitMeasure { Nrmse, Nmse, Mse };
static bool measureFromName(const std::string& name, FitMeasure& out);
// `goodnessOfFit(x, xref, measure)`. `value` comes back 1 x 1 for `Mse`
// and 1 x N for the other two, MATLAB's own shapes.
//
// The degenerate case is MATLAB's and is not a guard against dividing by
// zero: when a channel's reference is CONSTANT *and* the test data equals
// it -- both compared against `eps(||xref||_inf * Ns)` -- the ratio is
// taken to be 0 rather than 0/0. A constant reference that the model does
// NOT match still divides by zero and answers Inf, which is the honest
// answer and what MATLAB gives.
static bool goodnessOfFit(const ICoreMatrix& test, const ICoreMatrix& reference,
const FitMeasure& measure, ICoreMatrix& value,
std::string* whyNot = nullptr);
// `[yh, fit, x0] = compare(u, y, Ts, sys)` -- the console's numeric form of
// MATLAB's `compare(iddata(y, u, Ts), sys)`.
//
// ⚠ THE INITIAL STATE IS ESTIMATED, NOT ZERO, and a simulation from a zero
// state is a different fit. MATLAB's default `InitialCondition` is
// `'auto'`, which for a model with no noise component resolves to
// `'estimate'`: `x0` is the least-squares minimiser of
// `|| y - forced - Phi x0 ||` over the free-response basis
// `Phi(t, :) = C A^t`, exactly `findstates`' own rule. Measured on this
// row's twenty-sample fixture: 35.209 with the estimated state against
// 34.312 from zero.
//
// `initialState` is in the basis of `ICoreStateSpace::fromTransferFunction`
// when `sys` is a transfer function, so it is the one output of the three
// that is a property of the REALIZATION rather than of the model -- `yh`
// and `fit` are the same in every basis.
static bool compare(const ICoreStateSpace& system, const ICoreMatrix& input,
const ICoreMatrix& output, ICoreMatrix& simulated, double& fitPercent,
ICoreMatrix& initialState, std::string* whyNot = nullptr);
};
};ICoreNonlinearEquations.h#
src/ICoreBlocks/ICoreMath/Optimization/ICoreNonlinearEquations.h
ICoreNonlinearEquations#
ICoreNonlinearEquations.h:38 · class · nested Result · 1 declaration(s)
MATLAB's fsolve -- a system of nonlinear equations solved by the trust-region dogleg method (the toolbox board's T5.3).
class ICoreNonlinearEquations {
public:
// F(x1..xn) -> [f1..fm]. False means the call itself failed.
using VectorObjective =
std::function<bool(const std::vector<double>& x, std::vector<double>& F)>;
struct Result {
std::vector<double> x; // where the iteration stopped
std::vector<double> fvec; // F there
std::vector<double> jacobian; // m*n column-major, the finite-difference Jacobian there
int exitflag = 0; // MATLAB's integers: 1, 2, 3, -3, 0
int iterations = 0;
int funcCount = 0;
double firstOrderOptimality = 0.0; // norm(JAC'*F, inf)
};
// fsolve(F, x0) with MATLAB's default 'trust-region-dogleg'.
// `options` is a resolved `optimoptions(...)` set (T5.11); null, or one
// whose fields are all unset, runs MATLAB's own defaults, so the two are
// the same iteration. It is the LAST parameter because `whyNot` already
// had a default and every caller passes it positionally -- a parameter
// inserted before it would change existing calls silently.
static bool solve(const VectorObjective& F, const std::vector<double>& x0,
Result& out, std::string* whyNot = nullptr,
const ICoreOptimizationOptions::Values* options = nullptr);
// The forward-difference Jacobian, column by column, by the step rule
// above. `evaluations` is what it cost.
static bool finiteDifferenceJacobian(const VectorObjective& F, const std::vector<double>& x,
const std::vector<double>& Fx, size_t rows,
std::vector<double>& jacobian, int& evaluations,
std::string* whyNot = nullptr);
};
};ICoreNonlinearLeastSquares.h#
src/ICoreBlocks/ICoreMath/Optimization/ICoreNonlinearLeastSquares.h
ICoreNonlinearLeastSquares#
ICoreNonlinearLeastSquares.h:40 · class · nested Result · 1 declaration(s)
MATLAB's lsqnonlin and lsqcurvefit -- the nonlinear least-squares half of the Optimization Toolbox (the toolbox board's T5.4 and T5.5), by the DEFAULT algorithm, which is `trust-region-reflecti...
class ICoreNonlinearLeastSquares {
public:
// F(x1..xn) -> [f1..fm], the residual VECTOR (not its sum of squares --
// the trap for anyone reading the row's title too quickly). False means
// the call itself failed.
using VectorObjective =
std::function<bool(const std::vector<double>& x, std::vector<double>& F)>;
struct Result {
std::vector<double> x; // where the iteration stopped
std::vector<double> residual; // F there
double resnorm = 0.0; // sum(residual.^2), MATLAB's own second output
std::vector<double> jacobian; // m*n column-major
int exitflag = 0; // MATLAB's integers: 1, 2, 3, 0
int iterations = 0;
int funcCount = 0;
double firstOrderOptimality = 0.0;
};
// lsqnonlin(F, x0) with MATLAB's default trust-region-reflective and its
// default options (TolFun 1e-6, TolX 1e-6, MaxIter 400, MaxFunEvals 100*n).
// MATLAB requires at least as many equations as unknowns for this
// algorithm and says so by name; so does this.
// `options` is a resolved `optimoptions(...)` set (T5.11); null, or one
// whose fields are all unset, runs MATLAB's own defaults, so the two are
// the same iteration. It is the LAST parameter because `whyNot` already
// had a default and every caller passes it positionally -- a parameter
// inserted before it would change existing calls silently.
static bool solve(const VectorObjective& F, const std::vector<double>& x0,
Result& out, std::string* whyNot = nullptr,
const ICoreOptimizationOptions::Values* options = nullptr);
};
};ICoreNonlinearOptimization.h#
src/ICoreBlocks/ICoreMath/Optimization/ICoreNonlinearOptimization.h
ICoreNonlinearOptimization#
ICoreNonlinearOptimization.h:52 · class · nested Result · 1 declaration(s)
The Optimization Toolbox solvers that take a function HANDLE and search in n dimensions without constraints: MATLAB's fminunc and checkGradients (the toolbox board's T5.2 and T5.15).
class ICoreNonlinearOptimization {
public:
// f(x1..xn) -> fx. False means the call itself failed.
using ObjectiveN = std::function<bool(const std::vector<double>& x, double& fx)>;
// f(x1..xn) -> [fx, gradient], the two-output objective MATLAB reads under
// 'SpecifyObjectiveGradient' and the only shape `checkGradients` accepts.
// `gradient` is null when only the VALUE is wanted, which is how the
// finite-difference half asks -- a caller whose objective is expensive in
// its second output should skip it there, and MATLAB's own
// finitedifferences layer asks the same way. False means the call failed.
using ObjectiveWithGradient =
std::function<bool(const std::vector<double>& x, double& fx, std::vector<double>* gradient)>;
// What `fminunc` answers, in MATLAB's own vocabulary. `exitflag` carries
// MATLAB's integers unchanged -- 1 gradient below tolerance, 2 step below
// tolerance, 5 the line search could not reduce the objective further,
// 0 iteration or evaluation budget spent, -3 unbounded below the objective
// limit -- because they are an answer about the problem, not an error in
// the call.
struct Result {
std::vector<double> x; // the minimiser, in x0's own order
double fval = 0.0; // the objective there
std::vector<double> gradient; // n entries, the finite-difference gradient there
std::vector<double> hessian; // n*n column-major; empty unless asked for
int exitflag = 0;
int iterations = 0;
int funcCount = 0;
double stepsize = 0.0; // norm(deltaX, 2) of the last step
double lineSearchStepLength = 0.0; // output.lssteplength
double firstOrderOptimality = 0.0; // norm(grad, inf)
};
// fminunc(f, x0) with MATLAB's 'quasi-newton' algorithm and its default
// options (TolX 1e-6, TolFun 1e-6, MaxIter 400, MaxFunEvals 100*n,
// ObjectiveLimit -1e20). `wantHessian` is the six-output form: MATLAB
// computes the finite-difference Hessian only when it is asked for, with
// an eps^(1/4) step where the gradient uses sqrt(eps), and this does the
// same -- so asking for it changes `funcCount`, on both sides.
// `options` is an `optimoptions('fminunc', ...)` set that has already been
// resolved (T5.11); null, or a set whose fields are all unset, runs
// MATLAB's own defaults. `analytic` is the 'SpecifyObjectiveGradient'
// objective -- the two-output handle MATLAB reads when that option is on;
// it is used INSTEAD of the finite-difference gradient at every point,
// which is why turning the option on changes `funcCount` from
// (1 + n)*(1 + iterations) to 1 + iterations + the line search's trials
// (measured on R2026a: 66 against 22 on this row's own quartic) and moves
// the answer slightly -- an exact gradient is not the forward-difference
// one.
//
// ⚠ BOTH ARE THE LAST PARAMETERS RATHER THAN THE NATURAL ONES. `whyNot`
// has had a default since this function was written and every caller
// passes it positionally, so a parameter inserted before it would change
// the meaning of every existing call silently rather than failing to
// compile.
static bool unconstrainedMinimum(const ObjectiveN& f, const std::vector<double>& x0,
bool wantHessian, Result& out,
std::string* whyNot = nullptr,
const ICoreOptimizationOptions::Values* options = nullptr,
const ObjectiveWithGradient* analytic = nullptr);
// MATLAB's forward-difference gradient, the one every solver here uses:
// step h_i = sqrt(eps) * sign(x_i) * max(|x_i|, 1), with sign(0) = +1, and
// g_i = (f(x + h_i e_i) - f(x)) / h_i. `evaluations` is how many calls it
// cost, which the callers add to MATLAB's `funcCount`.
static bool finiteDifferenceGradient(const ObjectiveN& f, const std::vector<double>& x,
double fx, std::vector<double>& gradient,
int& evaluations, std::string* whyNot = nullptr);
// checkGradients(fun, x0) -- true when the gradient the objective returns
// agrees with the finite-difference one to `tolerance` (MATLAB's default
// 1e-6), relatively, at a point NEAR x0.
//
// ⚠ MATLAB perturbs x0 by `1e-3 * max(|x0|, 1) .* (1 - 2*rand(n, 1))` with
// its own RNG reset to seed 0, so the point tested is MATLAB's Mersenne
// Twister stream and not this one's (Z11). `maximumRelativeError` is
// therefore this console's number, not MATLAB's, and only the VERDICT
// crosses: a gradient that is right is right at either point, and one that
// is wrong by more than a rounding is wrong at either point. The caller
// says which of the two it publishes.
static bool gradientCheck(const ObjectiveWithGradient& f, const std::vector<double>& x0,
double tolerance, bool& valid, double& maximumRelativeError,
std::string* whyNot = nullptr);
};
};ICoreNonparametricIdentification.h#
src/ICoreBlocks/ICoreMath/Optimization/ICoreNonparametricIdentification.h
ICoreNonparametricIdentification#
ICoreNonparametricIdentification.h:33 · class · 1 declaration(s)
The NONPARAMETRIC estimators of the System Identification Toolbox -- MATLAB's covf, etfe and spa (toolbox board row T6.7).
class ICoreNonparametricIdentification {
public:
// ---- covf(z, M) --------------------------------------------------------
//
// The covariance function of the COLUMNS of `data`, biased, laid out the
// way MATLAB's own help block states it:
//
// R(i + (j - 1) * nz, k + 1) = (1 / N) * sum_t z_i(t) * z_j(t + k)
//
// so the result is nz^2 rows by M columns, k running 0 .. M-1, and every
// sum is divided by the FULL record length N rather than by the number of
// products in it (the biased estimate -- `xcorr(..., 'biased')`, which is
// what `@iddata/covf.m`'s own comment says this function is).
//
// ⚠ MATLAB REACHES THOSE NUMBERS THROUGH AN FFT and this reaches them by
// summing, for records small enough that its own `covf.m` would sum too
// (its `maxsize` trade-off). The two agree to 2.2e-16 on this row's
// record, measured before this file was written -- the sums are short and
// the FFT is exact for them. A row that needed the FFT path's exact
// rounding would have to say so; this one does not.
static bool covarianceFunction(const ICoreMatrix& data, const size_t& maxLag,
ICoreMatrix& out, std::string* whyNot = nullptr);
// ---- etfe(u, y, Ts, M, N) ---------------------------------------------
//
// The empirical transfer function estimate: the ratio of the output's
// Fourier transform to the input's, at `points` frequencies evenly spaced
// over (0, pi/Ts].
//
// `smoothing` is MATLAB's M -- 0 here means "not given", which is its `[]`
// and means no smoothing at all. A given M applies a Hamming lag window of
// about pi/M resolution, and MATLAB HALVES IT FIRST (`M = max(1, M/2)`,
// "this is to make better agreement with SPA" -- its comment), so
// `etfe(data, 8)` smooths with a 4-lag window and not an 8-lag one.
//
// `points` is MATLAB's N and defaults to 128. The transform length is
// `2 * ceil(Ncap / N) * N`, so a record shorter than the grid is ZERO
// PADDED to twice the grid and the estimate at 128 points from 40 samples
// is an interpolation of 40 samples' worth of information -- true of
// MATLAB's answer as much as of this one.
static bool empiricalTransferFunction(const ICoreMatrix& input, const ICoreMatrix& output,
const double& sampleTime, const size_t& smoothing,
const size_t& points,
std::vector<std::complex<double>>& response,
std::vector<double>& frequency,
std::string* whyNot = nullptr);
// ---- spa(u, y, Ts, M, w) ----------------------------------------------
//
// The Blackman-Tukey spectral estimate: covariance functions to `lagWindow`
// lags, weighted by a Hann lag window, transformed at each requested
// frequency, and divided --
//
// G(w) = PHI_yu(w) / PHI_u(w)
//
// -- with the noise spectrum PHI_v = Ts * |PHI_y - |PHI_yu|^2 / PHI_u|
// coming out of the same three transforms. `noiseSpectrum` is the
// `SpectrumData` of the `idfrd` MATLAB answers; it is returned here
// because the transforms behind it are already computed and dropping it
// would make the console's answer strictly poorer than MATLAB's.
//
// `lagWindow` 0 means "not given" -- see `defaultLagWindow`. An empty
// `frequencies` means MATLAB's own default grid, (1:128) * pi / 128 / Ts.
//
// ⚠ THE LAG WINDOW'S FIRST WEIGHT IS HALVED (`window(1) = window(1)/2` in
// `spa.m`) and the transform is then taken as `2 * real(...)`, which is
// the same thing as summing the negative lags -- but only because of that
// halving. An implementation that used the Hann weights as written and
// doubled would answer a spectrum too large by exactly R(0) everywhere.
static bool spectralAnalysis(const ICoreMatrix& input, const ICoreMatrix& output,
const double& sampleTime, const size_t& lagWindow,
const std::vector<double>& frequencies,
std::vector<std::complex<double>>& response,
std::vector<double>& frequency,
std::vector<double>& noiseSpectrum,
std::string* whyNot = nullptr);
// MATLAB's default lag window for `spa` on a record of `samples` samples.
//
// ⚠ IT IS NOT THE RULE `spa.m` APPEARS TO STATE, and this is the one
// measurement of this row that no reading of the source gives you. The
// file holds two rules, one for time-domain data and one for frequency-
// domain data:
//
// Mdef = max(min(30, floor(mean(Ncaps) / 10)), 2) % time domain
// Mdef = max(min(30, floor(mean(2*Ncaps) / 10)), 2) % frequency
//
// and R2026a takes the SECOND ONE for ordinary time-domain records.
// Measured on this machine over nine record lengths, reading the window
// back out of `spa(iddata(y, u, 1)).Report.WindowSize`: 12 -> 2, 20 -> 4,
// 25 -> 5, 40 -> 8, 60 -> 12, 100 -> 20, 150 -> 30, 200 -> 30, 320 -> 30.
// Every one of those is `min(30, floor(N / 5))` and none of them is
// `floor(N / 10)`, which would have answered 4 where MATLAB answers 8 on
// this board's own 40-sample record -- a spectral estimate half as
// resolved, differing in the FIRST digit (0.634 against 0.701 at the
// lowest frequency), and passing every test a session would think to
// write about it.
static size_t defaultLagWindow(const size_t& samples);
};
};ICoreOptimizationOptions.h#
src/ICoreBlocks/ICoreMath/Optimization/ICoreOptimizationOptions.h
ICoreOptimizationOptions#
ICoreOptimizationOptions.h:38 · class · nested Setting, Values · 8 declaration(s)
optimoptions, optimset and optimget -- the option SET the Optimization Toolbox solvers take, and the reader every solver on this console goes through to find out what it was asked for (toolbo...
class ICoreOptimizationOptions {
public:
// The solvers an option set can be made for. `None` is a legacy
// `optimset(...)` set, which names no solver: MATLAB's optimset is
// solver-agnostic and its 55 names are one flat list, with whichever
// solver receives the set deciding which of them it reads.
enum class Solver { Fminunc, Fsolve, Lsqnonlin, Lsqcurvefit, Lsqlin, Linprog,
Quadprog, Intlinprog, Fminimax, Fgoalattain, Fmincon,
Fminsearch, Fminbnd, Fzero, Lsqnonneg, None };
// One `name, value` pair as it was typed. A value is a NUMBER or a WORD;
// MATLAB's logical options (`SpecifyObjectiveGradient`) are numbers, and
// its `'on'`/`'off'` legacy spellings are words.
struct Setting {
std::string name;
bool isWord = false;
std::string word;
double number = 0.0;
};
// What every solver here reads out of a set.
//
// ⚠ AN UNSET FIELD IS NaN AND NOT A DEFAULT, and that is the whole reason
// the struct is shaped this way: the four solvers do NOT share defaults.
// `fminunc`'s TolX is 1e-6 and `fminbnd`'s is 1e-4; `fzero`'s is `eps`;
// `fminsearch`'s evaluation budget is `200*n` and the gradient solvers'
// is `100*n`. Filling one set of defaults in here would silently retune
// three of them, so each solver reads a field only when it was ASKED for
// and otherwise runs the number its own transcription carries. That is
// also what makes "no option set" and "optimoptions(solver)" the same
// run, which a regress case asserts.
struct Values {
double stepTolerance = notSet();
double functionTolerance = notSet();
double optimalityTolerance = notSet();
double objectiveLimit = notSet();
double maxIterations = notSet();
double maxFunctionEvaluations = notSet();
double constraintTolerance = notSet();
bool specifyObjectiveGradient = false;
bool specifyConstraintGradient = false;
};
// The "unset" marker, and the reader every solver uses: `given(v, d)` is
// `v` when the caller asked for it and `d` -- that solver's own default --
// when it did not.
static double notSet();
static bool wasGiven(const double& value);
static double given(const double& value, const double& fallback);
// 'fminunc' or '@fminunc', in any case. False for a word that names no
// solver at all, which is a different error from naming one optimoptions
// does not support.
static bool solverOf(const std::string& word, Solver& out);
static std::string upperCaseName(const Solver& solver);
static std::string lowerCaseName(const Solver& solver);
// `optimoptions(solver, name, value, ...)`. `existing` is the text of a
// set being updated (`optimoptions(opts, ...)`) or empty for the
// constructor form; when it is not empty its solver wins and `solver` is
// ignored.
static bool build(const Solver& solver, const std::string& existing,
const std::vector<Setting>& settings, std::string& text,
std::string* whyNot = nullptr);
// `optimset(name, value, ...)` and `optimset(oldopts, name, value, ...)`.
// The names stay in their LEGACY spelling, because that is the spelling
// `optimget` answers them under and the only one the four core solvers
// have ever had.
static bool buildLegacy(const std::string& existing, const std::vector<Setting>& settings,
std::string& text, std::string* whyNot = nullptr);
// `optimset(oldopts, newopts)` -- MATLAB's merge, where every option
// newopts SETS wins and every one it leaves unset keeps oldopts' value.
static bool mergeLegacy(const std::string& older, const std::string& newer,
std::string& text, std::string* whyNot = nullptr);
// `optimset(solver)` -- the option set that solver runs on, with its own
// defaults filled in. ⚠ TWO OF THEM ARE STRINGS: MATLAB answers
// `MaxIter = '200*numberofvariables'` for fminsearch, because the number
// cannot be written down until n is known -- and `resolve` reads a
// word-valued budget as "the solver's own rule" rather than as a number.
static bool defaultsFor(const Solver& solver, std::string& text,
std::string* whyNot = nullptr);
// `optimget(opts, name)`. `found` is false for a name the set does not
// carry -- MATLAB answers `[]` there, and the caller supplies its own
// default for the three-argument form. A name that is not an option at all
// is an error and answers false.
//
// ⚠ AN optimoptions SET IS REFUSED HERE, and that is MATLAB's own rule
// rather than a gap: `optimget(optimoptions('fminunc', ...), 'TolX')`
// answers "First argument must be an options structure created with
// OPTIMSET." there. The symmetry is inviting and it is not there -- this
// console answered the value until the parity suite said otherwise. A
// modern set is read by handing it to its SOLVER; MATLAB's other reader
// is property access, which C0.3 leaves this console without.
static bool read(const std::string& text, const std::string& named, bool& found,
Setting& value, std::string* whyNot = nullptr);
// Is this string an option set rather than an ordinary one? Used by the
// solvers to tell `fminunc(f, x0, opts)` from a string argument that is
// something else entirely. The DISPLAY of a set is this same text, which
// it already is as a `Kind::Str` -- there is no separate print form, for
// the reason the header comment gives.
static bool isOptionSet(const std::string& text);
// Read a set for `solver`. Answers false with the reason when the set
// names an option this console's algorithm cannot honour, or an option
// belonging to a DIFFERENT solver -- MATLAB refuses that too, and silently
// reading an fminunc set inside fsolve would honour the wrong defaults.
static bool resolve(const std::string& text, const Solver& solver, Values& out,
std::string* whyNot = nullptr);
// The algorithm this console runs for `solver`, in MATLAB's own spelling,
// or "" for a solver whose method is this console's own rather than one of
// MATLAB's named ones. `Algorithm` is honoured when it names this, and
// refused with this sentence otherwise.
static std::string implementedAlgorithm(const Solver& solver);
};
};ICoreParametricIdentification.h#
src/ICoreBlocks/ICoreMath/Optimization/ICoreParametricIdentification.h
ICoreParametricIdentification#
ICoreParametricIdentification.h:52 · class · 3 declaration(s)
The linear-regression estimators of the System Identification Toolbox -- MATLAB's arx, ar and ivar (toolbox row T6.3) -- and the order-search loops built on top of them, arxstruc, `selstruc...
class ICoreParametricIdentification {
public:
// MATLAB's five `ar` approaches. `ForwardBackward` is the DEFAULT -- not
// least squares, and not Burg, which is what a reader of the name usually
// expects -- and it fits the forward and the time-reversed equations
// together in one solve.
enum class ArApproach { ForwardBackward, LeastSquares, YuleWalker, Burg, GeometricLattice };
// How the data outside the measured interval is treated: `Pre` prepends n
// zeros, `Post` appends n, `PrePost` both. `None` is the default for every
// approach except `YuleWalker`, which ALWAYS uses `PrePost` -- MATLAB
// overrides a conflicting request with a warning -- and the padding is
// what makes Yule-Walker the autocorrelation method rather than the
// covariance one.
enum class ArWindow { None, Pre, Post, PrePost };
// `arx(u, y, Ts, [na nb nk])`. `denominator` comes back monic and na+1
// long, `numerator` nk+nb long with `nk` leading zeros -- MATLAB's own
// `m.A` and `m.B`.
//
// ⚠ THE FIRST SAMPLE USED IS NOT THE ONE THE MODEL ORDER SUGGESTS, and
// getting this wrong is not a rounding difference. The regression starts
// at t = max(na, nb - (nk == 0), 1) + 1 -- a rule that DOES NOT MENTION
// nk -- and any input regressor u(t - k) that reaches before the first
// sample is filled with ZERO rather than dropping its equation. The
// reading a careful implementer arrives at instead, "start at
// max(na, nb + nk - 1) + 1 so every regressor exists", is a different
// estimate: on the forty-sample record this row's corpus ships, at
// [na nb nk] = [1 2 3], it answers a1 = -0.443333 where MATLAB and this
// console answer -0.474552. `arx_time.m`'s `nmax` and its `I = jj > kl`
// guard are the two lines that say so.
//
// ⚠ The two rules AGREE whenever nk = 0 -- `nb - (nk == 0)` and
// `nb + nk - 1` are both `nb - 1` there -- so a case with no delay cannot
// tell them apart however carefully it is chosen. Only nk >= 1 with
// nb + nk - 1 above max(na, nb) separates them, which is why the corpus
// carries [1 2 3] and not just a sweep of nk.
//
// `lossFunction`, when asked for, is the mean squared residual of that
// same regression -- the number `arxstruc` scores an order set by.
static bool arx(const ICoreMatrix& input, const ICoreMatrix& output,
const size_t& na, const size_t& nb, const size_t& nk,
std::vector<double>& denominator, std::vector<double>& numerator,
double* lossFunction = nullptr, std::string* whyNot = nullptr);
// `ar(y, n[, approach[, window]])`. `denominator` comes back monic and
// n+1 long.
//
// `reflection` is MATLAB's second output and is filled only for `Burg` and
// `GeometricLattice`: row 1 the reflection coefficients with a leading 0,
// row 2 the loss function after each stage. It comes back EMPTY for the
// three approaches that are a solve rather than a recursion, which is what
// MATLAB leaves there too.
static bool ar(const ICoreMatrix& series, const size_t& order,
const ArApproach& approach, const ArWindow& window,
std::vector<double>& denominator, ICoreMatrix& reflection,
std::string* whyNot = nullptr);
// `ivar(y, na[, nc])` -- the instrumental-variable AR estimate. Four
// stages, and the third is the one nothing in the documentation mentions:
// a least-squares AR of order na + nc whitens the series, an ARX of that
// residual against the series gives a noise polynomial C, **C is reflected
// into the unit disc if it came out unstable** (`fstab`), and the
// instruments are the series filtered by 1/(C*C) after a shift of nc
// samples. The final solve is the IV normal equation
// `sum psi phi' theta = sum psi y`, which is NOT symmetric and is why the
// estimate is unbiased where `ar`'s least squares is not.
//
// `nc` defaults to `na` -- pass `na` for MATLAB's two-argument form.
static bool ivar(const ICoreMatrix& series, const size_t& na, const size_t& nc,
std::vector<double>& denominator, std::string* whyNot = nullptr);
// The spellings MATLAB takes on the call line, lower-cased by the caller:
// "fb", "ls", "yw", "burg", "gl" and "now", "prw", "pow", "ppw". Both
// answer false for a name they do not know, leaving `out` untouched.
static bool approachFromName(const std::string& name, ArApproach& out);
static bool windowFromName(const std::string& name, ArWindow& out);
// `arxstruc(ze, zv, NN)` -- row T6.11. One ARX fit per ROW of `orders`
// (each row `[na nb nk]`) on the estimation pair, scored on the validation
// pair, in MATLAB's own layout: `losses` is 4 x (rows + 1), the first row
// the loss per order set and the three below it that set's `[na; nb; nk]`.
// The EXTRA last column is a footer rather than a model:
// `losses(1, end)` is the number of validation samples and
// `losses(2, end)` the mean square of the validation output.
//
// ⚠ `arxstruc` DROPS ONE MORE SAMPLE THAN `arx` DOES. Its own start is
// `max(max(na) + 1, max(nb - (nk == 0)))` -- note the `+ 1`, which
// `arx_time.m` does not have -- so the loss it reports for an order set is
// not the residual of the model `arx` returns for that same set. Both
// rules are transcribed rather than unified, because unifying them would
// make one of the two rows disagree with MATLAB.
static bool arxstruc(const ICoreMatrix& estimationInput, const ICoreMatrix& estimationOutput,
const ICoreMatrix& validationInput, const ICoreMatrix& validationOutput,
const ICoreMatrix& orders, ICoreMatrix& losses,
std::string* whyNot = nullptr);
// `selstruc(V, c)` -- the numeric criteria only; the interactive
// `selstruc(V)` picker is Z2 and never reaches here. The penalty is
// `V * (1 + alpha * (na + nb) / N)` with `alpha` 0 for the plain minimum,
// 2 for `'aic'` and `log(N)` for `'mdl'` -- so this is the ORIGINAL
// Akaike form and not the `log` of anything, which is what the name leads
// a reader to write.
//
// `selectedOrder` is the winning `[na nb nk]` row; `criterionValues` is
// MATLAB's second output, the input matrix without its footer column and
// with the penalised losses replaced by their LOGARITHMS -- the one place
// a log does appear.
// `Alpha` takes the penalty weight from `alpha` -- MATLAB's numeric `c`,
// of which `selstruc(V, 0)` is the plain minimum -- while `Aic` and `Mdl`
// compute their own (2 and log(N)) from V's own sample count.
enum class SelectionCriterion { Alpha, Aic, Mdl };
static bool selstruc(const ICoreMatrix& losses, const SelectionCriterion& criterion,
const double& alpha, ICoreMatrix& selectedOrder,
ICoreMatrix& criterionValues, std::string* whyNot = nullptr);
static bool criterionFromName(const std::string& name, SelectionCriterion& out);
// `delayest(u, y, Ts[, na, nb, nkmin, nkmax])` -- the delay whose ARX fit
// has the smallest loss, scored by `arxstruc` with the SAME data as both
// the estimation and the validation set (which is what MATLAB passes) and
// chosen by `selstruc(V, 0)`. MATLAB's defaults are na = 2, nb = 2,
// nkmin = 0 and nkmax = nkmin + 40.
static bool delayest(const ICoreMatrix& input, const ICoreMatrix& output,
const size_t& na, const size_t& nb,
const int& nkMin, const int& nkMax,
int& delay, std::string* whyNot = nullptr);
};
};ICorePredictionErrorIdentification.h#
src/ICoreBlocks/ICoreMath/Optimization/ICorePredictionErrorIdentification.h
ICorePredictionErrorIdentification#
ICorePredictionErrorIdentification.h:59 · class · nested Orders, Model · 0 declaration(s)
The PREDICTION-ERROR estimators of the System Identification Toolbox -- MATLAB's armax, oe, bj and pem (toolbox row T6.4) -- over the general polynomial model A(q) y(t) = [B(q)/F(q)] u(t - ...
class ICorePredictionErrorIdentification {
public:
// Orders in MATLAB's own names. `nk` is the input delay in samples.
struct Orders {
size_t na = 0;
size_t nb = 0;
size_t nc = 0;
size_t nd = 0;
size_t nf = 0;
size_t nk = 0;
};
// MATLAB's A, B, C, D and F rows, in q^-1 and in MATLAB's own layout: A, C,
// D and F monic and one longer than their order, B exactly nk + nb long
// with nk leading zeros. A polynomial the structure does not have comes
// back as the scalar row {1} (and B as {0} when nb is zero), which is what
// `m.A`, `m.C`, `m.D` and `m.F` read on an idpoly that does not carry them.
struct Model {
std::vector<double> a{1.0};
std::vector<double> b;
std::vector<double> c{1.0};
std::vector<double> d{1.0};
std::vector<double> f{1.0};
// V at the answer -- e'e/N, the number MATLAB reports as Report.Fit.LossFcn.
double loss = 0.0;
// What the search did, for the row's note and for a refusal that needs
// to say why nothing better came back. `whyStop` is MATLAB's own
// sentence.
size_t iterations = 0;
size_t functionCount = 0;
std::string whyStop;
};
// `armax(u, y, Ts, [na nb nc nk])`, `oe(u, y, Ts, [nb nf nk])` and
// `bj(u, y, Ts, [nb nc nd nf nk])` all arrive here; which one it is is
// carried by the orders alone. Answers false with `whyNot` set when the
// record is too short for the orders, when a polynomial cannot be
// initialized, or when the search never left an infinite loss.
static bool estimate(const ICoreMatrix& input, const ICoreMatrix& output,
const Orders& orders, Model& out, std::string* whyNot = nullptr);
// MATLAB's `pem(data, init_sys)`: the same search from a GIVEN model
// instead of the initialisation above. `initial` supplies every polynomial
// and `orders` must agree with their lengths. MATLAB refuses a bare
// `pem(data)` without an initial model and so does the caller of this.
static bool refine(const ICoreMatrix& input, const ICoreMatrix& output,
const Model& initial, const Orders& orders,
Model& out, std::string* whyNot = nullptr);
// The initialisation on its own -- MATLAB's `inival_time`. Exposed because
// it is worth testing separately from the search: an initialisation that
// is one sample out lands the search in a different minimum, and the two
// failures look alike from the outside.
static bool initialize(const ICoreMatrix& input, const ICoreMatrix& output,
const Orders& orders, Model& out, std::string* whyNot = nullptr);
};
};ICoreRecursiveIdentification.h#
src/ICoreBlocks/ICoreMath/Optimization/ICoreRecursiveIdentification.h
ICoreRecursiveIdentification#
ICoreRecursiveIdentification.h:37 · class · 1 declaration(s)
The RECURSIVE estimators of the System Identification Toolbox that are still functions rather than objects -- MATLAB's rarx and rpem (toolbox board row T6.13).
class ICoreRecursiveIdentification {
public:
// How the gain is adapted at each step -- MATLAB's `adm` argument, and its
// own four spellings.
enum class Adaptation {
ForgettingFactor, // "ff": recursive least squares, `adg` is lambda
KalmanFilter, // "kf": `adg` is R1, the parameter drift covariance
NormalizedGradient, // "ng": `adg` is the gain, step divided by |phi|^2
UnnormalizedGradient// "ug": `adg` is the gain
};
static bool adaptationFromName(const std::string& name, Adaptation& out);
// rarx(z, [na nb nk], adm, adg).
//
// `data` is MATLAB's `z` and is OUTPUT FIRST: column 0 is y, the remaining
// columns are the inputs. (The estimators of T6.3 take `(u, y, ...)`; this
// one takes one matrix, in MATLAB's own order for it -- the argument-order
// trap T6.1's note names, and it is MATLAB's, not this console's.)
//
// `orders` is [na, nb..., nk...] with one nb and one nk per input, so its
// length is 2 * inputs + 1; a time series (no input column) takes [na]
// alone. Every nk must be one or more.
//
// `gain` is `adg`: a scalar for three of the four mechanisms, and for
// "kf" either a scalar (read as that multiple of the identity) or the full
// d x d drift covariance.
//
// `trajectory` comes back N x d -- one row per sample, parameters in
// MATLAB's order -- and `prediction` N x 1, holding `yhat(k)`, the value
// predicted for row k from the estimate held BEFORE that row was read.
//
// ⚠ THE INITIAL PARAMETER IS `eps`, NOT ZERO (`th0 = eps*ones(d,1)` in
// `rarx.m`), and P0 is `10000 * I`. The first is invisible in the answer
// and the second is not: a different P0 is a different trajectory for the
// whole record, not just for its first rows.
static bool arxTrajectory(const ICoreMatrix& data, const ICoreMatrix& orders,
const Adaptation& adaptation, const ICoreMatrix& gain,
ICoreMatrix& trajectory, ICoreMatrix& prediction,
ICoreMatrix& covariance, std::string* whyNot = nullptr);
// rpem(z, [na nb nc nd nf nk], adm, adg) -- the SAME recursion over the
// general model, and the reason `rarx` is not simply a special case of it
// in MATLAB either (they are two files, and this one is 254 lines to that
// one's 125):
//
// A(q) y(t) = [B(q)/F(q)] u(t) + [C(q)/D(q)] e(t)
//
// Zeroing the right orders gives every classical structure: ARX is
// [na nb 0 0 0 nk], ARMAX is [na nb nc 0 0 nk], output error is
// [0 nb 0 0 nf nk] and Box-Jenkins is the whole thing.
//
// ⚠ WHAT MAKES IT MORE THAN `rarx` WITH MORE PARAMETERS: the gain is
// formed from a GRADIENT vector psi and not from the regressor phi. Once
// the model has a C, D or F polynomial the prediction is no longer linear
// in the parameters, so the two vectors part company -- psi carries the
// filtered signals -- and the update uses psi where `rarx` uses phi. A
// reader who followed `rarx`'s shape and reached for phi here would get an
// answer that looks like a fit and is not the prediction-error one.
//
// ⚠ AND THE NOISE POLYNOMIALS ARE REFLECTED INTO THE UNIT DISC AT EVERY
// STEP (`fstab`), C and each F, before they are used and after they are
// written back into the parameter vector. It is not a safety net bolted
// on: an unstable C makes the next residual diverge, so the trajectory
// MATLAB answers is the one with the reflection in it and a run without it
// separates from MATLAB's within a few samples on any record that provokes
// one.
static bool predictionErrorTrajectory(const ICoreMatrix& data, const ICoreMatrix& orders,
const Adaptation& adaptation, const ICoreMatrix& gain,
ICoreMatrix& trajectory, ICoreMatrix& prediction,
ICoreMatrix& covariance, std::string* whyNot = nullptr);
};
};ICoreRecursiveLeastSquares.h#
src/ICoreBlocks/ICoreMath/Optimization/ICoreRecursiveLeastSquares.h
ICoreRecursiveLeastSquares#
ICoreRecursiveLeastSquares.h:18 · class · pImpl · 7 declaration(s)
Online (streaming) ARX transfer-function estimation via recursive least squares with a scalar exponential forgetting factor.
class ICoreRecursiveLeastSquares {
public:
// numOrder must be <= denOrder (causal/proper), matching ICoreTransferFunction's
// own causality requirement. forgettingFactor in (0, 1]; 1.0 = no forgetting.
ICoreRecursiveLeastSquares(const size_t& numOrder, const size_t& denOrder, const double& forgettingFactor = 1.0);
// Feeds one new (input, output) sample and updates the parameter estimate.
void update(const double& uk, const double& yk);
// Builds the currently-estimated transfer function (Ts = sampling time, must be > 0).
// CANONICAL: ICoreTransferFunction trims all-but-one leading zero coefficient, so
// the polynomials it hands back are as short as the estimate currently justifies —
// during warm-up, when theta is still all zeros, the numerator collapses to a single
// coefficient. Use it to display or to go on computing with; for anything that has a
// FIXED width to fill (a sized output port, a row of a collected matrix) use the two
// coefficient accessors below instead.
[[nodiscard]] ICoreTransferFunction getEstimatedTf(const double& Ts) const;
// The estimate in this class's own convention, always at full declared length:
// numerator descending [b0..b_numOrder] (numOrder+1 entries) and denominator
// descending monic [1, a1..a_denOrder] (denOrder+1 entries). Nothing is trimmed,
// so the lengths are a function of the configured orders alone and never of the
// values — which is what a fixed-size consumer needs, and what the generated
// code emits.
[[nodiscard]] std::vector<double> getNumeratorCoefficients() const;
[[nodiscard]] std::vector<double> getDenominatorCoefficients() const;
void reset();
// Declared, not implicit: the residue below is a unique_ptr to an Impl that
// is incomplete in this header, and only the .cpp can destroy one.
~ICoreRecursiveLeastSquares();
private:
class Impl; // the two-line residue; state lives here
std::unique_ptr<Impl> impl;
};ICoreScalarOptimization.h#
src/ICoreBlocks/ICoreMath/Optimization/ICoreScalarOptimization.h
ICoreScalarOptimization#
ICoreScalarOptimization.h:35 · class · 0 declaration(s)
Core MATLAB's three solvers over a function you pass in (the core-MATLAB board, C8.21-C8.23): fzero, fminbnd and fminsearch.
class ICoreScalarOptimization {
public:
// f(x) -> fx. False means the call itself failed.
using Objective = std::function<bool(double x, double& fx)>;
// f(x1..xn) -> fx, for the only one of the three that searches in n
// dimensions.
using ObjectiveN = std::function<bool(const std::vector<double>& x, double& fx)>;
// fzero(f, x0) -- search outward from x0 for a sign change, then solve
// (C8.21). MATLAB refuses rather than answering when the search leaves the
// reals or runs off to a non-finite value, and so does this.
// `options` is a resolved `optimset(...)` set (T5.11) -- the four core
// solvers take optimset's names and never optimoptions', which is
// MATLAB's own split and not this console's. Null, or a set whose fields
// are all unset, runs the default each of these three carries, and those
// are NOT the same number: fzero's TolX is `eps`, fminbnd's is 1e-4 and
// fminsearch's TolX and TolFun are both 1e-4 with a 200*n budget. It is
// the last parameter because `whyNot` already had a default.
static bool zero(const Objective& f, double x0, double& root,
std::string* whyNot = nullptr,
const ICoreOptimizationOptions::Values* options = nullptr);
// fzero(f, [a b]) -- a bracket the caller supplies. MATLAB requires the
// two ends to have OPPOSITE signs and says so by name when they do not;
// an end that is already zero is the answer.
static bool zeroInInterval(const Objective& f, double a, double b, double& root,
std::string* whyNot = nullptr,
const ICoreOptimizationOptions::Values* options = nullptr);
// fminbnd(f, a, b) (C8.23). `fval` is the value at the minimiser.
static bool minimumInInterval(const Objective& f, double a, double b,
double& minimiser, double& fval,
std::string* whyNot = nullptr,
const ICoreOptimizationOptions::Values* options = nullptr);
// fminsearch(f, x0) (C8.22). The starting point's LENGTH is the dimension,
// so a scalar x0 searches a line and a 3-vector searches a volume.
static bool minimumSearch(const ObjectiveN& f, const std::vector<double>& x0,
std::vector<double>& minimiser, double& fval,
std::string* whyNot = nullptr,
const ICoreOptimizationOptions::Values* options = nullptr);
};
};ICoreSubspaceIdentification.h#
src/ICoreBlocks/ICoreMath/Optimization/ICoreSubspaceIdentification.h
ICoreSubspaceIdentification#
ICoreSubspaceIdentification.h:60 · class · nested Options, Result · 1 declaration(s)
MATLAB's n4sid -- a state-space model of a chosen order estimated from input/output samples by the SUBSPACE method, transcribed from +idpack/@ssdata/n4sid_time.m (the toolbox board's T6.5).
class ICoreSubspaceIdentification {
public:
// MATLAB's `N4Weight`. `Auto` is `CVA` and that is source rather than a
// guess: `n4sid_time.m` line 115 reads `elseif strcmpi(n4wc,'auto'), n4wc =
// {'CVA'};`. `'SSARX'` and `'all'` are a second algorithm (`ssarx.m`, its
// own 340 lines over an ARX pre-estimate) and refuse by name rather than
// being silently answered by one of these two.
enum class Weight { Auto, CVA, MOESP };
struct Options {
Weight weight = Weight::Auto;
// MATLAB's `DisturbanceModel`: 'estimate' (the default, K is free) or
// 'none' (K fixed at zero). It is not a cosmetic flag -- it changes the
// automatic horizon rule, it makes the state matrix pass through
// `fstab`, and it changes B and D.
bool estimateDisturbance = true;
// MATLAB's `Feedthrough`: whether D is a free parameter. False by
// default for a model with states, and forced true for `nx = 0`.
bool feedthrough = false;
// MATLAB's `N4Horizon` -- {r, s1, s2}, the future and the two past
// horizons. Empty is 'auto' and is by far the common case.
std::vector<size_t> horizon;
};
// The model, in plain numbers rather than in a matrix type, which is the
// convention `ICorePredictionErrorIdentification::Model` set: this header
// is published as an API page and a caller that only wants the polynomials
// should not have to take ICoreMatrix with it. Every matrix here is ROW
// MAJOR and `order` is its side.
struct Result {
bool ok = false;
std::string message;
size_t order = 0;
// `nx * nx`, `nx`, `nx` and the scalar -- A, B, C, D of a SISO model.
// `nx = 0` leaves the first three empty and answers the static gain in
// `d` alone, which is what MATLAB's `n4sidZeroOrder` does.
std::vector<double> a, b, c;
double d = 0.0;
// The Kalman gain the subspace step implies, `nx` long. It is computed
// whether or not anything asks for it, because B, D and x0 are fitted
// to the error of the predictor it defines -- so a wrong K is a wrong
// B, not merely a missing noise model.
std::vector<double> k;
// MATLAB's second output: the initial state the fit used, `nx` long.
std::vector<double> x0;
// What `n4sid`'s estimation report calls N4Horizon and N4Weight.
std::vector<size_t> horizon;
std::string weight;
};
// `n4sid(u, y, Ts, nx)` in the console's own argument order -- input first,
// then the output, then the sample time, which is T6.1's convention for
// every estimator in this family and is required rather than optional here
// for the reason that row measured: in MATLAB an omitted sample time means
// two different things in one toolbox (`arx(u, y, [2 2 1])` answers a model
// with Ts = 1, `tfest(u, y, 2, 1)` a CONTINUOUS one with Ts = 0).
//
// ⚠ Ts SCALES NOTHING. Measured on R2026a: `n4sid(u, y, 2, 'Ts', 0.1)` and
// `n4sid(u, y, 2, 'Ts', 1)` answer the same A, B, C, D to the last bit and
// differ only in the sample time stamped on the model. The estimate is
// written in lags, so the sample time is a label on the answer and not an
// input to it. A Ts of zero is MATLAB's CONTINUOUS-time estimate, a
// different code path (`ctEstWithDTData`) that this class does not have and
// that refuses by name at the call site.
//
// u and y are vectors of the same length; nx is the model order, zero or
// more. False with `whyNot` set when the arguments do not conform or the
// data does not excite the model.
static bool n4sid(const ICoreMatrix& u, const ICoreMatrix& y, const double& Ts,
const size_t& nx, const Options& options, Result& out,
std::string* whyNot = nullptr);
// The name MATLAB spells a weight with, for the estimation report and for a
// refusal that has to say which weights exist.
static std::string weightName(const Weight& weight);
// False for a name that is not one of the two implemented; `whyNot` then
// carries the sentence that says which MATLAB has and which this console
// answers. Case-insensitive, as MATLAB's own option is.
static bool weightFromName(const std::string& name, Weight& out,
std::string* whyNot = nullptr);
};
};ICoreSystemIdentification.h#
src/ICoreBlocks/ICoreMath/Optimization/ICoreSystemIdentification.h
ICoreSystemIdentification#
ICoreSystemIdentification.h:20 · class · nested Result · 5 declaration(s)
Fits a discrete-time transfer function of a chosen order to logged input/output data.
class ICoreSystemIdentification {
public:
// The estimators, grouped by family. They differ in what they minimize and
// in how much they cost, not in the model they return.
enum class Method {
// --- Output-error / prediction-error, nonlinear ---
NonlinearLeastSquares = 0, // default: Levenberg-Marquardt on simulation error
Armax, // pseudo-linear regression with a C(q) noise model
// --- Linear-regression family (closed form or a few linear solves) ---
ArxLeastSquares, // equation-error, one least-squares solve
InstrumentalVariables, // ARX bias removed with model-simulated instruments
SteiglitzMcBride, // iteratively 1/A-prefiltered ARX; approaches the OE optimum
// --- Realization / subspace (no initial guess, no local minima) ---
SubspaceN4SID, // PO-MOESP style projection + SVD of the block-Hankel data
ImpulseRealizationEra, // least-squares FIR, then Eigensystem Realization (ERA)
RegularizedFirKernel, // TC/stable-spline regularized FIR (GCV tuned), then ERA
// --- Frequency domain, fitted to the empirical transfer function ---
FrequencyDomainLevy, // Levy's linearized rational fit
SanathananKoerner, // Levy with 1/|A| iterative reweighting
};
// Outcome of a fit. `ok` false means nothing usable came back and `message`
// says why; `ok` true with a non-empty `message` is a warning worth showing
// (a solver that hit its iteration cap, a method that could not honour the
// requested numerator order exactly, ...).
struct Result {
bool ok = false;
std::string message;
double fitPercent = 0.0; // NRMSE fit of the simulated output, 100 = exact
};
// u, y: input/output sample vectors of equal length (same length, same Ts).
// numOrder, denOrder: numerator/denominator degrees, numOrder <= denOrder.
// Ts: sample time, must be > 0.
// `outTf` receives the fitted model; the return value carries the status.
static Result fit(const ICoreMatrix& u, const ICoreMatrix& y,
const size_t& numOrder, const size_t& denOrder, const double& Ts,
const Method& method, ICoreTransferFunction& outTf);
// Convenience wrapper kept for existing callers (the console's tffit() --
// spelled tfest() until T6.2 gave that name to MATLAB's estimator -- and
// the offline IIR identification block): fits with `method` and returns just
// the model, logging any failure through ICoreRunDiagnosis as before.
static ICoreTransferFunction fitTransferFunction(const ICoreMatrix& u, const ICoreMatrix& y,
const size_t& numOrder, const size_t& denOrder,
const double& Ts,
const Method& method = Method::NonlinearLeastSquares);
// Drives tf's state-space realization with input sequence u (a vector) and
// returns the simulated output, one sample per entry of u, starting from a
// zero initial state. tf must be discrete (Ts > 0).
static ICoreMatrix simulate(const ICoreTransferFunction& tf, const ICoreMatrix& u);
///////////////////////////////
/// Method naming (UI)
//////////////////////////////
// Display names in presentation order; index 0 is the default method.
static const std::vector<std::string>& methodDisplayNames();
static std::string methodDisplayName(const Method& method);
// Returns false when `name` matches no method, leaving `out` untouched.
static bool methodFromDisplayName(const std::string& name, Method& out);
// One-line description of what the method does, for a tooltip or hint line.
static std::string methodSummary(const Method& method);
};
};