LXF82:PAW

Текущая версия на 16:47, 28 декабря 2008

• Метамодернизм в позднем творчестве В.Г. Сорокина
• ЛитРПГ - последняя отрыжка постмодерна
• "Ричард III и семиотика"
• 3D-визуализация обложки Ridero создаем обложку книги при работе над самиздатом.
• Архитектура метамодерна - говоря о современном искусстве, невозможно не поговорить об архитектуре. В данной статье будет отмечено несколько интересных принципов, характерных для построек "новой волны", столь притягательных и скандальных.
• Литература
• Метамодерн
• Рокер-Прометей против изначального зла в «Песне про советскую милицию» Вени Дркина, Автор: Нина Ищенко, к.ф.н, член Союза Писателей ЛНР - перепубликация из журнала "Топос".
• Как избавиться от комаров? Лучшие типы ловушек.
• Что делать если роблокс вылетает на windows
• Что делать, если ребенок смотрит порно?
• Почему собака прыгает на людей при встрече?
• Какое масло лить в Задний дифференциал (мост) Visco diff 38434AA050
• О чем может рассказать хвост вашей кошки?
• Верветки
• Отчетность бюджетных учреждений при закупках по Закону № 223-ФЗ
• Срок исковой давности как правильно рассчитать
• Дмитрий Патрушев минсельхоз будет ли преемником Путина
• Кто такой Владислав Поздняков? Что такое "Мужское Государство" и почему его признали экстремистским в России?
• Как правильно выбрать машинное масло в Димитровграде?
• Как стать богатым и знаменитым в России?
• Почему фильм "Пипец" (Kick-Ass) стал популярен по всему миру?
• Как стать мудрецом?
• Как правильно установить FreeBSD
• Как стать таким как Путин?
• Где лучше жить - в Димитровграде или в Ульяновске?
• Почему город Димитровград так называется?
• Что такое метамодерн?
• ВАЖНО! Временное ограничение движения автотранспортных средств в Димитровграде
• Тарифы на электроэнергию для майнеров предложено повысить

[править] PAW: приемы работы

ЧАСТЬ 2 Впечатлены возможностями PAW? Приготовьтесь к большему — сейчас Евгений Балдин продемонстрирует вам эту систему в деле!

Система анализа данных PAW или Physics Analysis Workstation для своей работы не требует доскональных знаний всех подсистем и команд. Но чтобы действовать эффективно, следует изучить логику построения команд и стандартные приемы. Это позволит в дальнейшем легко получать навыки для выполнения более сложных задач.

Если даже вы не планируете использовать PAW, в любом случае будет полезно присмотреться к этому инструменту, так как он сравнительно прост и основные концепции, необходимые для анализа данных, достаточно прозрачны для понимания. Создатели PAW действовали по принципу минимализма. Делалось только необходимое — никаких «рюшечек», зато просто. Жесткая структура команд дополнена возможностью писать свои функции и скрипты.

[править] Простейший анализ

Учиться лучше всего на примере реального анализа. Попробуем сделать нечто подобное.

Будем считать, что программа PAW уже запущена и мы находимся в рабочей директории. Вызов внешних команд обеспечивается с помощью инструкции SHELL (можно сократить до sh).

PAW > sh ls
ascii.png lkravg.dat PAW.metafile ee-ang.rz
th1.eps last.kumac last.kumacold sin.dat

Необходимо провести предварительный анализ данных, представленных в текстовом файле lkravg.dat:

1099279655 4119 0.8318 0.0014 1.13 5.99195
1099397693 4126 0.8404 0.0032 1.07 6.001685
…

Колонки соответствуют time_t — времени в секундах, номеру измерения, исследуемому значению, ошибке значения для текущего измерения и двум сторонним параметрам, от которых интересующее нас значение может зависеть. Задача: имея время и два сторонних параметра, попробовать предсказать исследуемое значение.

Анализировать можно и без модели явления, но чтобы правильно подготовить данные для исследования, необходимо ее иметь. Легенда для этих данных следующая: исследуемое значение представляет из себя степень «ухудшения» качества жидкого криптона (LKr — Liquid Kripton) в LKr-калориметре для регистрации энергии элементарных частиц. Качество вычисляется в относительных единицах по амплитуде сигнала от космических мюонов (эти частицы хорошо выделяются и есть всегда). Два параметра, от которых может зависеть качество: избыточное давление LKr и магнитное поле, в котором находится калориметр. Избыточное давление примерно линейно связано с температурой LKr, что влияет на амплитуду сигнала. В магнитном поле прямолинейная траектория мюона искажается, что тоже может влиять на амплитуду.

