Jupyter の ipynb ファイルの中身を見たい
2020.09.20投稿
必要性
Jupyter notebook 上でコードを書いていて、「あのとき作った Jupyter のコードをコピーしたい!」ということがよくある。その ipynb ファイルが、たまたま今開いている Jupyter のディレクトリの下にあればいいのだが、そうでない場合は、一度 Jupyter notebook を終了・再起動し所要の ipynb ファイルを開いてエディタにコピー、もとの作業に戻るため Jupyter notebook を再起動するという手間をかけていた。 なんとかならないかと思いながら調べていたら、効率よく他の ipynb ファイルの中身を確認する以下の3つの方法があることがわかったので、紹介しておく。 もしここに記載の方法を使ったことにより不都合なことが起きても、筆者は知りません。 自分の責任においてトライしてください。
VS code を使う
こちらのサイトに方法が記載されている。
私は VS Cod を使う方法を選択し使ってみた。
Mac の Get info の Open with: で、VS code をデフォルトにしておけば、ipynb ファイルをダブルクリックすることで開くことができる。
Jupyter lab を使う
Jupyter Lab なら複数のタブを開ける! そこでさっそくインストールして調子を見ているところです。 私のマシンではなんとなく動きがもっさりしている感じ。
複数のブラウザで Jupyter を開く
Jupyter notebook でも これ(http://localhost:8888/?token=xxxxxx......
)を url に入れてやると違うブラウザで Jupyter notebook を開ける。
例えば、私の場合、default は Safari であるが、terminal に出てくるこれ(http://localhost:8888/?token=xxxxxx......
)を Firefox や Chrome の url にいれることにより、それぞれのブラウザで notebook を開ける。
以 上
LG 27インチ4Kディスプレイ購入
9月7日注文、9月8日到着で、LGの27inch 4K display を購入しました。¥57,980也。ちょっと奮発しすぎたか。
4月より在宅勤務しています。会社に行くのは月に2〜3回ですかね。 しかし7月よりお客さんから設計の仕事をいただき、今はかなり忙しい。まあ色々検討してレポートを書いているのですが、ディスプレイ上の作業領域が「もっと欲しい!」ということで、サブディスプレイを購入したのでした。LGにした理由は、USB-Cポートがついていたからというのが大きいです。このポートがあると、ディスプレイとiMacをUSB-Cケーブル1本でつなぐことができます。 最初は1万円くらいのやつでいいかと思っていたのですが、色々Webを見ながら考えていたら、「目も悪くなってきたし、せっかくなら4Kにするか!」と思い立ったのでした。実際には27インチで4Kだと、100%のscaleでは文字が小さくて読めません。なので200%とかに拡大しています。「拡大したら4Kの意味ないじゃん」という向きもありますが、画面のドットそのものは精細なので、文字はとてもきれいです。 買った当初は縦置きにしようと思い、試しに縦置きで使ってみたのですが、上の方がとても遠く感じます。そこで通常の横置きに。図面などを見るにも視野が広くなり便利です。このディスプレイは横から縦へ回転させることができるので、縦置きにするのは、南北に長い地形図やGoogle earthを見るときだけとし、普段は横置きとすることにしました。
はじめは縦置きにしてみましたが、あまり使いよくない。
基本は横置きです。図面を見る時など視野が広くて便利。
Google Earthなどで縦に画像を表示させたい場合は便利。
とりあえず持っている機械類の集合写真。この写真はiPhoneで撮っているので、かわいそうに、iPhoneは写っていません。 正面が新しいディスプレイ、左側のWin機は会社支給のPC.
つなぎ方のイメージ図。iMac以外はHDMIでつないでいます。
以 上
Python 逆調整池運用解析修正版(Flood Routine)
逆調整池運用解析プログラムの修正版
- 関数FLOODの中の収束計算で、前回貯水池水位のみを既知としていたが、本来貯水池水位と放流量はペアなので、前回放流量も既知とした。プログラム作成時になぜ収束計算の中で前回流量を既知として確定していなかったのかは思い出せない。
- 関連出力は、データフレームにして一気に書き出すことにした。
以下最新版プログラム(2020.08.08)
# ============================================== # Flood Routine Analysis # ============================================== import sys import numpy as np import pandas as pd from scipy import interpolate def FLOOD(iOFC,fnameR1,fnameR2,fnameR3,ram,qrel): nr,rh,rv=INP_HV(fnameR1) # reservoir H-V nd,ti,q_in,ELini,Qoini=INP_INFLOW(fnameR2) # inflow time history no,elof,qof=INP_OUTFLOW(fnameR3) # outflow capacity pEL=np.zeros(nd,dtype=np.float64) # storage of reervoir water level pQo=np.zeros(nd,dtype=np.float64) # storage of outflow pDE=np.zeros(nd,dtype=np.float64) # storage of error pVO=np.zeros(nd,dtype=np.float64) # storage of reservoir volume pRA=np.zeros(nd,dtype=np.float64) # storage of ramp up rate pUD=np.zeros(nd,dtype=int) # storage of indicator od water level Vin=np.sum(q_in)*(ti[2]-ti[1])/2 # total volume of inflow Qimax=np.max(q_in) # Initial setting elv=ELini # elv=WLINI(q_in[0],elof,qof,no) # Initial reservoir water level EL=elv VOL=RET_V(nr,rh,rv,elv) # Initial reservoir volume q_out=Qoini qout0=Qoini print('q_in0={0:10.3f}'.format(q_in[0])) print('elv_ini={0:10.3f}'.format(EL)) print('q_ini={0:10.3f}'.format(q_out)) i=0 iUD=0 pUD[i]=iUD pEL[i]=EL pQo[i]=q_out pDE[i]=EL-elv pVO[i]=VOL pRA[i]=0 # Iterative calculation iUD=1 # Assume reservoir water level increase dh=0.0005 # Water level increment for searching equilibrium point eps=0.001 # Convergence criteria itmax=int(1.0/dh)*100 Vf=RET_V(nr,rh,rv,60.0) Vz=0 Vin2=0 Vot2=0 q1=q_out for i in range(0,nd-1): dt=ti[i+1]-ti[i] # Time interval Qin=0.5*(q_in[i]+q_in[i+1]) # Average inflow rate during dt hh=0.0 # Zero set of water level increment k=0 j=0 while 1: k=k+1 j=j+1; hh=hh+float(iUD)*dh # Update waler level increment elv1=EL elv2=EL+hh if iOFC==8: q2=OFC8(q1,qout0,elv2,elof,qof,no,ti[i],ram,qrel) Qout=0.5*(q1+q2) # Average outflow rate during dt dS=(Qin-Qout)*dt*3600.0 # Storage accumulated during dt R=VOL+dS # Total reservoir volume elv=RET_H(nr,rh,rv,R) # Reservoir water level at reservoir volume R if iUD==1 and j==10: # 10 times trial assuming reservoir water level increase if EL+hh>elv: # if reservoir water level is less than assumed water level, iUD=-1 # assume water level decrease hh=0.0 j=0 if iUD==-1 and j==10: # 10 times trial assuming reservoir water level decrease if EL+hh<elv: # if reservoir water level is greater than assumed water level, iUD=1 # assume water level decrease hh=0.0 j=0 if np.abs(EL+hh-elv)<eps: break # Judgement of convergence if itmax<k: break VOL=R # Cumulative volume EL=EL+hh # Elevation q_out=q2 # Outflow q1=q2 pUD[i+1]=iUD pEL[i+1]=EL pQo[i+1]=q_out pDE[i+1]=EL-elv pVO[i+1]=VOL pRA[i+1]=(pQo[i+1]-pQo[i])/(ti[i+1]-ti[i]) hmax=max(pEL)-ELini sys.stdout.write('Time: %10.3f\n'% np.max(ti)) sys.stdout.write('h : %10.3f\n'% hmax) sys.stdout.write('Qin : %10.3f %10.3f\n'% (np.min(q_in),np.max(q_in))) sys.stdout.write('Qout: %10.3f %10.3f\n'% (np.min(pQo),np.max(pQo))) sys.stdout.write('EL : %10.3f %10.3f\n'% (np.min(pEL),np.max(pEL))) sys.stdout.write('\n') # TWC of 80m downstream of Re-Regulating Dam fnameR='inp_twc.txt' data=np.loadtxt(fnameR,skiprows=0,usecols=(0,1)) qq=data[:,0] ee=data[:,1] ff=interpolate.interp1d(qq,ee) el_ds=ff(pQo) df = pd.DataFrame({ 'i' : np.arange(0,nd), # numper 'iUD': pUD, # indicator for water level 'ti' : ti, # time 'EL' : pEL, # water level of reservoir 'dEL': pDE, # error 'VOL': pVO, # reservoir volume 'Qi' : q_in, # inflow to reservoir 'Qo' : pQo, # outflow from reservoir 'RA' : pRA, # ramp up rate 'ELD': el_ds # downstream water level }) df = df.set_index('i') return df def WLINI(q,elof,qof,no): # To obtain initial water level at oitflow discharge q for i in range(0,no-1): if qof[i]<=q and q<=qof[i+1]: break if qof[no-1]<q: i=no-2 x1=qof[i] y1=elof[i] x2=qof[i+1] y2=elof[i+1] a=(y2-y1)/(x2-x1) b=(x2*y1-x1*y2)/(x2-x1) elv=a*q+b if q<0.0: elv=elof[0] return elv def OFC8(q1,qout0,elv2,elof,qof,no,tii,ram,qrel): Qlim=300.0 if 5.0<tii and tii<= 29.0 :qout0=qrel if 29.0<tii and tii<= 53.0 :qout0=qrel if 53.0<tii and tii<= 77.0 :qout0=qrel if 77.0<tii and tii<=101.0 :qout0=qrel if 101.0<tii: qout0=39.0 q2=OFC(elv2,elof,qof,no) rq2=q2 if q2<qout0: qout0=q2 q2=qout0 st=0.01 ts1=0.0 ts2=24.0 ts3=48.0 ts4=72.0 ts5=96.0 if ts1+st<=tii and tii<=ts1+5.0+st: if elv2<68.0: q2=q1+ram*st if rq2<q2: q2=rq2 if Qlim<=q2: q2=Qlim else: q2=rq2 if ts2+st<=tii and tii<=ts2+5.0+st: if elv2<68.0: q2=q1+ram*st if rq2<q2: q2=rq2 if Qlim<=q2: q2=Qlim else: q2=rq2 if ts3+st<=tii and tii<=ts3+5.0+st: if elv2<68.0: q2=q1+ram*st if rq2<q2: q2=rq2 if Qlim<=q2: q2=Qlim else: q2=rq2 if ts4+st<=tii and tii<=ts4+5.0+st: if elv2<68.0: q2=q1+ram*st if rq2<q2: q2=rq2 if Qlim<=q2: q2=Qlim else: q2=rq2 if ts5+st<=tii and tii<=ts5+5.0+st: if elv2<68.0: q2=q1+ram*st if rq2<q2: q2=rq2 if Qlim<=q2: q2=Qlim else: q2=rq2 return q2 def OFC(elv,elof,qof,no): # Free overflow discharge for i in range(0,no-1): if elof[i]<=elv and elv<=elof[i+1]: break if elof[no-1]<elv: i=no-2 x1=elof[i] y1=qof[i] x2=elof[i+1] y2=qof[i+1] a=(y2-y1)/(x2-x1) b=(x2*y1-x1*y2)/(x2-x1) q=a*elv+b if elv<=elof[0]: q=0.0 return q def RET_V(nr,rh,rv,elv): # To obtain reservoir cumulative volume from the water level for i in range(0,nr-1): if rh[i]<=elv and elv<=rh[i+1]: break if rh[nr-1]<elv: i=nr-2 x1=rv[i] y1=rh[i] x2=rv[i+1] y2=rh[i+1] a=(y2-y1)/(x2-x1) b=(x2*y1-x1*y2)/(x2-x1) v=(elv-b)/a return v def RET_H(nr,rh,rv,v): # To obtain reservoir water level from cumulative volume for i in range(0,nr-1): if rv[i]<=v and v<=rv[i+1]: break if rv[nr-1]<v: i=nr-2 x1=rv[i] y1=rh[i] x2=rv[i+1] y2=rh[i+1] a=(y2-y1)/(x2-x1) b=(x2*y1-x1*y2)/(x2-x1) elv=a*v+b return elv def INP_HV(fnameR1): # Input reservoir H-V data fin=open(fnameR1,'r') dat=fin.readlines() fin.close() n=len(dat) text=dat[0] text=text.strip() text=text.split() vcoef=float(text[0]) rh=np.zeros(n-1,dtype=np.float64) rv=np.zeros(n-1,dtype=np.float64) for i in range(0,n-1): text=dat[i+1] text=text.strip() text=text.split() rh[i]=float(text[0]) rv[i]=float(text[1])*vcoef nr=len(rh) return nr,rh,rv def INP_INFLOW(fnameR2): # Input time sequence of inflow fin=open(fnameR2,'r') dat=fin.readlines() fin.close() n=len(dat) text=dat[0] text=text.strip() text=text.split() ELini=float(text[0]) Qoini=float(text[1]) ti =np.zeros(n-1,dtype=np.float64) q_in=np.zeros(n-1,dtype=np.float64) for i in range(0,n-1): text=dat[i+1] text=text.strip() text=text.split() ti[i] =float(text[0]) q_in[i]=float(text[1]) nd=len(ti) return nd,ti,q_in,ELini,Qoini def INP_OUTFLOW(fnameR3): # Input outflow capacity (water level-discharge curve) fin=open(fnameR3,'r') dat=fin.readlines() fin.close() n=len(dat) elof=np.zeros(n,dtype=np.float64) qof =np.zeros(n,dtype=np.float64) for i in range(0,n): text=dat[i] text=text.strip() text=text.split() elof[i]=float(text[0]) qof[i] =float(text[1]) no=len(elof) return no,elof,qof def main(): lout=['out_150.csv','out_175.csv','out_200.csv','out_225.csv','out_250.csv'] lram=[150,175,200,225,250] lqrel=[89.0,87.1,82.9,81.5,80.7] for out,ram,qrel in zip(lout,lram,lqrel): iOFC=8 fnameR1='inp_RHV.txt' # H-V data of reservoir fnameR2='inp_c00.txt' # Inflow time history (Hydrograph) fnameR3='inp_OF8.txt' # Outflow capacity fnameW = out # Output file name df=FLOOD(iOFC,fnameR1,fnameR2,fnameR3,ram,qrel) df.to_csv(fnameW, sep=",") #============================ if __name__ == '__main__': main()
以 上
Relationship between critical depth and discharge coefficient for free overflow on dam spillway
In general, the discharge for free overflow on dam spillway can be expressed as follow:
On the other hand, a critical depth of rectangular cross section channel can be expressed as follow:
When it is assumed the water depth at the overflow crest of the dam is equal to the critical depth, the design head including approach velocity head can be expressed as follow:
Considering a unit length of overflow crest which means ,
From above, can be expressed as follow:
Therefore, discharge coefficient becomes a constant and it can be expressed as follow:
That is all.
Python:計算結果をエクセル出力し加工する(月平均値の算出とエクセル化)
能書き
報告書を作成しているとき、Pythonで計算した結果をWordに表形式で貼り付けたい場合がある。 これが、TeXやhtmlであれば、プロブラム内でタグを埋め込んだ形で標準出力し、コピー&ペーストでTeXもしくはhtml文書に貼り付けるのだが、相手がWordとなるとなかなか面倒である。PythonからWordを制御することができることも知っている(使ったことはない)が、そのためのみのコードを書くのも大変だし、後処理計算での活用などを考えるとエクセル出力しておき、それをwordに貼り付けるほうが、実用的である。
私がよく使う数表をwordに貼り付ける方法は以下の3つである。
(1) 計算結果をコンマ区切りで標準出力し、Wordにコピー&ペースト。その後はWordで[Insert]=>[table]=>[Convert Text to Table]=>[Separate test at Commas]=.[ok]の手順で、Word上で表に加工する方法。Word上で加工できるので仕上がりは結構綺麗にいく。
(2) 計算結果をエクセルに出力しフォーマットを整える。その後はエクセル上で範囲指定しクリップボードコピー、Word上で[Paste Special]=>[As Picture(PNG)]=>[ok]の手順で、画像としてWordに貼り付ける方法。画像としてwordに貼り付けるためWord上で加工できず、結果がぼやける場合もあり、最近は使わない。
(3)計算結果の表をhtml出力してブラウザで表示し、それをコピー&ペーストでwordに貼り付け、仕上げはword上で行う。ブラウザはsafariで行うのがやりやすいようである、
方法(3)は最近使い出したのだが、面倒なようでなかなか良い。 というのも計算結果は一回出力するのみでなく、次の処理で使う場合も多く、そうなると1ファイルに複数シートを配置できるエクセルで保存しておくのも後工程での使い勝手がよい。
エクセルからhtmlへの変換
エクセルからhtmlへの出力は、ネット上で便利なものが有り、以下のものを使わせてもらっている。
エクセルから変換したhtmlの表を貼り付けたものが下の表。二重線が消えているがhtml出力が最終成果ではないので気にしない。 必要あれば装飾をかければ良い。
Year | Jan | Feb | Mar | Apr | May | Jun | Jul | Aug | Sep | Oct | Nov | Dec | Ave |
2001 | 19.6 | 17.8 | 16.1 | 14.5 | 21.3 | 36.1 | 35.9 | 84.0 | 36.0 | 47.7 | 39.4 | 26.7 | 33.1 |
2002 | 22.1 | 20.4 | 18.5 | 16.7 | 24.1 | 19.9 | 43.1 | 36.7 | 44.8 | 32.9 | 29.6 | 23.8 | 27.8 |
2003 | 21.4 | 19.7 | 18.0 | 16.2 | 14.6 | 13.9 | 20.6 | 22.9 | 30.7 | 29.9 | 28.1 | 21.5 | 21.5 |
2004 | 19.2 | 17.6 | 15.9 | 14.3 | 21.1 | 37.7 | 47.0 | 48.7 | 110.6 | 47.1 | 32.0 | 24.2 | 36.2 |
2005 | 22.0 | 20.0 | 18.1 | 16.3 | 14.9 | 15.0 | 44.5 | 49.7 | 92.6 | 45.4 | 32.2 | 25.6 | 33.0 |
2006 | 22.3 | 20.5 | 18.6 | 16.7 | 15.2 | 16.1 | 32.7 | 65.6 | 104.4 | 62.3 | 34.7 | 25.6 | 36.3 |
2007 | 22.8 | 20.8 | 18.8 | 16.9 | 17.7 | 20.3 | 24.5 | 35.7 | 29.5 | 33.4 | 26.8 | 21.8 | 24.1 |
2008 | 20.0 | 18.3 | 16.6 | 15.0 | 13.7 | 25.0 | 56.6 | 83.9 | 74.4 | 41.0 | 36.4 | 28.9 | 35.9 |
2009 | 24.1 | 22.2 | 20.1 | 18.1 | 16.4 | 15.3 | 18.0 | 37.2 | 31.9 | 29.1 | 25.4 | 20.2 | 23.2 |
2010 | 18.6 | 17.0 | 15.3 | 13.8 | 12.4 | 11.8 | 34.5 | 114.9 | 89.1 | 45.6 | 32.9 | 24.4 | 36.0 |
2011 | 21.8 | 19.9 | 18.0 | 16.4 | 17.9 | 28.9 | 42.7 | 129.4 | 111.2 | 55.2 | 36.3 | 27.8 | 43.9 |
2012 | 24.9 | 22.7 | 20.4 | 18.4 | 17.8 | 25.9 | 41.6 | 60.0 | 72.2 | 42.4 | 33.8 | 26.5 | 33.9 |
2013 | 23.6 | 21.6 | 19.6 | 17.6 | 15.9 | 16.2 | 28.0 | 103.1 | 82.5 | 52.8 | 38.4 | 28.2 | 37.4 |
2014 | 23.7 | 21.8 | 19.8 | 17.8 | 16.0 | 18.3 | 57.8 | 56.8 | 58.2 | 34.4 | 28.8 | 23.1 | 31.5 |
2015 | 21.2 | 19.3 | 17.5 | 15.7 | 14.6 | 18.5 | 20.7 | 43.4 | 29.2 | 28.2 | 24.8 | 20.3 | 22.8 |
2016 | 18.8 | 17.1 | 15.5 | 14.0 | 12.6 | 13.6 | 39.1 | 48.7 | 58.3 | 35.8 | 29.9 | 23.9 | 27.3 |
2017 | 21.0 | 19.5 | 17.7 | 16.0 | 14.5 | 14.5 | 29.3 | 36.6 | 51.0 | 33.2 | 29.9 | 23.8 | 25.6 |
2018 | 21.1 | 19.6 | 17.8 | 16.1 | 14.8 | 33.8 | 40.6 | 58.9 | 41.9 | 64.7 | 35.2 | 26.6 | 32.7 |
2019 | 23.2 | 21.3 | 19.4 | 17.4 | 15.6 | ||||||||
2020 | |||||||||||||
Ave | 21.7 | 19.8 | 18.0 | 16.2 | 16.4 | 21.2 | 36.5 | 62.0 | 63.8 | 42.3 | 31.9 | 24.6 | 31.2 |
RGS1 (CA=3062km2) tank model monthly and annual discharge, unit: m3/s |
プログラム上の処理
Pythonプログラム上での出力したい表の処理は以下の手順で行う。
- 1.pandasのデータフレーム作成
- 2.一度データフレームをそのままエクセルに保存
- 3.保存したエクセルファイルを再度pythonで呼び出し、エクセうの書式指定、罫線などの可能を行い、再保存する。
データフレーム作成
df = pd.DataFrame({ 'Year' : lyy, 'Jan' : qqm[0:n,0], 'Feb' : qqm[0:n,1], 'Mar' : qqm[0:n,2], 'Apr' : qqm[0:n,3], 'May' : qqm[0:n,4], 'Jun' : qqm[0:n,5], 'Jul' : qqm[0:n,6], 'Aug' : qqm[0:n,7], 'Sep' : qqm[0:n,8], 'Oct' : qqm[0:n,9], 'Nov' : qqm[0:n,10], 'Dec' : qqm[0:n,11], 'Ave' : qqm[0:n,12] }) df = df.set_index('Year')
データフレーム保存
関連する内容を、1ファイルに、複数シートに分割して保存する。
fnameW='out_dis_month.xlsx' with pd.ExcelWriter(fnameW) as writer: df1_tank.to_excel(writer, sheet_name='RGS1_tank') df2_tank.to_excel(writer, sheet_name='RGS2_tank') df1_mlp.to_excel(writer, sheet_name='RGS1_mlp') df2_mlp.to_excel(writer, sheet_name='RGS2_mlp')
エクセルファイルの再読み込みと保存
この間に処理を記載する。
wb =openpyxl.load_workbook(fnameW) # .... # wb.save(fnameW)
フォーマット指定
# format for i in range(2, n+1): for j in range(2, m+1): ws.cell(row=i, column=j).number_format = '0.0'
罫線(前半は1重線、後半は2重線を指定)
from openpyxl.styles.borders import Border, Side # border line bcc = Border(top=Side(style='thin', color='000000'), bottom=Side(style='thin', color='000000'), left=Side(style='thin', color='000000'), right=Side(style='thin', color='000000') ) for i in range(1, n+1): for j in range(1, m+1): ws.cell(row=i, column=j).border = bcc bcc = Border(top=Side(style='double', color='000000'), bottom=Side(style='double', color='000000'), left=Side(style='double', color='000000'), right=Side(style='double', color='000000') ) for j in range(1, m+1): ws.cell(row=n, column=j).border = bcc for i in range(1, n+1): ws.cell(row=i, column=m).border = bcc
セルの結合(長い文字列を入力しセルを結合する)
ws.cell(row=n+1,column=1).value='RGS2 (CA=7860km2) tank model monthly and annual discharge, unit: m3/s' ss='A{0}:N{1}'.format(n+1,n+1); ws.merge_cells(ss)
セルの塗りつぶし
このプログラムには含まれていないが、塗りつぶしを行う場合は以下の要領で行う。
from openpyxl.styles import PatternFill # cell color fill = PatternFill(patternType='solid', fgColor='00ffff') for i in [6,11,16,21,26,31]: for j in range(3,17): ws.cell(row=i,column=j).fill = fill
主要計算部分の説明
例外処理を使っているので、その部分を解説。 このプログラムでは、1日でも欠測があれば、その月および年の平均は欠測としている。
main
で数表にする年と月を、リスト lyy
と lmm
で定義している。
lyy=['2001','2002','2003','2004','2005','2006','2007','2008','2009','2010', '2011','2012','2013','2014','2015','2016','2017','2018','2019','2020'] lmm=['01','02','03','04','05','06','07','08','09','10','11','12']
以下の関数で各月および各年の平均値を算出している。 引数の意味は以下の通り。
rf
: 日データを示すseries.lyyy
: リストlyy
と同じlmmm
: リストlmm
と同じ。
def cal_m(rf,lyyy,lmmm): qqm=np.zeros((len(lyyy)+1,len(lmmm)+1),dtype=np.float64) kda=np.zeros((len(lyyy)+1,len(lmmm)+1),dtype=np.float64) nac=0 for i,yy in enumerate(lyyy): for j,mm in enumerate(lmmm): ss=yy+'/'+mm # 処理する年/月を指定 try: # Series: rf[ss]に対する処理(合計およびデータ数カウント)を行う qqm[i,j]=rf[ss].sum() # sum in a month kda[i,j]=rf[ss].count() # availavle days nac=np.count_nonzero(np.isnan(rf[ss])) # unavailable days except KeyError: # もし[ss]に合致する年/月のデータがなければ配列にnp.nanを代入 qqm[i,j]=np.nan kda[i,j]=np.nan if 0<nac: # 仮に指定した年/月のデータがあってもnp.nanがあれば、すなわち月の全日のデータが揃っていなければnp.nanを代入 qqm[i,j]=np.nan kda[i,j]=np.nan qqm[i,-1]=np.sum(qqm[i,:]) # sum in a year 年間合計値 kda[i,-1]=np.sum(kda[i,:]) # sum in a year of available days 年間合計値 for j in range(len(lmmm)+1): qqm[-1,j]=np.nansum(qqm[:,j]) # nanを除外したデータの合計(各月合計値) kda[-1,j]=np.nansum(kda[:,j]) # nanを除外したデータ数の合計(各月合計値) qmean=qqm/kda # 月および年のデータ合計をデータ数で除すことにより平均値を算出 return qmean
プログラム例(全文)
プログラムは、ファイルから値を読み取り、上に示した表を作るもの。 実際のプログラムの出力はエクセル。 以下にPythonによるプログラム全文を示す。
import pandas as pd import numpy as np import datetime import openpyxl from openpyxl.styles.borders import Border, Side def formx(ws,n,m): # format for i in range(2, n+1): for j in range(2, m+1): ws.cell(row=i, column=j).number_format = '0.0' # border line bcc = Border(top=Side(style='thin', color='000000'), bottom=Side(style='thin', color='000000'), left=Side(style='thin', color='000000'), right=Side(style='thin', color='000000') ) for i in range(1, n+1): for j in range(1, m+1): ws.cell(row=i, column=j).border = bcc bcc = Border(top=Side(style='double', color='000000'), bottom=Side(style='double', color='000000'), left=Side(style='double', color='000000'), right=Side(style='double', color='000000') ) for j in range(1, m+1): ws.cell(row=n, column=j).border = bcc for i in range(1, n+1): ws.cell(row=i, column=m).border = bcc def wexcel(df1_tank,df2_tank,df1_mlp,df2_mlp): fnameW='out_dis_month.xlsx' with pd.ExcelWriter(fnameW) as writer: df1_tank.to_excel(writer, sheet_name='RGS1_tank') df2_tank.to_excel(writer, sheet_name='RGS2_tank') df1_mlp.to_excel(writer, sheet_name='RGS1_mlp') df2_mlp.to_excel(writer, sheet_name='RGS2_mlp') n=len(df1_tank.index)+1 m=len(df1_tank.columns)+1 wb =openpyxl.load_workbook(fnameW) # ws = wb.get_sheet_by_name('RGS1_tank') formx(ws,n,m) ws.cell(row=n+1,column=1).value='RGS1 (CA=3062km2) tank model monthly and annual discharge, unit: m3/s' ss='A{0}:N{1}'.format(n+1,n+1); ws.merge_cells(ss) # ws = wb.get_sheet_by_name('RGS2_tank') formx(ws,n,m) ws.cell(row=n+1,column=1).value='RGS2 (CA=7860km2) tank model monthly and annual discharge, unit: m3/s' ss='A{0}:N{1}'.format(n+1,n+1); ws.merge_cells(ss) # ws = wb.get_sheet_by_name('RGS1_mlp') formx(ws,n,m) ws.cell(row=n+1,column=1).value='RGS1 (CA=3062km2) MLP monthly and annual discharge, unit: m3/s' ss='A{0}:N{1}'.format(n+1,n+1); ws.merge_cells(ss) # ws = wb.get_sheet_by_name('RGS2_mlp') formx(ws,n,m) ws.cell(row=n+1,column=1).value='RGS2 (CA=7860km2) MLP monthly and annual discharge, unit: m3/s' ss='A{0}:N{1}'.format(n+1,n+1); ws.merge_cells(ss) # wb.save(fnameW) def rdata(fnameR): df=pd.read_csv(fnameR, header=0, index_col=0) # read excel data df.index = pd.to_datetime(df.index, format='%Y/%m/%d') return df def cal_m(rf,lyyy,lmmm): qqm=np.zeros((len(lyyy)+1,len(lmmm)+1),dtype=np.float64) kda=np.zeros((len(lyyy)+1,len(lmmm)+1),dtype=np.float64) nac=0 for i,yy in enumerate(lyyy): for j,mm in enumerate(lmmm): ss=yy+'/'+mm try: qqm[i,j]=rf[ss].sum() # sum in a month kda[i,j]=rf[ss].count() # availavle days nac=np.count_nonzero(np.isnan(rf[ss])) # unavailable days except KeyError: qqm[i,j]=np.nan kda[i,j]=np.nan if 0<nac: qqm[i,j]=np.nan kda[i,j]=np.nan qqm[i,-1]=np.sum(qqm[i,:]) # sum in a year kda[i,-1]=np.sum(kda[i,:]) # sum in a year of available days for j in range(len(lmmm)+1): qqm[-1,j]=np.nansum(qqm[:,j]) kda[-1,j]=np.nansum(kda[:,j]) qmean=qqm/kda return qmean def qq_mon(lyy,lmm,qq): qqm=cal_m(qq,lyy,lmm) lyy=lyy+['Ave'] n=len(lyy) df = pd.DataFrame({ 'Year' : lyy, 'Jan' : qqm[0:n,0], 'Feb' : qqm[0:n,1], 'Mar' : qqm[0:n,2], 'Apr' : qqm[0:n,3], 'May' : qqm[0:n,4], 'Jun' : qqm[0:n,5], 'Jul' : qqm[0:n,6], 'Aug' : qqm[0:n,7], 'Sep' : qqm[0:n,8], 'Oct' : qqm[0:n,9], 'Nov' : qqm[0:n,10], 'Dec' : qqm[0:n,11], 'Ave' : qqm[0:n,12] }) df = df.set_index('Year') return df def main(): lyy=['2001','2002','2003','2004','2005','2006','2007','2008','2009','2010', '2011','2012','2013','2014','2015','2016','2017','2018','2019','2020'] lmm=['01','02','03','04','05','06','07','08','09','10','11','12'] fname_list=[ 'df_rgs1_tank_result.csv', 'df_rgs2_tank_result.csv', 'df_rgs1_mlp_result.csv', 'df_rgs2_mlp_result.csv' ] for iii,fnameR in enumerate(fname_list): df0=rdata(fnameR) # daily discharge df0=df0['2001/01/01':'2019/05/31'] qq=pd.Series(df0['Q_pred'], index=df0.index) df=qq_mon(lyy,lmm,qq) print(fnameR) if iii==0: df1_tank=df if iii==1: df2_tank=df if iii==2: df1_mlp=df if iii==3: df2_mlp=df wexcel(df1_tank,df2_tank,df1_mlp,df2_mlp) #============== # Execution #============== if __name__ == '__main__': main()
以 上
GraphvizでMLPと通常の重回帰分析のイメージ図を作る
作例
GraphvizでMLPと通常重回帰分析のイメージ図(conceptual diagram あるいは schematic diagram かな)を作ってみた。 作例は以下の通り。
プログラミング環境は以下の通り。
Graphvizのインストール
まずはGraphvizのインストール。 brewとpipで2回行う。
brew install graphviz pip install graphviz
MLPのイメージ図作成プログラム
node の表示順を制御するため、gi.node('1','x[0]' ,pos='0,0!')
というように pos
を追記している。
from graphviz import Digraph g = Digraph(format='png') g.attr(rankdir='LR') g.attr(splines='false') g.attr(dpi='300') # Input subgraph gi=Digraph(name='cluster_i') gi.attr(label='Input') gi.attr(penwidth='0') gi.node('1','x[0]' ,pos='0,0!') gi.node('2','x[1]' ,pos='0,1!') gi.node('3','x[..]' ,pos='0,2!') gi.node('4','x[n]',pos='0,3!') # Hidden subgraph 1 gh1=Digraph(name='cluster_h1') gh1.attr(label='Hidden layer 1') gh1.attr(penwidth='0') gh1.node('5', 'h1[0]' ,pos='0,0!') gh1.node('6', 'h1[1]' ,pos='0,1!') gh1.node('7', 'h1[..]' ,pos='0,2!') gh1.node('8', 'h1[n]',pos='0,3!') # Hidden subgraph 2 gh2=Digraph(name='cluster_h2') gh2.attr(label='Hidden layer 2') gh2.attr(penwidth='0') gh2.node('9', 'h2[0]' ,pos='0,0!') gh2.node('10','h2[1]' ,pos='0,1!') gh2.node('11','h2[..]' ,pos='0,2!') gh2.node('12','h2[n]',pos='0,3!') # Output subgraph go=Digraph(name='cluster_o') go.attr(label='Output') go.attr(penwidth='0') go.node('13','y',pos='0,1.5!') g.subgraph(gi) g.subgraph(gh1) g.subgraph(gh2) g.subgraph(go) g.edge('1', '5') g.edge('1', '6') g.edge('1', '7') g.edge('1', '8') g.edge('2', '5') g.edge('2', '6') g.edge('2', '7') g.edge('2', '8') g.edge('3', '5') g.edge('3', '6') g.edge('3', '7') g.edge('3', '8') g.edge('4', '5') g.edge('4', '6') g.edge('4', '7') g.edge('4', '8') g.edge('5', '9') g.edge('5', '10') g.edge('5', '11') g.edge('5', '12') g.edge('6', '9') g.edge('6', '10') g.edge('6', '11') g.edge('6', '12') g.edge('7', '9') g.edge('7', '10') g.edge('7', '11') g.edge('7', '12') g.edge('8', '9') g.edge('8', '10') g.edge('8', '11') g.edge('8', '12') g.edge('9', '13') g.edge('10', '13') g.edge('11', '13') g.edge('12', '13') g.render('fig_3_mlp', view=True)
重回帰のイメージ図作成プログラム
edge(矢印線)に回帰係数を示す w[..]
を表示するため、g.edge('3', '5', label='w[..]')
というように label
を追記している。
from graphviz import Digraph g = Digraph(format='png') g.attr(rankdir='LR') g.attr(splines='false') g.attr(dpi='300') # Input subgraph gi=Digraph(name='cluster_i') gi.attr(label='Input') gi.attr(penwidth='0') gi.node('1','x[0]' ,pos='0,0!') gi.node('2','x[1]' ,pos='0,1!') gi.node('3','x[..]' ,pos='0,2!') gi.node('4','x[n]',pos='0,3!') # Output subgraph go=Digraph(name='cluster_o') go.attr(label='Output') go.attr(penwidth='0') go.node('5','y',pos='0,0!') g.subgraph(gi) g.subgraph(go) g.edge('1', '5', label='w[0]') g.edge('2', '5', label='w[1]') g.edge('3', '5', label='w[..]') g.edge('4', '5', label='w[n]') g.render('fig_3_reg', view=True)
以 上
Structural design formulas for penstocks embedded in rock
The structural design formulas for penstocks embedded in rock which are shown in 'Technical standards for gates and penstock' in Japan are described.
Allowable stresss of steel material
Welded joint efficiency
Location of welding | RT or UT carrying out for more than 5% of welded joint length |
RT or UT carrying out for less than 5% of welded joint length |
---|---|---|
Factory welding | 95% | 85% |
Site welding | 90% | 80% |
Required safety factor for buckling due to external pressure
Pipe shell
Stiffener
Design internal pressure
Definition of dimensions
Corrosion allowance
An allowance of more than or equal to 1.5mm shall be provided against corrossion and wear.
Minimum shell thickness
The minimum shell thickness used for the pressure lining part shall be more than or equal to those determined from following formula, if stiffeners are not used.
The minimum shell thickness shall not be less than 6mm even if the pipe diameter is small and stiffeners are used.
Stress calculation
Tensile stress in circumferential direction due to internal pressure
General
Internal pressure shared design by bedrock
Tensile stress in penstock axis direction
Stress due to restraint by stiffener
Stress due to temperature change
Stress due to Poisson's effect
Calculation for external presure
Without stiffener (E. Amstutz's formula)
With stiffener
Pipe shell proper (S. Timoshenko's formula)
Stiffener (E. Amstutz's formula)
That's all.