Coverage for encodermap/plot/plotting.py: 42%

1# -*- coding: utf-8 -*-

2# encodermap/plot/plotting.py

3################################################################################

4# Encodermap: A python library for dimensionality reduction.

8# Authors:

9# Kevin Sawade, Tobias Lemke

10#

11# Encodermap is free software: you can redistribute it and/or modify

12# it under the terms of the GNU Lesser General Public License as

13# published by the Free Software Foundation, either version 2.1

14# of the License, or (at your option) any later version.

15# This package is distributed in the hope that it will be useful to other

16# researches. IT DOES NOT COME WITH ANY WARRANTY WHATSOEVER; without even the

17# implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

18# See the GNU Lesser General Public License for more details.

19#

20# See <http://www.gnu.org/licenses/>.

21################################################################################

22"""Convenience functions for Plotting.

24Todo:

25 * Add interactive Plotting

26 * Find a way to use interactive plotting with less points but still cluster everything.

28"""

30##############################################################################

31# Imports

32##############################################################################

35from __future__ import annotations

37import os

38import shutil

39import subprocess

40import time

42import matplotlib as mpl

43import matplotlib.pyplot as plt

44import numpy as np

46from .._optional_imports import _optional_import

47from ..encodermap_tf1.misc import periodic_distance_np, sigmoid

48from ..misc.clustering import gen_dummy_traj, rmsd_centroid_of_cluster

50################################################################################

51# Optional Imports

52################################################################################

55md = _optional_import("mdtraj")

58##############################################################################

59# Globals

60##############################################################################

62__all__ = ["distance_histogram"]

64##############################################################################

65# Utilities

66##############################################################################

69def distance_histogram(

70 data, periodicity, sigmoid_parameters, axes=None, low_d_max=5, bins="auto"

71):

72 """

73 Plots the histogram of all pairwise distances in the data.

74 It also shows the sigmoid function and its normalized derivative.

76 Args:

77 data (np.ndarray): 2-dimensional numpy array. Columns should iterate over the dimensions of the datapoints,

78 i.e. the dimensionality of the data. The rows should iterate over datapoints.

79 periodicity (float): Periodicity of the data. Use float("inf") for non-periodic data.

80 sigmoid_parameters (tuple): Tuple of sketchmap sigmoid parameters in shape (sigma, a, b).

81 axes (Union[np.ndarray, None], optional): A numpy array of two matplotlib.axes objects or None. If None is provided,

82 the axes will be created. Defaults to None.

83 low_d_max (int, optional): Upper limit for plotting the low_d sigmoid. Defaults to 5.

84 bins (Union[str, int]. optional): Number of bins for histogram. Use 'auto' to let matplotlib decide how

85 many bins to use. Defaults to 'auto'.

88 """

89 vecs = periodic_distance_np(

90 np.expand_dims(data, axis=1), np.expand_dims(data, axis=0), periodicity

91 )

92 dists = np.linalg.norm(vecs, axis=2)

93 while True:

94 try:

95 dists = np.linalg.norm(dists, axis=2)

96 except np.AxisError:

97 break

98 dists = dists.reshape(-1)

100 if axes is None:

101 fig, axes = plt.subplots(2)

102 axe2 = axes[0].twinx()

103 counts, edges, patches = axe2.hist(

104 dists, bins=bins, density=True, edgecolor="black"

105 )

106 x = np.linspace(0, max(dists), 1000)

107

108 y = sigmoid(x, *sigmoid_parameters[:3])

109 edges_sig = sigmoid(edges, *sigmoid_parameters[:3])

110 dy = np.diff(y)

111 dy_norm = dy / max(dy)

112 axes[0].plot(x, y, color="C1", label="sigmoid")

113 axes[0].plot(x[:-1], dy_norm, color="C2", label="diff sigmoid")

114

115 axes[0].legend()

116 axes[0].set_xlabel("distance")

117 axes[0].set_ylim((0, 1))

118 axes[0].set_zorder(axe2.get_zorder() + 1)

119 axes[0].patch.set_visible(False)

120 axes[0].set_title("high-d")

121

122 x = np.linspace(0, low_d_max, 1000)