В первую очередь, необходимо прочитать данные из текстового файла. Для этого можно воспользоваться командой VECTOR/READ (help ve/re):

#чтение текстового файла в вектора
PAW > ve/re time,run,avg,avg_er,P,H lkravg.dat
# меняем тип маркера
PAW > set mtyp 2
# изобразим зависимость исследуемого значения от времени
PAW > ve/pl avg%time

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

Зависимость качества LKr от времени.

Сначала исключим временную зависимость. Воспользуемся для этого стандартной процедурой подгонки для векторов VECTOR/FIT X Y EY FUNC (help ve/fit), где X соответствует оси времени, Y — исследуемому значению, EY — ошибки определения значения, а FUNC — подгоночной функции. Вместо FUNC можно передать любую функцию, написанную на FORTRAN, но в нашем случае достаточно полинома первой степени. Для полинома в PAW есть сокращенное обозначение PN, где N - степень полинома.

#подгоняем временную зависимость прямой
PAW > ve/fit time avg avg_er P1
*********************************************************
*                                                       *
* Function minimization by SUBROUTINE HFITV             *
* Variable-metric method                                *
* ID = 0 CHOPT =                                        *
*                                                       *
*********************************************************
Convergence when estimated distance to minimum (EDM) .LT. 0.10E+01
EXT PARAMETER APPROXIMATE STEP FIRST
NO. NAME VALUE ERROR SIZE DERIVATIVE
1 P1 3.3620 0.50070E-03 0.85051E-02 -223.19
2 P2 -0.22939E-08 0.44872E-12 0.58032E-11 -0.27325E+12
CHISQUARE = 0.1680E+02 NPFIT = 1644

При выполнении процедуры подгонки PAW выдает стандартную «портянку». Основной интерес для пользователя представляют результаты подгонки. В случае полинома первой степени подгоняются два параметра P1 — свободная константа и P2 — коэффициент пропорциональности

Гистограмма ошибок измеряемого параметра.

Для исключения временной зависимости воспользуемся пакетом SIGMA (help sigma). С помощью команды sigma векторами можно манипулировать, как обычными переменными:

#создаем новый вектор уже без временной зависимости
PAW > sigma avg1=avg-3.3620+0.22939E-08*time

Прежде чем действовать дальше, оценим, какую точность в принципе можно ожидать. Вектору исследуемой величины avg соответствует вектор ошибок avg_er. Посмотрим, чему равны эти ошибки. Число измерений превосходит полторы тысячи, поэтому просмотреть все значения глазами не очень реально. Лучше создать гистограмму:

#создаем из значений вектора гистограмму N15
PAW > ve/pl avg_er 15
#включаем отображение статистики (число событий/среднее/разброс)
PAW > opt sta
#меняем цвет гистограммы на красный
PAW > set hcol 2
#вывод гистограммы N15 в диапазоне от 0. до 0.01
PAW > hi/plot 15(0.:0.01)

Из гистограммы видно, что большинство измерений имеют ошибку меньше 0.15·10-2-0.2·10-2. То есть, даже при самом удачном раскладе точность предсказания не будет выше этих 0.15 %-0.2 %.

Гистограмма ошибок измеряемого параметра

После того, как мы убрали основную (временную) зависимость, можно проверить, зависит ли изучаемая переменная от давления (P) и магнитного поля (H). Для этого воспользуемся способностью PAW создавать и отображать двумерные гистограмм:

#создаем двумерную гистограмму: давление от avg1
PAW > ve/pl P%avg1 10
#изобразим двумерную гистограмму 10 в виде поверхности
PAW > SURF 10 65 -30 1

Обратите внимание на разделитель % между векторами в команде VECTOR/PLOT. Значениям первого вектора противопоставляется ось ординат (ось Y), а значениям второго ось абсцисс (ось X).

Команда HISTOGRAMgfv/2D_PLOT/SURFACE [ID THETA PHI CHOPT] (help surf) позволяет отобразить двумерную гистограмму в виде поверхности. Здесь ID — номер гистограммы, THETA и PHI — углы поворота гистограммы и в сферической системе координат, CHOPT - опции изображения.

Сравнение двух гистограмм с помощью наложения.

Из картинки справа видно, что какая-то зависимость есть — избавимся от нее, как это было сделано с временной зависимостью. В результате подгонки сначала для давления, а потом для поля можно получить еще две формулы:

