2014年5月21日水曜日

bool

ブール型の演算について

どの言語もそうですが、整数型や浮動小数点型の演算はだいたい同じですが、ブール型の演算は少しずつやり方が異なっているような気がします。

pythonでの扱い方について記録しておきます。

a = True
b = False
a and b

とすると
False

となります。
a+b
とすると
 1
となり、
a+a
とすると
 2
になります。

つまり、ブール型はTrueが1、Falseが0の値になっていると考えて構わないと思います。

それで、今度はnumpy.arrayでの扱いです。
aa = numpy.array( [ a, a, b] )
とすると
array([ True,  True, False], dtype=bool)

となり、bool型として扱われていることがわかります。
で、今度は
aa and aa
とするとエラーになります。

実行した結果のエラーメッセージは下のようになっていました。
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
ValueError: The truth value of an array with more than one element is ambiguous.
 Use a.any() or a.all()

numpy.arrayは他の型ではそれぞれの要素に関して個別に演算するので、ブール型もそうかと思うとそうでもないようです。
で、今度は

aa = numpy.array( [True, True, False], dtype=bool)
bb = numpy.array( [False, True, False], dtype=bool)
aa + bb
aa * bb
とすると、それぞれ下のようになりました。

>>> aa+bb
array([ True,  True, False], dtype=bool)
>>> aa*bb
array([False,  True, False], dtype=bool)

つまり、要素個別に演算する場合は+と*でandとorの代用する必要があるということです。
andとor,notは単独のTrue, Falseといったブール値にしか使えません。


notで実現していた、反転をさせる場合はどうするかというと、
~aa
とします。
結果は
array([False, False,  True], dtype=bool)
となります。

ブールのarrayは結構有用です。
たとえば、

x = numpy.random.randint( 0, 10, 12 )
c = x>5
print x[c]

というように、条件にあったものだけ抽出する際に使います。
cは
array([ True, False,  True, False,  True, False,  True,  True,  True,
       False, False, False], dtype=bool)
といった風にbool型のarrayになっています。







2014年5月14日水曜日

フラクタル 第2回

以前、マンデルブロ集合の図を作成しました。
今回pythonはあまり関係ないです。


マンデルブロ集合の画像を線分データのベクター画像(svg)として変換しました。
PDF


そのベクターデータを元に彫刻してみました。

工作関連は別なブログに載せてありますので、詳しくはそちらを見て下さい。


2014年5月12日月曜日

フラクタル

pythonで複素数の練習を兼ねて
フラクタル図形を作ってみました。
マンデルブロ集合というものです。


フラクタルが何かについては専門書または他のサイトなどで調べてみてください。
マンデルブロ集合の定義は簡単な複素数漸化式で記述されます。
つまりマンデルブロ図形というのはある複素数cについて、

z(n+1) = z(n)*z(n) + c
z(0) = 0

とした場合、n→∞で発散するか否かを色分けして複素平面に表示したものです。

複素数平面の点をランダムに選び、n回漸化式を解いてz(n)を作ります。
その時の、
log |z(n)|
の値を画素の輝度値としました。
forループを2重に使うと時間がかかるので、numpy.arrayで一気に計算します。

m = 10000
c = numpy.random.uniform( -2.5,2.5,m ) + 1j * numpy.random.uniform( -2.5,2.5,m )
z = numpy.zeros( m, dtype=numpy.complex )
for i in range(16):
     z = z**2 + c

と言う感じです。

opencvを使って画像化します。
しかし、案外時間が掛かるので、zの漸化式計算と、画像化を分けて処理しました。

Mandelbrot set

30Mポイントのサンプリングをして、そこから画像化しましたが、埋め切れていません。画像を拡大してみるとポツポツと黒い点が見えます。
実は細くつながっています。が、データ点数(サンプリング)で埋め切れていないのと画像化の解像度が荒いため、つながっているようには見えません。






2014年5月9日金曜日

複素数

python numpyでの複素数の扱い方についてです。


複素数は型としてはnumpy.complexと定義されています。
これを利用してみます。

つまり、

numpy.array( [3, 2, 1], dtype=numpy.complex )