123 y = sigmoid(x, *sigmoid_parameters[3:])

124 dy = np.diff(y)

125 dy_norm = dy / max(dy)

126 idx = np.argmin(

127 np.abs(np.expand_dims(edges_sig, axis=1) - np.expand_dims(y, axis=0)), axis=1

128 )

129 edges_x = x[idx]

130

131 axes[1].plot(x, y, color="C1", label="sigmoid")

132

133 axes[1].legend()

134 axes[1].set_xlabel("distance")

135 axes[1].set_ylim((0, 1))

136 axes[1].set_title("low-d")

137 for i in range(len(edges)):

138 if edges_x[i] != edges_x[-1]:

139 axes[1].annotate(

140 "",

141 xy=(edges[i], 0),

142 xytext=(edges_x[i], 0),

143 xycoords=axes[0].transData,

144 textcoords=axes[1].transData,

145 arrowprops=dict(facecolor="black", arrowstyle="-", clip_on=False),

146 )

147 axes[0].figure.tight_layout()

148 return axes[0], axe2, axes[1]

149

150

151def _zoomingBoxManual(ax1, ax2, color="red", linewidth=2, roiKwargs={}, arrowKwargs={}):

152 """Fakes a zoom effect between two mpl.axes.Axes.

153

154 Uses mpl.patches.ConnectionPatch and mpl.patches.Rectangle

155 to make it seem like ax2 is a zoomed in version of ax1.

156 Instead of defining the coordinates of the zooming rectangle

157 The axes limits of ax2 are used.

158

159 Args:

160 ax1 (plt.axes): The axes with the zoomed-out data.

161 ax2 (plt.axes): The second axes with the zoomed-in data.

162 color (str): The color of the zoom effect. Is passed into mpl,

163 thus can be str, or tuple, ... Defaults to 'red'

164 linewidth (int): The linewidth. Defaults to 2.

165 roiKwargs (dict): Keyworded arguments for the rectangle.

166 Defaults to {}.

167 arrowKwargs (dict): Keyworded arguments for the arrow.

168 Defaults to {}.

169

170 """

171 limits = np.array([*ax2.get_xlim(), *ax2.get_ylim()])

172 roi = limits

173 roiKwargs = dict(

174 dict(

175 [

176 ("fill", False),

177 ("linestyle", "dashed"),

178 ("color", color),

179 ("linewidth", linewidth),

180 ]

181 ),

182 **roiKwargs,

183 )

184 ax1.add_patch(

185 mpl.patches.Rectangle(

186 [roi[0], roi[2]], roi[1] - roi[0], roi[3] - roi[2], **roiKwargs

187 )

188 )

189 arrowKwargs = dict(

190 dict([("arrowstyle", "-"), ("color", color), ("linewidth", linewidth)]),

191 **arrowKwargs,

192 )

193 corners = np.vstack([limits[[0, 1, 1, 0]], limits[[2, 2, 3, 3]]]).T

194 con1 = mpl.patches.ConnectionPatch(

195 xyA=corners[0],

196 xyB=corners[1],

197 coordsA="data",

198 coordsB="data",

199 axesA=ax2,

200 axesB=ax1,

201 )

202 con1.set_color([0, 0, 0])

203 ax2.add_artist(con1)

204 con1.set_linewidth(2)

205 con2 = mpl.patches.ConnectionPatch(

206 xyA=corners[3],

207 xyB=corners[2],

208 coordsA="data",

209 coordsB="data",

210 axesA=ax2,

211 axesB=ax1,

212 )

213 con2.set_color([0, 0, 0])

214 ax2.add_artist(con2)

215 con2.set_linewidth(2)

216

217

218def render_vmd(

219 filepath,

220 rotation=[0, 0, 0],

221 scale=1,

222 script_location="auto",

223 image_location="auto",

224 debug=False,

225 image_name="",

226 drawframes=False,

227 ssupdate=True,

228 renderer="tachyon",

229 additional_spheres=[],

230 additional_lines=[],

231 surf=None,

232 custom_script=None,

233):

