Add Rsplit_surface and clean-stream.py

3 files changed, 254 insertions, 0 deletions

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 <fedor.chervinskii@gmail.com>
+#   2013 Thomas White <taw@physics.org>
+#
+# 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 <http://www.gnu.org/licenses/>.
+
+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 <point group>"
+	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 <fedor.chervinskii@gmail.com>
+#
+# 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 <http://www.gnu.org/licenses/>.
+
+
+"""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 <fedor.chervinskii@gmail.com>
+#
+# 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 <http://www.gnu.org/licenses/>.
+
+
+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)