とすると

 array([ 3.+0.j,  2.+0.j,  1.+0.j])

というふうに複素数配列となります。


複素数の演算もできます
a = numpy.array( [3, 2, 1], dtype=numpy.complex )
a ** 2

とすると2乗します。
実行結果は

 array([ 9.+0.j,  4.+0.j,  1.+0.j])

となります。

指数計算はそのまま掛けて計算されます。
つまり、python風な数式で表すと
z=x+jy
とすると
z**2 = (x+jy) * (x+jy) = x**2 - y**2 + j 2*x*y
となります。


複素数には色々と実数にはなかった計算方式が出てくるので注意が必要です。
実数部は
numpy.real(z)
虚数部は
numpy.imag(z)
とします。
これらの関数はzがarrayでもOKです。その場合、出力もarrayになります。

複素共役は
numpy.conj(z)
とします。
複素共役が何かわからない場合は検索してみて下さい。
誰かがわかりやすく説明してくれていると思います。

絶対値はnumpy.abs(z)
で計算できます。


複素数の乱数発生は簡単にはできませんでした。
一様分布の場合
z = numpy.random.uniform( -5, 5, 10 ) + 1j * numpy.random.uniform( -5, 5, 10 )
としました。
結果は

array([-4.48658518-4.19169279j,  2.39542899+4.87030019j,
        4.03412469-1.00727465j,  1.72313031+0.99345025j,
        0.68053779-0.385337j  ,  2.08333535-0.32240792j,
       -0.65638447-0.27364433j,  1.36415791-1.89435328j,
        0.19786534-4.51657661j,  2.44853396+3.04850488j])

となりました。
もっとスマートな方法があれば、是非ご教示下さい。





2014年5月7日水曜日

svgwriteのエラー

pythonでSVGファイルを生成する必要があったので、
svgwriteというモジュールをインストールしてみました。

インストールは簡単にできます。
pip install svgwrite
です。



サンプルを実行してみようとすると
エラー…

import下だけなのにエラーが出ました。
エラーメッセージを貼り付けると下のようにxml関連らしいです。

>>> import svgwrite
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "c:\Python26\Lib\site-packages\svgwrite\__init__.py", line 51, in <module
>
    from svgwrite.drawing import Drawing
  File "c:\Python26\Lib\site-packages\svgwrite\drawing.py", line 29, in <module>

    from svgwrite.container import SVG, Defs
  File "c:\Python26\Lib\site-packages\svgwrite\container.py", line 27, in <modul
e>
    from svgwrite.base import BaseElement
  File "c:\Python26\Lib\site-packages\svgwrite\base.py", line 12, in <module>
    from svgwrite.etree import etree
  File "c:\Python26\Lib\site-packages\svgwrite\etree.py", line 33, in <module>
    original_serialize_xml = etree._serialize_xml
AttributeError: 'module' object has no attribute '_serialize_xml'


下のサイトを参考に少し修正しました。
https://bitbucket.org/mozman/svgwrite/pull-request/3/fix-etreepy/diff


pythonインストールフォルダに
Lib/site-packages
というフォルダがあります。
この下にいわゆるライブラリがインストールされます。

インストールされていれば、
この下にsvgwriteというフォルダが有ると思います。
その中にetree.py
というファイルがあるので、

33行目
original_serialize_xml = etree._serialize_xml

を削除して、最終行に移動します。
現状のエラーでは循環呼び出しをしてしまったらしいです。

これで、ひとまずimportは出来ました。

使い方などはまた次の機会にレポートします。



2014年4月25日金曜日

arrayデータの16進パターン

numpyのarrayにはいくつかのメソッドがあります。

覚えておくとたまに便利になる機能がいくつかあります。
まず、tostringです。
動作確認しつつまとめてみます。

tostringメソッドは文字通り、arrayのデータをstringデータとして出力します。
数値を文字列に変換するわけではありません。
どういうときに使用するかというと、場合それぞれですが、
たとえば、シリアル通信やTCP/IP通信などで文字列として受け取ったデータを、数値として扱いたい場合などに使えます。
最近では、画像データの変換として使いました。