234 """Render pdb file with combination of vmd, tachyon and image magick.

235

236 This function creates a standardised vmd tcl/tk script and writes it

237 to disk. Then vmd is called with the subprocess package and used to

238 create a tachyon input file. Tachyon is then called to render the image

239 with ambient occlusion and soft lighting. The output is a high quality

240 targa (.tga) image, which is converted to png using image magick.

241

242 Args:

243 filepath (str): Location of the pdb file which should be rendered.

244 rotation ([x_rot, y_rot, z_rot], optional): List of rotation values. Defaults to [0, 0, 0].

245 scale (float, optional): By how much the structure should be scaled. Defaults to 1.

246 script_location (str, optional): Where to save the script. Script will be removed

247 after finish nonehteless. Defaults to 'auto' and writes into cwd.

248 image_location (str, optional): Where to render the images file to. Will be

249 deleted nonetheless. Don't give an extension for this. Defaults to 'auto' and

250 writes into cwd.

251 debug (bool, optional): Print debug info. Defaults to False.

252 image_name (str, optional): This string will be used to save the image to after it has

253 been rendered and converted to png. This will not be deleted. Defaults to ''.

254 drawframes (bool, optional): If a trajectory is loaded, this will render all frames in it.

255 Defaults to False.

256 ssupdate (bool, optional): Updates the secondary structure for every frame. Normally

257 vmd uses the secondary structure of the first frame. Setting this to True calcs

258 the sec struct for every frame. Defaults to True.

259 renderer (str, optional): Which renderer to use.

260 * 'tachyon' uses the external Tachyon rendered. So vmd -> .dat -> .tga -> .png.

261 * 'snapshot' uses the vmd internal snapshot renderer.

262 Defaults to 'tachyon'.

263 additional_spheres (list, optional): Draw spheres around two subunits to make

264 them visually more distinct. Takes a list of lists. Each list in the main

265 list should contain 4 values [x, y, z, r] (r for radius). Defaults to [].

266 additional_lines (list, optional): A list of additional lines that should be added to the

267 script. Please refert to the vmd manual for further info. Defaults to [].

268 surf (Union[str, None], optional): A string defining the surface renderer. Can either be

269 'quicksurf' or 'surf'. If None is provided, the surface won't be rendered (falls back

270 to cartoon representation). Defaults to None.

271 custom_script (Union[str, None], optional): Provide a completely custom script as this option.

272 The render commands will still be appended to this script. If None is provided, the

273 default script will be used.

274

275 See also:

276 See this nice webpage about rendering publication worthy images with vmd.

277 https://www.ks.uiuc.edu/Research/vmd/minitutorials/tachyonao/

278

279 Returns:

280 image (np.ndarray): This array contains the raw pixel data.

281 Can be used with matplotlib to have a quick view of the image.

282

283 Examples:

284 >>> pdb_file = '/path/to/pdbfile.pdb'

285 >>> image = render_vmd(pdb_file, scale=2)

286 >>> plt.imshow(image)

287

288 """

289 if "." in image_location: 289 ↛ 290line 289 didn't jump to line 290, because the condition on line 289 was never true

290 raise Exception(

291 "The argument image_location does not take a file extension, because the name is used for a .dat, .tga and .png file."

292 )

293

294 # add a shebang to the script

295 # script = '#!/home/soft/bin/vmd\n\n'

296

297 # print debug hello world

298 script = 'puts "Hello World"\n'

299

300 # if a list of files is provided we iterate over them

301 if isinstance(filepath, list): 301 ↛ 302line 301 didn't jump to line 302, because the condition on line 301 was never true

302 for i, file in enumerate(filepath):

303 # load molecule and change representation

304 script += f"mol new {file}\n"

305 if surf is None:

306 script += f"mol modstyle 0 {i} newcartoon 0.3 50\n"

307 script += f"mol modcolor 0 {i} structure\n"

308 elif surf == "quicksurf":

309 script += f"mol modstyle 0 {i} quicksurf 0.6 0.7 0.7 Medium\n"

310 else:

311 script += f"mol modstyle 0 {i} {surf}\n"

312 script += f"mol modmaterial 0 {i} AOChalky\n"

313 if drawframes and md.load(file).n_frames > 1:

314 if renderer == "STL":

315 import warnings

316

