start dftplots examples

.gitignore

@@ -1,6 +1,7 @@
 env/*
 *.pyc
 *.swp
+.DS_Store
 output/*
 ccsgp_get_started/.changes
 ccsgp_get_started.egg-info

ccsgp_get_started/dftplots/README.md

@@ -0,0 +1,3 @@
+### DFTplots
+
+scripts to generate plots for DFT benchmarking paper 

ccsgp_get_started/dftplots/gp_bindenergy.py

@@ -0,0 +1,138 @@
+import logging, argparse, os, sys
+import numpy as np
+from collections import OrderedDict
+from ..ccsgp.ccsgp import make_panel
+from ..examples.utils import getWorkDirs
+from ..ccsgp.utils import getOpts, colorscale
+from ..ccsgp.config import default_colors
+
+def gp_bindenergy(guest):
+  """example for plotting from a text file via numpy.loadtxt
+
+  1. prepare input/output directories
+  2. load the data into an OrderedDict() [adjust axes units]
+  3. call ccsgp.make_panel with data from 2
+
+  Also see::
+
+    $ python -m ccsgp_get_started.dftplots.gp_bindenergy -h
+
+  for help on the command line options.
+
+  :param guest: guest molecule
+  :type guest: str
+  :ivar inDir: input directory according to package structure and guest
+  :ivar outDir: output directory according to package structure
+  :ivar data: OrderedDict with datasets to plot as separate keys
+  :ivar file: data input file for specific guest, format: [x y] OR [x y dx dy]
+  :ivar funct: functional, filename stem of input file
+  :ivar file_url: absolute url to input file
+  :ivar nfiles: number of files
+  :ivar dx: bar width estimated for specific amount of datasets
+  """
+  # prepare color arrays for dataset
+  my_color_set = [default_colors[3]] + [default_colors[i] for i in range(0, 3)]
+  my_color_array = []
+  for icol,color in enumerate(my_color_set):
+      my_color_array.append(colorscale(color[-7:-1], 0.8))
+      if icol == 0: continue
+      my_color_array.append(colorscale(color[-7:-1], 1.4))
+  # prepare input/output directories
+  inDir, outDir = getWorkDirs()
+  inDir = os.path.join(inDir, guest)
+  if not os.path.exists(inDir): # catch missing guest
+    return "guest %s doesn't exist" % guest
+  # prepare data
+  dpt_dict = OrderedDict()
+  if 'CO2' in guest: key = '(a)'
+  elif 'CH4' in guest: key = '(b)'
+  else: key = '(a)'
+  dpt_dict[key] = [[], [], []] #dpt
+  # exp-data
+  exp_filename = os.path.join(inDir,'experiments.dat')
+  use_exp_data = os.path.exists(exp_filename)
+  isM = guest.startswith('M-')
+  if use_exp_data:
+    exp_data = np.loadtxt(open(exp_filename, 'rb')) # load exp_data
+    if len(exp_data.shape) < 2:
+        exp_data = np.array([exp_data])
+    if isM: exp_data[:,1] /= 10.
+  yfake = -5 if use_exp_data else -40
+  if isM: yfake = 1
+  dpt_dict[key][0].append(np.array([[1,yfake]]))
+  dpt_dict[key][1].append(
+    'with boxes lt -1 lc {}'.format(default_colors[-10]) if use_exp_data else \
+    'with lines lc rgb "white"'
+  )
+  dpt_dict[key][2].append('experiments' if use_exp_data else ' ')
+  # sim-data
+  files = [
+      os.path.join(inDir, funct+ext)
+      for funct,ext in map(os.path.splitext, os.listdir(inDir))
+      if funct != 'experiments' and ext == '.dat'
+  ]
+  nfiles = len(files)
+  dx = 0.7/nfiles
+  for idx,file_url in enumerate(files):
+      funct = os.path.splitext(os.path.basename(file_url))[0]
+      data_import = np.loadtxt(open(file_url, 'rb')) # load data
+      data_import[:,0] += (idx - (nfiles-1)/2.) * dx
+      if isM: data_import[:,1] /= 10.
+      dpt_dict[key][0].append(data_import)
+      dpt_dict[key][1].append('with boxes lt -1 lc %s' % my_color_array[idx])
+      dpt_dict[key][2].append(funct)
+  # relative differences
+  if guest != 'M-O_H2O_relax' and use_exp_data:
+    dpt_dict[' '] = [[], [], []] #dpt
+    for idx,d in enumerate(dpt_dict[key][0][1:]):
+        if not isM and idx == 4: continue
+        try:
+          exp_data_mod = exp_data if idx != 0 else exp_data[(0,6),:]
+        except:
+          exp_data_mod = np.array([exp_data[0,:]])
+        diff = np.copy(d)
+        if len(diff) != len(exp_data_mod):
+          for i,dp in enumerate(diff):
+            if i >= len(exp_data_mod) or round(dp[0]) != exp_data_mod[i,0]:
+              exp_data_mod = np.insert(exp_data_mod, i, [round(dp[0]), dp[1]], 0)
+        diff[:,1] /= exp_data_mod[:,1]
+        diff[:,1] -= 1.
+        dpt_dict[' '][0].append(diff)
+        dpt_dict[' '][1].append('with boxes lt -1 lc %s' % my_color_array[idx])
+        dpt_dict[' '][2].append(dpt_dict[key][2][idx+1])
+  logging.debug(dpt_dict) # shown if --log flag given on command line
+  # generate plot using ccsgp.make_panel
+  exp_data_gpcalls = [
+    'arrow {} from {},{} to {},{} lw 4 lc {} nohead front'.format(
+      i+5, dp[0]-nfiles*dx/2., dp[1], dp[0]+nfiles*dx/2., dp[1], default_colors[-10]
+    ) for i,dp in enumerate(exp_data)
+  ] if use_exp_data else []
+  guest_split = guest.split('_')
+  ylabel = '{} distance for {} binding (nm)'.format(guest_split[0], guest_split[1]) if isM \
+    else 'electronic binding energy (kJ/mol)'
+  make_panel(
+      dpt_dict = dpt_dict,
+      gpcalls = [
+          'boxwidth {} absolute'.format(dx),
+          'style fill solid 1.0 border lt -1',
+          'xtics ("Mg"1, "Mn"2, "Fe"3, "Co"4, "Ni"5, "Cu"6, "Zn"7)',
+      ] + exp_data_gpcalls,
+      name = os.path.join(outDir, guest), yreverse = (not isM),
+      key = [ 'at graph 1.02, 1.17', 'maxrows 2', 'width -1.1', 'nobox' ],
+      ylabel = ylabel, xlabel = 'metals', rmargin = 0.99,
+      tmargin = 0.93 if guest != 'M-O_H2O_relax' and use_exp_data else 0.88,
+      size='7in,5.5in' if guest != 'M-O_H2O_relax' and use_exp_data else '4in,5.5in',
+      debug = True, layout = '1x2' if guest != 'M-O_H2O_relax' and use_exp_data else '1x1'
+  )
+  return 'done'
+
+if __name__ == '__main__':
+  parser = argparse.ArgumentParser()
+  parser.add_argument("guest", help="guest = input subdir with txt files")
+  parser.add_argument("--log", help="show log output", action="store_true")
+  args = parser.parse_args()
+  loglevel = 'DEBUG' if args.log else 'WARNING'
+  logging.basicConfig(
+    format='%(message)s', level=getattr(logging, loglevel)
+  )
+  print gp_bindenergy(args.guest)

ccsgp_get_started/dftplots/gp_interaction.py

@@ -0,0 +1,90 @@
+import logging, argparse, os, sys
+import numpy as np
+from collections import OrderedDict
+from ..ccsgp.ccsgp import make_panel
+from ..examples.utils import getWorkDirs
+from ..ccsgp.utils import getOpts, colorscale
+from ..ccsgp.config import default_colors
+from string import ascii_lowercase
+
+def gp_interaction():
+  """example for plotting from a text file via numpy.loadtxt
+
+  1. prepare input/output directories
+  2. load the data into an OrderedDict() [adjust axes units]
+  3. call ccsgp.make_panel with data from 2
+
+  Below is an output image for country initial T and the 4 most populated
+  countries for this initial (click to enlarge). Also see::
+
+    $ python -m ccsgp_get_started.examples.gp_datdir -h
+
+  for help on the command line options.
+
+  .. image:: pics/T.png
+     :width: 450 px
+
+  .. image:: pics/U.png
+     :width: 450 px
+
+  :param initial: country initial
+  :type initial: str
+  :param topN: number of most populated countries to plot
+  :type topN: int
+  :ivar inDir: input directory according to package structure and initial
+  :ivar outDir: output directory according to package structure
+  :ivar data: OrderedDict with datasets to plot as separate keys
+  :ivar file: data input file for specific country, format: [x y] OR [x y dx dy]
+  :ivar country: country, filename stem of input file
+  :ivar file_url: absolute url to input file
+  :ivar nSets: number of datasets
+  """
+  # prepare color arrays for dataset
+  my_color_set = [default_colors[i] for i in [0,2]]
+  my_color_array = []
+  for icol, color in enumerate(my_color_set):
+      my_color_array.append(colorscale(color[-7:-1], 0.8))
+      my_color_array.append(colorscale(color[-7:-1], 1.4))
+  # prepare input/output directories
+  inDir, outDir = getWorkDirs()
+  dpt_dict = OrderedDict()
+  for fold_idx,fold in enumerate(['totalE', 'Ex', 'Ec_LDA', 'Ec_vdW']):
+    energy = os.path.splitext(fold)[0]
+    key = '('+list(ascii_lowercase)[fold_idx]+')'
+    # prepare data
+    dpt_dict[key] = [[], [], []] #dpt
+    for index,funct in enumerate(os.listdir(os.path.join(inDir, energy))):
+      funct_name = os.path.splitext(funct)[0]
+      file_url = os.path.join(inDir, energy, funct)
+      data_import = np.loadtxt(open(file_url, 'rb')) #load data
+      data_import[:,0] /= 10.
+      dpt_dict[key][0].append(data_import)
+      dpt_dict[key][1].append('with linespoints lw 4 pt 18 ps 1.5 lt 2 lc %s' %
+                  my_color_array[index])
+      dpt_dict[key][2].append(funct_name)
+  #print dpt_dict
+  logging.debug(dpt_dict) # shown if --log flag given on command line
+  # generate plot using ccsgp.make_plot
+  make_panel(
+    dpt_dict = dpt_dict,
+    name = os.path.join(outDir, 'interaction'),
+    xr = [1.8, 4.53], yr = [-55, 20],
+    key = ['bottom right', 'maxrows 2', 'width -1.1', 'nobox' ],
+    ylabel = 'DFT binding energy (kJ/mol)',
+    xlabel = 'Mg-O distance, r (nm)', rmargin = 0.98, size='7in,7.5in',
+    lines = {'x=0': 'lc {} lw 4 lt 1'.format(default_colors[-10])},
+    debug = False, layout = '2x2', key_subplot_id = 2
+  )
+  return 'done'
+
+if __name__ == '__main__':
+  parser = argparse.ArgumentParser()
+  #parser.add_argument("initial", help="country initial = input subdir with txt files")
+  #parser.add_argument("topN", help="number of most populated countries to plot")
+  parser.add_argument("--log", help="show log output", action="store_true")
+  args = parser.parse_args()
+  loglevel = 'DEBUG' if args.log else 'WARNING'
+  logging.basicConfig(
+    format='%(message)s', level=getattr(logging, loglevel)
+  )
+  print gp_interaction()