Coverage for dragonfly_sim/utils/dv_sensitivity.py: 88%
165 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-17 03:37 +0000
1"""Total-derivative norms per individual design parameter.
3Motivation
4----------
5``PODBasis.construct_basis``'s ``local_weighted`` mode weights snapshots by a plain
6Euclidean distance across the snapshot parameter vector. When design variables of
7different physical character share that vector -- control point motions in metres
8alongside angle of attack in radians -- an equal distance in different dimensions
9implies a wildly different change in the state, so the weighting ends up dominated
10by whichever variable happens to have the loosest bound. The intended fix is a
11diagonal weighting matrix on the distance, one weight per design parameter, sized
12by how strongly the solution responds to that parameter.
14Adjoint-based optimization never forms ``du/dx``, so the stand-in used here is the
15total derivative of the OUTPUT quantities, ``dF/dx``, which the gradient evaluation
16already computes. This module extracts it and reports, for each individual design
17parameter ``j``,
19 s_j = || [dF_1/dp_j, ..., dF_m/dp_j] ||_2
21over all optimization outputs F (objective and constraints).
23Where the numbers come from
24---------------------------
25``SimulatorBase.compute_optimization_derivatives()`` returns ``df`` of shape
26``(1, n_dv)`` and ``dc`` of shape ``(n_c, n_dv)`` -- the total derivatives of the
27objective and constraints with respect to the flattened design variable vector, in
28PHYSICAL (unscaled) design variable units. modopt applies its own
29``o_scaler``/``x_scaler`` only afterwards, so what the simulator hands back is
30exactly what is wanted here. ``attach_total_derivative_probe`` wraps that method on
31the simulator INSTANCE, so nothing in ``csdl_alpha`` or ``modopt`` is modified and
32no additional solves are performed -- it reads the dict that is already being
33returned to the optimizer.
35``sim.dv_meta`` / ``sim.c_meta`` / ``sim.o_meta`` (built by
36``csdl_alpha.backends.simulator.build_opt_metadata``) map each csdl Variable to its
37``l_ind:u_ind`` slice of the flat vector, which is what lets every column be
38labelled with its design variable and its index inside that variable.
40Units: ``param_scales``
41-----------------------
42If the POD parameter vector is not the design variable vector verbatim,
43``param_scales[j]`` is the factor mapping one to the other,
44``p_j = param_scales[j] * x_j``, so the reported sensitivities are divided by it:
46 dF/dp_j = (dF/dx_j) / param_scales[j]
48That puts the printed numbers directly in the units the POD distance is measured
49in. In DRAGonFLy the snapshot parameter vector (``DG_windtunnel_model.
50_parameter_vector``: the flat shape parameters, then alpha) is the design variable
51vector, so ``param_scales`` is normally left out.
53Dependencies
54------------
55Deliberately numpy-only in its numeric core -- no csdl, petsc or dolfinx --
56following the same rule as ``ffd_dv_utils``, so it can be exercised without an MPI
57environment. ``extract_metadata`` is the single function that touches csdl objects,
58and it only reads plain attributes (``.name``, ``.shape``, ``.size``). Callers pass
59their own ``printer`` (typically ``PETSc.Sys.Print``) to keep output on rank 0.
60"""
62import numpy as np
63from typing import NamedTuple
66class DVSensitivityMeta(NamedTuple):
67 """Labels for the rows and columns of the stacked total-derivative matrix.
69 param_labels : n_dv strings, one per individual design parameter, in the
70 order of the flattened design variable vector
71 group_slices : [(name, l_ind, u_ind)] per design VARIABLE, same order
72 output_labels: n_out strings, objective first then constraints, matching the
73 row order of ``stack_total_derivatives``
74 """
75 param_labels: list
76 group_slices: list
77 output_labels: list
79 @property
80 def n_dv(self):
81 return len(self.param_labels)
83 @property
84 def n_out(self):
85 return len(self.output_labels)
88def _entry_labels(name, shape, size):
89 """One label per entry of a variable, flattened in C order.
91 C order is not a free choice: ``build_opt_metadata`` fills the flat vector
92 with ``var.value.flatten()``, and ``ffd_dv_utils.build_cp_motion_dv`` builds
93 the matching bounds the same way, so entry ``[i, j, k]`` of a design variable
94 is at flat index ``ravel_multi_index((i, j, k), shape)``.
95 """
96 if size == 1:
97 return [name]
98 shape = tuple(int(s) for s in shape)
99 return ["{}[{}]".format(name, ",".join(str(i) for i in np.unravel_index(k, shape)))
100 for k in range(size)]
103def extract_metadata(sim):
104 """Build a DVSensitivityMeta from a csdl_alpha simulator.
106 The ONLY function here that touches csdl objects. Everything downstream works
107 on numpy arrays and plain python lists.
109 Variables created by ``CustomOperation.create_output`` carry no name (it does
110 not call ``add_name``), so unnamed outputs fall back to positional labels. Give
111 the objective a name with ``add_name`` at the call site if the printout should
112 identify it.
113 """
114 group_slices = []
115 param_labels = []
116 for k, var in enumerate(sim.dv_meta):
117 name = var.name if var.name is not None else "dv{}".format(k)
118 l_ind = sim.dv_meta[var]['l_ind']
119 u_ind = sim.dv_meta[var]['u_ind']
120 group_slices.append((name, l_ind, u_ind))
121 param_labels.extend(_entry_labels(name, var.shape, u_ind - l_ind))
123 # Sorted by l_ind rather than trusting dict order: the row order of the
124 # stacked matrix below is defined by these index ranges (see
125 # SimulatorBase._assemble_jacs), not by insertion order.
126 output_labels = []
127 for var in sorted(sim.o_meta, key=lambda v: sim.o_meta[v]['l_ind']):
128 name = var.name if var.name is not None else "objective"
129 output_labels.append("{} (objective)".format(name))
131 for k, var in enumerate(sorted(sim.c_meta, key=lambda v: sim.c_meta[v]['l_ind'])):
132 name = var.name if var.name is not None else "constraint{}".format(k)
133 size = sim.c_meta[var]['u_ind'] - sim.c_meta[var]['l_ind']
134 output_labels.extend(_entry_labels(name, var.shape, size))
136 return DVSensitivityMeta(param_labels=param_labels,
137 group_slices=group_slices,
138 output_labels=output_labels)
141def stack_total_derivatives(df, dc):
142 """Stack the objective gradient and constraint Jacobian into one (n_out, n_dv).
144 ``df`` is ``(1, n_dv)`` and ``dc`` is ``(n_c, n_dv)``, exactly as
145 ``compute_optimization_derivatives`` returns them; either may be None (an
146 unconstrained problem, or a feasibility problem with no objective). Returns
147 None when both are.
148 """
149 blocks = []
150 if df is not None:
151 df = np.atleast_2d(np.asarray(df, dtype=np.double))
152 blocks.append(df)
153 if dc is not None:
154 dc = np.atleast_2d(np.asarray(dc, dtype=np.double))
155 blocks.append(dc)
157 if len(blocks) == 0:
158 return None
159 if len(blocks) == 2 and blocks[0].shape[1] != blocks[1].shape[1]:
160 raise ValueError(
161 "objective gradient and constraint Jacobian disagree on the number of "
162 "design parameters: {} vs {}".format(blocks[0].shape[1], blocks[1].shape[1]))
164 return np.concatenate(blocks, axis=0)
167def total_derivative_norms(J, meta=None, param_scales=None, output_values=None):
168 """Column and block norms of the stacked total-derivative matrix.
170 Parameters
171 ----------
172 J : (n_out, n_dv) array from stack_total_derivatives
173 meta : DVSensitivityMeta, optional -- needed for the labelled block norms
174 param_scales : (n_dv,) array, optional. p_j = param_scales[j]*x_j; J is
175 divided by it so every reported quantity is per unit of the POD
176 parameter, not per unit of the design variable. See the module docstring.
177 output_values : (n_out,) array, optional -- the values of the outputs
178 themselves, used for the magnitude-normalized variant.
180 Returns
181 -------
182 dict with keys 'J' (the scaled matrix), 'column_norms',
183 'column_norms_relative' (None without output_values) and 'block_norms'
184 (None without meta).
185 """
186 J = np.asarray(J, dtype=np.double)
187 if J.ndim != 2:
188 raise ValueError("J must be 2-D (n_out, n_dv), got shape {}".format(J.shape))
189 n_out, n_dv = J.shape
191 if param_scales is not None:
192 param_scales = np.asarray(param_scales, dtype=np.double)
193 if param_scales.shape != (n_dv,):
194 raise ValueError("param_scales must have shape ({},), got {}"
195 .format(n_dv, param_scales.shape))
196 if np.any(param_scales == 0.):
197 raise ValueError("param_scales contains a zero; it is a per-parameter "
198 "unit conversion and must be nonzero")
199 J = J / param_scales[None, :]
201 column_norms = np.linalg.norm(J, ord=2, axis=0)
203 # The outputs are not commensurate -- in the airfoil problem D is O(1e-2) and
204 # L is O(4e-1) -- so the raw stacked norm above is dominated by whichever
205 # output happens to be largest, and any weight derived from it inherits that
206 # bias. Dividing each row by its own output value gives the relative-response
207 # variant, which is the scale-free one. Outputs that are exactly zero are left
208 # unnormalized rather than blowing the whole vector up to inf.
209 column_norms_relative = None
210 if output_values is not None:
211 output_values = np.asarray(output_values, dtype=np.double).flatten()
212 if output_values.shape != (n_out,):
213 raise ValueError("output_values must have shape ({},), got {}"
214 .format(n_out, output_values.shape))
215 denom = np.where(np.abs(output_values) > 0., np.abs(output_values), 1.)
216 column_norms_relative = np.linalg.norm(J / denom[:, None], ord=2, axis=0)
218 block_norms = None
219 if meta is not None:
220 if meta.n_dv != n_dv or meta.n_out != n_out:
221 raise ValueError(
222 "metadata describes a ({}, {}) problem but J is {}"
223 .format(meta.n_out, meta.n_dv, J.shape))
224 block_norms = {}
225 for i, output_label in enumerate(meta.output_labels):
226 for group_name, l_ind, u_ind in meta.group_slices:
227 block_norms[(output_label, group_name)] = float(
228 np.linalg.norm(J[i, l_ind:u_ind], ord=2))
230 return {'J': J,
231 'column_norms': column_norms,
232 'column_norms_relative': column_norms_relative,
233 'block_norms': block_norms}
236def format_total_derivative_norms(norms, meta, header=None, max_entries=64):
237 """Render the norms as a printable multi-line table.
239 The per-group summary and the per-output x per-group block norms are always
240 shown; the full per-parameter listing only when there are at most
241 ``max_entries`` design parameters, since a large FFD block would otherwise
242 bury the rest of the log.
243 """
244 column_norms = norms['column_norms']
245 relative = norms['column_norms_relative']
246 block_norms = norms['block_norms']
248 lines = ["=" * 78]
249 lines.append("Total-derivative norms per design parameter"
250 + ("" if header is None else " ({})".format(header)))
251 lines.append("outputs stacked: " + ", ".join(meta.output_labels))
252 lines.append("=" * 78)
254 lines.append("per design-variable group, ||dF/dp||_2 over all outputs:")
255 lines.append(" {:<24s} {:>5s} {:>12s} {:>12s} {:>12s}".format(
256 "group", "n", "min", "median", "max"))
257 for group_name, l_ind, u_ind in meta.group_slices:
258 block = column_norms[l_ind:u_ind]
259 lines.append(" {:<24s} {:>5d} {:>12.6e} {:>12.6e} {:>12.6e}".format(
260 group_name, u_ind - l_ind, block.min(), np.median(block), block.max()))
262 if block_norms is not None:
263 lines.append("")
264 lines.append("per output x group, ||dF_i/dp_group||_2:")
265 group_names = [group_name for group_name, _, _ in meta.group_slices]
266 lines.append(" {:<30s}".format("output")
267 + "".join("{:>15s}".format(name) for name in group_names))
268 for output_label in meta.output_labels:
269 lines.append(" {:<30s}".format(output_label)
270 + "".join("{:>15.6e}".format(block_norms[(output_label, name)])
271 for name in group_names))
273 if meta.n_dv <= max_entries:
274 lines.append("")
275 header_line = " {:>5s} {:<28s} {:>14s}".format("idx", "parameter", "||dF/dp||")
276 if relative is not None:
277 header_line += " {:>14s}".format("||dF/dp||/|F|")
278 lines.append(header_line)
279 for j, label in enumerate(meta.param_labels):
280 line = " {:>5d} {:<28s} {:>14.6e}".format(j, label, column_norms[j])
281 if relative is not None:
282 line += " {:>14.6e}".format(relative[j])
283 lines.append(line)
284 else:
285 lines.append("")
286 lines.append(" ({} design parameters; per-parameter listing suppressed, "
287 "max_entries={})".format(meta.n_dv, max_entries))
289 lines.append("=" * 78)
290 return "\n".join(lines)
293def print_total_derivative_norms(deriv_dict, meta, printer=print, param_scales=None,
294 header=None, max_entries=64):
295 """Compute and print the norms for one ``compute_optimization_derivatives`` result.
297 Returns the norms dict, or None when the result carried no derivatives at all.
298 """
299 J = stack_total_derivatives(deriv_dict.get('df'), deriv_dict.get('dc'))
300 if J is None:
301 printer("Total-derivative probe: no objective gradient or constraint "
302 "Jacobian in this result; nothing to report.")
303 return None
305 output_values = _stack_output_values(deriv_dict, J.shape[0])
306 norms = total_derivative_norms(J, meta=meta, param_scales=param_scales,
307 output_values=output_values)
308 printer(format_total_derivative_norms(norms, meta, header=header,
309 max_entries=max_entries))
310 return norms
313def _stack_output_values(deriv_dict, n_out):
314 """The objective and constraint VALUES, in the row order of the stacked J.
316 Returns None if they are absent or do not add up to n_out, in which case the
317 magnitude-normalized column norms are simply skipped -- they are a diagnostic,
318 not worth failing the probe over.
319 """
320 blocks = []
321 for key in ('f', 'c'):
322 value = deriv_dict.get(key)
323 if value is not None:
324 blocks.append(np.asarray(value, dtype=np.double).flatten())
325 if len(blocks) == 0:
326 return None
327 stacked = np.concatenate(blocks)
328 return stacked if stacked.shape == (n_out,) else None
331def attach_total_derivative_probe(sim, printer=print, param_scales=None,
332 record_to=None,
333 record_attr='total_derivative_norms_per_iteration',
334 eval_idx_getter=None, on_norms=None,
335 max_entries=64):
336 """Print the per-parameter total-derivative norms on every gradient evaluation.
338 Wraps ``compute_optimization_derivatives`` on the simulator INSTANCE (an
339 instance attribute shadowing the class method), which is what keeps this out of
340 csdl_alpha and modopt entirely. modopt guards that call with its own warm-``x``
341 check, so it fires once per new design point and costs no additional solves --
342 the wrapper only reads the dict already being returned.
344 Parameters
345 ----------
346 sim : a csdl_alpha SimulatorBase set up as an optimization problem
347 printer : callable taking one string; pass PETSc.Sys.Print under MPI. The
348 derivatives are replicated across ranks (the custom operations allreduce
349 their jacvec products), so rank-0 printing loses nothing.
350 param_scales : (n_dv,) design-variable-to-parameter unit conversion; see the
351 module docstring. Its length is checked against the design variable count
352 HERE, at attach time, rather than at the first gradient evaluation -- that
353 is the guard that catches a design variable being added without the
354 snapshot parameter vector being updated together.
355 record_to : optional object (a DataStore) with a list attribute named
356 ``record_attr``, appended with ``(eval_idx, column_norms)`` per call.
357 eval_idx_getter : callable returning the evaluation index these derivatives
358 belong to. Without one, the probe counts its own calls.
359 on_norms : optional callable receiving the full norms dict after each
360 gradient evaluation. This is how a consumer that wants to ACT on the
361 sensitivities gets them -- see
362 ReducedOrderModel.set_distance_sensitivities, which turns them
363 into the diagonal weights for the local_weighted POD distance. Called
364 inside the same guard as the reporting, so a consumer that raises
365 degrades the diagnostic rather than the optimization.
367 Returns
368 -------
369 callable that detaches the probe.
370 """
371 if getattr(sim, '_total_derivative_probe_attached', False):
372 raise ValueError("a total-derivative probe is already attached to this "
373 "simulator; detach it before attaching another")
375 meta = extract_metadata(sim)
377 if param_scales is not None:
378 param_scales = np.asarray(param_scales, dtype=np.double).flatten()
379 if param_scales.shape != (meta.n_dv,):
380 raise ValueError(
381 "param_scales has {} entries but the problem has {} design "
382 "parameters ({})".format(param_scales.shape[0], meta.n_dv,
383 ", ".join(name for name, _, _ in meta.group_slices)))
385 if record_to is not None and not hasattr(record_to, record_attr):
386 raise AttributeError("record_to has no attribute '{}' to append to"
387 .format(record_attr))
389 original = sim.compute_optimization_derivatives
390 call_count = [0]
392 def wrapped(*args, **kwargs):
393 deriv_dict = original(*args, **kwargs)
395 # The probe is a diagnostic bolted onto the optimizer's critical path;
396 # it must never be the reason an optimization dies. Report and carry on.
397 try:
398 eval_idx = eval_idx_getter() if eval_idx_getter is not None else call_count[0]
399 norms = print_total_derivative_norms(
400 deriv_dict, meta, printer=printer, param_scales=param_scales,
401 header="eval_idx {}".format(eval_idx), max_entries=max_entries)
402 if norms is not None and record_to is not None:
403 getattr(record_to, record_attr).append(
404 (eval_idx, norms['column_norms'].copy()))
405 if norms is not None and on_norms is not None:
406 on_norms(norms)
407 except Exception as exc:
408 printer("WARNING: total-derivative probe failed: {!r}".format(exc))
410 call_count[0] += 1
411 return deriv_dict
413 sim.compute_optimization_derivatives = wrapped
414 sim._total_derivative_probe_attached = True
416 def detach():
417 # Deleting the instance attribute re-exposes the class method underneath.
418 try:
419 del sim.compute_optimization_derivatives
420 except AttributeError:
421 pass
422 sim._total_derivative_probe_attached = False
424 return detach