From 5514263bb9bf0584b636215714e8381f7e6de21e Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 26 Jan 2018 09:29:47 +0100 Subject: Add scripts/plot-pr and scripts/plot-contourmap --- scripts/plot-contourmap | 153 +++++++++++++++++++++++++++++++++++++ scripts/plot-pr | 196 ++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 349 insertions(+) create mode 100755 scripts/plot-contourmap create mode 100755 scripts/plot-pr (limited to 'scripts') diff --git a/scripts/plot-contourmap b/scripts/plot-contourmap new file mode 100755 index 00000000..053a6be9 --- /dev/null +++ b/scripts/plot-contourmap @@ -0,0 +1,153 @@ +#!/usr/bin/env python + +import numpy as np +import matplotlib +import matplotlib.cm as cm +import matplotlib.mlab as mlab +import matplotlib.pyplot as plt +import sys +import fnmatch +import os +from matplotlib.widgets import RadioButtons,Button + +im=None +cnt=None +centre_marker=None + +def next_click(w): + global crystal + crystal += 20 + print("Crystal %i" % crystal) + update_graph() + + +def prev_click(w): + global crystal + crystal -= 20 + print("Crystal %i" % crystal) + update_graph() + + +def iteration_click(label): + global cycle + cycle = label.split(" ")[1].split("\n")[0] + update_graph() + + +def variable_click(label): + global varpair + varpair = label + update_graph() + + +def update_graph(): + + global im, cnt, centre_marker + + filename="pr-logs/grid-crystal%i-cycle%s-%s.dat" % (crystal,cycle,varpair) + print filename + + with open(filename, "r") as f: + + l = f.readline().split(None,3) + + min1 = float(l[0]) + max1 = float(l[1]) + cur1 = float(l[2]) + label1 = l[3] + + l = f.readline().split(None,3) + min2 = float(l[0]) + max2 = float(l[1]) + cur2 = float(l[2]) + label2 = l[3] + + extent = (min1, max1, min2, max2) + Z = np.loadtxt(fname=filename, ndmin=2, delimiter=" ", skiprows=2) + + ax.set_xlim([min1,max1]) + ax.set_ylim([min2,max2]) + + if not im: + im = ax.imshow(Z, interpolation='none', origin='lower', + cmap=cm.gray, extent=extent, aspect='auto') + im.autoscale() + else: + im.set_data(Z) + im.set_extent(extent) + im.autoscale() + + levels = np.arange(0.0, 1.6, 0.1) + + if cnt: + for coll in cnt.collections: + coll.remove() + + cnt = ax.contour(Z, levels, origin='lower', linewidths=1, extent=extent) + + ax.set_title(filename) + ax.set_xlabel(label1) + ax.set_ylabel(label2) + plt.flag() + + if centre_marker: + centre_marker.remove() + centre_marker, = plt.plot(cur1, cur2, 'bH', color='r') + + fig.canvas.draw() + +fig = plt.figure(figsize=(10,5)) +fig.subplots_adjust(left=0.05, bottom=0.05, right=0.70, top=0.95) + +# Find out what there is to plot +crystals = [] +cycles = [] +varpairs = [] +for file in os.listdir("pr-logs"): + if not fnmatch.fnmatch(file, "grid-*.dat"): + continue + sp = file.rstrip(".dat").split("-") + crystals.append(int(sp[1].lstrip("crystal"))) + cycles.append(sp[2].lstrip("cycle")) + varpairs.append(sp[3]+"-"+sp[4]) + +crystals = sorted(set(crystals)) +cycles = sorted(set(cycles)) +varpairs = sorted(set(varpairs)) + +crystal = crystals[0] +cycle = cycles[0] +varpair = varpairs[0] + +# Iteration selector +ax = plt.axes([0.75, 0.55, 0.2, 0.40], facecolor="lightgoldenrodyellow") +iterations = ["Iteration "+str(f) for f in cycles] +if iterations[0] != "Iteration 0": + print("Whoops, couldn't find iteration 0!") + sys.exit(1) +iterations[0] = "Iteration 0\n(initial scaling only)" +if iterations[len(iterations)-1] != "Iteration F": + print("Whoops, couldn't find final iteration!") +else: + iterations[len(iterations)-1] = "Iteration F\n(after final merge)" +iteration = RadioButtons(ax, iterations) +iteration.on_clicked(iteration_click) + +# Variable selector +ax = plt.axes([0.75, 0.20, 0.2, 0.30], facecolor="lightgoldenrodyellow") +variable = RadioButtons(ax, varpairs) +variable.on_clicked(variable_click) + +# Crystal selector +ax = plt.axes([0.75, 0.08, 0.20, 0.06]) +crystal_prev = Button(ax, "Previous crystal") +crystal_prev.on_clicked(prev_click) +ax = plt.axes([0.75, 0.01, 0.20, 0.06]) +crystal_next = Button(ax, "Next crystal") +crystal_next.on_clicked(next_click) + +ax = plt.axes([0.1, 0.1, 0.6, 0.8]) +update_graph() +plt.colorbar(im, shrink=0.8) + +plt.show() diff --git a/scripts/plot-pr b/scripts/plot-pr new file mode 100755 index 00000000..25427cae --- /dev/null +++ b/scripts/plot-pr @@ -0,0 +1,196 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# +# Plot post-refinement results +# +# Copyright © 2017-2018 Deutsches Elektronen-Synchrotron DESY, +# a research centre of the Helmholtz Association. +# +# Author: +# 2017-2018 Thomas White +# + +import sys +import matplotlib +import numpy as np +from matplotlib import pyplot as plt +from operator import truediv, itemgetter +from matplotlib.widgets import RadioButtons,SpanSelector,Button + +crystal=0 +resmin=0 +resmax=100e9 +inum="0" +first=True + +def update_graph(): + + global inum, resmin, resmax, first + + # Update pobs/pcalc scatter plot + xlim = pgraphE.get_xlim() + ylim = pgraphE.get_ylim() + pgraphE.cla() + pgraphE.scatter([x[0] for x in pgraph if x[2]==inum if x[3]resmin], + [x[1] for x in pgraph if x[2]==inum if x[3]resmin], marker="+") + pgraphE.grid(b=True, which='major', color='k', linestyle='--') + if not first: + pgraphE.set_xlim(xlim) + pgraphE.set_ylim(ylim) + pgraphE.set_xlabel("pcalc") + pgraphE.set_ylabel("pobs") + pgraphE.set_title("Observed and calculated partialities") + + # Update spectrum plot + pobsE.set_data([x[0] for x in specgraph if x[3]==inum if x[4]resmin], + [x[2] for x in specgraph if x[3]==inum if x[4]resmin]) + pcalcE.set_data([x[0] for x in specgraph if x[3]==inum if x[4]resmin], + [x[1] for x in specgraph if x[3]==inum if x[4]resmin]) + specgraphE.set_title("Spectrum plot for crystal %i" % crystal) + + first = False + plt.draw() + + +def iteration_click(label): + global inum + inum = label.split(" ")[1].split("\n")[0] + update_graph() + + +def resolutionchange(rmin,rmax): + global resmin,resmax + resmin = rmin + resmax = rmax + update_graph() + + +def read_spectrum(filename): + + specgraph = [] + + try: + fh = open(filename, 'r') + except IOError: + print("Failed to read "+filename) + raise IOError + + fline = fh.readline() + fline = fh.readline() # Ignore header (two lines) + + while True: + fline = fh.readline() + if not fline: + break + kval = float(fline.split()[0]) + res = float(fline.split()[1]) + pcalc = float(fline.split()[2]) + pobs = float(fline.split()[3]) + inum = fline.split()[4] + specgraph.append([kval,pcalc,pobs,inum,res]) + + return sorted(specgraph, key=itemgetter(0)) + + +def next_click(w): + global specgraph, crystal + crystal += 20 + specgraph = read_spectrum("pr-logs/specgraph-crystal%i.dat" % crystal) + update_graph() + + +def prev_click(w): + global specgraph, crystal + crystal -= 20 + specgraph = read_spectrum("pr-logs/specgraph-crystal%i.dat" % crystal) + update_graph() + + +def set_sensible_axes(w): + pgraphE.set_xlim([-0.1,1.1]) + pgraphE.set_ylim([-0.2,2.0]) + specgraphE.set_ylim([-5.0,5.0]) + plt.draw() + + +# Read the pgraph file +pgraph = [] + +try: + fh = open("pr-logs/pgraph.dat", 'r') +except IOError: + print("Failed to read "+filename) + raise IOError + +fline = fh.readline() # Ignore header + +while True: + fline = fh.readline() + if not fline: + break + res = float(fline.split()[4]) + Ipart = float(fline.split()[5]) + pc = float(fline.split()[6]) + po = float(fline.split()[7]) + inum = fline.split()[8] + pgraph.append([pc,po,inum,res,Ipart]) + + +# Read the spectrum file +specgraph = read_spectrum("pr-logs/specgraph-crystal%i.dat" % crystal) + +fig = plt.figure(figsize=(15,6)) +fig.subplots_adjust(left=0.05, right=0.65) + +# pobs/pcalc +pgraphE = fig.add_subplot(121) + +# Spectrum +specgraphE = fig.add_subplot(122) +pobsE = specgraphE.plot([0,1],[0,1],label="pobs")[0] +pcalcE = specgraphE.plot([0,1],[0,1],label="pcalc")[0] +specgraphE.set_xlabel("khalf / m^-1") +specgraphE.set_ylabel("Partiality") +specgraphE.grid(b=True, which='major', color='k', linestyle='--') +specgraphE.legend() + +# Iteration selector +ax = plt.axes([0.75, 0.50, 0.2, 0.40], facecolor="lightgoldenrodyellow") +iterations = sorted(set(["Iteration "+str(x[2]) for x in pgraph])) +if iterations[0] != "Iteration 0": + print("Whoops, couldn't find iteration 0!") + sys.exit(1) +iterations[0] = "Iteration 0\n(initial scaling only)" +if iterations[len(iterations)-1] != "Iteration F": + print("Whoops, couldn't find final iteration!") +else: + iterations[len(iterations)-1] = "Iteration F\n(after final merge)" +iteration = RadioButtons(ax, iterations) +iteration.on_clicked(iteration_click) +iteration.set_active(len(iterations)-1) + +# Resolution selector +ax = plt.axes([0.75, 0.25, 0.2, 0.15], facecolor="lightgoldenrodyellow") +ax.hist([x[3] for x in pgraph],weights=[x[4] for x in pgraph],bins=100) +ax.set_xlabel("Resolution / m^-1") +ax.set_ylabel("Frequency") +ax.set_title("Resolution selector") +resolution = SpanSelector(ax, resolutionchange, 'horizontal', span_stays=True, useblit=True) + +# Crystal switcher +ax = plt.axes([0.75, 0.10, 0.09, 0.05]) +crystal_prev = Button(ax, "Previous crystal") +crystal_prev.on_clicked(prev_click) +ax = plt.axes([0.86, 0.10, 0.09, 0.05]) +crystal_next = Button(ax, "Next crystal") +crystal_next.on_clicked(next_click) + +# Set sensible axes +ax = plt.axes([0.75, 0.01, 0.09, 0.05]) +sensible_axes = Button(ax, "Set sensible axes") +sensible_axes.on_clicked(set_sensible_axes) + +specgraphE.relim() +specgraphE.autoscale_view(True,True,True) + +plt.show() -- cgit v1.2.3