Предположим, у нас есть большое количество практически идентичных дифрактограмм. Допустим, снимались они на одном и том же дифрактометре (стало быть, имеем одну и ту же базовую линию), на одной и той же длине волны (то есть имеем приблизительно одно и то же положение рефлексов). Также, пусть на картине отсутствует перекрытие рефлексов (к нему надо подходить всё-таки творчески и не доверять целиком компьютеру).
Задача — обработать потоково несколько дифракционных картин сразу, результаты аппроксимации (положение центра, площадь под кривой, ширина на полувысоте) требуется записать в общий файл.
Команды в 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
На этом я заканчиваю описание скрипта, существенно упрощающего мне жизнь. Пожелания и предложения приветствуются! Далее попробую использовать положения рефлексов для определения параметров элементарной ячейки структуры.
Читать далее…
