Учебная страница курса биоинформатики,
год поступления 2013
Исследуйте качество восстановления функции ЭП от одной переменной в зависимости от того, какие и сколько гармоник ряда Фурье используются для ее восстановления
Описание задания
Также см. пояснение
Задание включает следующие этапы.
- Создание модельной функции ЭП в одномерной элементарной ячейке
- Расчет параметров сигнала, моделирующих экспериментальные данные: амплитуды и фазы
- Восстановление функции ЭП по модельным ("экспериментальным") данным
- Оценка качества восстановления функции ЭП
Этапы 3 - 4 предстоит повторить несколько раз.
Cкрипты для выполнения задания
- Скрипты проверены для python 2.7. Под python 3.0 не проверял.
- Для запуска дома необходимы библиотеки numpy и matplotlib. Скачайте и установите.
- Скопировать к себе обязательно надо файл fourier.py. Все скрипты используют этот модуль.
Запуск без параметров python <имя скрипта> выдает информацию о том, что делает скрипт.
- Запуск с параметром -h сообщает какие есть параметры и как их задавать.
1. Задание функции
Модель для компьютерного эксперимента. На отрезке [0,30] (ангстремы) расположены две молекулы. Атомы в молекуле связаны ковалентно и находятся на расстоянии 1-1.5 анстрем друг от друга. Молекулы расположены на расстоянии 3-5 ангстрем (водородная связь или гидрофобное взаимодействие между ними). Всего 5-7 атомов (2+3 или 2+3+2; возможны другие варианты). Электронные плотности (ЭП) атомов описываются гауссовой кривой. Максимум ЭП в центре атома приблизительно пропорционален числу электронов в атоме. Берите разные атомы.
Функция задается на интервале [0,30]. Моделируется в 1D график электронной плотности молекулы. 30 следует понимать как 30 ангстрем. Функция имеет вид нескольких гауссовых кривых с центром в разных точках. Иногда - на расстоянии 1-1.5 ангстрем - модель ковалентно связанных атомов; иногда 3-5 ангстрем.
- Формат файла с функцией:
#X Y [строка с # в первой строке игнорируется]
0.0 1.01
0.01 1.07
... ....
29.97 1.00
Функцию можно создать с помощью скрипта compile-func.py. Гауссова функция определяется числами lambda,beta,gamma по формуле: gauss = lambda*exp(-(beta^2)*(X-gamma)^2). Здесь a^2=a*a. Пример задания параметров скрипта:
-g 30,3,3+40,3,4.3+2,3.5,6.5+30,3,7.5 (сумма четырех гауссовых функций: их max в точках 3,4.3,6.5,7.5, высота 30,40,2,30 соответственно; ширина "колокола" порядка 1; третий - низкий - пик моделирует атом водорода);
- ширина колокола должна примерно соответствовать диаметру атома - порядка 1 с небольшим ангстрема (электронная плотность резко спадает на таком расстоянии от центра);
- высота зависит от числа электронов в атоме, поэтому может различаться.
2. Расчет амплитуд и фаз сигналов, моделирующих экспериментальные данные
Амплитуды и фазы рассчитываются по входной функции ЭП. При моделировании экспериментальных данных учитывается, что в эксперименте, во-первых, определяются амплитуды не для всех сигналов; во-вторых, интенсивности сигналов (следовательно, и амплитуды) определяются с ошибкой; в-третьих, фазы определятся для всех измененных сигналов, но тоже с ошибкой.
Коэффициенты Фурье рассчитываются с помощью скрипта func2fourier.py
Входной файл – функция, полученная как выходной файл скрипта compile-func.py .
- Выходной файл имеет формат:
<номер гармоники> <амплитуда> <фаза>
Добавление гауссового шума к амплитудам (параметр -F <число>) и фазам (-P <число>) искажает все вычисленные амплитуды и фазы. Пример: -F 20 (шум 20%) приводит к тому, что к каждой амплитуде прибавляется случайное число, распределенное нормально с параметрами: среднее = 0, среднее квадратичное отклонение (сигма)=0.2*F. Аналогично действует параметр -P <число>
Отбор гармоник можно сделать одним из способов:
- Скопировать файл с коэффициентами Фурье и убрать лишние строки.
- Пометить лишние строки знаком "#" в первой позиции строки (такие строки при чтении игнорируются).
Запустить скрипт fourier-filter.py , который отфильтрует нужные строки (разбирайтесь сами - help'ы есть).
3. Восстановление функции ЭП по амплитудам и фазам части сигналов
Восстановление функции по отобранным гармоникам выполняется скриптом fourier2func.py. Кроме графика, выдается также файл с исходной и восстановленной функцией. Для отключения графика исходной функции служит параметр `-s' (если захотите экспериментировать с коллегой).
Используйте возможности графического окна для выделения нужной части графика и его сохранения.
4. Как сравнить восстановленную функцию с исходной
Отличное восстановление – по графику восстановленной функции можно определить положение максимума всех гауссовых слагаемых функции ("атомов")
Хорошее восстановление – можно угадать положение всех максимумов, зная число слагаемых ("атомов"), хотя на восстановленной функции максимумы от атомов не отличимы от шума
Среднее восстановление – положение каких-то атомов определить по восстановленной функции нельзя, других - можно
Плохое восстановление – положение атомов определить не представляется возможным; можно только предсказать примерный размер "молекулы"
Можете привлечь коллегу: покажите график восстановленной функции, не показывая графика исходной, и попросите определить положение атомов.
Cкрипты лежат также на P: в директории y11/Term_7/Fourier