317 warnings.warn(

318 "Rendering multiple frames with STL may lead to undesired results. Instead of yielding the union of all single-frame surfaces, you will get a mishmash of all surfaces with intersection faces etc."

319 )

320 script += f"mol drawframes 0 {i} 0:1:999\n"

321 else:

322 # load molecule and change representation

323 script += f"mol new {filepath}\n"

324 if surf is None: 324 ↛ 327line 324 didn't jump to line 327, because the condition on line 324 was never false

325 script += "mol modstyle 0 0 newcartoon 0.3 50\n"

326 script += "mol modcolor 0 0 structure\n"

327 elif surf == "quicksurf":

328 script += "mol modstyle 0 0 quicksurf 0.6 0.7 0.7 Medium\n"

329 else:

330 script += f"mol modstyle 0 0 {surf}\n"

331 script += "mol modmaterial 0 0 AOChalky\n"

332 if drawframes: 332 ↛ 335line 332 didn't jump to line 335, because the condition on line 332 was never false

333 script += "mol drawframes 0 0 0:1:999\n"

334

335 if ssupdate: 335 ↛ 359line 335 didn't jump to line 359, because the condition on line 335 was never false

336 print(

337 "\033[93m"

338 + "For the ssupdate function to work encodermap/vmd/sscache.tcl will be sourced within vmd. If no Error is thrown the file is present."

339 + "\033[0m"

340 )

341 sscache_location = (

342 os.path.split(os.path.split(os.path.split(__file__)[0])[0])[0]

343 + "/vmd/sscache.tcl"

344 )

345 if not os.path.isfile(sscache_location): 345 ↛ 346line 345 didn't jump to line 346, because the condition on line 345 was never true

346 raise FileNotFoundError(

347 f"The sscache.tcl script is not here: {sscache_location}. Please put it there."

348 )

349 script += f"source {sscache_location}\n"

350 script += "start_sscache 0\n"

351 # script += "proc update_secondary_structure_assigment { args } {"

352 # script += " foreach molid [molinfo list] {"

353 # script += " mol ssrecalc $molid"

354 # script += " }"

355 # script += "}"

356 # script += "trace variable vmd_frame(0) w update_secondary_structure_assigment"

357

358 # change some parameters to make a nice image

359 script += "color Display Background white\n"

360 script += "color Axes Labels black\n"

361 script += "display depthcue off\n"

362 script += "display ambientocclusion on\n"

363 script += "display aoambient 1.0\n"

364 script += "display aodirect 0.3\n"

365 script += "display antialias on\n"

366 # script += 'display resize 2000 2000\n'

367 script += "axes location off\n"

368

369 # scale and rotate

370 script += f"rotate x by {rotation[0]}\n"

371 script += f"rotate y by {rotation[1]}\n"

372 script += f"rotate z by {rotation[2]}\n"

373 script += f"scale by {scale}\n"

374

375 # define image location

376 if image_location == "auto": 376 ↛ 380line 376 didn't jump to line 380, because the condition on line 376 was never false

377 image_location = os.getcwd() + "/vmdscene"

378

379 # add spheres

380 if np.any(additional_spheres): 380 ↛ 381line 380 didn't jump to line 381, because the condition on line 380 was never true

381 for _, color in zip(additional_spheres, ["grey", "iceblue"]):

382 x, y, z, r = np.round(_, 2)

383 script += f"draw color {color}\n"

384 script += f"draw sphere {{ {x} {y} {z} }} radius {r} resolution 25\n"

385 script += "draw material Transparent\n"

386

387 # add additional lines

388 if additional_lines: 388 ↛ 389line 388 didn't jump to line 389, because the condition on line 388 was never true

389 for line in additional_lines:

390 script += line + "\n"

391

392 if custom_script is not None: 392 ↛ 393line 392 didn't jump to line 393, because the condition on line 392 was never true

393 script = custom_script

394

395 # render command. Alternatively, I can use external Tachyon, which makes better images

396 if renderer == "tachyon": 396 ↛ 398line 396 didn't jump to line 398, because the condition on line 396 was never false

397 script += f"render Tachyon {image_location}.dat\n"

398 elif renderer == "snapshot":

399 script += "render aasamples TachyonInternal 6\n"

