Пишем скрипт для Fityk

Предположим, у нас есть большое количество практически идентичных дифрактограмм. Допустим, снимались они на одном и том же дифрактометре (стало быть, имеем одну и ту же базовую линию), на одной и той же длине волны (то есть имеем приблизительно одно и то же положение рефлексов). Также, пусть на картине отсутствует перекрытие рефлексов (к нему надо подходить всё-таки творчески и не доверять целиком компьютеру).

Задача — обработать потоково несколько дифракционных картин сразу, результаты аппроксимации (положение центра, площадь под кривой, ширина на полувысоте) требуется записать в общий файл.

Команды в 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+ %имя_функции для получения полной информации о пике. Итак, полный цикл обработки одной дифракционной картины должен быть таким:

1. Открываем файл с дифрактограммой
2. Удаляем базовую линию
3. Определяем активную область
4. «Угадываем» положение пика, описанного определенной функцией
5. Производим уточнение
6. Записываем параметры пика в файл
7. Удаляем аппроксимирующую функцию
8. Переходим к другой активной области и возвращаемся к пункту 3.
9. Открываем новый файл с дифрактограммой

Для одного файла на языке 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

На этом я заканчиваю описание скрипта, существенно упрощающего мне жизнь. Пожелания и предложения приветствуются! Далее попробую использовать положения рефлексов для определения параметров элементарной ячейки структуры.

Копируете статью - поставьте ссылку на оригинал!