From b94386e57353dbd6aba74e8c605647091dc85403 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 16 Sep 2013 10:16:02 +0200 Subject: Add Rsplit_surface and clean-stream.py --- scripts/Rsplit_surface | 62 +++++++++++++++++++++++++++ scripts/Rsplit_surface.py | 105 ++++++++++++++++++++++++++++++++++++++++++++++ scripts/clean-stream.py | 87 ++++++++++++++++++++++++++++++++++++++ 3 files changed, 254 insertions(+) create mode 100755 scripts/Rsplit_surface create mode 100644 scripts/Rsplit_surface.py create mode 100644 scripts/clean-stream.py (limited to 'scripts') diff --git a/scripts/Rsplit_surface b/scripts/Rsplit_surface new file mode 100755 index 00000000..3cd22d67 --- /dev/null +++ b/scripts/Rsplit_surface @@ -0,0 +1,62 @@ +#!/bin/sh + +# Rsplit_surface +# +# Plot Rsplit as a contour map +# +# Copyright © 2013 Fedor Chervinskii +# Copyright © 2013 Deutsches Elektronen-Synchrotron DESY, +# a research centre of the Helmholtz Association. +# +# Authors: +# 2013 Fedor Chervinskii +# 2013 Thomas White +# +# This file is part of CrystFEL. +# +# CrystFEL is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# CrystFEL is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with CrystFEL. If not, see . + +filename=$1 +pdbfilename=$2 +pg=$3 + +if [ x$filename = "x" -o x$pdbfilename = "x" -o x$pg = "x" ]; then + echo "Syntax: Rsplit_surface my.stream my.pdb " + exit 1 +fi + +python clean-stream.py $filename clean.stream +filename=clean.stream + +LIMIT=3 + +for ((a=1; a <= LIMIT ; a++)) +do + let string=$a + ./alternate-stream $filename ${filename/./1.} ${filename/./2.} + if (($a > 1)); then + rm $filename + fi + process_hkl -i ${filename/./1.} -o ${filename/.stream/1.hkl} -y $pg + process_hkl -i ${filename/./2.} -o ${filename/.stream/2.hkl} -y $pg + compare_hkl ${filename/.stream/1.hkl} ${filename/.stream/2.hkl} -y $pg -p ${pdbfilename} --fom=rsplit --shell-file=rsplit$string.dat + rm ${filename/./2.} + rm ${filename/.stream/1.hkl} + rm ${filename/.stream/2.hkl} + filename=${filename/.stream/1.stream} +done + +rm $filename + +python Rsplit_surface.py $1 diff --git a/scripts/Rsplit_surface.py b/scripts/Rsplit_surface.py new file mode 100644 index 00000000..ae4e01ac --- /dev/null +++ b/scripts/Rsplit_surface.py @@ -0,0 +1,105 @@ +#!/usr/bin/python +# coding=utf-8 + +from __future__ import division +from array import array +from numpy import * + +import matplotlib + +# Rsplit_surface.py +# +# Plot Rsplit as a contour map +# +# Copyright © 2013 Fedor Chervinskii +# +# Authors: +# 2013 Fedor Chervinskii +# +# This file is part of CrystFEL. +# +# CrystFEL is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# CrystFEL is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with CrystFEL. If not, see . + + +"""there could be problems with dependencies, in this case, """ +#matplotlib.use('PS') + +import matplotlib.pyplot as plt +from matplotlib import cm +from matplotlib.mlab import griddata +from mpl_toolkits.axes_grid1 import host_subplot +import mpl_toolkits.axisartist as AA +import sys + +plt.rc('text', usetex=True) +plt.rc('font', family='serif') + +data_name = sys.argv[1] + +Nfile = open('N.dat', 'r') +line = array(Nfile.readlines())[0] +N0 = eval(line.split()[0]) +steps = arange(1, 10) + +verts = [] + +for i in steps : + + filename='rsplit'+str(i)+'.dat' + + file = open(filename, 'r') + + lines = array(file.readlines()) + N = lines.size + + for j in range(1, N) : + numbers = lines[j].split() + if 0 < eval(numbers[1]) <= 100 : verts.append((eval(numbers[0]),i,eval(numbers[1]))) + +x, y, z = zip(*verts) + +fig = plt.figure() + +ax = host_subplot(111, axes_class=AA.Axes) + +X = arange(1, 5, 0.05) +Y = arange(0, 10, 0.1) +X, Y = meshgrid(X, Y) +Z = griddata(x,y,z,X,Y) + +cm = ax.pcolormesh(X, Y, Z, cmap = plt.get_cmap('spectral'), vmin=0, vmax=100) + +plt.title(r"$R_{split}$ surface for %s data $N_0$ = %d" % (data_name, N0), fontsize=16) + +ax.set_xlabel(r"1/d", fontsize=16) +ax.set_ylabel(r"Number of patterns", fontsize=16) +ypoints = [50000,40000,30000,20000,10000,9000,8000,7000,6000,5000,4000,3000,2000,1000,900,800,700,600,500,400,300,200,100] +temp = [] +for ypoint in ypoints : temp.append(N0/ypoint) +ax.set_yticks(log2(temp)+1) +labels = [] +for ypoint in ypoints : + if ypoint in [50000,10000,5000,1000,500] : labels.append(str(ypoint)) + else : labels.append("") + +ax.set_yticklabels(labels) + +ax.set_ylim(log2(temp[0])+1,log2(temp[-1])+1) + +fig.colorbar(cm, shrink=0.5, aspect=10) + +plt.savefig(data_name[:-7] + '.png') +#plt.savefig(data_name[:-7] + '.ps') +plt.show() + diff --git a/scripts/clean-stream.py b/scripts/clean-stream.py new file mode 100644 index 00000000..0982f455 --- /dev/null +++ b/scripts/clean-stream.py @@ -0,0 +1,87 @@ +#!/usr/bin/python +# coding=utf-8 + +# clean-stream.py +# +# Remove non-indexed frames from a stream +# +# Copyright © 2013 Fedor Chervinskii +# +# Authors: +# 2013 Fedor Chervinskii +# +# This file is part of CrystFEL. +# +# CrystFEL is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# CrystFEL is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with CrystFEL. If not, see . + + +from __future__ import division +from itertools import islice +import sys +import re + + +infile_1 = open (sys.argv[1],'r') +infile_2 = open (sys.argv[1],'r') +outfile = open (sys.argv[2],'w') +Nfile = open ('N.dat','w') + +suited = False +counter1 = -1 +counter2 = 0 +start_index = 0 +n_patt = 0 +num_peaks = 0 +indexed = False +num_suited = 0 + +line = infile_1.readline() + +reg0 = re.compile('----- Begin chunk') +reg1 = re.compile('----- End chunk') +reg6 = re.compile('^indexed_by') +reg7 = re.compile('none') + +while (line != ''): + + counter1 += 1 + + if (reg6.match(line)): + indexed = ( False if reg7.search(line) else True) + + if (reg0.match(line)): + if (start_index == 0) : + while (counter2 < counter1) : + outline = infile_2.readline() + outfile.write(outline) + counter2 += 1 + start_index = counter1 + + if (reg1.match(line)): + n_patt += 1 + if indexed : + suited = True + if suited : + num_suited += 1 + while (counter2 <= counter1) : + outline = infile_2.readline() + if suited : + outfile.write(outline) + counter2 += 1 + suited = False + + line = infile_1.readline() + +print '%d suited of %d patterns have been extracted and saved as %s' % (num_suited, n_patt, sys.argv[2]) +Nfile.write('%d' % num_suited) -- cgit v1.2.3