madman.helpers.ase.workflows.converge

  1r"""Convergence."""
  2
  3from collections.abc import Callable, Generator, Mapping, Sequence
  4from numbers import Integral, Real
  5from typing import Literal, Optional, TypedDict
  6
  7import numpy as np
  8import seaborn as sns
  9from ase.parallel import parprint
 10from matplotlib import pyplot as plt
 11from sympy import Expr, Function, Symbol, latex, sympify
 12
 13from madman.utilities import gen_regular_grid
 14
 15
 16sns.set_theme(
 17    context="talk",
 18    style="white",
 19    rc={"figure.titlesize": "medium", "axes.formatter.useoffset": False},
 20)
 21
 22
 23class ConvergenceParameter:
 24    r"""Convergence parameter."""
 25
 26    def __init__(
 27        self,
 28        keyword: str,
 29        *,
 30        value_min: Real,
 31        value_max: Optional[Real] = None,
 32        value_spc: Optional[Real] = None,
 33        objective: Optional[Callable[[Real], np.ndarray[Real] | Real]] = None,
 34        direction: Optional[Real] = None,
 35        criterion: Optional[
 36            Literal[
 37                "slope-neg",
 38                "slope-abs",
 39                "1st-min",
 40                "slope-abs-avg",
 41                "slope-abs-max",
 42            ]
 43        ] = None,
 44        threshold: Optional[Real] = None,
 45        stability: Optional[Integral] = None,
 46        symb: Optional[Symbol | str] = None,
 47        unit: Optional[Expr | str] = None,
 48        obj_symb: Optional[Sequence[Symbol | str, ...] | Symbol | str] = None,
 49        obj_unit: Optional[Sequence[Expr | str, ...] | Expr | str] = None,
 50    ) -> None:
 51        r"""Initialize ConvergenceParameter object.
 52
 53        Args:
 54            keyword: Parameter keyword.
 55            value_min: Parameter minimum value.
 56            value_max: Parameter maximum value.
 57            value_spc: Parameter value spacing.
 58            objective: Univariate objective function.
 59            direction: Convergence direction.
 60            criterion: Convergence criterion.
 61            threshold: Convergence threshold.
 62            stability: Convergence stability.
 63            symb: Parameter symbol for plots.
 64            unit: Parameter unit for plots.
 65            obj_symb: Objective function symbol(s) for plots.
 66            obj_unit: Objective function unit(s) for plots.
 67
 68        Raises:
 69            TypeError: If keyword is invalid.
 70            TypeError: If objective is non-callable.
 71            TypeError: If direction is non-real.
 72            ValueError: If direction is zero.
 73            ValueError: If criterion is unknown.
 74            TypeError: If threshold is non-real.
 75            ValueError: If threshold is non-positive.
 76            TypeError: If stability is non-integer.
 77            ValueError: If stability is invalid.
 78            TypeError: If symbol is invalid.
 79            TypeError: If unit is invalid.
 80            TypeError: If objective symbol is invalid.
 81            TypeError: If objective unit is invalid.
 82        """
 83        if not isinstance(keyword, str):
 84            raise TypeError("Invalid keyword!")
 85        self._keyword = keyword
 86
 87        self._value_min = float(value_min)
 88        self._value_max = float(value_max) if value_max is not None else None
 89        self._value_spc = float(value_spc) if value_spc is not None else None
 90
 91        value_seq = (
 92            gen_regular_grid(self._value_min, self._value_max, self._value_spc)
 93            if self._value_max is not None and self._value_spc is not None
 94            else [self._value_min]
 95        )
 96        self._value_seq = list(value_seq)
 97
 98        if objective is not None and not callable(objective):
 99            raise TypeError("Non-callable objective!")