400 script += f"render TachyonInternal {image_location}.tga\n"

401 elif renderer == "STL":

402 script += "axes location off\n"

403 script += f"render STL {image_location}.stl\n"

404 elif renderer == "Wavefront":

405 script += "axes location off\n"

406 script += f"render Wavefront {image_location}.obj\n"

407 else:

408 raise NotImplementedError(

409 "Other renderers than tachyon and snaphsot currently not supported."

410 )

411

412 # list molecules and quit

413 script += "mol list\n"

414 script += "quit"

415

416 if debug: 416 ↛ 417line 416 didn't jump to line 417, because the condition on line 416 was never true

417 print(script)

418

419 # write the script

420 if script_location == "auto": 420 ↛ 422line 420 didn't jump to line 422, because the condition on line 420 was never false

421 script_location = os.getcwd() + "/vmd_script.tcl"

422 with open(script_location, "w") as f:

423 f.write(script)

424

425 # call vmd -e script

426 cmd = f"vmd -e {script_location} -dispdev none"

427 if debug: 427 ↛ 428line 427 didn't jump to line 428, because the condition on line 427 was never true

428 print(cmd)

429 proc = subprocess.Popen(

430 cmd,

431 stdin=subprocess.PIPE,

432 stdout=subprocess.PIPE,

433 stderr=subprocess.PIPE,

434 shell=True,

435 )

436 (stdout, stderr) = proc.communicate()

437 if debug: 437 ↛ 438line 437 didn't jump to line 438, because the condition on line 437 was never true

438 print(stdout.decode("utf-8"))

439 print("\n")

440 print(stderr.decode("utf-8"))

441

442 # check if image has been written

443 if renderer == "tachyon": 443 ↛ 448line 443 didn't jump to line 448, because the condition on line 443 was never false

444 assert os.path.isfile(

445 f"{image_location}.dat"

446 ), "Tachyon datafile not generated by renderer"

447 else:

448 assert os.path.isfile(

449 f"{image_location}.tga"

450 ), f"Snapshot image not created. {stderr.decode()} {stdout.decode()}"

451

452 time.sleep(2)

453 assert os.path.isfile(

454 f"{image_location}.tga"

455 ), f"Tachyon datafile not generated by renderer. Here's the script:\n\n{script}\n\n"

456

457 if renderer == "tachyon":

458 # call Tachyon and render

459 cmd = f"/usr/bin/tachyon -aasamples 12 {image_location}.dat -res 2000 2000 -fullshade -format TARGA -o {image_location}.tga"

460 if debug:

461 print(cmd)

462 proc = subprocess.Popen(

463 cmd,

464 stdin=subprocess.PIPE,

465 stdout=subprocess.PIPE,

466 stderr=subprocess.PIPE,

467 shell=True,

468 )

469 (stdout, stderr) = proc.communicate()

470 if debug:

471 print(stdout.decode("utf-8"))

472 print("\n")

473 print(stderr.decode("utf-8"))

474

475 # check if image has been written

476 assert os.path.isfile(

477 f"{image_location}.tga"

478 ), "Tachyon renderer did not render image"

479

480 if renderer == "STL":

481 if image_name:

482 shutil.copyfile(f"{image_location}.stl", image_name)

483 import trimesh

484

485 mesh = trimesh.load(f"{image_location}.stl")

486 os.remove(f"{image_location}.stl")

487 return mesh

488

489 if renderer == "Wavefront":

490 if image_name:

491 shutil.copyfile(f"{image_location}.obj", image_name)

492 shutil.copyfile(f"{image_location}.mtl", image_name.replace(".obj", ".mtl"))

493 print(

494 f"Find the rendered images at {image_name} and {image_name.replace('.obj', '.mtl')}."

495 )

496 return None

497

498 # convert to png

499 cmd = f"/usr/bin/convert {image_location}.tga {image_location}.png"

500 if debug:

501 print(cmd)

502 proc = subprocess.Popen(

503 cmd,

504 stdin=subprocess.PIPE,

505 stdout=subprocess.PIPE,

506 stderr=subprocess.PIPE,

507 shell=True,

508 )

509 (stdout, stderr) = proc.communicate()

510 if debug:

