Update kdeplot

doc/auxiliary-tools.rst

@@ -210,31 +210,40 @@ Averaged phonon-phonon interaction :math:`P_{\mathbf{q}j}` (in eV^2)
 feature is briefly introduced below since it may be useful. Scipy is
 needed to use this script.
 
-This script plots distribution of phonon modes in the
-frequency-lifetime plane. Its density is estimated using Gaussian-KDE
-using `scipy
+This script draws density of phonon modes in the frequency-lifetime
+plane. Its density is estimated using Gaussian-KDE using `scipy
 <https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.gaussian_kde.html>`_.
+Then (frequency, lifetime)-data points are superimposed on the density
+plot.
 
 ``kdeplot`` reads a result of the thermal conductivity calculation as
 the first argument::
 
-   % kdeplot --nbins=200 kappa-m191919.hdf5
+   % kdeplot kappa-m111111.hdf5
+
+This calculation takes some time from minutes to hours depending on
+mesh numbers and nbins. Therefore it is recommended to start with
+smaller mesh and gradually to increase mesh numbers and nbins up to
+satisfaction.
 
-(This may take long time. On a relatively fast machine, it took 20 mins.)
 After finishing the calculation, the plot is saved in
-``lifetime.png``. The black dots show the phonon modes. The density is
-estimated from these dots. The drawing area is automatically set to
-make look good, e.g, higher lifetime side is not drawn if the density
-is negligible.
+``lifetime.png``. The black dots show original data points. The
+drawing area is automatically set to make the look good, where its
+higher lifetime side is not drawn if all density beyond a lifetime
+value is smaller than some percentage (default: 10%) of the maximum
+density.
+
+The following plot is drawn with a
+19x19x19 mesh and nbins=200 and the ``Si-PBEsol`` example is used to
+generate the data.
 
 .. |ikde| image:: Si-kdeplot.png
         :width: 50%
 
 |ikde|
 
-(This plot is based on the ``Si-PBEsol`` example.)
 
-Option
+Options
 ~~~~~~~
 
 ``--temperature``
@@ -252,18 +261,60 @@ setting temperature range and step.
 ^^^^^^^^^^^^
 
 This option controls the resolution of the density plot. The default
-value is 100.
+value is 100. With larger nbins, the resolution of the plot becomes
+better, but the computation will take more time.
+
+::
+
+   % kdeplot --nbins=200 kappa-m111111.hdf5
 
 ``--cutoff``, ``--fmax``
 ^^^^^^^^^^^^^^^^^^^^^^^^^^
 
 The option ``--cutoff`` (``--fmax``) sets the maximum value of
-lifetime (frequency) to be included as data points. Normally increasing
-this value from the chosen value without specifying this option
-does nothing since automatic control of drawing area cut high lifetime
-(frequency) side if the density is low.
+lifetime (frequency) to be included as data points **before**
+Gaussian-KDE. Normally increasing this value from the chosen value
+without specifying this option does nothing since automatic control of
+drawing area cuts high lifetime (frequency) side if the density is low.
+
+``--ymax``
+^^^^^^^^^^^
+
+Maximum value of drawing region of lifetime (y) axis. This switches
+off automatic determination of drawing area, therefore as a side
+effect, the computation will be roughly twice faster.
+
+::
+
+   % kdeplot --ymax=60 kappa-m111111.hdf5
+
+``--cmap``
+^^^^^^^^^^^
+
+**New** (``kdeplot`` is a stand-alone script. The latest version is
+download at
+https://github.com/atztogo/phono3py/blob/develop/scripts/kdeplot.)
 
 
+Color map to be used for the density plot. It's given by the name
+presented at the matplotlib documentation,
+https://matplotlib.org/users/colormaps.html.
 
+::
 
+   % kdeplot --cmap="OrRd" kappa-m111111.hdf5
 
+``--min_density``
+^^^^^^^^^^^^^^^^^^
+
+**New** (``kdeplot`` is a stand-alone script. The latest version is
+download at
+https://github.com/atztogo/phono3py/blob/develop/scripts/kdeplot.)
+
+The density threshold is specified by the ratio with respect to
+maximum density. The default value is 0.1. When ``--ymax`` is
+specified together, this option is ignored.
+
+::
+
+   % kdeplot --min_density=0.01 kappa-m111111.hdf5

scripts/kdeplot

@@ -35,7 +35,7 @@ def collect_data(gamma, weights, frequencies, t_index, cutoff, max_freq):
 
     return x, y
 
-def run_KDE(x, y, nbins, y_max=None):
+def run_KDE(x, y, nbins, y_max=None, min_density=0.1):
     """Running Gaussian-KDE by scipy
     """
 
@@ -59,7 +59,7 @@ def run_KDE(x, y, nbins, y_max=None):
         for i, r_zi in enumerate((zi.T)[::-1]):
             if indices:
                 indices.append(nbins - i - 1)
-            elif np.max(r_zi) > zi_max / 10:
+            elif np.max(r_zi) > zi_max * min_density:
                 indices = [nbins - i - 1]
         short_nbinds = len(indices)
         
@@ -72,7 +72,8 @@ def run_KDE(x, y, nbins, y_max=None):
 
     return xi, yi, zi, short_nbinds
 
-def plot(plt, xi, yi, zi, x, y, short_nbinds, nbins, y_max=None):    
+def plot(plt, xi, yi, zi, x, y, short_nbinds, nbins,
+         y_max=None, cmap='viridis'):    
     #
     # Plotting
     #
@@ -88,7 +89,7 @@ def plot(plt, xi, yi, zi, x, y, short_nbinds, nbins, y_max=None):
             y_cut.append(_y)
     
     fig = plt.figure()
-    plt.pcolormesh(xi[:,:nbins], yi[:,:nbins], zi[:,:nbins])
+    plt.pcolormesh(xi[:,:nbins], yi[:,:nbins], zi[:,:nbins], cmap=cmap)
     plt.colorbar()
     
     plt.scatter(x_cut, y_cut, s=5, c='k', marker='.', linewidth=0)
@@ -108,15 +109,18 @@ def plot(plt, xi, yi, zi, x, y, short_nbinds, nbins, y_max=None):
 # Arg-parser
 parser = argparse.ArgumentParser(
     description="Plot property density with gaussian KDE")
-parser.add_argument("--fmax", help="Max frequency to plot",
-                    type=float, default=None)
+parser.add_argument("--cmap", dest="cmap", default="viridis",
+                    help="Matplotlib cmap")
 parser.add_argument("--cutoff",
                     help=("Property (y-axis) below this value is included in "
                           "data before running Gaussian-KDE"),
                     type=float, default=None)
-parser.add_argument("--ymax",
-                    help="Set maximum y of draw area",
+parser.add_argument("--fmax", help="Max frequency to plot",
                     type=float, default=None)
+parser.add_argument("--min_density",
+                    help=("Minimum density ratio with respect to maximum "
+                          "density used to determine drawing region"),
+                    type=float, default=0.1)
 parser.add_argument('--nbins', type=int, default=100,
                     help=("Number of bins in which data are assigned, "
                           "i.e., determining resolution of plot")),
@@ -125,6 +129,9 @@ parser.add_argument('--temperature', type=float, default=300.0,
                     dest='temperature',
                     help='Temperature to output data at')
 parser.add_argument("--title", dest="title", default=None, help="Plot title")
+parser.add_argument("--ymax",
+                    help="Set maximum y of draw area",
+                    type=float, default=None)
 parser.add_argument('filenames', nargs='*')
 args = parser.parse_args()
 
@@ -168,6 +175,9 @@ for i, t in enumerate(temperatures):
         t_index = i
         break
 
+print("Temperature at which lifetime density is drawn: %7.3f"
+      % temperatures[t_index])
+
 #
 # Set data
 #
@@ -189,8 +199,11 @@ if args.nu:
 for gamma, s in zip(gammas, symbols):
     x, y = collect_data(gamma, weights, frequencies,
                         t_index, args.cutoff, args.fmax)
-    xi, yi, zi, short_nbinds = run_KDE(x, y, args.nbins, y_max=args.ymax)
-    fig = plot(plt, xi, yi, zi, x, y, short_nbinds, args.nbins, y_max=args.ymax)
+    xi, yi, zi, short_nbinds = run_KDE(x, y, args.nbins,
+                                       y_max=args.ymax,
+                                       min_density=args.min_density)
+    fig = plot(plt, xi, yi, zi, x, y, short_nbinds, args.nbins,
+               y_max=args.ymax, cmap=args.cmap)
     if s:
         fig.savefig("lifetime-%s.png" % s)
     else: