Восстановление структуры по "экспериментальным данным РСА" в 1D модели



1. Задание функции.



Модель для компьютерного эксперимента (параметры взяты из примера variant7.txt)

На отрезке [0,30] Å расположены молекулы. Атомы в молекулах связаны ковалентно и находятся на расстоянии 1-1.5 Å друг от друга. Молекулы расположены на расстоянии 3-5 Å (предполагаются водородная связь или гидрофобное взаимодействие). Всего 6 атомов. Электронные плотности (ЭП) атомов описываются гауссовой кривой, максимум которой в центре атома приблизительно пропорционален числу электронов в атоме. Гауссова функция определяется числами λ (высота пика, соотвествует числу электронов), β (ширина пика), γ (максимум пика) по формуле:
gauss = λ*exp(-(β^2)*(X-γ)^2).

Использованная команда:
python compile-func.py -g 10,2,5+10,2,7+30.0,2,12+30.0,2,16+30.0,2,22+30,2,24
Выдача: func.txt

По рис. 1 видно, что полученная функция соотвествует 2 двухатомнынм молекулам (крайние пики) и 2 отдельным атомам (пики в центре).

Рисунок 1. Смоделированная функция ЭП.




2. Расчет амплитуд и фаз сигналов, моделирующих экспериментальные данные.



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

Использованная команда при разложении функции в ряд Фурье:
python func2fourier.py -i func.txt -o range.txt
Выдача: range.txt, содержит 499 гармоник, для каждой из которых восстановлены амплитуда (F) и фаза (phi).



3. Восстановление функции ЭП по амплитудам и фазам части сигналов.



Необходимые технические определения:
Отличное восстановление – по графику восстановленной функции можно определить положение максимума всех гауссовых слагаемых функции (наших атомов).
Хорошее восстановление – можно угадать положение всех максимумов, зная число слагаемых (атомов), хотя на восстановленной функции максимумы от атомов не отличимы от шума.
Cреднее восстановление – положение каких-то атомов определить по восстановленной функции нельзя, других - можно.
Плохое восстановление – положение атомов определить не представляется возможным; можно только предсказать примерный размер молекулы.
1) Использование полного набора гармоник дает идеальное восстановление.
Полный набор гармоник ряда Фурье - такой, для которого известны все гармоники с номерами 0, 1, 2, ..., n.
Разрешение полного набора гармоник - период гармоники с номером n, т.е. с наибольшим номером.
Период гармоники - расстояние между соседними максимумами синусоиды; его также называют длиной волны этой гармоники, хотя никакой физической волны нет.


Использовання команда:
python fourier2func.py -f func.txt -i range.txt
Выдача: two_func.txt с функцией и ее восстановленным вариантом.


Рисунок 2. Восстановленная по всем гармоникам функция ЭП.



2) Для хорошего восстановления исходной функции ЭП необходимо определить минимальное число используемых гармоник. Так, fourier-filter.py позволяет отобрать некоторые гармоники и учитывать только их при восстановлении (в моем случае менялся параметр n, 1≤n≤40).

Использованные команды:
python fourier-filter.py -r 0-n -i range.txt -o range_n.txt
python fourier2func.py -f func.txt -i range_n.txt -o range_filn.txt0
Выдача - аналогичные выше приведенному текстовые файлы и графики восстановленной функции.

Как можно наблюдать на рис. 3-8, при n=40 восстановление получается практически идеальным, поэтому далее для вычислений будет использоваться 41 гармоника.


Рисунок 3. ЭП, восстановленная по первым 2 гармоникам (n=1).



Рисунок 4. ЭП, восстановленная по первым 6 гармоникам (n=5).



Рисунок 5. ЭП, восстановленная по первым 11 гармоникам (n=10).



Рисунок 6. ЭП, восстановленная по первой 21 гармонике (n=20).



Рисунок 7. ЭП, восстановленная по первой 31 гармонике (n=30).



Рисунок 8. ЭП, восстановленная по первой 41 гармонике (n=40).



3) Поскольку в реальном эксперименте присутствуют ошибки в измерении интенсивности сигналов (амплитуд) и фаз, нужно учесть этот факт в модели, что обеспечивается за счет парметров гауссовского шума: -F для амплитуд (число а), -P для фаз (число b).

Использованные команды:
python func2fourier.py -F a -P b -i func.txt -o output1.txt
python fourier-filter.py -r 0-40 -i output1.txt -o output2.txt
python fourier2func.py -f func.txt -i output2.txt -o output3.txt
Рисунок 9. ЭП, восстановленная по данным с зашумленной амплитудой (F=10).



Рисунок 10. ЭП, восстановленная по данным с зашумленной амплитудой (F=15).



Рисунок 11. ЭП, восстановленная по данным с зашумленной фазой (F=10).



Рисунок 12. ЭП, восстановленная по данным с зашумленной фазой (F=15).



Рисунок 13. ЭП, восстановленная по данным с зашумленной амплитудой и фазой (F=10, P=10).



Рисунок 14. ЭП, восстановленная по данным с зашумленной амплитудой и фазой (F=15, P=15).



Рисунок 15. ЭП, восстановленная по данным с зашумленной амплитудой и фазой (F=20, P=20).



Как можно заметить (рис. 9-15), зашумление в 10% ни по амплитуде, ни по фазе не дает существенного искажения. При повышении уровня шума до 15% положения молекул и атомов все еще однозначно определяются. Совместный шум в 10% по фазе и амплитуде не вносит особого вклада, 15% - более ощутим, а начиная с шума в 20% становится сложно выявить наименьшие из истинных пиков. Другие комбинации параметров на рисунках приведены не были, однако при повторении итераций я обнаружила, что фаза вносит чуть больший вклад, нежели интенсивность, в искажение функции ЭП.

4) Теперь посмотрим, как функция восстанавливается без первых n гармоник.

Использованные команды:
python fourier-filter.py -r n-40 -i range.txt -o rann_40.txt
python fourier2func.py -f func.txt -i rann_40.txt -o n_40.txt
Рисунок 16. ЭП, восстановленная по всем гармоникам, кроме первой (1-40).



Рисунок 17. ЭП, восстановленная по всем гармоникам, кроме первых двух (2-40).



Рисунок 18. ЭП, восстановленная по всем гармоникам, кроме первых десяти (10-40).



Рисунок 19. ЭП, восстановленная по гармоникам 0-25, 30-40.



Рисунок 20. ЭП, восстановленная по гармоникам 0-40 и 50.



Когда удаляется первая гармоника, то из функции вычитается константа (рис. 16). При удалении дальнейших гармоник из функции вычитаются суммы синусоид так, что удаление небольшого числа первых гармоник достаточно сильно влияет на качество восстановленной функции. Хотя даже в этом случае можно легко определить позиции атомов (рис. 17). Функция, восстановленная после удаления десяти первых очень сильно отличается от истинной ЭП - пики и шумы примерно на одном уровне и не отличимы друг от друга (рис. 18). При удалении 5 гармоник из середины набора никаких резких ухудшений восстановления функции не наблюдается (рис. 19), при этом удаление одного и того же числа гармоник из разных частей набора неэквивалентно. При добавлении дополнительной 50-ой гармоники решение не становится лучше, оно так же близко к истинному, как и для набора гармоник 0-40.

С учетом ранее введенных определений касательно качества восстановления ЭП по результатам моделирования составлена данная таблица.