Предположим, у нас есть большое количество практически идентичных
дифрактограмм. Допустим, снимались они на одном и том же
дифрактометре (стало быть, имеем одну и ту же базовую линию), на одной
и той же длине волны (то есть имеем приблизительно одно и то же
положение рефлексов). Также, пусть на картине отсутствует перекрытие
рефлексов (к нему надо подходить всё-таки творчески и не доверять
целиком компьютеру).
Команды в fityk
Команды, которые мы собираемся давать описаны в руководстве программы,
однако программа настолько удобна, что при работе в графическом режиме
автоматически пишет в эхо-области какую команду вы выполнили. Стало
быть, если вы открыли дифракционную картину, то увидите что-то типа:
@0 < '/home/user/path/to/profile.dat'
Более общие способы загрузки файлов, типа использования определенных
колонок, стандартного отклонения и прочие способы работы с данными
смотрите
в официальной документации. Но
общий синтаксис таков:
слот < имя_файла [:xcol:ycol:scol:block ] [доп. опции …]
Если вы удалите базовую линию, то в эхо области выведется сообщение типа
(разумеется, в зависимости от того, какой вид базовой линии вы
выбрали, результат будет отличаться):
Y = y - spline[29.9167, 8.58396, 70.1091, 7.88817](x)
Для выделения активной области служит команда:
a = (x_min < x < x_max)
Для того, чтобы установить нулевое приближение профильной линии пика,
служит команда (в данном случае для функции PseudoVoigt):
%p = guess PseudoVoigt
Как я говорил в предыдущем посте, программа достаточно хорошо
«угадывает» начальные профильные параметры.
Уточнение профиля начинается по команде
fit
Экспортирование данных производится командой
info слот (выражение, ...) > имя_файла
Я использую команду info+ %имя_функции для получения полной
информации о пике.
Итак, полный цикл обработки одной дифракционной картины должен быть
таким:
- Открываем файл с дифрактограммой
- Удаляем базовую линию
- Определяем активную область
- «Угадываем» положение пика, описанного определенной функцией
- Производим уточнение
- Записываем параметры пика в файл
- Удаляем аппроксимирующую функцию
- Переходим к другой активной области и возвращаемся к пункту 3.
- Открываем новый файл с дифрактограммой
Для одного файла на языке cfityk это будет выглядеть так:
@0 < '/home/user/path/to/profile.dat'
Y = y - spline[29.9167, 8.58396, 70.1091, 7.88817](x)
a = (30<x<34)
%p = guess Voigt
fit
info+ %p > ./peaks
delete %p
a = (38<x<41)
%p = guess Voigt
fit
info+ %p >> ./peaks
delete %p
a = (44<x<48)
%p = guess Voigt
fit
info+ %p >> ./peaks
delete %p
a = (56<x<59)
%p = guess Voigt
fit
info+ %p >> ./peaks
delete %p
a = (65<x<70)
%p = guess Voigt
fit
info+ %p >> ./peaks
delete %p
Чтобы сделать это для нескольких файлов, надо написать скрипт.
Примерно год назад, я уже его написал, одну часть на Perl, а другую на
bash с применением sed и awk. Недавно я на него взглянул и
ужаснулся. Не то, чтобы было непонятно, но как-то некрасиво. С Perl
я только лишь познакомился, и, вероятно, сделал всё совсем не так, как
это следует делать. Недавно стал изучать Python, поэтому на нём и
написал скрипт, который будет потоково обрабатывать сразу несколько
файлов и распихивать все выходные данные по файлам.
Пишем скрипт на Python
Для начала определимся с концепцией. Допустим у нас в директории n
дифрактограмм, которые надо обработать совершенно одинаковым
способом. Я предлагаю запихнуть информацию о базовой линии и активных
областях в один файл, его я назвал base_active.list, у меня он
выглядит так:
Y = y - spline[29.9167, 8.58396, 70.1091, 7.88817](x) #baseline
30 33
38 41
44 48
55 59
65 69
То есть, в первой строчке та самая формула, описывающая базовую линию,
а далее x_min x_max, фигурирующие в определении активной области.
Мы должны указать какой функцией аппроксимировать все пики. Я решил
это сделать в интерактивном режиме. То есть после запуска скрипта
появляется сообщение какую функцию выбрать - у меня всё повешено на
цифры.
Итак, вводная часть выглядит так:
#!/usr/bin/env python
#-*-coding: utf-8 -*-
import sys, math, os, re
try:
os.rename('outfile.list', 'outfile.list.bak')
except: pass
try:
os.rename('center.out', 'center.out.bak')
except: pass
outfile = open('outfile.list', 'w')
center_outfile = open('center.out', 'w')
try:
input_file = sys.argv[1]
except:
print "Usage:", sys.argv[0], "input files"
sys.exit(1)
try:
active = open('base_active.list', 'r')
except:
print "Need base_active.list file of baseline formula and active areas"
sys.exit(1)
input_FUNCTION = input("Enter type of function:\n\
1 - Pearson7A\n\
2 - Pearson7\n\
3 - PseudoVoigt\n\
4 - Voigt\n\
5 - Lorentzian\n\
6 - Gaussian\n\
")
if input_FUNCTION == 1:
FUNCTION = 'Pearson7A'
elif input_FUNCTION == 2:
FUNCTION = 'Pearson7'
elif input_FUNCTION == 3:
FUNCTION = 'PseudoVoigt'
elif input_FUNCTION == 4:
FUNCTION = 'Voigt'
elif input_FUNCTION == 5:
FUNCTION = 'Lorentzian'
elif input_FUNCTION == 6:
FUNCTION = 'Gaussian'
else:
print "Inappropriate value"
sys.exit(1)
Первым делом делаем бэкап старых выходных файлов.
try:
os.rename('outfile.list', 'outfile.list.bak')
except: pass
try:
os.rename('center.out', 'center.out.bak')
except: pass
Я сохранял один файл outfile.list, в котором по порядку идут файлы и
параметры их пиков:
20.dat Pearson7A
N center area FWHM shape
0 32.4379 330.0 0.159 1.606
1 40.0144 356.4 0.121 0.481
2 46.5596 177.8 0.169 0.863
3 57.8951 175.0 0.191 0.832
4 67.9423 189.1 0.228 0.624
diffrac.dat Pearson7A
N center area FWHM shape
0 31.9074 435.0 0.164 1.571
1 39.3616 1876.0 0.146 0.504
2 45.7733 190.3 0.150 0.712
3 56.9166 223.9 0.175 0.717
4 66.7720 698.2 0.164 0.485
и файл center.out в котором строки содержат имя файла и
положения центров пиков по порядку:
20 32.4379 40.0144 46.5596 57.8951 67.9423
diffrac 31.9074 39.3616 45.7733 56.9166 66.7720
Затем открываю на запись файлы outfile.list и center.out:
outfile = open('outfile.list', 'a')
center_outfile = open('center.out', 'a')
Проверяем есть ли хотя бы один передаваемый скрипту параметр:
try:
input_file = sys.argv[1]
except:
print "Usage:", sys.argv[0], "input files"
sys.exit(1)
sys.srgv[0] соответствует самому скрипту.
Аналогично проверяем есть ли в данной директории файл с базовой линией
и активными зонами base_active.list.
Ну и наконец, задаём в интерактивном режиме функцию профиля.
Далее, я решил код оформить в виде двух функций: формирующей входной
файл для cfityk и выкусывающей из выходного файла нужные параметры.
Первая функция выглядит следующим образом:
def create_cfityk_input(active):
base = active.readline()
myfile = ""
while True:
line = active.readline()
if line == "":
break
x = line.split()
x1 = x[0]
x2 = x[1]
myfile += "a = ("+ x1+ " < x < "+ x2+ ")\n"
myfile += "%p = guess" + FUNCTION + "\n"
myfile += "fit\n"
myfile += "info+%p >> ./peaks\n"
myfile += "delete %p\n\n"
active.close()
return myfile, base
Сначала читается первая строчка - базовая линия. Затем, в цикле, все
оставшиеся строчки до конца файла. Элементы в этих строчках разделены
пробельными символами, поэтому легко выкусываем x1 и x2 в активной
области. И формируем тело файла. Параметры пиком пишем в файл
'peaks'. Функция возвращает основное тело файла и строку с базовой
линией. Возвращаемое значение я присвоил кортежу:
(fit_section, base) = create_cfityk_input(active)
Далее необходимо скормить полученный файл программе cfityk. Это я
реализовал во второй функции:
def fit_one_file(input_file):
try:
os.remove('peaks')
except:
pass
myfirstlines = '@0<' + input_file + "\n" + base + "\n"
cfityk_input = open('cfityk_input', 'w')
cfityk_input.write(myfirstlines)
cfityk_input.write(fit_section)
cfityk_input.close()
os.system("cfityk < cfityk_input") # запускаем fityk
os.remove("cfityk_input") # удаляю входной файл для fityk - больше не нужен
Сначала удаляем старый файл 'peaks' (если есть), в который мы пишем
параметры рефлексов. Затем формируем первые две линии во входном
файле, содержащие собственно имя файла и базовую линию. Наконец,
запускаем cfityk и удаляем входной файл.
После прохода по каждому входному файлу формируется файл peaks,
имеющий вид:
%p = Pearson7($_1, $_2, $_3, $_4)
Pearson7(height, center, hwhm, shape=2) = height/(1+((x-center)/hwhm)^2*(2^(1/shape)-1))^shape
height = $_1 = ~1951.83 = 1951.83 [auto]
center = $_2 = ~31.9079 = 31.9079 [auto]
hwhm = $_3 = ~0.0920368 = 0.0920368 [auto]
shape = $_4 = ~2.59149 = 2.59149 [auto]
FWHM: 0.184074
Area: 421.841
%p = Pearson7($_5, $_6, $_7, $_8)
Pearson7(height, center, hwhm, shape=2) = height/(1+((x-center)/hwhm)^2*(2^(1/shape)-1))^shape
height = $_5 = ~153.244 = 153.244 [auto]
center = $_6 = ~39.3615 = 39.3615 [auto]
hwhm = $_7 = ~0.090969 = 0.090969 [auto]
shape = $_8 = ~1.1861 = 1.1861 [auto]
FWHM: 0.181938
Area: 39.7215
…
Чтобы вырезать отсюда значения положения центра, полуширины, площади
под кривой и параметра shape (который присутствует только для функций
Pearson7A, Pearson7, PseudoVoigt, Voigt, представляющих различные
комбинации гаусса и лоренца) можно воспользоваться регулярными
выражениями.
Ценную информацию о регулярных выражениях в Python можно найти на
INTUIT.RU.
Для определения положения центра сгодятся выражения вида:
pattern_center = r"center.+=\s(.*)\[auto\]"
Аналогично будет для определения shape.
В зависимости от выбранной функции будут различные выражения для
вырезания площади под кривой.
Для Voigt, PseudoVoigt, Pearson7, Gaussian, Lorentzian:
pattern_area = r"Area:\s(.*)"
Для Pearson7A:
pattern_area = r"area.+=\s(.*)\[auto\]"
И для ширины на полувысоте:
pattern_fwhm = r"FWHM:\s(.*)"
Собственно, и всё, остальное чисто технические моменты.
Полный скрипт:
#!/usr/bin/env python
#-*-coding: utf-8 -*-
import sys, math, os, re
try:
os.rename('outfile.list', 'outfile.list.bak')
except: pass
try:
os.rename('center.out', 'center.out.bak')
except: pass
outfile = open('outfile.list', 'w')
center_outfile = open('center.out', 'w')
try:
input_file = sys.argv[1]
except:
print "Usage:", sys.argv[0], "input files"
sys.exit(1)
try:
active = open('base_active.list', 'r')
except:
print "Need base_active.list file of baseline formula and active areas"
sys.exit(1)
input_FUNCTION = input("Enter type of function:\n\
1 - Pearson7A\n\
2 - Pearson7\n\
3 - PseudoVoigt\n\
4 - Voigt\n\
5 - Lorentzian\n\
6 - Gaussian\n\
")
if input_FUNCTION == 1:
FUNCTION = 'Pearson7A'
elif input_FUNCTION == 2:
FUNCTION = 'Pearson7'
elif input_FUNCTION == 3:
FUNCTION = 'PseudoVoigt'
elif input_FUNCTION == 4:
FUNCTION = 'Voigt'
elif input_FUNCTION == 5:
FUNCTION = 'Lorentzian'
elif input_FUNCTION == 6:
FUNCTION = 'Gaussian'
else:
print "Inappropriate value"
sys.exit(1)
def create_cfityk_input(active):
base = active.readline()
myfile = ""
while True:
line = active.readline()
if line == "":
break
x = line.split()
x1 = x[0]
x2 = x[1]
myfile += "a = ("+ x1+ " < x < "+ x2+ ")\n"
myfile += "%p = guess" + FUNCTION + "\n"
myfile += "fit\n"
myfile += "info+%p >> ./peaks\n"
myfile += "delete %p\n\n"
active.close()
return myfile, base
def fit_one_file(input_file):
try:
os.remove('peaks')
except:
pass
myfirstlines = '@0<' + input_file + "\n" + base + "\n"
cfityk_input = open('cfityk_input', 'w')
cfityk_input.write(myfirstlines)
cfityk_input.write(fit_section)
cfityk_input.close()
os.system("cfityk < cfityk_input") # запускаем fityk
os.remove("cfityk_input") # удаляю входной файл для fityk - больше не нужен
# регулярные выражения для выделения положения центра,
# площади под кривой и ширины на полувысоте
pattern_center = r"center.+=\s(.*)\[auto\]"
center_re = re.compile(pattern_center)
if FUNCTION == 'Voigt' or FUNCTION == 'PseudoVoigt' \
or FUNCTION == 'Pearson7' or FUNCTION == 'Gaussian' \
or FUNCTION == 'Lorentzian':
pattern_area = r"Area:\s(.*)"
elif FUNCTION == 'Pearson7A':
pattern_area = r"area.+=\s(.*)\[auto\]"
area_re = re.compile(pattern_area)
pattern_fwhm = r"FWHM:\s(.*)"
fwhm_re = re.compile(pattern_fwhm)
if FUNCTION == 'Pearson7A' or FUNCTION == 'Pearson7' \
or FUNCTION == 'PseudoVoigt' or FUNCTION == 'Voigt':
complex_functions = 1
pattern_shape = r"shape.+=\s(.*)\[auto\]"
shape_re = re.compile(pattern_shape)
# открываю вывод cfityk
buff = ""
peaks_params = open('peaks', 'r')
while True:
line = peaks_params.readline()
if line == "":
break
# создаем буфер из выходного файла
buff += line
peaks_params.close()
param_center = center_re.findall(buff)
param_area = area_re.findall(buff)
param_fwhm = fwhm_re.findall(buff)
if complex_functions == 1:
# функции Pearson7A, Pearson7, PseudoVoigt, Voigt
param_shape = shape_re.findall(buff)
outfile.write('%s\t %s\n' % (input_file, FUNCTION) )
outfile.write("N\tcenter\tarea\tFWHM\tshape\n")
center, area, fwhm, shape = [], [], [], []
for i in range(len(param_center)):
center.append(float(param_center[i]))
try:
area.append(float(param_area[i]))
except:
pass
fwhm.append(float(param_fwhm[i]))
if complex_functions ==1:
shape.append(float(param_shape[i]))
outfile.write('%g\t%0.4f\t%0.1f\t%0.3f\t%0.3f\n' \
% (i, center[i], area[i], fwhm[i], shape[i]) )
outfile.write("\n")
pattern_name = r'(\S+)\..*$'
name_re = re.compile(pattern_name)
nameWithoutExtension = name_re.findall(input_file)[0]
center_outfile.write('%s\t' % nameWithoutExtension)
for i in range(len(center)-1):
center_outfile.write('%0.4f\t' % center[i])
center_outfile.write('%0.4f\n' % center[-1])
(fit_section, base) = create_cfityk_input(active)
for i in range(1, len(sys.argv)):
""" фит всех входных файлов, кроме собственно скрипта """
fit_one_file(sys.argv[i])
outfile.close()
center_outfile.close()
Запускаем командой:
python script_fityk.py file1.dat file2.dat file3.dat
На выходе получим файлы outfile.list формата:
имя_файла профильная_функция
номер_пика положение_центра площадь_под_кривой ширина_на_полувысоте форма(если есть)
и center.out:
имя_файла положения_центров_пиков
Меня больше всего интересуют положения центров, поэтому я их и выделил
в отдельный файл, хотя ничто не мешает сделать подобное для площадей
пиков и полуширин. Можно это же сделать всё и средствами командной
оболочки bash, обработав файл outfile.list.
Делаем скрипт исполняемым
Для того, чтобы запускать скрипт из любой директории, нужно сделать
его исполняемым:
chmod u+x script_fityk.py
и положить его куда-нибудь, где хранятся в домашней директории
бинарники. У меня это ~/bin. Эта директория должна быть в $PATH,
для этого в файле ~/.bashrc должны быть строки типа:
export PATH=$PATH:$HOME/bin
то есть к уже существующей переменной окружения $PATH добавляем ещё
одну директорию. После этого обновите оболочку посредством
перезапуска терминала или
source ~/.bashrc
На этом я заканчиваю описание скрипта, существенно упрощающего мне
жизнь. Пожелания и предложения приветствуются!
Далее попробую использовать положения рефлексов для определения
параметров элементарной ячейки структуры.
