Source Code for Module pipe_control.opendx

  1  ############################################################################### 
  2  #                                                                             # 
  3  # Copyright (C) 2003-2014 Edward d'Auvergne                                   # 
  4  # Copyright (C) 2014 Troels E. Linnet                                         # 
  5  #                                                                             # 
  6  # This file is part of the program relax (http://www.nmr-relax.com).          # 
  7  #                                                                             # 
  8  # This program is free software: you can redistribute it and/or modify        # 
  9  # it under the terms of the GNU General Public License as published by        # 
 10  # the Free Software Foundation, either version 3 of the License, or           # 
 11  # (at your option) any later version.                                         # 
 12  #                                                                             # 
 13  # This program is distributed in the hope that it will be useful,             # 
 14  # but WITHOUT ANY WARRANTY; without even the implied warranty of              # 
 15  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the               # 
 16  # GNU General Public License for more details.                                # 
 17  #                                                                             # 
 18  # You should have received a copy of the GNU General Public License           # 
 19  # along with this program.  If not, see <http://www.gnu.org/licenses/>.       # 
 20  #                                                                             # 
 21  ############################################################################### 
 22   
 23  # Module docstring. 
 24  """Module containing the base class for the OpenDX space mapping classes.""" 
 25   
 26   
 27  # Python module imports. 
 28  from copy import deepcopy 
 29  from numpy import float64, array, zeros 
 30  from time import asctime, localtime 
 31   
 32  # relax module imports. 
 33  from lib.errors import RelaxError 
 34  from lib.io import open_write_file, write_data 
 35  from extern.numpy_future import percentile 
 36  from lib.software.opendx.files import write_config, write_general, write_point, write_program 
 37  from pipe_control import value 
 38  from specific_analyses.api import return_api 
 39   
 40   
 41 -def map(params=None, map_type='Iso3D', spin_id=None, inc=20, lower=None, upper=None, axis_incs=10, file_prefix="map", dir="dx", point=None, point_file="point", chi_surface=None, create_par_file=False): 
 42      """Map the space corresponding to the spin identifier and create the OpenDX files. 
 43   
 44      @keyword params:            The list of model parameters to map. 
 45      @type params:               list of str 
 46      @keyword map_type:          The type of map to create.  The available options are: 
 47                                      - 'Iso3D', a 3D isosurface visualisation of the space. 
 48      @type map_type:             str 
 49      @keyword spin_id:           The spin identification string. 
 50      @type spin_id:              str 
 51      @keyword inc:               The resolution of the plot.  This is the number of increments per dimension. 
 52      @type inc:                  int 
 53      @keyword lower:             The lower bounds of the space to map.  If supplied, this should be a list of floats, its length equal to the number of parameters in the model. 
 54      @type lower:                None or list of float 
 55      @keyword upper:             The upper bounds of the space to map.  If supplied, this should be a list of floats, its length equal to the number of parameters in the model. 
 56      @type upper:                None or list of float 
 57      @keyword axis_incs:         The number of tick marks to display in the OpenDX plot in each dimension. 
 58      @type axis_incs:            int 
 59      @keyword file_prefix:       The file prefix for all the created files. 
 60      @type file_prefix:          str 
 61      @keyword dir:               The directory to place the files into. 
 62      @type dir:                  str or None 
 63      @keyword point:             If supplied, a red sphere will be placed at these coordinates. 
 64      @type point:                None or list of float 
 65      @keyword point_file:        The file prefix for the point output files. 
 66      @type point_file:           str or None 
 67      @keyword create_par_file:   Whether to create a file with parameters and associated chi2 value. 
 68      @type create_par_file:      bool 
 69      """ 
 70   
 71      # Check the args. 
 72      if inc <= 1: 
 73          raise RelaxError("The increment value needs to be greater than 1.") 
 74      if axis_incs <= 1: 
 75          raise RelaxError("The axis increment value needs to be greater than 1.") 
 76   
 77      # Space type. 
 78      if map_type.lower() == "iso3d": 
 79          if len(params) != 3: 
 80              raise RelaxError("The 3D isosurface map requires a 3 parameter model.") 
 81   
 82          # Create the map. 
 83          Map(params, spin_id, inc, lower, upper, axis_incs, file_prefix, dir, point, point_file, chi_surface, create_par_file) 
 84      else: 
 85          raise RelaxError("The map type '" + map_type + "' is not supported.") 
 86   
 87   
 88   
 89 -class Map: 
 90      """The space mapping base class.""" 
 91   
 92 -    def __init__(self, params, spin_id, inc, lower, upper, axis_incs, file_prefix, dir, point, point_file, chi_surface, create_par_file): 
 93          """Map the space upon class instantiation.""" 
 94   
 95          # Initialise. 
 96          ############# 
 97   
 98          # Function arguments. 
 99          self.params = params 
