Pyscoで高速化

Pyscoってのを使うと早くなるといううわさがあるのでメモ
import psyco
psyco.full()
と冒頭に入れるらしい。

Super bound unboundについて

superという関数について調べる。


class super(object)
| super(type) -> unbound super object
| super(type, obj) -> bound super object; requires isinstance(obj, type)
| super(type, type2) -> bound super object; requires issubclass(type2, type)
| Typical use to call a cooperative superclass method:
| class C(B):
| def meth(self, arg):
| super(C, self).meth(arg)

super(type)っての要するにtypeの基底クラスを与えるのか？
では、
super(type,obj) ->bound super object.ってどういうこと？


しらべてみるとboundってのはメソッドがどっかのインスタンスに属している。unboundは属していないってことらしい。

boundメソッドは、そのまま関数として使える。unboundメソッドはそのまま関数として使おうとするとだめ。

次のコードを例にとって考えてみる。

class Bassclass:

    def __init__(self,str1):

        self.a = str1+"bass class"

    def func1(self):

        return self.a



print Bassclass().func1()

C++だと、コンストラクタにselfなんて渡さなかった（たぶん）けど、pythonだとselfを渡さないとだめ。また、何も引数がいらなさそうなfunc1()もselfを渡すような形で定義してやら無いとダメ。
上のコードだとBassclassっていうスコープの中だからいちいちself.aなんて書かなくてもいいのかと思ってたけどそうじゃないみたい。


疑問点が１つ。
class Bassclass(object)
という文があるが、これを
class Bassclass()
にするとなんか知らんけどだめ。

新スタイルクラスとクラシックスタイルクラスってのがあるらしい。
新スタイルクラスではobjectの派生クラスとして定義しないとダメらしい。
んで、新スタイルクラスではスタティックメソッドとクラスメソッドが使えるそうな。
スタティックメソッドは、インスタンスなしで使用できるメソッド。普通のクラスに属さないメソッドみたいなもんか。
クラスメソッドは、クラスを渡すメソッド。こっちは全く分からん。

class Bassclass(object):

    def __init__(self,str1):

        self.a = str1+"bass class"



    def func1(self):

        return self.a



class Derivedclass(Bassclass):

    def __init__(self, str1, str2):

        super(self.__class__, self).__init__(str1)

        self.b = str2

    def func1(self):

        return self.a

    def func2(self):

        return self.b



d = Derivedclass("ouou",12)

print d.func1()

print d.func2()

Gnuplot-py のサンプル(マンデルブロー集合) Gnuplot-py sample code (for Mandelbrot set)

Gnuplot-pyライブラリをインストールすると、どこかにdemo.pyができているのでこいつを参考にしてコードを書いた。

import psyco

psyco.full()

import Gnuplot, Gnuplot.funcutils

from scipy import *



STOP_VALUE = 4



#CX_AREA=[-1.495,-1.475]

#CY_AREA=[-0.02,0.045]

CX_AREA=[0,0.4]

CY_AREA=[-0.65,-0.35]

SIZE = 300

YSTEP = (CY_AREA[1]-CY_AREA[0])/SIZE



XSTEP = (CX_AREA[1]-CX_AREA[0])/SIZE



def mandelbrot_seq(x,y):

    mandelbrot_seq = 0.+0.j

    c=x+y*1j

    for n in range(0,100):

        mandelbrot_seq=mandelbrot_seq**2+c

        if abs(mandelbrot_seq)>STOP_VALUE:#avoid overflow(mandelbrot_seq is too large number)

            break

        else:

            pass

    return abs(mandelbrot_seq**2+c-mandelbrot_seq)





g = Gnuplot.Gnuplot(debug=1)

#g('set parametric')

#g('set data style lines')

g('set hidden')

g('set pm3d map')

#g('set contour base')

#g.title('An example of a surface plot')

g.xlabel('Re axis')

g.ylabel('Im axis')



x = arange(CX_AREA[0],CX_AREA[1],XSTEP)