511 print(stdout.decode("utf-8"))

512 print("\n")

513 print(stderr.decode("utf-8"))

514

515 # read image

516 image = plt.imread(f"{image_location}.png")

517

518 # write image if name has been provided

519 if image_name:

520 if os.path.isabs(image_name):

521 shutil.copyfile(f"{image_location}.png", image_name)

522 else:

523 shutil.copyfile(f"{image_location}.png", os.getcwd() + f"/{image_name}")

524

525 # remove temporary files

526 if renderer == "tachyon":

527 os.remove(f"{image_location}.dat")

528 os.remove(f"{image_location}.tga")

529 os.remove(f"{image_location}.png")

530 # os.remove(f'{script_location}')

531

532 # return matplotlib image object

533 return image

534

535

536def render_movie(path, scatter_data, dummy_traj):

537 pass

538

539

540def plot_cluster(

541 trajs, pdb_path, png_path, cluster_no=None, col="user_selected_points"

542):

543 from mpl_toolkits.axes_grid1 import make_axes_locatable

544

545 fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(ncols=2, nrows=2)

546 fig.set_size_inches(20, 20)

547

548 if cluster_no is None: 548 ↛ 549line 548 didn't jump to line 549, because the condition on line 548 was never true

549 cluster_no = trajs.CVs[col].max()

550

551 # prepare ax1 to make the two side histograms

552 divider = make_axes_locatable(ax4)

553 axHistx = divider.append_axes("top", size=1.2, pad=0.1) # , sharex=ax1)

554 axHisty = divider.append_axes("right", size=1.2, pad=0.1) # , sharey=ax1)

555

556 # some data management

557 data = trajs.lowd

558 where = np.where(trajs.CVs[col] == cluster_no)

559 not_where = np.where(trajs.CVs[col] != cluster_no)

560 x = data[:, 0]

561 y = data[:, 1]

562

563 # scatter everything grey and cluster blue

564 ax1.scatter(*data[where].T)

565 ax1.scatter(*data[not_where].T, c="grey", s=5)

566 ax1.set_xlabel("x in a.u.")

567 ax1.set_ylabel("y in a.u.")

568 ax1.set_title(f"Scatter of low-dimensional data")

569

570 # density

571 bin_density = 46

572 log_density = True

573

574 # ax2 gets hexbin density

575 # x_bins = np.linspace(x.min(), x.max(), bin_density)

576 # y_bins = np.linspace(y.min(), y.max(), bin_density)

577 H, xedges, yedges = np.histogram2d(x=x, y=y, bins=bin_density)

578 extent = [xedges[0], xedges[-1], yedges[0], yedges[-1]]

579 xcenters = np.mean(np.vstack([xedges[0:-1], xedges[1:]]), axis=0)

580 ycenters = np.mean(np.vstack([yedges[0:-1], yedges[1:]]), axis=0)

581 X, Y = np.meshgrid(xcenters, ycenters)

582 extent = [xedges[0], xedges[-1], yedges[0], yedges[-1]]

583 if log_density: 583 ↛ 587line 583 didn't jump to line 587, because the condition on line 583 was never false

584 with np.errstate(divide="ignore"): # ignore division by zero error

585 F = np.log(H)

586 else:

587 F = H

588 # mappable = ax2.hexbin(x=X.ravel(), y=Y.ravel(), C=F.T.ravel(), cmap=plt.cm.turbo_r, extent=extent,

589 # norm=mpl.colors.PowerNorm(1), gridsize=bin_density +1)

590 cmap = plt.get_cmap("turbo").with_extremes(under="w")

591 mappable = ax2.contourf(

592 X, Y, H.T, cmap=cmap, levels=np.linspace(0.001, H.max(), 20)

593 )

594 ax2.set_xlabel("x in a.u.")

595 ax2.set_ylabel("y in a.u.")

596 ax2.set_title("Log density of points")

597

598 # colorbar for ax2

599 # colorbar

600 # use the axes divider method to add colorbar

601 ax_divider = make_axes_locatable(ax2)

602 # add colorbaraxis to work with ticks and whatnot

603 cax = ax_divider.append_axes("right", size="7%", pad="2%")

