Add energy analyze plot

nifty5/__init__.py

@@ -73,7 +73,7 @@ from .minimization.energy_adapter import EnergyAdapter
 from .minimization.metric_gaussian_kl import MetricGaussianKL
 
 from .sugar import *
-from .plot import Plot
+from .plot import Plot, energy_history_analysis
 
 from .library.smooth_linear_amplitude import (
     SLAmplitude, LinearSLAmplitude, CepstrumOperator)

nifty5/plot.py

@@ -17,6 +17,7 @@
 
 import os
 
+import matplotlib.pyplot as plt
 import numpy as np
 
 from . import dobj
@@ -523,3 +524,56 @@ class Plot(object):
             _plot(self._plots[i], ax, **self._kwargs[i])
         fig.tight_layout()
         _makeplot(kwargs.pop("name", None))
+
+
+def energy_history_analysis(fname_outer, fname_inner, fname_sampling,
+                            fname=None, figsize=[30, 10]):
+    """Visualizes the behaviour of the various minimizers during a
+    minimization procedure.
+
+    Parameters
+    ----------
+    fname_outer : string
+        File name of file which was written by an :class:`IterationController`
+        and contains the energy history of the Newton minimizer.
+    fname_inner : string
+        File name of file which was written by an :class:`IterationController`
+        and contains the energy history of the inverter inside of the Newton
+        minimizer.
+    fname_sampling : string
+        File name of file which was written by an :class:`IterationController`
+        and contains the energy history of the inverter which controls the
+        sampling.
+    fname : string
+        File name of the output plot. None lets matplotlib open an interactive
+        window and no file is written. Default is None.
+    figsize : tuple of float
+        figsize of output plot (see `matplotlib.pyplot.subplots`). Default is
+        [30, 10].
+    """
+    # FIXME Visualize convergence criteria as well
+    tsa, esa, _ = np.loadtxt(fname_sampling, delimiter=' ').T
+    tou, eou, _ = np.loadtxt(fname_outer, delimiter=' ').T
+    tin, ein, _ = np.loadtxt(fname_inner, delimiter=' ').T
+    t0 = np.min([tsa, tou, tin])
+    tsa = (tsa-t0)/3600
+    tou = (tou-t0)/3600
+    tin = (tin-t0)/3600
+
+    fig, (ax0, ax1) = plt.subplots(nrows=2, sharex=True, figsize=figsize)
+    p1 = ax0.scatter(tin, ein, marker='x', c='r', label='Newton inverter',
+                     alpha=.3)
+    p2 = ax0.scatter(tsa, esa, marker='x', label='Sampling inverter', alpha=.3)
+    ax0.set_ylim([1.5*np.min(esa), np.max(esa) - 0.1*np.min(esa)])
+    ax0.set_ylabel('Conjugate gradient energy')
+    ax0.legend([p1, p2], [p1.get_label(), p2.get_label()])
+    ax1.scatter(tou, eou, marker='>', c='g')
+    ax1.set_ylabel('Newton energy')
+    ax1.set_yscale('log')
+    ax1.set_xlabel("Time [h]")
+    plt.tight_layout()
+    if fname is None:
+        plt.savefig(fname)
+    else:
+        plt.show()
+    plt.close()