例を挙げます。
import numpy
x = numpy.arange(-3,3)
とするとxの中身は
array([-3, -2, -1,  0,  1,  2])
となります。
そこで、x
s = x.tostring()
とすると、sの中身は
'\xfd\xff\xff\xff\xfe\xff\xff\xff\xff\xff\xff\xff\x00\x00\x00\x00\x01\x00\x00\x00\x02\x00\x00\x00'
となり、16進のバイナリパターンが文字列として出力されています。
つまり、arangeで作成した段階で、4[byte](32bit)のint型として扱われ、リトルエンディアン方式で保管されているようです。

これを文字列として保存して、arrayとして戻すためには、
numpy.frombuffer( s, dtype=int )
とします。結果は、
array([-3, -2, -1,  0,  1,  2])
となります。
浮動小数点の場合はdtypeをfloatに指定します。
その他、32bitや64bit、または整数型か浮動小数点型か、符号有無などの細かい設定が必要な場合は、numpyの下に定義されている型を利用します。
たとえば

符号無整数
numpy.uint8
numpy.uint16
numpy.uint32
numpy.uint64

符号付整数
numpy.int8
numpy.int16
numpy.int32
numpy.int64

浮動小数点
numpy.float16
numpy.float32
...

複素数
numpy.complex
numpy.complex64
...

などです。



ちなみに、数値データをそれぞれ文字列に変換する場合は、いくつか方法がありますが、
a.astype( str )
とすればできます。
listとして出力するなら
map( str, a )
としてもよいかもしれません。
今度、それぞれの方法について処理時間などの比較をしてみます。



例で使用したnumpyは1.8.1 (python 2.6)を使用しています。




2014年4月24日木曜日

数理統計学・データサイエンスとpython


pythonを使って統計アルゴリズムを記述することができます。
現在、統計アルゴリズムをコーディングする上でいくつかの比較的簡便な方法が開発されてきました。統計言語Sを基にしたR言語もその一つです。

ここでは最近急速に発展してきたpythonを使った方法について考察してみます。
最近の流れとしては大規模データ、いわゆるビッグデータの分析に注目が集まっているようです。
pythonはこれまでにも書いてきたようにスクリプト言語で基本インタプリタですので、書き方によっては計算時間がネックになってしまうことがあります。
さらには、もともと統計分析を目的として開発されたSまたはR言語のような機能の豊富さには負けます。
しかし、pythonはそれを上回るほどの開発効率やソースコードの簡潔さが特徴で、ライブラリも多岐にわたるため応用範囲が広いという利点があります。

科学技術計算で使うためのライブラリとしては、有名なものだけでも

  • Numpy
  • Scipy
  • Sklearn
  • Pandas
  • Jubatus

などがあります。

Numpyについてはこれまでも何回か説明してきましたが、統計に限らず数値演算の基本的な機能を提供します。
Numpyのもっとも重要な機能はArrayだと思います。
Arrayというのはつまり配列ですが、python標準のlistより、ベクトル的な計算を意識したもので、これを効果的に使うことにより、効率的な計算が可能です。
NumpyにはMatrixという型もあります。
Matrixを使うと行列演算が簡単にかけるので便利ではありますが、
Arrayの方が使い勝手が良いと思います。


ScipyはNumpyではカバーしきれなかった科学技術寄りの機能を実現しています。
たとえば、統計アルゴリズムの例でいうと、乱数発生はNumpyの機能が使えます。つまり、
numpy.random.randn(100)
などとすると100個の正規乱数がnumpy.array形式で生成されます。
しかし、正規分布の確率密度関数のグラフが書きたいとなれば、
scipy.statsの機能が使えます。

import numpy
import scipy.stats
p = scipy.stats.norm.pdf( numpy.arange(-3,3,0.1) )

とすれば、pに-3から3までの正規密度関数の値が入ります。
正規分布のほかにもcauchy分布やgammaやvon misesなどの分布もあります。
統計の論文に出てきそうな標準的確率分布はあらかたあります。


残りのライブラリの説明はまた今度します。
興味があれば、それぞれググってみてください。