604 # define colorbar norm. I like to work with values between 0 and 1

605 # initialize colormap

606 cb = plt.colorbar(mappable, cax=cax)

607 cax.set_ylabel("Number of points")

608

609 # cluster on ax4

610 # x hist

611 spines = [k for k in axHistx.spines.values()]

612 spines[1].set_linewidth(0)

613 spines[3].set_linewidth(0)

614 axHistx.set_xticks([])

615 H, edges, patches = axHistx.hist(data[:, 0][where], bins=50)

616 axHistx.set_ylabel("count")

617 axHistx.set_title("Scatter of Cluster")

618

619 # y hist

620 spines = [k for k in axHisty.spines.values()]

621 spines[1].set_linewidth(0)

622 spines[3].set_linewidth(0)

623 axHisty.set_yticks([])

624 H, edges, patches = axHisty.hist(

625 data[:, 1][where], bins=50, orientation="horizontal"

626 )

627 axHisty.set_xlabel("count")

628

629 # scatter data

630 ax4.scatter(x=data[where, 0], y=data[where, 1])

631 spines = [k for k in ax4.spines.values()]

632 spines[3].set_linewidth(0)

633 spines[1].set_linewidth(0)

634 ax4.set_xlabel("x in a.u.")

635 ax4.set_ylabel("y in a.u.")

636

637 # annotate rms

638 rms = np.sqrt(

639 (1 / len(data[where]))

640 * np.sum(

641 (data[where, 0] - np.mean(data[where, 0])) ** 2

642 + (data[where, 1] - np.mean(data[where, 1])) ** 2

643 )

644 )

645 text = f"RMS = {np.round(rms, decimals=5)}"

646 ax4.text(0.05, 0.95, text, transform=ax1.transAxes)

647

648 # annotate geometric center

649 centroid = [np.mean(x[where]), np.mean(y[where])]

650 ax4.scatter(*centroid, s=50, c="C1")

651 ax4.annotate(

652 "geom. center",

653 xy=centroid,

654 xycoords="data",

655 xytext=(0.95, 0.95),

656 textcoords="axes fraction",

657 arrowprops=dict(facecolor="black", shrink=0.05, fc="w", ec="k", lw=2),

658 horizontalalignment="right",

659 verticalalignment="top",

660 color="C1",

661 )

662

663 # annotate rmsd center

664 # view, dummy_traj = gen_dummy_traj(trajs, cluster_no, max_frames=100, col=col)

665 # index, distances, centroid = rmsd_centroid_of_cluster(dummy_traj, parallel=False)

666 # idx = np.round(np.linspace(0, len(where) - 1, 100)).astype(int)

667 # where = where[idx]

668 # centroid = data[where[0][::5][index]]

669 # ax4.scatter(*centroid, s=50, c='C2')

670 # ax4.annotate('rmsd center',

671 # xy=centroid, xycoords='data',

672 # xytext=(0.95, 0.85), textcoords='axes fraction',

673 # arrowprops=dict(facecolor='black', shrink=0.05, fc="w", ec="k", lw=2),

674 # horizontalalignment='right', verticalalignment='top', color='C2')

675

676 # make vmd snapshot

677 try:

678 image = render_vmd(

679 pdb_path, drawframes=True, renderer="tachyon", debug=False, scale=1.5

680 )

681 ax3.imshow(image)

682 [k.set_linewidth(0) for k in ax3.spines.values()]

683 ax3.set_xticks([])

684 ax3.set_yticks([])

685 ax3.set_title("Image of cluster")

686 except:

687 ax3.annotate("VMD Rendering not possible", (0.5, 0.5))

688 pass

689

690 # # calculate distances between rmsd centroid and all other points

691 # distances = scipy.spatial.distance.cdist(centroid.reshape(1, 2), np.stack([x, y]).T)

692 # H, edges, patches = ax3.hist(distances.flatten(), color='C1')

693 # ax3.set_title("Distances to rmsd centroid.")

694 # ax3.set_xlabel("Distance in a.u.")

695 # ax3.set_ylabel("Count")

696

697 plt.suptitle(f"Cluster {cluster_no}")

698 plt.savefig(png_path, transparent=False)

699 plt.close(fig)