Function FormInequalityPart

include/cppoptlib/function_penalty.h:155–194

// to gradients automatically via `MaxZeroExpression`.
template <typename TScalar, DifferentiabilityMode Mode, int TDim>
FunctionExpr<TScalar, Mode, TDim> FormInequalityPart(
    const ConstrainedOptimizationProblem<TScalar, Mode, TDim>& prob,
    const LagrangeMultiplierState<TScalar>& mult_state,
    const PenaltyState<TScalar>& pen_state) {
  FunctionExpr<TScalar, Mode, TDim> inequalityPart =
      ConstExpression<TScalar, Mode, TDim>(0);
  const TScalar penalty = pen_state.penalty;
  // A zero or negative penalty in the PHR formulation is ill-defined
  // (we would divide by zero in the constant offset).  The outer
  // loop's auto-scaling guarantees `rho > 0` before this function is
  // called; this guard returns the zero constant if a user bypasses
  // auto-scaling with `rho = 0`.
  if (penalty <= TScalar{0}) {
    return inequalityPart;
  }
  for (size_t j = 0; j < prob.inequality_constraints.size(); j++) {
    const TScalar mu = mult_state.inequality_multipliers[j];
    const auto& g = prob.inequality_constraints[j];
    // argument(x) = mu - penalty * g(x).  At the inner solver's
    // tape level this is a scalar - scalar*Function composition
    // built out of expression templates; `MaxZeroExpression` then
    // clamps the result at zero and carries the piecewise gradient.
    auto argument = mu - penalty * g;
    auto positive_part = MaxZeroExpression<decltype(argument)>(argument);
    // half_inv_rho = 1 / (2 rho).  Folded as a scalar multiplier so
    // the expression template still produces a single weighted
    // squared term per constraint.
    const TScalar half_inv_rho = TScalar{1} / (TScalar{2} * penalty);
    inequalityPart =
        inequalityPart + half_inv_rho * (positive_part * positive_part);
    // Subtract the constant offset `mu^2 / (2 rho)`.  A constant in
    // the composite does not affect the inner solver's minimum but
    // makes the reported composite value match the closed-form PHR
    // expression (useful for diagnostics and tests).
    const TScalar constant_offset = mu * mu * half_inv_rho;
    inequalityPart =
        inequalityPart - ConstExpression<TScalar, Mode, TDim>(constant_offset);
  }
  return inequalityPart;
}

// --- Forming the Penalty-only composite (no Lagrangian part) ---
//