y = arange(CY_AREA[0],CY_AREA[1],YSTEP)



g.splot(Gnuplot.funcutils.compute_GridData(x,y, mandelbrot_seq, binary=0))

g('set pm3d')をg('set pm3d map')に変更すると底面のマップのみが表示されるのでこっちのほうが見やすい。

マンデルブロー集合の描画 Rendering Mandelbrot set

pythonの練習でマンデルブロー集合を描画するコードを書いた。

Fig:実行結果　The result of execution

from scipy import *

from pylab import *

STOP_VALUE = 4

CX_AREA=[-1,0.5]

CY_AREA=[-1,1]

STEP = 0.01

def judge_conv_mandelbrot_seq(c):

    mandelbrot_seq = 0.+0.j

    for n in range(0,50):

        mandelbrot_seq=mandelbrot_seq**2+c

        if (abs(real(mandelbrot_seq))>STOP_VALUE)or(abs(imag(mandelbrot_seq))>STOP_VALUE):

            break

#    if abs(mandelbrot_seq)    if abs(mandelbrot_seq)>STOP_VALUE:

        return 0

    else:

        return 1

    

mandelbrot_set=[]

for cx in arange(CX_AREA[0],CX_AREA[1],STEP):

    for cy in arange(CY_AREA[0],CY_AREA[1],STEP):

        if judge_conv_mandelbrot_seq(complex(cx,cy)):

             mandelbrot_set.append(complex(cx,cy))

        else:

            pass

plot(real(mandelbrot_set),imag(mandelbrot_set),',')



show()

Gnuplot-pyを使った3次元グラフの書き方

3D plot をしたい。
matplotlibにplot3Dという関数があってこいつで3Dグラフをレンダリングできるらしいのだが、バージョン違いで上手く動かない。
そこでGnuplotを使う。

Sample Code

import Gnuplot

import numpy



gp = Gnuplot.Gnuplot(debug=1)

x=[0,1,0]

y=[0,1,0]

y2=[0,1,0]



d3=Gnuplot.Data(x,y,y2)

gp.splot(d3)

こんな感じでいける。
matplotlibと違って、いったんGnuplot.Data()でなんか知らないけど型を変換しないといけないらしい。
まだ良く分かっていないのでとりあえずメモ。

PythonでRC回路の過渡解析をしてみる その2 数値計算 by Scipy編

前回建てた微分方程式と初期条件
Vc'=1/(RC)(Vin-Vc)

Vin(0) = 0.0
Vc(0)= Vcc

前回、手計算で解いた結果
Vc(t)=Vcc×exp{-1/(RC)×t}

でも、定数変化法で解くので手間が掛かる。


手計算ではなく数値計算する方法。

Pythonの数値計算モジュールScipyを使う。図１に実行結果を示す。

図１：実行結果

from scipy import *

from pylab import *

deriv = lambda y,t : array([y[1],-y[0]])

# Initial parameters

R=10**3

C=10**-7

Vcc=5.





def Vin(t):

    return 0.



def Vcdash(Vc,t):

    return (Vin(t)-Vc[0])/(R*C)



start=0

end=10**-3

numsteps=10000

time=linspace(start,end,numsteps)



from scipy import integrate

Vc0=Vcc

Vc=integrate.odeint(Vcdash,Vc0,time)



#Plot Setting

figure(1,figsize=(6,4)).subplots_adjust(bottom=0.11)

xlabel('t')

grid(True)

hold(True)



plot(time,Vc)

#Setting Window title and Save Figure

title('Transient Analysis of RC Circuit')

savefig('transient analysis.png',dpi=72)



show()

integrate.odeint()関数には、微分方程式と初期値と時間に相当する行列を渡すと計算してくれます。
微分方程式はこの場合Vcdash()です。
なんで微分方程式が
def Vcdash(Vc,t):
return (Vin(t)-Vc[0])/(R*C)
って表現されるのか？その理由を以下に示します。

＜リストを用いて表現される理由＞
微分方程式を数値解析的に解くには差分方程式だったか漸化式だったかを用いる。（たしか）
差分方程式@wikipedia