100          self.spin_id = spin_id 
101          self.n = len(params) 
102          self.inc = inc 
103          self.axis_incs = axis_incs 
104          self.file_prefix = file_prefix 
105          self.dir = dir 
106          self.point_file = point_file 
107   
108          # Define nested listed, which holds parameter values and chi2 value. 
109          self.par_chi2_vals = [] 
110   
111          # The specific analysis API object. 
112          self.api = return_api() 
113   
114          # Points. 
115          if point != None: 
116              # Check if list is a nested list of lists. 
117              point_list = [] 
118              if isinstance(point[0], float): 
119                  point_list.append(array(point, float64)) 
120                  self.point = point_list 
121                  self.num_points = 1 
122              else: 
123                  for i in range(len(point)): 
124                      point_list.append(array(point[i], float64)) 
125                  self.point = point_list 
126                  self.num_points = i + 1 
127          else: 
128              self.num_points = 0 
129   
130          # Get the default map bounds. 
131          self.bounds = zeros((self.n, 2), float64) 
132          for i in range(self.n): 
133              # Get the bounds for the parameter i. 
134              bounds = self.api.map_bounds(self.params[i], self.spin_id) 
135   
136              # No bounds found. 
137              if not bounds: 
138                  raise RelaxError("No bounds for the parameter " + repr(self.params[i]) + " could be determined.") 
139   
140              # Assign the bounds to the global data structure. 
141              self.bounds[i] = bounds 
142   
143          # Lower bounds. 
144          if lower != None: 
145              self.bounds[:, 0] = array(lower, float64) 
146   
147          # Upper bounds. 
148          if upper != None: 
149              self.bounds[:, 1] = array(upper, float64) 
150   
151          # Setup the step sizes. 
152          self.step_size = zeros(self.n, float64) 
153          self.step_size = (self.bounds[:, 1] - self.bounds[:, 0]) / self.inc 
154   
155   
156          # Create all the OpenDX data and files. 
157          ####################################### 
158   
159          # Get the date. 
160          self.get_date() 
161   
162          # Create the strings associated with the map axes. 
163          self.map_axes() 
164   
165          # Generate the map. 
166          self.create_map() 
167   
168          ## Generate the file with parameters and associated chi2 value. 
169          if create_par_file: 
170              self.create_par_chi2(file_prefix=self.file_prefix, par_chi2_vals=self.par_chi2_vals) 
171   
172              # Generate the matplotlib script file, to plot surfaces 
173              self.matplotlib_surface_plot() 
174   
175          ## Generate the file with parameters and associated chi2 value for the points send to dx. 
176          if self.num_points >= 1 and create_par_file: 
177              # Calculate the parameter and associated chi2 values for the points. 
178              par_chi2_vals = self.calc_point_par_chi2() 
179   
180              ## Generate the file with parameters and associated chi2 value. 
181              self.create_par_chi2(file_prefix=self.point_file, par_chi2_vals=par_chi2_vals) 
182   
183          # Default the chi2 surface values, for Innermost, Inner, Middle and Outer Isosurface. 
184          if chi_surface == None: 
185              all_chi2 = array(self.all_chi, float64) 
186              innermost = percentile(all_chi2, 10) 
187              inner = percentile(all_chi2, 20) 
188              middle = percentile(all_chi2, 50) 
189              outer = percentile(all_chi2, 90) 
190              chi_surface = [innermost, inner, middle, outer] 
191   
192          # Create the OpenDX .net program file. 
193          write_program(file_prefix=self.file_prefix, point_file=self.point_file, dir=self.dir, inc=self.inc, N=self.n, num_points=self.num_points, labels=self.labels, tick_locations=self.tick_locations, tick_values=self.tick_values, date=self.date, chi_surface = chi_surface) 
194   
195          # Create the OpenDX .cfg program configuration file. 
196          write_config(file_prefix=self.file_prefix, dir=self.dir, date=self.date) 
197   
198          # Create the OpenDX .general file. 
199          write_general(file_prefix=self.file_prefix, dir=self.dir, inc=self.inc) 
200   
201          # Create the OpenDX .general and data files for the given point. 
202          if self.num_points >= 1: 
203              write_point(file_prefix=self.point_file, dir=self.dir, inc=self.inc, point=self.point, num_points=self.num_points, bounds=self.bounds, N=self.n) 
204   
205   
206 -    def calc_point_par_chi2(self): 
207          """Function for chi2 value for the points.""" 
208   
209          # Print out. 
210          print("\nCalculate chi2 value for the point parameters.") 
211   
212          # Define nested listed, which holds parameter values and chi2 value. 
213          par_chi2_vals = [] 
214   
215          # Loop over the points. 
216          for i in range(self.num_points): 
217              i_point = self.point[i] 
218   
219              # Set the parameter values. 
220              if self.spin_id: 
221                  value.set(val=i_point, param=self.params, spin_id=self.spin_id, force=True) 
222              else: 
223                  value.set(val=i_point, param=self.params, force=True) 
224   
225              # Calculate the function values. 
226              if self.spin_id: 
227                  self.api.calculate(spin_id=self.spin_id, verbosity=0) 
228              else: 
229                  self.api.calculate(verbosity=0) 
230   
231              # Get the minimisation statistics for the model. 
232              if self.spin_id: 
233                  k, n, chi2 = self.api.model_statistics(spin_id=self.spin_id) 
234              else: 
235                  k, n, chi2 = self.api.model_statistics(model_info=0) 
236   
237              # Assign value to nested list. 
238              par_chi2_vals.append([i, i_point[0], i_point[1], i_point[2], chi2]) 
239   
240          # Return list 
241          return par_chi2_vals 
242   
243   
244 -    def create_map(self): 
245          """Function for creating the map.""" 
246   
247          # Print out. 
248          print("\nCreating the map.") 
249   
250          # Open the file. 
251          map_file = open_write_file(file_name=self.file_prefix, dir=self.dir, force=True) 
252   
253          # Generate and write the text of the map. 
254          self.map_3D_text(map_file) 
255   
256          # Close the file. 
257          map_file.close() 
258   
259   
260 -    def create_par_chi2(self, file_prefix, par_chi2_vals): 
261          """Function for creating file with parameters and the chi2 value.""" 
262   
263          # Print out. 
264          print("\nCreating the file with parameters and the chi2 value.") 
265   
266          # Open the file. 
267          par_file = open_write_file(file_name=file_prefix+'.par', dir=self.dir, force=True) 
268   
269          # Copy the nested list to sort it. 
270          par_chi2_vals_sort = deepcopy(par_chi2_vals) 
271   
272          # Then sort the value. 
273          par_chi2_vals_sort.sort(key=lambda values: values[4]) 
274   
275          # Collect the data structure, which is a list of list of strings. 
276          data = [] 
277          for i, line in enumerate(par_chi2_vals): 
278              line_sort = par_chi2_vals_sort[i] 
279   
280              # Convert values to strings. 
281              line_str = ["%3.5f"%j for j in line] 
282              line_sort_str = ["%3.5f"%j for j in line_sort] 
283   
284              # Convert the index from float to index. 
285              line_str[0] = "%i" % line[0] 
286              line_sort_str[0] = "%i" % line_sort[0] 
287   
288              # Merge the two lists and append to data. 
289              data_list = line_str + line_sort_str 
290              data.append(data_list) 
291   
292          # Make the headings. 
293          headings = ['i'] + self.params + ['chi2'] 
294          headings += headings 
295   
296          # Add "_sort" to headings. 
297          headings[5] = headings[5] + "_sort" 
298          headings[6] = headings[6] + "_sort" 
299          headings[7] = headings[7] + "_sort" 
300          headings[8] = headings[8] + "_sort" 
301          headings[9] = headings[9] + "_sort" 
302   
303          # Write the parameters and chi2 values to file. 
304          write_data(out=par_file, headings=headings, data=data) 
305   
306          # Close the file. 
307          par_file.close() 
308   
309   
310 -    def get_date(self): 
311          """Function for creating a date string.""" 
312   
313          self.date = asctime(localtime()) 
314   
315   
316 -    def map_3D_text(self, map_file): 
317          """Function for creating the text of a 3D map.""" 
318   
319          # Initialise. 
320          values = zeros(3, float64) 
321          percent = 0.0 
322          percent_inc = 100.0 / (self.inc + 1.0)**(self.n - 1.0) 
323          print("%-10s%8.3f%-1s" % ("Progress:", percent, "%")) 
324   
325          # Collect all chi2, to help finding a reasobale chi level for the Innermost, Inner, Middle and Outer Isosurface. 
326          all_chi = [] 
327   
328          # Fix the diffusion tensor. 
329          unfix = False 
330          if hasattr(cdp, 'diff_tensor') and not cdp.diff_tensor.fixed: 
331              cdp.diff_tensor.fixed = True 
332              unfix = True 
333   
334          # Initial value of the first parameter. 
335          values[0] = self.bounds[0, 0] 
336   
337          # Define counter 
338          counter = 0 
339   
340          # Loop over the first parameter. 
341          for i in range((self.inc + 1)): 
342              # Initial value of the second parameter. 
343              values[1] = self.bounds[1, 0] 
344   
345              # Loop over the second parameter. 
346              for j in range((self.inc + 1)): 
347                  # Initial value of the third parameter. 
348                  values[2] = self.bounds[2, 0] 
349   
350                  # Loop over the third parameter. 
351                  for k in range((self.inc + 1)): 
352                      # Set the parameter values. 
353                      if self.spin_id: 
354                          value.set(val=values, param=self.params, spin_id=self.spin_id, verbosity=0, force=True) 
355                      else: 
356                          value.set(val=values, param=self.params, verbosity=0, force=True) 
357   
358                      # Calculate the function values. 
359                      if self.spin_id: 
360                          self.api.calculate(spin_id=self.spin_id, verbosity=0) 
361                      else: 
362                          self.api.calculate(verbosity=0) 
363   
364                      # Get the minimisation statistics for the model. 
365                      if self.spin_id: 
366                          k, n, chi2 = self.api.model_statistics(spin_id=self.spin_id) 
367                      else: 
368                          k, n, chi2 = self.api.model_statistics(model_info=0) 
369   
370                      # Set maximum value to 1e20 to stop the OpenDX server connection from breaking. 
371                      if chi2 > 1e20: 
372                          map_file.write("%30f\n" % 1e20) 
373                      else: 
374                          map_file.write("%30f\n" % chi2) 
375   
376                          # Save all values of chi2. To help find reasonale level for the Innermost, Inner, Middle and Outer Isosurface. 
377                          all_chi.append(chi2) 
378   
379                      # Assign value to nested list. 
380                      self.par_chi2_vals.append([counter, values[0], values[1], values[2], chi2]) 
381   
382                      # Add to counter. 
383                      counter += 1 
384   
385                      # Increment the value of the third parameter. 
386                      values[2] = values[2] + self.step_size[2] 
387   
388                  # Progress incrementation and printout. 
389                  percent = percent + percent_inc 
390                  print("%-10s%8.3f%-8s%-8g" % ("Progress:", percent, "%,  " + repr(values) + ",  f(x): ", chi2)) 
391   
392                  # Increment the value of the second parameter. 
393                  values[1] = values[1] + self.step_size[1] 
394   
395              # Increment the value of the first parameter. 
396              values[0] = values[0] + self.step_size[0] 
397   
398          # Unfix the diffusion tensor. 
399          if unfix: 
400              cdp.diff_tensor.fixed = False 
401   
402          # Save all chi2 values. 
403          self.all_chi = all_chi 
404   
405 -    def map_axes(self): 
406          """Function for creating labels, tick locations, and tick values for an OpenDX map.""" 
407   
408          # Initialise. 
409          self.labels = "{" 
410          self.tick_locations = [] 
411          self.tick_values = [] 
412          loc_inc = float(self.inc) / float(self.axis_incs) 
413   
414          # Loop over the parameters 
415          for i in range(self.n): 
416              # Parameter conversion factors. 
417              factor = self.api.return_conversion_factor(self.params[i]) 
418   
419              # Parameter units. 
420              units = self.api.return_units(self.params[i]) 
421   
422              # Labels. 
423              if units: 
424                  self.labels = self.labels + "\"" + self.params[i] + " (" + units + ")\"" 
425              else: 
426                  self.labels = self.labels + "\"" + self.params[i] + "\"" 
427   
428              if i < self.n - 1: 
429                  self.labels = self.labels + " " 
430              else: 
431                  self.labels = self.labels + "}" 
432   
433              # Tick values. 
434              vals = self.bounds[i, 0] / factor 
435              val_inc = (self.bounds[i, 1] - self.bounds[i, 0]) / (self.axis_incs * factor) 
436   
437              string = "" 
438              for j in range(self.axis_incs + 1): 
439                  string = string + "\"" + "%.2f" % vals + "\" " 
440                  vals = vals + val_inc 
441              self.tick_values.append("{" + string + "}") 
442   
443              # Tick locations. 
444              string = "" 
445              val = 0.0 
446              for j in range(self.axis_incs + 1): 
447                  string = string + " " + repr(val) 
448                  val = val + loc_inc 
449              self.tick_locations.append("{" + string + " }") 
450   
451   
452 -    def matplotlib_surface_plot(self): 
453          """Function to write matplotlib script to plot surfaces of parameters.""" 
454   
455          # Add ".par" to file prefix 
456          mapfile_name = '"%s.par"' % self.file_prefix 
457   
458          # If point file_file is different from None 
459          if self.point_file != None: 
460              pointfile_name = '"%s.par"' % self.point_file 
461          else: 
462              pointfile_name = "None" 
463   
464          # Open the file. 
465          plot_file = open_write_file(file_name=self.file_prefix+'.py', dir=self.dir, force=True) 
466   
467          matplotlib_file = [ 
468              'from copy import deepcopy'+"\n", 
469              'import numpy as np'+"\n", 
470              'import scipy.interpolate'+"\n", 
471              'from numpy.ma import masked_where'+"\n", 
472              ''+"\n", 
473              'from mpl_toolkits.mplot3d import axes3d'+"\n", 
474              'import matplotlib.pyplot as plt'+"\n", 
475              'from matplotlib import cm'+"\n", 
476              ''+"\n", 
477              '# Open file and get header.'+"\n", 
478              'mapfile_name = %s'%mapfile_name+"\n", 
479              'pointfile_name = %s'%pointfile_name+"\n", 
480              ''+"\n", 
481              'mapfile = open(mapfile_name, "r")'+"\n", 
482              'lines = mapfile.readlines()'+"\n", 
483              'mapfile.close()'+"\n", 
484              'header = lines[0].split()[1:]'+"\n", 
485              ''+"\n", 
486              '# Prepare the dtype for reading file.'+"\n", 
487              'dtype_str = "i8,f8,f8,f8,f8,i8,f8,f8,f8,f8"'+"\n", 
488              ''+"\n", 
489              'print("Fileheader is: %s"%header)'+"\n", 
490              'print("Value types are: %s"%dtype_str)'+"\n", 
491              ''+"\n", 
492              '# Load the data.'+"\n", 
493              'data = np.genfromtxt(fname=mapfile_name, dtype=dtype_str, names=header)'+"\n", 
494              ''+"\n", 
495              '# Load the point data'+"\n", 
496              'if pointfile_name:'+"\n", 
497              '    # Load the point data.'+"\n", 
498              '    data_p = np.genfromtxt(fname=pointfile_name, dtype=dtype_str, names=header)'+"\n", 
499              '    '+"\n", 
500              ''+"\n", 
501              '# Define where to cut the data, as the minimum.'+"\n", 
502              'header_min = header[6:10]'+"\n", 
503              ''+"\n", 
504              '# Define to cut at min map point.'+"\n", 
505              'map_min_par0 = data[header_min[0]][0]'+"\n", 
506              'map_min_par1 = data[header_min[1]][0]'+"\n", 
507              'map_min_par2 = data[header_min[2]][0]'+"\n", 
508              'map_min_chi2 = data[header_min[3]][0]'+"\n", 
509              ''+"\n", 
510              '# Now get the headers for the data.'+"\n", 
511              'header_val = header[1:5]'+"\n", 
512              ''+"\n", 
513              '# Define which 2D maps to create, as a list of 2 parameters, and at which third parameter to cut the values.'+"\n", 
514              'maps_xy = [header_val[0], header_val[1], header_val[2], map_min_par2]'+"\n", 
515              'maps_xz = [header_val[0], header_val[2], header_val[1], map_min_par1]'+"\n", 
516              'maps_yz = [header_val[1], header_val[2], header_val[0], map_min_par0]'+"\n", 
517              ''+"\n", 
518              'maps = [maps_xy, maps_xz, maps_yz]'+"\n", 
519              ''+"\n", 
520              '# Nr of columns is number of maps.'+"\n", 
521              'nr_cols = 1'+"\n", 
522              '# Nr of rows, is 2, for 3d projection and imshow'+"\n", 
523              'nr_rows = 2'+"\n", 
524              ''+"\n", 
525              '# Loop over the maps:'+"\n", 
526              'for x_par, y_par, z_par, z_cut in maps:'+"\n", 
527              '    # Define figure'+"\n", 
528              '    fig = plt.figure()'+"\n", 
529              ''+"\n", 
530              '    # Define c_par'+"\n", 
531              '    c_par = header_val[3]'+"\n", 
532              ''+"\n", 
533              '    # Now get the values for the map.'+"\n", 
534              '    map_x = data[x_par]'+"\n", 
535              '    map_y = data[y_par]'+"\n", 
536              '    map_z = data[z_par]'+"\n", 
537              '    map_c = data[c_par]'+"\n", 
538              ''+"\n", 
539              '    # Now define which map to create.'+"\n", 
540              '    mask_xy = masked_where(map_z == z_cut, map_z)'+"\n", 
541              '    map_mask_x = map_x[mask_xy.mask]'+"\n", 
542              '    map_mask_y = map_y[mask_xy.mask]'+"\n", 
543              '    map_mask_c = map_c[mask_xy.mask]'+"\n", 
544              ''+"\n", 
545              '    # Define min and max values.'+"\n", 
546              '    map_mask_x_min = map_mask_x.min()'+"\n", 
547              '    map_mask_x_max = map_mask_x.max()'+"\n", 
548              '    map_mask_y_min = map_mask_y.min()'+"\n", 
549              '    map_mask_y_max = map_mask_y.max()'+"\n", 
550              '    map_mask_c_min = map_mask_c.min()'+"\n", 
551              '    map_mask_c_max = map_mask_c.max()'+"\n", 
552              ''+"\n", 
553              '    # Set up a regular grid of interpolation points'+"\n", 
554              '    int_points = 300'+"\n", 
555              '    xi, yi = np.linspace(map_mask_x_min, map_mask_x_max, int_points), np.linspace(map_mask_y_min, map_mask_y_max, int_points)'+"\n", 
556              '    xi, yi = np.meshgrid(xi, yi)'+"\n", 
557              ''+"\n", 
558              '    # Interpolate to create grid'+"\n", 
559              '    ci = scipy.interpolate.griddata((map_mask_x, map_mask_y), map_mask_c, (xi, yi), method="linear")'+"\n", 
560              ''+"\n", 
561              '    # Set which x, y, z to plot'+"\n", 
562              '    x_p = xi'+"\n", 
563              '    y_p = yi'+"\n", 
564              '    c_p = deepcopy(ci)'+"\n", 
565              ''+"\n", 
566              '    # Cut map at a certain height.'+"\n", 
567              '    # First get index os largest values'+"\n", 
568              '    #z_max = map_mask_c_max'+"\n", 
569              '    z_max = map_mask_c_min + 0.5*map_mask_c_min'+"\n", 
570              '    ci_mask = masked_where(ci >= z_max, ci)'+"\n", 
571              ''+"\n", 
572              '    # Replace with 0.0'+"\n", 
573              '    c_p[ci_mask.mask] = 0.0'+"\n", 
574              '    # Find new max'+"\n", 
575              '    new_max = np.max(c_p)'+"\n", 
576              ''+"\n", 
577              '    # Insert values in array.'+"\n", 
578              '    c_p[ci_mask.mask] = new_max'+"\n", 
579              ''+"\n", 
580              '    # Define min.'+"\n", 
581              '    z_min = map_mask_c_min - 0.5*map_mask_c_min'+"\n", 
582              ''+"\n", 
583              '    # Create figure and plot'+"\n", 
584              '    ax = fig.add_subplot(nr_rows, nr_cols, 1, projection="3d")'+"\n", 
585              '    ax.plot_surface(x_p, y_p, c_p, rstride=8, cstride=8, alpha=0.3)'+"\n", 
586              ''+"\n", 
587              '    # Possible add scatter points for map.'+"\n", 
588              '    #ax.scatter(map_x, map_y, map_c, c="b", marker="o", s=5)'+"\n", 
589              ''+"\n", 
590              '    # One could also make the mesh just from the values, but this require much memory.'+"\n", 
591              '    ##ax.scatter(x_p, y_p, c_p, c="y", marker="o", s=5)'+"\n", 
592              ''+"\n", 
593              '    # Add contour levels on sides.'+"\n", 
594              '    ax.contour(x_p, y_p, c_p, zdir="z", offset=z_min, cmap=cm.coolwarm)'+"\n", 
595              '    ax.contour(x_p, y_p, c_p, zdir="x", offset=map_mask_x_min, cmap=cm.coolwarm)'+"\n", 
596              '    ax.contour(x_p, y_p, c_p, zdir="y", offset=map_mask_y_min, cmap=cm.coolwarm)'+"\n", 
597              ''+"\n", 
598              '    # Add scatter values, for 5 lowest values.'+"\n", 
599              '    x_par_min = x_par + "_sort"'+"\n", 
600              '    y_par_min = y_par + "_sort"'+"\n", 
601              '    c_par_min = c_par + "_sort"'+"\n", 
602              '    mp_x = data[x_par_min][0:5]'+"\n", 
603              '    mp_y = data[y_par_min][0:5]'+"\n", 
604              '    mp_c = data[c_par_min][0:5]'+"\n", 
605              '    ax.scatter(mp_x[0], mp_y[0], mp_c[0], c="r", marker="o", s=200)'+"\n", 
606              '    ax.scatter(mp_x[1:], mp_y[1:], mp_c[1:], c="g", marker="o", s=100)'+"\n", 
607              ''+"\n", 
608              '    # Add points from file, as the closest point in map.'+"\n", 
609              '    if pointfile_name:'+"\n", 
610              '        if data_p[x_par].ndim == 0:'+"\n", 
611              '            points_x = np.asarray([data_p[x_par]])'+"\n", 
612              '            points_y = np.asarray([data_p[y_par]])'+"\n", 
613              '        else:'+"\n", 
614              '            points_x = data_p[x_par]'+"\n", 
615              '            points_y = data_p[y_par]'+"\n", 
616              ''+"\n", 
617              '        # Normalize, by division of largest number of map'+"\n", 
618              '        points_x_norm = points_x / map_mask_x_max'+"\n", 
619              '        points_y_norm = points_y / map_mask_y_max'+"\n", 
620              '        map_mask_x_norm = map_mask_x / map_mask_x_max'+"\n", 
621              '        map_mask_y_norm = map_mask_y / map_mask_y_max'+"\n", 
622              ''+"\n", 
623              '        p_x = []'+"\n", 
624              '        p_y = []'+"\n", 
625              '        p_c = []'+"\n", 
626              '        # Now calculate the Euclidean distance in the space, to find map point best represents the point.'+"\n", 
627              '        for i, point_x_norm in enumerate(points_x_norm):'+"\n", 
628              '            point_y_norm = points_y_norm[i]'+"\n", 
629              ''+"\n", 
630              '            # Get the distance.'+"\n", 
631              '            dist = np.sqrt( (map_mask_x_norm - point_x_norm)**2 + (map_mask_y_norm - point_y_norm)**2)'+"\n", 
632              ''+"\n", 
633              '            # Return the indices of the minimum values along an axis.'+"\n", 
634              '            min_index = np.argmin(dist)'+"\n", 
635              '            p_x.append(map_mask_x[min_index])'+"\n", 
636              '            p_y.append(map_mask_y[min_index])'+"\n", 
637              '            p_c.append(map_mask_c[min_index])'+"\n", 
638              ''+"\n", 
639              '        # Convert to numpy array'+"\n", 
640              '        p_x = np.asarray(p_x)'+"\n", 
641              '        p_y = np.asarray(p_y)'+"\n", 
642              '        p_c = np.asarray(p_c)'+"\n", 
643              ''+"\n", 
644              '        # Plot points'+"\n", 
645              '        ax.scatter(p_x, p_y, p_c, c="m", marker="o", s=50)'+"\n", 
646              ''+"\n", 
647              ''+"\n", 
648              '    # Set label'+"\n", 
649              '    ax.set_xlabel("%s"%x_par)'+"\n", 
650              '    ax.set_ylabel("%s"%y_par)'+"\n", 
651              '    ax.set_zlabel("%s"%c_par)'+"\n", 
652              ''+"\n", 
653              ''+"\n", 
654              '    # Set limits'+"\n", 
655              '    ax.set_zlim(z_min, z_max)'+"\n", 
656              ''+"\n", 
657              ''+"\n", 
658              '    # Create figure and plot'+"\n", 
659              '    ax = fig.add_subplot(nr_rows, nr_cols, 2)'+"\n", 
660              '    fig_imshow = ax.imshow(ci, vmin=map_mask_c_min, vmax=map_mask_c_max, origin="lower", extent=[map_mask_x_min, map_mask_x_max, map_mask_y_min, map_mask_y_max])'+"\n", 
661              ''+"\n", 
662              '    # Add scatter values, for 5 lowest values.'+"\n", 
663              '    ax.scatter(mp_x[0], mp_y[0], c=mp_c[0], marker="o", s=200)'+"\n", 
664              '    ax.scatter(mp_x[1:], mp_y[1:], c="g", marker="o", s=100)'+"\n", 
665              ''+"\n", 
666              '    # Also add point to this map.'+"\n", 
667              '    if pointfile_name:'+"\n", 
668              '        # Plot points'+"\n", 
669              '        ax.scatter(p_x, p_y, c="m", marker="o", s=50)'+"\n", 
670              ''+"\n", 
671              '    # Set label'+"\n", 
672              '    ax.set_xlabel("%s"%x_par)'+"\n", 
673              '    ax.set_ylabel("%s"%y_par)'+"\n", 
674              ''+"\n", 
675              '    # Add colorbar.'+"\n", 
676              '    fig.subplots_adjust(right=0.8)'+"\n", 
677              '    cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.3])'+"\n", 
678              '    fig.colorbar(fig_imshow, cax=cbar_ax)'+"\n", 
679              ''+"\n", 
680              '# Show plot.'+"\n", 
681              'plt.show()'+"\n", 
682              ''+"\n", 
683          ] 
684   
685          # Loop over the lines and write. 
686          for line in matplotlib_file: 
687              plot_file.write(line) 
688   
689          # Close the file. 
690          plot_file.close() 
691