#поправка на давление
PAW > sigma avg2=avg1-0.10901+0.99336E-01*P
#поправка на магнитное поле
PAW > sigma avg3=avg2+0.10844E-01-0.18258E-02*H
#сравниваем две гистограммы до и после поправок
#на давление и магнитное поле
#переключаем цвет на красный
PAW > set hcol 2
#рисуем гистограмму по значениям вектора
PAW > ve/plot avg3
#переключаем цвет на синий
PAW > set hcol 4
#рисуем вторую гистограмму поверх ‘s’ - superimpose
PAW > ve/plot avg1 ! ‘s’

Обратите внимание на команду SET (help GRAPHICS/SET). Эта инструкция позволяет менять параметры графического представления данных. Если запустить ее без опций, то будет выдан список переменных, которые устанавливаются с помощью этой команды.

На рисунке справа представлены две гистограммы. Красная — конечный результат, а синяя — то, что осталось после исключения временной зависимости. Очевидно, что после учета давления и поля разброс уменьшился. Несимметричность гистограмм указывает на то, что временная зависимость на самом деле нелинейная. В реальности при подгонке использовалась экспоненциальная зависимость плюс некоторая константа. Выбор подгоночной функции был продиктован физической моделью явления, но по большому счету на этом временном интервале она не сильно отличается от обычной прямой.

В результате всех действий была получена зависимость:

AVG=3.46+2.29•10-9-9*time-9.9*10-2-2P+1.8*10-3-3H

Сравнение двух гистограмм. Подгонка функцией Гаусса.

Эта функция позволяет предсказывать исследуемое значение в зависимости от времени, давления и магнитного поля, но какова точность этого предсказания? Для ее оценки достаточно получить распределение разницы значений между экспериментальными точками и данной зависимостью, а затем взять его ширину. Можно просто оценить ширину распределения «на глазок», посмотрев на обсуждаемый рисунок, а можно получить ее из известной функции, более-менее описывающей это распределение, например, из функции Гаусса:

#создаем из значений вектора avg1 гистограмму N11
PAW > ve/plot avg1 11
#подгоняем гистограмму N11 в диапазоне (-0.007,0.007)
#функцией Гаусса (G - Gauss)
PAW > hi/fit 11(-0.007:0.007) G
…
EXT PARAMETER …
NO. NAME VALUE ERROR …
1 Constant 51.608 1.9494 …
2 Mean 0.39292E-03 0.20076E-03 …
3 Sigma 0.51728E-02 0.27301E-03 …
…
#создаем из значений вектора avg3 гистограмму N13
PAW > ve/plot avg3 13
PAW > hi/fit 13(-0.004:0.004) G
…
EXT PARAMETER …
NO. NAME VALUE ERROR …
1 Constant 75.655 3.0477 …
2 Mean 0.20614E-03 0.11857E-03 …
3 Sigma 0.29684E-02 0.16655E-03 …
…
#делим графическое окно на две зоны по оси ординат (help zone)
PAW > zone 1 2
#меняем цвет гистограммы на синий
PAW > set hcol 4
#рисуем гистограмму в диапазоне (-0.009,0.009)
PAW > hi/pl 11(-0.009:0.009)
#меняем цвет гистограммы на красный
PAW > set hcol 2
PAW > hi/pl 13(-0.009:0.009)

Из всех значений в данном случае интересно только значение ширины распределения, или SIGMA^[1].

Если учитывать только временную зависимость, то точность предсказания будет примерно 0.52 %, если же учесть давление и магнитное поле, то точность улучшится до 0.30 %, что существенно хуже идеальных 0.15-0.2 %. Очевидно, что остались еще какие-то неучтенные систематики.

Оценим, какую систематику удалось выбрать, учтя зависимость от давления и магнитного поля. Для этого воспользуемся системой COMIS (help comis):

PAW > comis
PAW
CS> type sqrt(0.52**2-0.30**2)
MND> end
*T SQRT(0.52**2-0.30**2) = 0.4247352
PAW
CS> end

В общем, неплохо. Кстати, более тщательный анализ никаких кардинальных улучшений не дал — удалось только уменьшить «хвосты» и сделать итоговое распределение более симметричным за счет выбора более сложной подгоночной функции.

