Стрелки в Pylab

Pylab представляет собой очень мощную библиотеку для графопостроения. При проведении расчётов и анализа данных в Python, при помощи pylab не представляет труда вывести результат на график. Однако просто графики зачастую не являются информативными для неподготовленного читателя. В этом случае, помогают подсказки, сделанные на графиках, в частности, стрелки с указаниями. Кроме того, стрелки являются незаменимыми при построении графиков с несколькими y-осями. В этом случае весьма желательно указывать с помощью стрелок к какой оси относится тот или иной график.

Pylab.arrow

Казалось бы, при всей мощи, от pylab следовало бы ожидать удобной возможности нанесения такого рода деталей на графики. Но, надо признать, что здесь всё не очень очевидно. К примеру, имеется функция pylab.arrow, принимающая начальные координаты (x,y) и (dx,dy). Казалось бы, нет ничего проще.

import pylab as pl
import numpy as np
pl.figure(figsize=(4,4))
pl.title('pylab.arrow, square ratio')
z = np.arange(10)
pl.plot(z,z,'ro')
pl.arrow(2,2,3,4, head_width=0.5, head_length=1)
pl.show()

Однако, если нарисовать график в осях, существенно отличающихся по масштабу, то вид будет ужасный. Например:

x = np.arange(0.1, 1000, 0.1)
y = np.log(x)

pl.figure(figsize=(6,4))
pl.title('pylab.arrow, different axis scale')
pl.plot(x, y, 'b-')
pl.arrow(400,2,200,3, head_width=0.5, head_length=1)
pl.arrow(200,2,300,4, head_width=0.5, head_length=100)
pl.ylim(ymin = 0)
pl.show()

Мало того, что форма стрелки искажённая, так ещё и надо вручную задавать размеры наконечника head_length и head_width под каждый масштаб.

Pylab.annotate

Однако оказалось, что есть более приличный способ отрисовки стрелок, pl.annotate. На самом деле основное предназначение pl.annotate состоит в том, чтобы отмечать стрелками заданные точки, и одновременно указывать подписи. Например, таким образом:

Для этого был использован следующий код:

pl.figure(figsize=(6,4))
pl.plot(x, y, 'b-')
pl.title('pylab.annotate example')
pl.annotate('Label', xy = (500, np.log(500)), xytext=(300,1),\
                arrowprops=dict(facecolor='red', width=0.5,\
                                    headwidth=10, shrink=0.05),\
                ha='center', va='baseline', fontsize='large')

В качестве аргументов передаются

• строка с надписью
• координаты точки, куда указывает стрелка
• координаты надписи
• свойства стрелки в ввиде словаря
• остальные свойства, касающиеся свойств текста

Все свойства текста описаны в мануале по pylab.annotate, посмотреть который можно в интерактивном режиме питона. Например:

$ ipython
> import pylab as pl
> help(pl.annotate)

Что касается свойств стрелки, то их можно посмотреть так:

$ ipython
> import matplotlib
> help(matplotlib.lines.Line2D)

Из основных характеристик стрелок, на которые следует обратить внимание, отмечу следующие:

• facecolor — цвет наконечника
• width — толщина линии в пунктах
• headwidth — ширина наконечника в пунктах
• frac — доля длины, которую занимает наконечник
• shrink — отступ, от заданной точки xy до конца наконечника (чтобы наконечник не втыкался точно в заданную точку, а оставил некоторое пространство). Измеряется в долях длины.

Что касается свойств текста, то я бы в первую очередь отметил следующие:

• ha (horizontalalignment) — способ выравнивания по горизонтали ('center', 'right', 'left')
• va (verticalalignment) — соответственно, по вертикали ('center', 'top', 'bottom', 'baseline')
• rotation — поворот текста (угол в градусах, 'vertical', 'horizontal')
• fontsize — размер шрифта (можно указать в пунктах, а можно и словами 'small', 'medium', 'large', 'x-large', 'xx-large')

Как видно, при использовании pylab.annotate форма стрелки получилась адекватной. А если вместо надписи оставить пустую строку, то просто получится стрелка от xytext до xy (при условии, что shrink = 0)

Лично мне весьма неудобно выписывать все свойства arrowprops всякий раз, когда надо нарисовать стрелку, поэтому я нацарапал простейшую функцию, включающую в себя обязательные аргументы, такие как строка с надписью и координаты начальной и конечной точек, а также необязательные, касающиеся характеристик стрелок и текста:

def my_arrow(label, xy_from, xy_to, color='blue', shrink=0.05, width=0.5, headwidth=10):
    pl.annotate(label, xy = xy_to, xytext=xy_from,\
                    arrowprops=dict(facecolor=color, width=width,\
                                        shrink=shrink, frac=0.1),\
                    ha='center', va='baseline', fontsize='large')