100        self._objective = objective
101
102        if len(self._value_seq) > 1:
103            if not isinstance(direction, Real):
104                raise TypeError("Non-real direction!")
105            if direction == 0.0:
106                raise ValueError("Zero direction!")
107            self._direction = int(direction / abs(direction))
108
109            if self._direction < 0.0:
110                self._value_seq = sorted(self._value_seq, reverse=True)
111
112            if criterion not in {
113                "slope-neg",
114                "slope-abs",
115                "1st-min",
116                "slope-abs-avg",
117                "slope-abs-max",
118            }:
119                raise ValueError("Unknown criterion!")
120            self._criterion = criterion
121
122            if not isinstance(threshold, Real):
123                raise TypeError("Non-real threshold!")
124            if threshold <= 0.0:
125                raise ValueError("Non-positive threshold!")
126            self._threshold = float(threshold)
127
128            if not isinstance(stability, Integral):
129                raise TypeError("Non-integer stability!")
130            if stability < 1 or stability > len(self._value_seq) - 1:
131                raise ValueError("Invalid stability!")
132            self._stability = int(stability)
133
134        else:
135            self._direction = None
136            self._criterion = None
137            self._threshold = None
138            self._stability = None
139
140        if symb is None or symb == "":
141            symb = f"\\text{{{self._keyword}}}"
142        if isinstance(symb, str):
143            self._symb = Symbol(symb)
144        elif isinstance(symb, Symbol):
145            self._symb = symb
146        else:
147            raise TypeError("Invalid symbol!")
148
149        if unit is None or unit == "":
150            unit = "1"
151        if isinstance(unit, str):
152            unit = sympify(unit, rational=True)
153            for s in unit.free_symbols:
154                s_new = Symbol(f"\\mathrm{{{s.name}}}")
155                unit = unit.replace(s, s_new)
156            self._unit = unit
157        elif isinstance(unit, Expr):
158            self._unit = unit
159        else:
160            raise TypeError("Invalid unit!")
161
162        if obj_symb is None or obj_symb == "":
163            obj_symb = r"\text{objective}"
164        if isinstance(obj_symb, (Symbol, str)):
165            obj_symb = [obj_symb]
166        if not isinstance(obj_symb, Sequence) or not all(
167            isinstance(s, (Symbol, str)) for s in obj_symb
168        ):
169            raise TypeError("Invalid objective symbol!")
170
171        self._obj_symb = []
172        for s in obj_symb:
173            if s == "":
174                s = r"\text{objective}"
175            if isinstance(s, str):
176                s = Symbol(s)
177            self._obj_symb.append(s)
178
179        if obj_unit is None or obj_unit == "":
180            obj_unit = "1"
181        if isinstance(obj_unit, (Expr, str)):
182            obj_unit = [obj_unit]
183        if not isinstance(obj_unit, Sequence) or not all(
184            isinstance(s, (Expr, str)) for s in obj_unit
185        ):
186            raise TypeError("Invalid objective unit!")
187
188        self._obj_unit = []
189        for u in obj_unit:
190            if u == "":
191                u = "1"
192            if isinstance(u, str):
193                u = sympify(u, rational=True)
194                for s in u.free_symbols:
195                    s_new = Symbol(f"\\mathrm{{{s.name}}}")
196                    u = u.replace(s, s_new)
197            self._obj_unit.append(u)
198
199        self._obj_value_seq = [None] * len(self._value_seq)
200        self._obj_slope_seq = [None] * len(self._value_seq)
201        self._value_opt = np.nan
202        self._converged = False
203
204    @property
205    def keyword(self) -> str:
206        r"""Parameter keyword."""
207        return self._keyword
208
209    @property
210    def value_min(self) -> float:
211        r"""Parameter minimum value."""
212        return self._value_min
213
214    @property
215    def value_max(self) -> float:
216        r"""Parameter maximum value."""
217        return self._value_max
218
219    @property
220    def value_spc(self) -> float:
221        r"""Parameter target spacing."""
222        return self._value_spc
223
224    @property
225    def objective(self) -> Optional[Callable[[Real], np.ndarray[Real] | Real]]:
226        r"""Objective function."""
227        return self._objective
228
229    @property
230    def direction(self) -> Optional[int]:
231        r"""Convergence direction."""
232        return self._direction
233
234    @property
235    def criterion(self) -> Optional[str]:
236        r"""Convergence criterion."""
237        return self._criterion
238
239    @property
240    def threshold(self) -> Optional[float]:
241        r"""Convergence threshold."""
242        return self._threshold
243
244    @property
245    def stability(self) -> Optional[int]:
246        r"""Convergence stability."""
247        return self._stability
248
249    @property
250    def symb(self) -> Symbol:
251        r"""Parameter plot symbol."""
252        return self._symb
253
254    @property
255    def unit(self) -> Expr:
256        r"""Parameter plot unit."""
257        return self._unit
258
259    @property
260    def obj_symb(self) -> list[Symbol, ...]:
261        r"""Objective plot symbol(s)."""
262        return self._obj_symb
263
264    @property
265    def obj_unit(self) -> list[Expr, ...]:
266        r"""Objective plot unit(s)."""
267        return self._obj_unit
268
269    @property
270    def value_seq(self) -> list[float, ...]:
271        r"""Parameter values considered."""
272        return self._value_seq
273
274    @property
275    def obj_value_seq(self) -> list[np.ndarray[Real] | Real, ...]:
276        r"""Objective function values."""
277        return self._obj_value_seq
278
279    @property
280    def obj_slope_seq(self) -> list[np.ndarray[Real] | Real, ...]:
281        r"""Objective slope sequence."""
282        return self._obj_slope_seq
283
284    @property
285    def value_opt(self) -> float:
286        r"""Parameter optimal value."""
287        return self._value_opt
288
289    @property
290    def converged(self) -> float:
291        r"""True if optimal parameter value is converged."""
292        return self._converged
293
294    def converge(self) -> None:
295        r"""Converge parameter.
296
297        Returns:
298            True if convergence succeeded.
299
300        Raises:
301            ValueError: If objective function is undefined.
302            ValueError: If criterion is incompatible.
303        """
304        self.clear()
305
306        if self._objective is None:
307            raise ValueError("Undefined objective function!")
308
309        if len(self._value_seq) == 1:
310            self._value_opt = self._value_seq[0]
311            parprint(f"{self._keyword} fixed to {self._value_opt}...")
312            self._converged = True
313            return None
314
315        for i, value in enumerate(self._value_seq):
316            obj_value = self._objective(value)
317            if (
318                np.ndim(obj_value) == 0
319                and self._criterion not in {"slope-neg", "slope-abs", "1st-min"}
320            ) or (
321                np.ndim(obj_value) >= 1
322                and self._criterion not in {"slope-abs-avg", "slope-abs-max"}
323            ):
324                raise ValueError(
325                    f"`{self.keyword}` convergence criterion incompatible"
326                    + " with dimensionality of objective function values!"
327                )
328            self._obj_value_seq[i] = obj_value
329
330            if i > 0:
331                obj_diff = obj_value - self._obj_value_seq[i - 1]
332                diff = value - self._value_seq[i - 1]
333                obj_slope = obj_diff / diff
334                self._obj_slope_seq[i - 1] = obj_slope
335
336            if i >= self._stability:
337                i_stab = i - self._stability
338
339                if self._criterion == "slope-neg":
340                    is_converged = [
341                        obj_slope < 0.0 and abs(obj_slope) < self._threshold
342                        for obj_slope in self._obj_slope_seq[i_stab:i]
343                    ]
344                elif self._criterion == "slope-abs":
345                    is_converged = [
346                        abs(obj_slope) < self._threshold
347                        for obj_slope in self._obj_slope_seq[i_stab:i]
348                    ]
349                elif self._criterion == "1st-min":
350                    is_converged = [
351                        self._obj_value_seq[i_stab] < obj_value
352                        or abs(self._obj_value_seq[i_stab] - obj_value)
353                        < self._threshold
354                        for obj_value in self._obj_value_seq[i_stab + 1 : i + 1]
355                    ]
356                elif self._criterion == "slope-abs-avg":
357                    is_converged = [
358                        np.mean(np.abs(obj_slope)) < self._threshold
359                        for obj_slope in self._obj_slope_seq[i_stab:i]
360                    ]
361                elif self._criterion == "slope-abs-max":
362                    is_converged = [
363                        np.amax(np.abs(obj_slope)) < self._threshold
364                        for obj_slope in self._obj_slope_seq[i_stab:i]
365                    ]
366                else:
367                    is_converged = [False]
368
369                if all(is_converged):
370                    self._value_opt = self._value_seq[i_stab]
371                    parprint(
372                        f"{self._keyword} converged to {self._value_opt}..."
373                    )
374                    self._converged = True
375                    return None
376
377        self._value_opt = self._value_seq[-1]
378        parprint(f"{self._keyword} did not converge!")
379        return None
380
381    def plot(
382        self, qnty: Literal["obj-value", "obj-slope"]
383    ) -> Generator[plt.Axes, None, None]:
384        r"""Plot quantity vs convergence parameter.
385
386        Args:
387            qnty: Quantity to plot.
388
389        Yields:
390            Plot axes.
391
392        Raises:
393            ValueError: If quantity is unknown.
394        """
395        if qnty not in {"obj-value", "obj-slope"}:
396            raise ValueError("Unknown quantity!")
397
398        obj_dim = np.ndim(self._obj_value_seq[0])
399        if obj_dim > 2:
400            parprint("Cannot plot objective with values dimension > 2!")
401            return None
402
403        if obj_dim == 2:
404            n_plots = self._obj_value_seq[0].shape[0]
405        else:
406            n_plots = 1
407
408        if obj_dim == 0:
409            obj_value_seq = [
410                obj_value if obj_value is not None else np.nan
411                for obj_value in self._obj_value_seq
412            ]
413            obj_slope_seq = [
414                obj_slope if obj_slope is not None else np.nan
415                for obj_slope in self._obj_slope_seq
416            ]
417        else:
418            obj_value_seq = [
419                (
420                    np.reshape(obj_value, [n_plots, -1])
421                    if obj_value is not None
422                    else None
423                )
424                for obj_value in self._obj_value_seq
425            ]
426            obj_slope_seq = [
427                (
428                    np.reshape(obj_slope, [n_plots, -1])
429                    if obj_slope is not None
430                    else None
431                )
432                for obj_slope in self._obj_slope_seq
433            ]
434
435        label = f"${latex(self._symb)}$ / ${latex(self._unit)}$"
436
437        if n_plots > 1 and len(self._obj_symb) == 1:
438            obj_symbs = []
439            for i in range(n_plots):
440                s = self._obj_symb[0]
441                s = s.replace(s, Symbol(s.name + f"[{i + 1}]"))
442                obj_symbs.append(s)
443        else:
444            obj_symbs = self._obj_symb
445
446        if n_plots > 1 and len(self._obj_unit) == 1:
447            obj_units = n_plots * self._obj_unit
448        else:
449            obj_units = self._obj_unit
450
451        delta = Function(r"\Delta")
452        obj_value_labels = [
453            f"${latex(obj_symb)}$ / ${latex(obj_unit)}$"
454            for obj_symb, obj_unit in zip(obj_symbs, obj_units)
455        ]
456        obj_slope_labels = [
457            f"${latex(delta(obj_symb) / delta(self._symb))}$"
458            + " / "
459            + f"${latex(obj_unit / self._unit)}$"
460            for obj_symb, obj_unit in zip(obj_symbs, obj_units)
461        ]
462
463        for i in range(n_plots):
464            _, ax = plt.subplots(tight_layout=True)
465
466            if qnty == "obj-value":
467                ax.set_ylabel(obj_value_labels[i])
468            else:
469                ax.set_ylabel(obj_slope_labels[i])
470
471            if obj_dim == 0:
472                ax.set_xlabel(label)
473                ax.axvline(self._value_opt, linestyle="--", color="k")
474                ax.plot(
475                    self._value_seq,
476                    obj_value_seq if qnty == "obj-value" else obj_slope_seq,
477                    "o-",
478                    color="k",
479                )
480            else:
481                if qnty == "obj-value":
482                    for value, obj_value in zip(self._value_seq, obj_value_seq):
483                        if obj_value is not None:
484                            ax.plot(
485                                obj_value[i],
486                                "-" if value == self._value_opt else "--",
487                                color="k" if value == self._value_opt else None,
488                                label=f"{value:.2g}",
489                            )
490                else:
491                    ax.set_xlabel(label)
492                    obj_slope_abs_avg_seq = [
493                        (
494                            np.abs(obj_slope[i]).mean()
495                            if obj_slope is not None
496                            else np.nan
497                        )
498                        for obj_slope in obj_slope_seq
499                    ]
500                    obj_slope_abs_max_seq = [
501                        (
502                            np.abs(obj_slope[i]).max()
503                            if obj_slope is not None
504                            else np.nan
505                        )
506                        for obj_slope in obj_slope_seq
507                    ]
508                    ax.plot(
509                        self._value_seq,
510                        obj_slope_abs_avg_seq,
511                        "o-",
512                        color="k",
513                        label=r"$\text{avg}$",
514                    )
515                    ax.plot(
516                        self._value_seq,
517                        obj_slope_abs_max_seq,
518                        "v-",
519                        color="k",
520                        label=r"$\text{max}$",
521                    )
522
523                ax.legend(
524                    loc="upper right",
525                    bbox_to_anchor=(1.0, 1.0),
526                    edgecolor="k",
527                    title_fontsize="x-small",
528                    fontsize="x-small",
529                    title=label if qnty == "obj-value" else None,
530                )
531
532            if qnty == "obj-value" and self._criterion == "1st-min":
533                if not np.isnan(self._value_opt):
534                    i_opt = self._value_seq.index(self._value_opt)
535                    ax.axhspan(
536                        obj_value_seq[i_opt] - self._threshold,
537                        obj_value_seq[i_opt] + self._threshold,
538                        color="k",
539                        alpha=0.5,
540                    )
541
542            if qnty == "obj-slope" and self._criterion == "slope-neg":
543                ax.axhspan(-self._threshold, 0.0, color="k", alpha=0.5)
544
545            if qnty == "obj-slope" and self._criterion == "slope-abs":
546                threshold = self._threshold
547                ax.axhspan(-threshold, threshold, color="k", alpha=0.5)
548
549            if qnty == "obj-slope" and self._criterion in {
550                "slope-abs-avg",
551                "slope-abs-max",
552            }:
553                ax.axhspan(0.0, self._threshold, color="k", alpha=0.5)
554
555            yield ax
556
557    def clear(self) -> None:
558        r"""Clear convergence results."""
559        self._obj_value_seq = [None] * len(self._value_seq)
560        self._obj_slope_seq = [None] * len(self._value_seq)
561        self._value_opt = np.nan
562        self._converged = False
563
564
565class MultivariateConvergenceParameterSettings(TypedDict, total=False):
566    r"""Multivariate convergence parameter settings."""
567
568    value_min: Real
569    r"""Parameter minimum value."""
570
571    value_max: Optional[Real]
572    r"""Parameter maximum value."""
573
574    value_spc: Optional[Real]
575    r"""Paramater value spacing."""
576
577    direction: Optional[Real]
578    r"""Convergence direction."""
579
580    criterion: Optional[
581        Literal[
582            "slope-neg",
583            "slope-abs",
584            "1st-min",
585            "slope-abs-avg",
586            "slope-abs-max",
587        ]
588    ]
589    r"""Convergence criterion."""
590
591    threshold: Optional[Real]
592    r"""Convergence threshold."""
593
594    stability: Optional[Integral]
595    r"""Convergence stability."""
596
597    symb: Optional[Symbol | str]
598    r"""Parameter symbol for plots."""
599
600    unit: Optional[Expr | str]
601    r"""Parameter unit for plots."""
602
603
604class MultivariateConvergence:
605    r"""Multivariate convergence."""
606
607    def __init__(
608        self,
609        objective: Callable[[Real, ...], np.ndarray[Real] | Real],
610        params: Mapping[str, MultivariateConvergenceParameterSettings],
611        *,
612        crop: bool = False,
613        req_sc: bool = False,
614        niter_max: Integral = 1,
615        obj_symb: Optional[Sequence[str, ...] | str] = None,
616        obj_unit: Optional[Sequence[str, ...] | str] = None,
617    ) -> None:
618        r"""Initialize MultivariateConvergence object.
619
620        Args:
621            objective: Multivariate objective function; caching is recommended.
622            params: Convergence parameters, specified as a mapping of the
623                multivariate objective function argument keywords into the
624                corresponding convergence settings.
625                The settings are a dictionary of the ConvergenceParameter
626                argument keywords 'value_min', 'value_max', 'value_spc',
627                'direction', 'criterion', 'threshold', 'stability', 'symb' and
628                'unit', into the desired values.
629                The remaining ConvergenceParameter arguments are 'objective',
630                'obj_symb' and 'obj_unit'.
631                'objective' is set to the restriction of the multivariate
632                objective function to the univariate objective function of the
633                current parameter for the current value of the other parameters.
634                'obj_symb' and 'obj_unit' are inherited.
635            crop: If True, parameter values preceding the current optimal value
636                are removed at each iteration step.
637            req_sc: If True, self-consistency is required to achieve
638                convergence; self-consistency means that the final optimal
639                parameter values are equal to the final working parameter
640                values.
641            niter_max: Maximum number of iterations.
642            obj_symb: Objective function symbol(s) for plots.
643            obj_unit: Objective function unit(s) for plots.
644
645        Raises:
646            TypeError: If objective is non-callable.
647            TypeError: If `crop` is not boolean.
648            TypeError: If `req_sc` is not boolean.
649            TypeError: If maximum number of iterations is non-integer.
650            ValueError: If maximum number of iterations is below 1.
651        """
652        if not callable(objective):
653            raise TypeError("Non-callable objective!")
654        self._objective = objective
655
656        if not isinstance(crop, bool):
657            raise TypeError("Non-boolean `crop` parameter!")
658        self._crop = crop
659
660        if not isinstance(req_sc, bool):
661            raise TypeError("Non-boolean `req_sc` parameter!")
662        self._req_sc = req_sc
663
664        if not isinstance(niter_max, Integral):
665            raise TypeError("Non-integer maximum number of iterations!")
666        if niter_max < 1:
667            raise ValueError("Maximum number of iterations below 1!")
668        self._niter_max = int(niter_max)
669
670        self._params_seq = [{}] + [None] * (self._niter_max - 1)
671        self._wrk_pt_seq = [{}] + [None] * (self._niter_max - 1)
672        self._opt_pt_seq = [None] * self._niter_max
673        self._values_opt = {}
674
675        for keyword, settings in params.items():
676            settings.update({"obj_symb": obj_symb, "obj_unit": obj_unit})
677            param = ConvergenceParameter(keyword, objective=None, **settings)
678            self._params_seq[0].update({keyword: param})
679            self._wrk_pt_seq[0].update({keyword: param.value_seq[0]})
680            self._values_opt.update({keyword: np.nan})
681
682        for param in self._params_seq[0].values():
683            param._objective = type(self).ret_obj_wrt_param(
684                self._objective,
685                keyword=param.keyword,
686                wrk_pt=self._wrk_pt_seq[0],
687            )
688        self._converged = False
689
690    @property
691    def objective(self) -> Callable[[Real, ...], np.ndarray[Real] | Real]:
692        r"""Multivariate objective function."""
693        return self._objective
694
695    @property
696    def crop(self) -> bool:
697        r"""Crop parameter values at each iteration."""
698        return self._crop
699
700    @property
701    def req_sc(self) -> bool:
702        r"""Require self-consistency."""
703        return self._req_sc
704
705    @property
706    def niter_max(self) -> int:
707        r"""Maximum number of iterations."""
708        return self._niter_max
709
710    @property
711    def params_seq(self) -> list[dict[str, ConvergenceParameter], ...]:
712        r"""Convergence parameters sequence."""
713        return self._params_seq
714
715    @property
716    def wrk_pt_seq(self) -> list[dict[str, Real], ...]:
717        r"""Working point sequence."""
718        return self._wrk_pt_seq
719
720    @property
721    def opt_pt_seq(self) -> list[dict[str, Real], ...]:
722        r"""Optimal point sequence."""
723        return self._opt_pt_seq
724
725    @property
726    def values_opt(self) -> dict[str, Real]:
727        r"""Optimal parameter values."""
728        return self._values_opt
729
730    @property
731    def converged(self) -> bool:
732        r"""True if converged."""
733        return self._converged
734
735    def run(self) -> bool:
736        r"""Run multivariate convergence.
737
738        Returns:
739            True if convergence succeeded.
740        """
741        self.clear()
742
743        iter_i = 0
744        while iter_i < self._niter_max:
745            if iter_i > 0:
746                params = {}
747                wrk_pt = self._opt_pt_seq[iter_i - 1].copy()
748                for param in self._params_seq[iter_i - 1].values():
749                    if param.direction is None:
750                        params.update({param.keyword: param})
751                    else:
752                        value_min = (
753                            wrk_pt[param.keyword]
754                            if self._crop is True and param.direction > 0
755                            else param.value_min
756                        )
757                        value_max = (
758                            wrk_pt[param.keyword]
759                            if self._crop is True and param.direction < 0
760                            else param.value_max
761                        )
762                        param_obj = type(self).ret_obj_wrt_param(
763                            self._objective,
764                            keyword=param.keyword,
765                            wrk_pt=wrk_pt,
766                        )
767
768                        params.update(
769                            {
770                                param.keyword: ConvergenceParameter(
771                                    param.keyword,
772                                    value_min=value_min,
773                                    value_max=value_max,
774                                    value_spc=param.value_spc,
775                                    objective=param_obj,
776                                    direction=param.direction,
777                                    criterion=param.criterion,
778                                    threshold=param.threshold,
779                                    stability=param.stability,
780                                    symb=param.symb,
781                                    unit=param.unit,
782                                    obj_symb=param.obj_symb,
783                                    obj_unit=param.obj_unit,
784                                )
785                            }
786                        )
787                self._wrk_pt_seq[iter_i] = wrk_pt
788                self._params_seq[iter_i] = params
789
790            opt_pt = {}
791            for param in self._params_seq[iter_i].values():
792                param.converge()
793                if param.converged is False:
794                    parprint(
795                        f"{param.keyword} convergence failed!"
796                        + " Fixing to best value..."
797                    )
798                opt_pt.update({param.keyword: param.value_opt})
799            self._opt_pt_seq[iter_i] = opt_pt
800
801            if self._opt_pt_seq[iter_i] == self._wrk_pt_seq[iter_i]:
802                parprint("Self-consistency achieved...")
803                if all(
804                    param.converged is True
805                    for param in self._params_seq[iter_i].values()
806                ):
807                    parprint(
808                        "Multivariate convergence achieved..."
809                    )
810                    self._converged = True
811                self._values_opt = self._opt_pt_seq[iter_i]
812                return None
813            iter_i += 1
814
815        if self._req_sc is True:
816            parprint("Self-consistent multivariate convergence failed!")
817        else:
818            parprint("Non-self-consistent multivariate convergence performed...")
819        self._values_opt = self._opt_pt_seq[-1]
820        return None
821
822    def plot(
823        self,
824        qnty: Literal["obj-value", "obj-slope"],
825    ) -> dict[str, list[Generator[plt.Axes, None, None], ...]]:
826        r"""Plot quantity vs convergence parameters.
827
828        Args:
829            qnty: Quantity to plot.
830
831        Returns:
832            Plot axes.
833        """
834        plots = {keyword: [] for keyword in self._values_opt}
835
836        for params, wrk_pt in zip(self._params_seq, self._wrk_pt_seq):
837            if params is None:
838                break
839            for keyword, param in params.items():
840                wrk_pt_txt = [
841                    f"${latex(p.symb)}={latex(wrk_pt[k]*p.unit)}$"
842                    for k, p in params.items()
843                    if k != keyword
844                ]
845                wrk_pt_txt = ", ".join(wrk_pt_txt)
846
847                def param_plots(param, wrk_pt_txt):
848                    for ax in param.plot(qnty):
849                        ax.set_title(wrk_pt_txt, fontsize="x-small")
850                        yield ax
851
852                plots[keyword].append(param_plots(param, wrk_pt_txt))
853
854        return plots
855
856    def clear(self) -> None:
857        r"""Clear convergence results."""
858        i = 1
859        while i < len(self._params_seq):
860            self._params_seq[i] = None
861            self._wrk_pt_seq[i] = None
862            self._opt_pt_seq[i] = None
863            i += 1
864        for param in self._params_seq[0].values():
865            param.clear()
866        self._opt_pt_seq[0] = None
867        self._values_opt = {keyword: np.nan for keyword in self._values_opt}
868        self._converged = False
869
870    @staticmethod
871    def ret_obj_wrt_param(
872        objective: Callable[[Real, ...], np.ndarray[Real] | Real],
873        *,
874        keyword: str,
875        wrk_pt: Mapping[str, Real],
876    ) -> Callable[[Real], np.ndarray[Real] | Real]:
877        r"""Return univariate objective from multivariate objective.
878
879        Note:
880            The univariate objective w.r.t. to a parameter of the multivariate
881            objective is obtained by fixing all the other parameters to the
882            values specified by the working point.
883
884        Args:
885            objective: Multivariate objective function.
886            keyword: Univariate objective parameter keyword.
887            wrk_pt: Working point specified as a mapping of keywords of the
888                multivariate objective parameters into the corresponding
889                working point values.
890
891        Returns:
892            Univariate objective.
893        """
894
895        def obj_wrt_param(value):
896            kwargs = wrk_pt.copy()
897            kwargs.update({keyword: value})
898            return objective(**kwargs)
899
900        return obj_wrt_param