解きたい微分方程式

Vc'=1/(RC)(Vin-Vc)

を前進差分で差分方程式にしてやると、
Vc'=(Vc(t+h)-Vc(t))/h=1/(RC)(Vin-Vc(t)) --(1)
ここで、hを数値計算するメッシュの分解能として考えると（＝メッシュがhずつ刻まれていると考えると）、
Vc(t) =V[n]
Vc(t+h)=v[n+1]
とおける。
すると(1)式は

(Vc[t+1]+V[t])/h=(Vin-Vc[t])/RC
V[t+1]=Vin/RC-Vc[t](h-1/RC)

となり、ここでhを0に近づければ、（計算機では０にはできないが）
Vc[t+1]≒1/(RC)(Vin-Vc[t]) --(2)


(1),(2)を見比べると、Vc'はVc[t+1]に、VcはVc[t]に置き換わっただけということが分かる。ということで、微分方程式はリストを使って表現できる。
2階微分方程式の場合、例えば、

y''=ay'+by+c

を解きたい場合には、

y'=x

とおくと

x'=ax+by+c

のように、1階常微分方程式になりscipy.integrate.odeint()関数で扱える形式になる。
ただし、今度は、yの初期値だけでなくy'=xの初期値も必要となってくるので注意が必要。
参考Scipy CookBook




ついでに、このプログラムのVin(t)という関数は電源電圧を表現している、これをsin(ωt)を返す関数にした例↓。



図２:Vin(t)=Vcc*sin(10**5*t)とした場合の実行結果

PythonでRC回路の過渡解析をしてみる その１ 手計算

↓の回路を対象に計算。最終的にはVcの式が知りたい。



まずは微分方程式を建てる。

Vc=1/C∫Idt --(1)
I=Vr/R --(2)
Vin-Vc=Vr --(3)
この３つの式をまとめる。
Vcをみたいので、VcとVinの関係式に持ち込むという方針で変形する。
(3)を(2)に代入すると、
I=(Vin-Vc)/R --(3')
(3')を(1)に代入すると
Vc=1/C∫(Vin+Vc)/Rdt
Rは定数なので移動して
Vc=1/(RC)∫(Vin-Vc)dt --(微分方程式)


<初期条件を建てる>
次に初期値を設定。
設定としては、Vinの電圧がVcc→0.0に落ちたということにします。あとは、コンデンサはフルチャージされていることとする。

Vin(-0) = Vcc
Vin(+0) = 0.0
Vc(0)= Vcc（コンデンサには電荷が完全チャージされている状態を想定。）

上で導出した微分方程式を並べる。

Vc=1/(RC)∫(Vin-Vc)dt --(微分方程式)
話が前後しますが、微分方程式を微分形式に直します。両辺を微分すると、

Vc'=1/(RC)(Vin-Vc) --(4)

微分したら積分定数が消える気がするが、積分定数とはすなわち初期条件のことなので（たぶん）ＯＫ。

<手計算で解く>
この微分方程式は、非斉次方程式ですので、定数変化法で解きます。

まずは、斉次方程式の一般解を求める。
Vc'=1/(RC)(Vin-Vc)
を斉次化するには
非斉次項Vinを
Vin=0とおけばよいので
Vc'=-1/(RC)Vc --(7)

次に、(5)に初期条件を代入してやると、
Vc(0)=Vcc=A×exp(0)
∴A=Vcc
よって
斉次方程式の一般解＝Vcc×exp{-1/(RC)×t} --(8)

(8)より
非斉次方程式の一般解＝B(t)×Vcc×exp{-1/(RC)×t}
となっているはずなので、B(t)を求める
初期条件を代入すると、
B(t)=1
となる。よって、

Vc(t)=Vcc×exp{-1/(RC)×t}
となる。
この式は、t=0でVc=Vccとなり、t>0の領域で指数関数的に減少するカーブを描く。
これはコンデンサの放電時の波形。