Package pipe_control :: Module opendx
[hide private]
[frames] | no frames]

Source Code for Module pipe_control.opendx

  1  ############################################################################### 
  2  #                                                                             # 
  3  # Copyright (C) 2002-2005,2007-2008,2012-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