Теперь нет ничего проще, нарисовать стрелку:

pl.figure(figsize=(6,4))
pl.title('pylab.annotate')
pl.plot(x, y, 'b-')
my_arrow('Label', (100,1), (600, 5))
my_arrow('', (300,1), (800, 5))
# красные точки для понимания каким образом позиционируется текст
pl.plot([100,600],[1,5], 'ro')
pl.ylim(ymin = 0)

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

Стоит заметить, что этим вид стрелок не ограничивается. Можно также рисовать стрелки с рюшечками, стоит только в arrowprops упомянуть про arrowstyle.

pl.figure(figsize=(6,4))
pl.title('Fancy arrow')
pl.plot(x, y, 'b-')
pl.annotate('Fancy arrow', xy = (500, np.log(500)), xytext=(300,1),\
                arrowprops=dict(facecolor='red',arrowstyle='fancy',\
                           connectionstyle='angle3, angleA=0, angleB=90'),\
                fontsize='x-large')
pl.ylim(ymin = 0)

Для более детального ознакомления следует посмотреть в примеры mpl_examples/pylab_examples/annotation_demo2.py

Читать далее…

Однотипные срезы массивов в python

Известно, что в python есть такая замечательная вещь, как срезы массивов (списков, кортежей и т.д.). По сути, в качестве начала и конца среза передаются положения элементов в массиве. Например,