Быстрый анализ позволяет оценить, к какой точности имеет смысл стремиться. Это очень важно, так как сложность получения большей точности увеличивается от требуемой точности существенно нелинейным образом. Фактически на только что изложенный анализ по разным причинам ушел примерно месяц реального времени. Правда, основные проблемы были вовсе не технические. В частности очень много времени ушло на осознание, что исследуемая величина зависит от времени — это не казалось очевидным.

[править] Ntuple

Разность между экспериментом и предсказания в зависимости от наложенного условия на величину ошибки.

То, что только что было сделано с помощью векторов можно проделать с помощью ntuple. Для этого надо сначала создать ntuple (NTUPLE/CREATE), а затем считать в него текстовый файл (NTUPLE/READ).

# создаем ntuple c ID=1
PAW > nt/cre 1 ‘LKr quality’ 6 ! ! time run avg er P H
#читаем в ntuple с ID 1 текстовый файл
PAW > nt/read 1 lkravg.dat
==> 1644 events have been read
PAW > nt/print 1
*****************************************************************
* NTUPLE ID= 1 ENTRIES= 1644 LKr quality
*****************************************************************
* Var numb * Name * Lower * Upper *
*****************************************************************
* 1 * TIME * 0.109828E+10 * 0.113892E+10 *
* 2 * RUN * 0.406400E+04 * 0.709400E+04 *
* 3 * AVG * 0.742800E+00 * 0.846400E+00 *
* 4 * ER * 0.700000E-03 * 0.201000E-01 *
* 5 * P * 0.104400E+01 * 0.116000E+01 *
* 6 * H * 0.000000E+00 * 0.704971E+01 *
*****************************************************************
#включаем отображение статистики
PAW > opt stat
#отрисовываем разницу между экспериментом и предсказанием
#с различным ограничением на величину ошибки er
PAW > set hcol 1
PAW > nt/pl 1.avg-(3.456-2.29E-9*time-9.9E-2*P+1.8E-3*H)
PAW > set hcol 2
PAW > nt/pl 1.avg-(3.456-2.29E-9*time-9.9E-2*P+1.8E-3*H) er<0.0015 ! ! ! ‘s’
PAW > set hcol 3
PAW > nt/pl 1.avg-(3.456-2.29E-9*time-9.9E-2*P+1.8E-3*H) er>0.0015 ! ! ! ‘s’

Преимущество при использовании ntuple заключается в том, что интерактивно можно накладывать условия для фильтрации данных. Как правило, ntuple создаются с помощью внешних программ, а PAW используется уже для интерактивного анализа. В стандартной документации к PAW очень подробно описано, как это делается.

[править] Гистограммы

Подгонка гистограммы с помощью функции Гаусса и пользовательской функции.

Гистограммы — это базовые объекты PAW. Непосредственно перед отображением данные почти всегда преобразуются в одно- или двумерную гистограмму. На гистограмму можно просто смотреть, а можно подгонять какой-либо теоретической зависимостью (HISTOGRAM/FIT).

#поиск rz-файлов в базовой директории (файл создан внешней программой)
PAW > sh ls *.rz
ee-ang.rz
#чтение файла
#известно что в файле есть ntuple c ID=1
PAW > hi/fil 1 ee-ang.rz
#создаем гистограмм N20 из переменной E1 с некоторыми
условиями
PAW > nt/plot 1.E1 f1=11&&f2=-11&&E1<2 ! ! ! ! 20
#подгоняем гистограмму распределением Гаусса
PAW > hi/fit 20 G

Из рисунка видно, что наблюдаемое распределение функцией Гаусса не подгоняется. Надо что-то делать. Очевидно, что теоретическая функция должна быть как минимум несимметрична. Для этой цели подойдет так называемый логарифмический гаусс. Для подгонки надо создать файл loggaus.for примерно следующего содержания:

C Файл loggaus.for
real function loggaus(x)
C С помощью этого common-блока PAW получает доступ
C к параметрам функции
common/PAWPAR/PAR(4)
sqrtln4=1.177410022515475
A=PAR(1)
pike=PAR(2)
sigma=PAR(3)*pike/100.
assim=PAR(4)
loggaus=0.
if (abs(assim).le.1.E-6) then
assim=sign(1.E-6,assim)
endif
if (sigma.le.0.) goto 10
xx=1.+sinh(assim*sqrtln4)/sqrtln4*(x-pike)/sigma
if (xx<1.E-07) goto 10
loggaus=A*exp(-((log(xx)/assim)**2+assim**2)/2.)
10 continue
end

Функция зависит от четырех параметров: A — амплитуда, pike - местоположение пика, sigma — ширины распределения в процентах, assim — асимметрии.

