logscale option added (#258)

* logscale option added * imports improved missing directory creation added improved documentation * better plots * few issues in Fluka plots generation fixed documentation improved * unnecessary files reverted
DataMedSci · Oct 12, 2017 · 5602023 · 5602023
1 parent 05f4be4
commit 5602023
Show file tree

Hide file tree

Showing 21 changed files with 112 additions and 29 deletions.
diff --git a/.gitattributes b/.gitattributes
diff --git a/docs/21.txt.png b/docs/21.txt.png
diff --git a/docs/ex_cyl.png b/docs/ex_cyl.png
diff --git a/docs/ex_cyl.txt.png b/docs/ex_cyl.txt.png
diff --git a/docs/ex_yzmsh.png b/docs/ex_yzmsh.png
diff --git a/docs/ex_yzmsh.txt.png b/docs/ex_yzmsh.txt.png
diff --git a/docs/ex_yzmsh.txt_error.png b/docs/ex_yzmsh.txt_error.png
diff --git a/docs/ex_yzmsh2.txt.png b/docs/ex_yzmsh2.txt.png
diff --git a/docs/ex_yzmsh_error.png b/docs/ex_yzmsh_error.png
diff --git a/docs/ex_yzmsh_grey.png b/docs/ex_yzmsh_grey.png
diff --git a/docs/ex_yzmsh_logz.png b/docs/ex_yzmsh_logz.png
diff --git a/docs/ex_zmsh_logxy.png b/docs/ex_zmsh_logxy.png
diff --git a/docs/ex_zmsh_logy.png b/docs/ex_zmsh_logy.png
diff --git a/docs/gnuplot_converter.rst b/docs/gnuplot_converter.rst
@@ -20,15 +20,18 @@ As expected we will get :bash:`proton0001_fort.dat` file. Now a gnuplot script c
 
     convertmc gnuplot proton0001_fort.21
 
-We will get a script :bash:`proton0001_forttxt.plot` with following content::
+We will get a script :bash:`proton0001_fort.plot` with following content::
 
     set term png
     set output "proton0001_fort.png"
-    set title ""
+    set title "dose_z_1"
+    set xlabel ""
+    set ylabel "dose_z_1 []"
+    max(x,y) = (x > y) ? x : y
     plot  './proton0001_fort.dat' u 1:2 w l lt 1 lw 2 lc -1 title 'mean value'
 
 After running :bash:`gnuplot proton0001_fort.plot` we will get a plot
 
-.. figure:: proton0001_forttxt.png
+.. figure:: proton0001_fort.png
     :scale: 80 %
-    :alt: sample file proton0001_forttxt.png generated with gnuplot converter
+    :alt: sample file proton0001_fort.png generated with gnuplot converter
diff --git a/docs/image_converter.rst b/docs/image_converter.rst
@@ -12,31 +12,62 @@ Conversion is done using standard command::
 
     convertmc image --many "*.bdo"
 
-After converting data with 1-D scoring grid, following plot can be generated
+After converting data with 1-D scoring grid, following plot can be generated:
 
-.. figure:: ex_cyl.txt.png
+.. figure:: ex_cyl.png
     :scale: 80 %
-    :alt: sample file ex_cyl.txt.png generated with image converter
+    :alt: sample file ex_cyl.png generated with image converter
 
 Data containing 2-D scoring grid are visualised as heatmap with color denoting scored value.
 
-.. figure:: ex_yzmsh.txt.png
+.. figure:: ex_yzmsh.png
     :scale: 80 %
-    :alt: sample file ex_yzmsh.txt.png generated with image converter
+    :alt: sample file ex_yzmsh.png generated with image converter
 
 
 Options
 -------
 
+Logarithmic scale
+^^^^^^^^^^^^^^^^^
+
+User can also set logscale on one or more axis in the plots using `--log` option.
+
+An example plot with logarithmic scale on Y axis::
+
+    convertmc image --many "*.bdo" --log y
+
+.. figure:: ex_zmsh_logy.png
+    :scale: 80 %
+    :alt: sample file ex_zmsh_logy.png generated with image converter
+
+Scale can be also change on two axis at once::
+
+    convertmc image --many "*.bdo" --log x y
+
+.. figure:: ex_zmsh_logxy.png
+    :scale: 80 %
+    :alt: sample file ex_zmsh_logxy.png generated with image converter
+
+
+An example plot with 2-D heatmap and logarithmic scale on Z (color) axis::
+
+    convertmc image --many "*.bdo" --log z
+
+.. figure:: ex_yzmsh_logz.png
+    :scale: 80 %
+    :alt: sample file ex_yzmsh_logz.png generated with image converter
+
+
 Colormap
 ^^^^^^^^
 
 When generating heatmaps it is also possible to specify colormap. List of available colormaps is
 available here: http://matplotlib.org/users/colormaps.html. By default colormap called `gnuplot2` is used.
-An example plot obtained with other colormap (`Grays`) can be obtained with following command::
+An example plot obtained with other colormap (`Greys`) can be obtained with following command::
 
     convertmc image --many "*.bdo" --colormap Greys
 
-.. figure:: ex_yzmsh2.txt.png
+.. figure:: ex_yzmsh_grey.png
     :scale: 80 %
-    :alt: sample file ex_yzmsh2.txt.png generated with image converter
+    :alt: sample file ex_yzmsh_grey.png generated with image converter
diff --git a/docs/proton0001_fort.png b/docs/proton0001_fort.png
diff --git a/docs/proton0001_forttxt.png b/docs/proton0001_forttxt.png
diff --git a/pymchelper/detector.py b/pymchelper/detector.py
@@ -343,7 +343,8 @@ def merge_list(input_file_list,
         first.error /= np.sqrt(first.counter)  # np.sqrt() always returns np.float64
 
     if output_file is None:
-        output_file = input_file_list[0][:-4]
+        # strip file extension from the first item on input file list
+        output_file = os.path.splitext(input_file_list[0])[0]
 
     output_dir = os.path.dirname(output_file)
     if output_dir:  # output directory has been found, output_file is not a plain file in current dir

diff --git a/pymchelper/readers/fluka.py b/pymchelper/readers/fluka.py
@@ -59,4 +59,4 @@ def read(self, detector, nscale=1):
         # set units : detector.units are [x,y,z,v,data,detector_title]
         detector.units = [""] * 9
 
-        detector.title = usr.title.decode('ascii')
+        detector.title = usr.detector[0].name.decode('ascii')
diff --git a/pymchelper/run.py b/pymchelper/run.py
@@ -7,7 +7,7 @@
 import argparse
 
 from pymchelper.detector import merge_list, merge_many, Converters, ErrorEstimate
-from pymchelper.writers.plots import ImageWriter
+from pymchelper.writers.plots import ImageWriter, PlotAxis
 
 logger = logging.getLogger(__name__)
 
@@ -67,6 +67,13 @@ def main(args=sys.argv[1:]):
 
     parser_image = subparsers.add_parser(Converters.image.name, help='converts to PNG images')
     add_default_options(parser_image)
+    axis_names = [x.name for x in PlotAxis]
+    parser_image.add_argument('-l', '--log',
+                              help='set logscale for plot axis',
+                              nargs='+',
+                              choices=axis_names,
+                              default={},
+                              type=str)
     parser_image.add_argument("--colormap",
                               help='image color map, see http://matplotlib.org/users/colormaps.html '
                                    'for list of possible options (default: ' + ImageWriter.default_colormap + ')',
@@ -120,7 +127,8 @@ def main(args=sys.argv[1:]):
         if parsed_args.output is not None:
             output_dir = os.path.dirname(parsed_args.output)
             if output_dir and not os.path.exists(output_dir):
-                raise IOError("Directory {}/ does not exist.".format(output_dir))
+                logger.warning("Directory {:s} does not exist, creating.".format(output_dir))
+                os.makedirs(output_dir)
 
         # TODO add filename discovery
         files = sorted(glob.glob(parsed_args.input))

diff --git a/pymchelper/writers/plots.py b/pymchelper/writers/plots.py
@@ -1,11 +1,18 @@
-import os
 import logging
+import os
+from enum import IntEnum
 
 import numpy as np
 
 logger = logging.getLogger(__name__)
 
 
+class PlotAxis(IntEnum):
+    x = 1
+    y = 2
+    z = 3
+
+
 class PlotDataWriter:
     def __init__(self, filename, options):
         self.filename = filename
@@ -97,7 +104,7 @@ def write(self, detector):
             y_axis_number = detector.axis_data(1, plotting_order=True).number
             ylabel = detector.units[y_axis_number]
 
-            # for 2-D plots writte additional awk script to convert data
+            # for 2-D plots write additional awk script to convert data
             # as described in gnuplot faq: http://www.gnuplot.info/faq/faq.html#x1-320003.9
             with open(self.awk_script_filename, 'w') as script_file:
                 logger.info("Writing: " + self.awk_script_filename)
@@ -124,6 +131,7 @@ def __init__(self, filename, options):
         if not self.plot_filename.endswith(".png"):
             self.plot_filename += ".png"
         self.colormap = options.colormap
+        self.axis_with_logscale = {PlotAxis[name] for name in options.log}
 
     default_colormap = 'gnuplot2'
 
@@ -135,8 +143,20 @@ def _save_2d_error_plot(self, detector, xlist, ylist, elist):
         import matplotlib
         matplotlib.use('Agg')
         import matplotlib.pyplot as plt
+        from matplotlib import colors
+
+        # configure logscale on X and Y axis (both for positive and negative numbers)
+        if PlotAxis.x in self.axis_with_logscale:
+            plt.xscale('symlog')
+        if PlotAxis.y in self.axis_with_logscale:
+            plt.yscale('symlog')
 
-        plt.pcolormesh(xlist, ylist, elist.clip(0.0), cmap=self.colormap)
+        if PlotAxis.z in self.axis_with_logscale:
+            norm = colors.LogNorm(vmin=elist[elist > 0].min(), vmax=elist.max())
+        else:
+            norm = colors.Normalize(vmin=elist.min(), vmax=elist.max())
+
+        plt.pcolormesh(xlist, ylist, elist.clip(0.0), cmap=self.colormap, norm=norm)
         y_axis_number = detector.axis_data(1, plotting_order=True).number
         y_axis_name = detector.units[6 + y_axis_number]
         plt.ylabel(self._make_label(detector.units[y_axis_number], y_axis_name))
@@ -147,9 +167,14 @@ def _save_2d_error_plot(self, detector, xlist, ylist, elist):
         plt.close()
 
     def write(self, detector):
-        import matplotlib
-        matplotlib.use('Agg')
-        import matplotlib.pyplot as plt
+        try:
+            import matplotlib
+            matplotlib.use('Agg')
+            import matplotlib.pyplot as plt
+            from matplotlib import colors
+        except ImportError:
+            logger.error("Matplotlib not installed, output won't be generated")
+            return
 
         # skip plotting 0-D and 3-D data
         if detector.dimension in (0, 3):
@@ -172,14 +197,21 @@ def write(self, detector):
         plt.xlabel(self._make_label(detector.units[x_axis_number], x_axis_name))
         xlist = list(detector.axis_values(0, plotting_order=True))  # make list of values from generator
 
+        # configure logscale on X and Y axis (both for positive and negative numbers)
+        if PlotAxis.x in self.axis_with_logscale:
+            plt.xscale('symlog')
+
+        if PlotAxis.y in self.axis_with_logscale:
+            plt.yscale('symlog')
+
         # 1-D plotting
         if detector.dimension == 1:
 
             # add optional error area
             if np.any(detector.error):
                 plt.fill_between(xlist,
                                  (data - error).clip(0.0),
-                                 (data + error).clip(0.0, 1.05 * (detector.v.max())),
+                                 (data + error).clip(0.0, 1.05 * data.max()),
                                  alpha=0.2, edgecolor='#CC4F1B', facecolor='#FF9848', antialiased=True)
             plt.ylabel(self._make_label(detector.units[4], detector.title))
             plt.plot(xlist, data)
@@ -194,11 +226,6 @@ def write(self, detector):
             ylist = np.asarray(ylist).reshape(shape_tuple)
             zlist = data.reshape(shape_tuple)
 
-            # add error plot if error data present
-            if np.any(detector.error):
-                elist = error.reshape(shape_tuple)
-                self._save_2d_error_plot(detector, xlist, ylist, elist)
-
             x_axis_number = detector.axis_data(0, plotting_order=True).number
             x_axis_name = detector.units[6 + x_axis_number]
             plt.xlabel(self._make_label(detector.units[x_axis_number], x_axis_name))
@@ -207,8 +234,21 @@ def write(self, detector):
             y_axis_name = detector.units[6 + y_axis_number]
             plt.ylabel(self._make_label(detector.units[y_axis_number], y_axis_name))
 
-            plt.pcolormesh(xlist, ylist, zlist, cmap=self.colormap)
+            # configure logscale on Z axis
+            if PlotAxis.z in self.axis_with_logscale:
+                norm = colors.LogNorm(vmin=data[data > 0].min(), vmax=data.max())
+            else:
+                norm = colors.Normalize(vmin=data.min(), vmax=data.max())
+
+            plt.pcolormesh(xlist, ylist, zlist, cmap=self.colormap, norm=norm)
+
             cbar = plt.colorbar()
             cbar.set_label(detector.units[4], rotation=270, verticalalignment='bottom')
+
         plt.savefig(self.plot_filename)
         plt.close()
+
+        # add 2-D error plot if error data present
+        if detector.dimension == 2 and np.any(detector.error):
+            elist = error.reshape(shape_tuple)
+            self._save_2d_error_plot(detector, xlist, ylist, elist)