x = range(10)
print x
> [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
print x[1:3]
> [1, 2]
print x[0:5]
> [0, 1, 2, 3, 4]

Видно, что по сути отсчёт ведётся не по элементам, а по промежуткам между элементами, при этом на выходе даются все элементы, попадающие в заданный интервал промежутков. Особый интерес представляют отрицательные индексы, являющиеся по сути обратным отсчётом промежутков с конца массива.

print x[1:-1]
> [1, 2, 3, 4, 5, 6, 7, 8]
print x[0:-2]
> [0, 1, 2, 3, 4, 5, 6, 7]

Однако, что делать, если, скажем, имеется набор экспериментальных данных, которые мы хотим обрезать, основываясь на значениях, нежели на положении элементов в массивах? Допустим, есть список x,

from math import *
import numpy as np
x = np.arange(0, 2*pi, 0.1)

из которого надо вырезать 1.2<x<2.7.

Традиционный способ

В принципе, это можно осуществить перебором по всем элементам списка:

begin,end = 0,0
x0, x1 = 1.2, 2.7
for i in range(len(x)):
    if begin == 0:
        if x[i] > x0:
            begin = i
        else: pass
    else: pass
    if int(end) == 0:
        if x[i] > x1:
            end = i
        else: pass
    else: pass
print begin, end
print x[begin:end]

В результате получим:

> 12 28
[ 1.2  1.3  1.4  1.5  1.6  1.7  1.8  1.9  2.   2.1  2.2  2.3  2.4  2.5  2.6 2.7]

В принципе. то, чего и ожидали. Теперь можно этот цикл for оформить внутри функции, которая будет принимать в качестве аргументов список и границы.

def conventional_cut(x, (x0,x1)):
    begin,end = 0,0
    for i in range(len(x)):
        if begin == 0:
            if x[i] > x0:
                begin = i
            else: pass
        else: pass
        if end == 0:
            if x[i] > x1:
                end = i
            else: pass
        else: pass
    return x[begin:end]
print conventional_cut(x, (0.9, 1.4))

> [ 1.   1.1  1.2  1.3]

Несколько массивов

С практической точки зрения удобно резать сразу два списка: агрумент и функцию. Допустим, отснимали большой спектр y от x, из которого интерес представляет только часть. Таким образом, необходимо откусить часть от x, а также соответствующую часть y.

Сделать это легко, слегка модернизировав функцию:

def cut_xy((x,y), (x0,x1)):
    begin,end = 0,0
    for i in range(len(x)):
        if begin == 0:
            if x[i] > x0:
                begin = i
            else: pass
        else: pass
        if end == 0:
            if x[i] > x1:
                end = i
            else: pass
        else: pass
    return x[begin:end], y[begin:end]

То есть мы будем использовать те же индексы для списка y, что и для x. Приведу пример, для большей наглядности, графический:

from math import *
import numpy as np
import pylab as pl

x = np.arange(0, 4*pi, 0.01)
y = np.sin(x)/x
y = y + 0.1*(np.random.rand(len(y))-0.5)
# представим, что x,y - экспериментальные данные

xnew, ynew = cut_xy((x,y), (0.75*pi, 1.5*pi))
pl.plot(x, y, 'r-')
pl.plot(xnew, ynew, 'b-')

pl.xlabel("x")
pl.ylabel("Signal")
pl.show()

В результате получим то, что и следовало ожидать.

Выборочное соответствующее вырезание двух массивов

Разумеется, никто не запрещает резать не два, а три и более массивов одновременно, немного видоизменив функцию cut_xy.

Использование numpy.nonzero

Традиционный метод работает, однако может показаться корявым. По-видимому, более элегантное решение может быть получено с использованием функции nonzero из библиотеки numpy.

Пример работы numpy.nonzero:

from math import *
import numpy as np
x = np.arange(0, 2*pi, 0.1)
print x[np.nonzero(x<2.7)]
print x[np.nonzero(x>1.4)]

> [ 0.   0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9  1.   1.1  1.2  1.3  1.4
  1.5  1.6  1.7  1.8  1.9  2.   2.1  2.2  2.3  2.4  2.5  2.6]
> [ 1.4  1.5  1.6  1.7  1.8  1.9  2.   2.1  2.2  2.3  2.4  2.5  2.6  2.7  2.8
  2.9  3.   3.1  3.2  3.3  3.4  3.5  3.6  3.7  3.8  3.9  4.   4.1  4.2  4.3
  4.4  4.5  4.6  4.7  4.8  4.9  5.   5.1  5.2  5.3  5.4  5.5  5.6  5.7  5.8
  5.9  6.   6.1  6.2]

В скобках np.nonzero() указывается булево выражение, при этом nonzero возвращает индексы ненулевых элементов. Выражение x<2.7 является булевым и имеет следующий вид:

[ True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True False False False False False False False False False
 False False False False False False False False False False False False
 False False False False False False False False False False False False
 False False False]

То есть, если x[i]>2.7, то результат True, в противном случае False. Стало быть, передавая x<2.7 в качестве аргумента функции nonzero, и скармливая результат исходному массиву x в качестве индексов среза, получаем обрезанный x, где каждый элемент x не превышает 2.7.

Трудность заключается в том, что не удаётся в качестве аргумента nonzero затолкать сразу несколько условий. Таким образом, если надо откусить массив с двух краёв, то необходимо сначала откусить массив слева, а потом результат — справа:

X = x[np.nonzero(x>1.4)]
X = X[np.nonzero(X<2.7)]
print X
> [ 1.4  1.5  1.6  1.7  1.8  1.9  2.   2.1  2.2  2.3  2.4  2.5  2.6]

Замечу, что при использовании nonzero левая граница промежутка включается в искомый массив.

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

X = x[np.nonzero(x>1.4)]
X = X[np.nonzero(X<2.7)]
Y = y[np.nonzero(x>1.4)]
Y = Y[np.nonzero(X<2.7)]

print "X =", X
print "Y =", Y

На выходе:

X = [ 1.4  1.5  1.6  1.7  1.8  1.9  2.   2.1  2.2  2.3  2.4  2.5  2.6]
Y = [ 0.66167563  0.62029339  0.59024345  0.57468639  0.5094329   0.5281753
  0.4542492   0.39952115  0.37289637  0.36049412  0.266617    0.19706189
  0.247443  ]

Можно аналогичным образом затолкать это в функцию:

def cut((x,y),(x0,x1)):
    after = np.nonzero(x>x0)
    X = x[after]
    before = np.nonzero(X<x1)
    X = X[before]
    Y = y[after]
    Y = Y[before]
    return X,Y

которая принимает в качестве аргументов исходные (x,y) и интервалы x.

Графический пример:

from math import *
import numpy as np
import pylab as pl

x = np.arange(0, 2*pi, 0.1)
y = np.sin(x)/x
y = y + 0.1*(np.random.rand(len(y))-0.5)
# представим, что x,y - экспериментальные данные

pl.plot(x, y, 'r-')
X,Y = cut((x,y),(4,5))
pl.plot(X,Y,'g-')

pl.xlabel("x [channels]")
pl.ylabel("Signal [counts]")
pl.show()

Результат аналогичной процедуры при использовании numpy.nonzero

Читать далее…

Уточнение данных полиномом n-степени в Python (scipy, polyfitw)

Задача - произвести подгонку данных полиномом n-степени. Известно, что каждое, экспериментально измеренное значение, представляет собой два числа: собственно значение и погрешность измерения. Таким образом, желательно произвести подгонку с учётом весовых коэффициентов каждой точки. Ведь может оказаться так, что в эксперименте есть точки, измеренные с меньшей точностью и степенью достоверности, но которые не желательно просто так выбрасывать, но и нельзя учитывать наравне с другими, более достоверными.

В качестве инструмента выбираем Python, благодаря расширяемости этого языка.

Модули в Python

Scipy

В Python для численных операций имеются модули SciPy и NumPy. Однако, несмотря на всю мощь этих модулей, я не нашёл в них возможности производить уточнение именно с учетом весовых коэффициентов. В scipy содержится функция polyfit, которая принимает только значения аргумента, функции и степени полинома, но не учитывает погрешности.

Модуль CARSMath

На сайте The Consortium for Advanced Radiation Sources есть разработанный Марком Риверсом модуль CARSMath. Он содержит функцию polyfitw, которая помимо аргумента, функции и степени полинома принимает также и весовые коэффициенты (именно весовые коэффициенты, а не погрешность, то есть 1/погрешность).

Установка модуля CARSMath

Скачайте модуль отсюда. Бросьте тарбол, куда-нибудь, например, в ~/python/python_epics.

Распакуйте:

tar xvf python_epics.tar

Нас интересует файл CARSMath.py.

Теперь, чтобы использовать этот модуль, надо его бросать в рабочую директорию с вашим скриптом, что очень неудобно. Поэтому надо воспользоваться переменной окружения $PYTHONPATH. Для этого отредактируйте файл ~/.bashrc. Необходимо добавить следующую строчку:

export PYTHONPATH="$HOME/distrib/python/lib"

Естественно, директория должна быть указана та, где лежит скачанный модуль CARSMath (я его перекинул в ~/distrib/python/lib, там же у меня лежат и другие модули, не распространяемые с дистрибутивом).

После этого инициализируем оболочку заново:

source ~/.bashrc

Если вы запускаете иксы не через startx, а через менеджер входа в систему, то вам надо узнать каким образом в этом случае инициализировать переменные окружения. Ключевое слово .XClients, если я ничего не путаю.

Теперь можно, находясь в любой директории, импортировать модуль CARSMath стандартным образом:

import CARSMath

Сравнение двух способов подгонки

Сравним два способа подгонки экспериментальных данных полиномом n-й степени:

• с помощью polyfit из scipy
• с помощью polyfitw из CARSMath

Для этого создадим массив точек, сделаем несколько «выбросов» с большими погрешностями и аппроксимируем полиномом (в одном случае без учёта погрешностей, а во втором - с учётом).

Ниже проведён скрипт, реализующий эти действия.

#!/usr/bin/env python
#-*-coding: utf-8 -*-

import scipy
from numpy import arange
import pylab as pl
from CARSMath import polyfitw

x = arange(0, 2, 0.01)
y = 0.5 * x**3 - 2 * x**2 + x
error = 0.005 * scipy.rand(len(x))

# генерим "выпавшие" точки с большими погрешностями
for i in range(30, 60):
    y[i] = y[i] - 2* abs(scipy.rand(1))
    error[i] = 1 * abs(scipy.rand(1))

# находим коэффициенты полинома без учёта погрешностей
fit_coeffs = scipy.polyfit(x, y, 3)
# структура коэффициентов: (x**(max), x**(max-1), x**(max-2), ...)

def fit_func(x, fit_coeffs):
    """
    фит без учёта весов

    """
    fit_func = 0
    # максимальная степень полинома.
    # В данном случае 3 = 4 - 1
    max_degree = len(fit_coeffs)-1

    for i in range(max_degree):
        fit_func += fit_coeffs[i]* x**(max_degree-i)
    return fit_func


# коэффициенты полинома с учетом весов 1/error
# i-й коэффициент соответствует i-й степени полинома
# здесь полная степень полинома = 3
weight_coeffs = polyfitw(x, y, 1/error, 3)

def weight_func(x, weight_coeffs):
    """
    фит с учётом весов

    """
    weight_func = 0

    # обратите внимание на другую структуру вектора weight_coeffs
    # (x**0, x**1, x**2, x**3)
    for i in range(len(weight_coeffs)):
        weight_func += weight_coeffs[i] * x**i
    return weight_func


pl.errorbar(x, y, error, fmt= 'bo')
# график уточнения без "весов"
pl.plot(x, fit_func(x, fit_coeffs), 'r-', \
        linewidth=2, label='no weights')
# график уточнения с "весами"
pl.plot(x, weight_func(x, weight_coeffs), 'g-', \
        linewidth=2, label='with weights')

pl.title('Comparison between two fitting methods')
pl.legend()
pl.show()

Здесь я создал функции fit_func, weight_func, возвращающие значения подгоночной функции, хотя в таком простом примере этого можно было и не делать.

В результате получим следующую картину:

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

Читать далее…