#создаем вектор параметров с начальными значениями
PAW > ve/cre par(4) r 25. 1.4 5. 0.
# создаем вектор с минимально допустимыми значениями (на глазок)
PAW > ve/cre pmin(4) r 10. 1.3 1. -1.
# создаем вектор с максимально допустимыми значениями
PAW > ve/cre pmax(4) r 40. 1.5 10. 1.
# создаем вектор, для ошибок подгонки
PAW > ve/cre err(4) r
# просим PAW подогнать распределение теоретической функцией
# опции подгонки:
# B - учитывать минимально/максимально допустимые значения
# M - перейти интерактивную сессию Minuit
# M обычно не используется, так как действия PAW при подгонке
# по умолчанию вполне разумны
PAW > hi/fit 20 loggaus.for “BM” 4 par ! pmin pmax err
…
# задаем имена параметров
Minuit > name 1 A
Minuit > name 2 pike
Minuit > name 3 sigma
Minuit > name 4 assim
# задаем метод минимизации (migrad, обычно, самый подходящий)
Minuit > migrad
# просим попробовать улучшить подгонку (дольше, но чуть точнее)
Minuit > improve
…
Minuit > exit
…
#смотрим результаты подгонки
PAW > ve/print par
PAR(1) = 35.4279
PAR(2) = 1.4809
PAR(3) = 4.30618
PAR(4) = -0.999999
#смотрим ошибки
PAW > ve/print err
ERR(1) = 3.61113
ERR(2) = 0.00484378
ERR(3) = 0.339507
ERR(4) = 0.174973
#нарисовать гистограмму 20 еще раз
# e - рисовать статистические ошибки в бинах
PAW > hi/plot 20 e

При подгонке этого распределения основной интерес представляло его ширина: =(4.3±0.3)%. Важно не только значение подгонки, но и оценка ошибки. Например, результат для sigma отличается от того, что должно быть в идеале больше чем на десять ошибок — можно сделать вывод, что есть какая-то, даже не проблема, а «плюха».

[править] Функции

Пример одномерной функции.

Функции также являются базовым объектом для PAW. Для отрисовки одномерных функций используется команда FUNCTION/PLOT.

#включить логарифмический масштаб для оси Y
PAW > opt logy
#нарисовать одномерную функцию
PAW > fun/plot (sin(x)/x)**2+0.1 -10 10
#вернуться к линейному масштабу для оси Y
PAW > opt liny

Обратите внимание на инструкцию opt (help GRAPHICS/OPTION). Эта инструкция по своим функциям схожа с командой set, но в отличии от нее отвечает за организацию представления данных, а не за графическое оформление.

Множество Мандельброта. Опции CONT3 и SURF4, соответственно.

Работу с двумерными функциями продемонстрируем на классическом фрактальном изображении имени Мандельброта. Создадим код на FORTRAN:

C Из официальной документации к PAW
C Файл mandel.for
real function MANDEL(XP)
dimension XP(2)
data NMAX/30/
x=XP(1)
y=XP(2)
xx=0.
yy=0.
do n=1,NMAX
tt=xx*xx-yy*yy+x
yy=2.*xx*yy+y
xx=tt
if (4..lt.xx*xx+yy*yy) goto 1
enddo
1 MANDEL=FLOAT(n)/FLOAT(NMAX)
end

В случае двумерных функций проблема отображения стоит гораздо острее, чем у одномерных. Двумерные функции для отображения преобразуются в гистограммы (help fun2)

# По результатам вычисления mandel.for создаем гистограмму 10
PAW > fun2 10 mandel.for 100 −2.4 .8 100 −1.2 1.2 ‘ ‘
# Выводим гистограмму 10 как контур
PAW > hi/pl 10 cont3
# Выводим гистограмму 10 как поверхность
PAW > hi/pl 10 surf4

Если разрешение не удовлетворяет, то можно создать гистограмму не 100×100, как в примере, а 1000×1000.

[править] Заключение

Объять необъятное совершенно не реально, особенно при лимите на объем текста. Официальная документация содержит около 500 страниц, причем, один алфавитный указатель занимает 17 страниц. PAW - матерый программный продукт, которому уже двадцать лет. Этому пакету есть достойный приемник, правда, не лишенный недостатков, но об этом в следующий раз.

1. ↑ В случае гауссоподобного распределения в диапазоне (<x>- ,<x>+) лежит примерно 68 % от всех событий. <x> — среднее значение или MEAN, соответствует SIGMA